# source:git/Singular/LIB/signcond.lib@edfd36

spielwiese
Last change on this file since edfd36 was edfd36, checked in by Hans Schönemann <hannes@…>, 18 years ago
• Property mode set to `100644`
File size: 10.1 KB
Line
1// \$Id: signcond.lib,v 1.3 2005-05-06 11:17:33 Singular Exp \$
2// E. Tobis  12.Nov.2004
3// last change 16. Apr. 2005 (G.-M. Greuel)
4///////////////////////////////////////////////////////////////////////////////
5category="Symbolic-numerical solving"
6info="
7LIBRARY: signcond.lib Routines for computing realizable sign conditions
8AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
9
10OVERVIEW:  Routines to determine the number of solutions of a multivariate
11           polynomial system which satisfy a given sign configuration.
12           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
13           Geometry\", Springer, 2003.
14
15PROCEDURES:
16  signcnd(P,I)   The sign conditions realized by polynomials of P on a V(I)
17  psigncnd(P,l)  Pretty prints the output of signcnd (l)
18  firstoct(I)    The number of elements of V(I) with every coordinate > 0
19
20KEYWORDS: real roots,sign conditions
21";
22
23LIB "mrrcount.lib";
24LIB "linalg.lib";
25///////////////////////////////////////////////////////////////////////////////
26
27proc firstoct(ideal I)
28"USAGE:    firstoct(i); i ideal
29RETURN:   number: the number of points of V(i) lying in the first octant
30ASSUME:   i is a Groebner basis
32EXAMPLE:  example firstoct; shows an example"
33{
34  ideal firstoctant;
35  int j;
36  list result;
37  int n;
38
39  if (isparam(I)) {
40    ERROR("This procedure cannot operate with parametric arguments");
41  }
42
43  for (j = nvars(basering);j > 0;j--) {
44    firstoctant = firstoctant + var(j);
45  }
46
47  result = signcnd(firstoctant,I);
48
49  list fst;
50  for (j = nvars(basering);j > 0;j--) {
51    fst[j] = 1;
52  }
53
54  n = isIn(fst,result[1]);
55
56  if (n != -1) {
57    return (result[2][n]);
58  } else {
59    return (0);
60  }
61}
62example
63{
64  echo = 2;
65  ring r = 0,(x,y),dp;
66  ideal i = (x-2)*(x+3)*x,y*(y-1);
67  firstoct(i);
68}
69///////////////////////////////////////////////////////////////////////////////
70
71proc signcnd(ideal P,ideal I)
72"USAGE:     signcnd(P,I); ideal P,I
73RETURN:    list: the sign conditions realized by the polynomials of P on V(I).
74           The output of signcnd is a list of two lists. Both lists have the
75           same length. That length is the number of sign conditions realized
76           by the polynomials of P on the set V(i).
77           Each element of the first list indicates a sign condition of the
78           polynomials of P.
79           Each element of the second list indicates how many elements of V(I)
80           give rise to the sign condition expressed by the same position on
81           the first list.
82           See the example for further explanation of the output.
83ASSUME:    I is a Groebner basis
84NOTE:      The procedure psigncnd performs some pretty printing of this output
86EXAMPLE:   example signcnd; shows an example"
87{
88  ideal B;
89
90  // Cumulative stuff
91  matrix M;
92  matrix SQs;
93  matrix C;
94  list Signs;
95  list Exponents;
96
97  // Used to store the precalculated SQs
98  list SQvalues;
99  list SQpositions;
100
101  int i;
102
103  // Variables for each step
104  matrix Mi;
105  matrix M3x3[3][3];
106  matrix M3x3inv[3][3]; // Constant matrices
107  matrix c[3][1];
108  matrix sq[3][1];
109  int j;
110  list exponentsi;
111  list signi;
112  int numberOfNonZero;
113
114  if (isparam(P) || isparam(I)) {
115    ERROR("This procedure cannot operate with parametric arguments");
116  }
117
118  M3x3 = matrix(1,3,3);
119  M3x3 = 1,1,1,0,1,-1,0,1,1; // The 3x3 matrix
120  M3x3inv = inverse(M3x3);
121
122  // First, we compute sturmquery(1,V(I))
123  I = groebner(I);
124  B = qbase(I);
125  sq[1,1] = sturmquery(1,B,I); // Number of real roots in V(I)
126  SQvalues = SQvalues + list(sq[1,1]);
127  SQpositions = SQpositions + list(1);
128
129  // We initialize the cumulative variables
130  M = matrix(1,1,1);
131  Exponents = list(list());
132  Signs = list(list());
133
134  i = 1;
135
136  while (i <= size(P)) { // for each poly in P
137
138    sq[2,1] = sturmquery(P[i],B,I);
139    sq[3,1] = sturmquery(P[i]^2,B,I);
140
141
142    c = M3x3inv*sq;
143
144    // We have to eliminate the 0 elements in c
145    exponentsi = list();
146    signi = list();
147
148
149    // We determine the list of signs which correspond to a nonzero
150    // number of roots
151    numberOfNonZero = 3;
152
153    if (c[1,1] != 0) {
154      signi = list(0);
155    } else {
156      numberOfNonZero--;
157    }
158
159    if (c[2,1] != 0) {
160      signi = signi + list(1);
161    } else {
162      numberOfNonZero--;
163    }
164
165    if (c[3,1] != 0) {
166      signi = signi + list(-1);
167    } else {
168      numberOfNonZero--;
169    }
170
171    // We now determine the little matrix we'll work with,
172    // and the list of exponents
173    if (numberOfNonZero == 3) {
174      Mi = M3x3;
175      exponentsi = list(0,1,2);
176    } else {if (numberOfNonZero == 2) {
177      Mi = matrix(1,2,2);
178      Mi[1,2] = 1;
179      if (c[1,1] != 0 && c[2,1] != 0) { // 0,1
180        Mi[2,1] = 0;
181        Mi[2,2] = 1;
182      } else {if (c[1,1] != 0 && c[3,1] != 0) { // 0,-1
183        Mi[2,1] = 0;
184        Mi[2,2] = -1;
185      } else { // 1,-1
186        Mi[2,1] = 1;
187        Mi[2,2] = -1;
188      }}
189      exponentsi = list(0,1);
190    } else {if (numberOfNonZero == 1) {
191      Mi = matrix(1,1,1);
192      exponentsi = list(0);
193    }}}
194
195    // We store the Sturm Queries we'll need later
196    if (numberOfNonZero == 2) {
197      SQvalues = SQvalues + list(sq[2,1]);
198      SQpositions = SQpositions + list(size(Exponents)+1);
199    } else {if (numberOfNonZero == 3) {
200      SQvalues = SQvalues + list(sq[2,1],sq[3,1]);
201      SQpositions = SQpositions + list(size(Exponents)+1,size(Exponents)*2+1);
202    }}
203
204    // Now, we accumulate information
205    M = tensor(Mi,M);
206    Signs = expprod(Signs,signi);
207    Exponents = expprod(Exponents,exponentsi);
208
209    i++;
210  }
211
212  // At this point, we have the cumulative matrix,
213  // the vector of exponents and the matching sign conditions.
214  // We have to solve the big linear system to finish.
215
216  M = inverse(M);
217
218  // We have to compute the constants vector (the Sturm Queries)
219
220  SQs = matrix(1,size(Exponents),1);
221
222  j = 1; // We'll iterate over the presaved SQs
223
224  for (i = 1;i <= size(Exponents);i++) {
225    if (j <= size(SQvalues)) {
226      if (SQpositions[j] == i) {
227        SQs[i,1] = SQvalues[j];
228        j++;
229      } else {
230      SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
231      }
232    } else {
233        SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
234    }
235  }
236
237  C = M*SQs;
238
239  list result;
240  result[2] = list();
241  result[1] = list();
242
243  // We have to filter the 0 elements of C
244  for (i = 1;i <= size(Signs);i++) {
245    if (C[i,1] != 0) {
246      result[1] = result[1] + list(Signs[i]);
247      result[2] = result[2] + list(C[i,1]);
248    }
249  }
250
251  return (result);
252}
253example
254{ echo = 2;
255  ring r = 0,(x,y),dp;
256  ideal i = (x-2)*(x+3)*x,y*(y-1);
257  ideal P = x,y;
258  list l = signcnd(P,i);
259
260  size(l[1]);     // = the number of sign conditions of P on V(i)
261
262  //Each element of l[1] indicates a sign condition of the polynomials of P.
263  //The following means P[1] > 0, P[2] = 0:
264  l[1][2];
265
266  //Each element of l[2] indicates how many elements of V(I) give rise to
267  //the sign condition expressed by the same position on the first list.
268  //The following means that exactly 1 element of V(I) gives rise to the
269  //condition P[1] > 0, P[2] = 0:
270  l[2][2];
271}
272///////////////////////////////////////////////////////////////////////////////
273
274proc psigncnd(ideal P,list l)
275"USAGE:     psigncnd(P,l); ideal P, list l
276RETURN:    list: a formatted version of l
278EXAMPLE:   example psigncnd; shows an example"
279{
280  string s;
281  int n = size(l[1]);
282  int i;
283
284  for (i = 1;i <= n;i++) {
285    s = s + string(l[2][i]) + " elements of V(I) satisfy " + psign(P,l[1][i])
286        + sprintf("%n",12);
287  }
288  return(s);
289}
290example
291{
292  echo = 2;
293  ring r = 0,(x,y),dp;
294  ideal i = (x-2)*(x+3)*x,(y-1)*(y+2)*(y+4);
295  ideal P = x,y;
296  list l = signcnd(P,i);
297  psigncnd(P,l);
298}
299///////////////////////////////////////////////////////////////////////////////
300
301static proc psign(ideal P,list s)
302{
303  int i;
304  int n = size(P);
305  string output;
306
307  output = "{P[1]";
308
309  if (s[1] == -1) {
310    output = output + " < 0";
311  };
312  if (s[1] == 0) {
313    output = output + " = 0";
314  };
315  if (s[1] == 1) {
316    output = output + " > 0";
317  };
318
319  for (i = 2;i <= n;i++) {
320    output = output + ",";
321    output = output + "P[" + string(i) + "]";
322    if (s[i] == -1) {
323      output = output + " < 0";
324    };
325    if (s[i] == 0) {
326      output = output + " = 0";
327    };
328    if (s[i] == 1) {
329      output = output + " > 0";
330    };
331
332  }
333  output = output + "}";
334  return (output);
335}
336///////////////////////////////////////////////////////////////////////////////
337
338static proc isIn(list a,list b) //a is a list. b is a list of lists
339{
340  int i,j;
341  int found;
342
343  found = 0;
344  i = 1;
345  while (i <= size(b) && !found) {
346    j = 1;
347    found = 1;
348    if (size(a) != size(b[i])) {
349      found = 0;
350    } else {
351      while(j <= size(a)) {
352        found = found && a[j] == b[i][j];
353        j++;
354      }
355    }
356    i++;
357  }
358
359  if (found) {
360    return (i-1);
361  } else {
362    return (-1);
363  }
364}
365///////////////////////////////////////////////////////////////////////////////
366
367static proc expprod(list A,list B) // Computes the product of the list of lists A and the list B.
368{
369  int i,j;
370  list result;
371  int la,lb;
372
373  if (size(A) == 0) {
374    A = list(list());
375  }
376
377  la = size(A);
378  lb = size(B);
379
380  result[la*lb] = 0;
381
382
383  for (i = 0;i < lb;i++) {
384    for (j = 0;j < la;j++) {
385      result[i*la+j+1] = A[j+1] + list(B[i+1]);
386    }
387  }
388
389  return (result);
390}
391///////////////////////////////////////////////////////////////////////////////
392
393static proc initlist(int n) // Returns an n-element list of 0s.
394{
395  list l;
396  int i;
397  l[n] = 0;
398  for (i = 1;i < n;i++) {
399    l[i] = 0;
400  }
401  return(l);
402}
403///////////////////////////////////////////////////////////////////////////////
404
405static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate
406{
407  int i;
408  int n;
409  poly result;
410
411  n = size(exp);
412  result = 1;
413
414  for (i = 1;i <= n; i++) {
415    result = result * (P[i]^exp[i]);
416  }
417  return (result);
418}
419///////////////////////////////////////////////////////////////////////////////
420
421static proc incexp(list exp)
422{
423  int k;
424
425  k = 1;
426
427  while (exp[k] == 2) { // We assume exp is not the last exponent (i.e. 2,...,2)
428    exp[k] = 0;
429    k++;
430  }
431
432  // exp[k] < 2
433  exp[k] = exp[k] + 1;
434
435  return (exp);
436}
437///////////////////////////////////////////////////////////////////////////////
438
439
Note: See TracBrowser for help on using the repository browser.