source: git/Singular/LIB/signcond.lib @ 380a17b

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