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

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