1 | // $Id: rootsmr.lib,v 1.4 2008-03-18 15:51:41 Singular Exp $ |
---|
2 | // E. Tobis 12.Nov.2004, April 2004 |
---|
3 | // last change 7. May 2005 (G.-M. Greuel) |
---|
4 | /////////////////////////////////////////////////////////////////////////////// |
---|
5 | category="Teaching" |
---|
6 | info=" |
---|
7 | LIBRARY: rootsmr.lib Counting the number of real roots of polynomial systems |
---|
8 | AUTHOR: Enrique A. Tobis, etobis@dc.uba.ar |
---|
9 | |
---|
10 | OVERVIEW: Routines for counting the number of real roots of a multivariate |
---|
11 | polynomial system. Two methods are implemented: deterministic |
---|
12 | computation of the number of roots, via the signature of a certain |
---|
13 | bilinear form (nrRootsDeterm); and a rational univariate projection, |
---|
14 | using a pseudorandom polynomial (nrRootsProbab). It also includes a |
---|
15 | command to verify the correctness of the pseudorandom answer. |
---|
16 | References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic |
---|
17 | Geometry\", Springer, 2003. |
---|
18 | |
---|
19 | PROCEDURES: |
---|
20 | nrRootsProbab(I) Number of real roots of 0-dim ideal (probabilistic) |
---|
21 | nrRootsDeterm(I) Number of real roots of 0-dim ideal (deterministic) |
---|
22 | symsignature(m) Signature of the symmetric matrix m |
---|
23 | sturmquery(h,B,I) Sturm query of h on V(I) |
---|
24 | matbil(h,B,I) Matrix of the bilinear form on R/I associated to h |
---|
25 | matmult(f,B,I) Matrix of multiplication by f (m_f) on R/I in the basis B |
---|
26 | tracemult(f,B,I) Trace of m_f (B is an ordered basis of R/I) |
---|
27 | coords(f,B,I) Coordinates of f in the ordered basis B |
---|
28 | randcharpoly(B,I,n) Pseudorandom charpoly of univ. projection, n optional |
---|
29 | verify(p,B,i) Verifies the result of randcharpoly |
---|
30 | randlinpoly(n) Pseudorandom linear polynomial, n optional |
---|
31 | powersums(f,B,I) Powersums of the roots of a char polynomial |
---|
32 | symmfunc(S) Symmetric functions from the powersums S |
---|
33 | univarpoly(l) Polynomial with coefficients from l |
---|
34 | qbase(i) Like kbase, but the monomials are ordered |
---|
35 | |
---|
36 | KEYWORDS: real roots, univariate projection |
---|
37 | "; |
---|
38 | /////////////////////////////////////////////////////////////////// |
---|
39 | LIB "linalg.lib"; // We use charpoly |
---|
40 | LIB "rootsur.lib"; // We use varsigns |
---|
41 | |
---|
42 | proc nrRootsProbab(ideal I, list #) |
---|
43 | "USAGE: nrRootsProbab(I,[n]); ideal I, int n |
---|
44 | RETURN: int: the number of real roots of the ideal I by a probabilistic |
---|
45 | algorithm |
---|
46 | ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed |
---|
47 | by using std. If I is already a Groebner basis (i.e. if |
---|
48 | attrib(I,"isSB"); returns 1) then this Groebner basis will be |
---|
49 | used, hence it must be one w.r.t. (any) global ordering. This may |
---|
50 | be useful if the ideal is known to be a Groebner basis or if it |
---|
51 | can be computed faster by a different method. |
---|
52 | NOTE: If n<10 is given, n is the number of digits being used for |
---|
53 | constructing a random characteristic polynomial, a bigger n is |
---|
54 | more safe but slower (default: n=5). |
---|
55 | If printlevel>0 the number of complex solutions is displayed |
---|
56 | (default: printlevel=0). |
---|
57 | SEE ALSO: nrroots, nrRootsDeterm, randcharpoly, solve |
---|
58 | EXAMPLE: example nrRootsProbab; shows an example" |
---|
59 | { |
---|
60 | //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). |
---|
61 | //Then the algorithm needs: |
---|
62 | //1 std(I) and ~n NF computations (of randcharpoly wrt I) |
---|
63 | |
---|
64 | if (isparam(I)) { |
---|
65 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
66 | } |
---|
67 | int pr = printlevel-voice+2; |
---|
68 | int v; |
---|
69 | int n=5; |
---|
70 | if (size(#) == 1) { |
---|
71 | n=#[1]; |
---|
72 | } |
---|
73 | if (attrib(I,"isSB")!=1) { |
---|
74 | I = std(I); |
---|
75 | } |
---|
76 | |
---|
77 | ideal b = qbase(I); |
---|
78 | v = size(b); |
---|
79 | if (v == 0) { |
---|
80 | ERROR("ideal is not 0-dimensional"); |
---|
81 | } |
---|
82 | dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); |
---|
83 | |
---|
84 | poly p = randcharpoly(b,I,n); |
---|
85 | |
---|
86 | return (nrroots(p)); |
---|
87 | } |
---|
88 | |
---|
89 | example |
---|
90 | { |
---|
91 | echo = 2; |
---|
92 | ring r = 0,(x,y,z),lp; |
---|
93 | ideal i = (x-1)*(x-2),(y-1)^3*(x-y),(z-1)*(z-2)*(z-3)^2; |
---|
94 | nrRootsProbab(i); //no of real roots (using internally std) |
---|
95 | |
---|
96 | i = groebner(i); //using the hilbert driven GB computation |
---|
97 | int pr = printlevel; |
---|
98 | printlevel = 2; |
---|
99 | nrRootsProbab(i); |
---|
100 | printlevel = pr; |
---|
101 | } |
---|
102 | /////////////////////////////////////////////////////////////////////////////// |
---|
103 | |
---|
104 | proc nrRootsDeterm(ideal I) |
---|
105 | "USAGE: nrRootsDeterm(I); ideal I |
---|
106 | RETURN: int: the number of real roots of the ideal I by a deterministic |
---|
107 | algorithm |
---|
108 | ASSUME: If I is not a Groebner basis, then a Groebner basis will be computed |
---|
109 | by using std. If I is already a Groebner basis (i.e. if |
---|
110 | attrib(I,"isSB"); returns 1) then this Groebner basis will be |
---|
111 | used, hence it must be one w.r.t. (any) global ordering. This may |
---|
112 | be useful if the ideal is known to be a Groebner basis or if it |
---|
113 | can be computed faster by a different method. |
---|
114 | NOTE: If printlevel>0 the number of complex solutions is displayed |
---|
115 | (default: printlevel=0). The procedure nrRootsProbab is usually faster. |
---|
116 | SEE ALSO: nrroots, nrRootsProbab, sturmquery, solve |
---|
117 | EXAMPLE: example nrRootsDeterm; shows an example" |
---|
118 | { |
---|
119 | //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)). |
---|
120 | //Then the algotithm needs: |
---|
121 | //1 std(I) and (1/2)n*(n+1)^2 ~ 1/2n^3 NF computations (of monomials wrt I) |
---|
122 | |
---|
123 | if (isparam(I)) { |
---|
124 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
125 | } |
---|
126 | int pr = printlevel-voice+2; |
---|
127 | int v; |
---|
128 | |
---|
129 | if (attrib(I,"isSB")!=1) { |
---|
130 | I = std(I); |
---|
131 | } |
---|
132 | |
---|
133 | ideal b = qbase(I); |
---|
134 | v = size(b); |
---|
135 | if (v == 0) { |
---|
136 | ERROR("ideal is not 0-dimensional"); |
---|
137 | } |
---|
138 | dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity"); |
---|
139 | |
---|
140 | return (sturmquery(1,b,I)); |
---|
141 | } |
---|
142 | |
---|
143 | example |
---|
144 | { |
---|
145 | echo = 2; |
---|
146 | ring r = 0,(x,y,z),lp; |
---|
147 | ideal I = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; |
---|
148 | nrRootsDeterm(I); //no of real roots (using internally std) |
---|
149 | |
---|
150 | I = groebner(I); //using the hilbert driven GB computation |
---|
151 | int pr = printlevel; |
---|
152 | printlevel = 2; |
---|
153 | nrRootsDeterm(I); |
---|
154 | printlevel = pr; |
---|
155 | } |
---|
156 | /////////////////////////////////////////////////////////////////////////////// |
---|
157 | |
---|
158 | proc symsignature(matrix m) |
---|
159 | "USAGE: symsignature(m); m matrix. m must be symmetric. |
---|
160 | RETURN: int: the signature of m |
---|
161 | SEE ALSO: matbil,sturmquery |
---|
162 | EXAMPLE: example symsignature; shows an example" |
---|
163 | { |
---|
164 | int positive, negative, i, j; |
---|
165 | list l; |
---|
166 | poly variable; |
---|
167 | |
---|
168 | if (isparam(m)) { |
---|
169 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
170 | } |
---|
171 | |
---|
172 | if (!isSquare(m)) { |
---|
173 | ERROR ("m must be a square matrix"); |
---|
174 | } |
---|
175 | |
---|
176 | // We check whether m is symmetric |
---|
177 | for (i = 1;i <= nrows(m);i++) { |
---|
178 | for (j = i;j <= nrows(m);j++) { |
---|
179 | if (m[i,j] != m[j,i]) { |
---|
180 | ERROR ("m must be a symmetric matrix"); |
---|
181 | } |
---|
182 | } |
---|
183 | } |
---|
184 | |
---|
185 | poly f = charpoly(m); // Uses the last variable of the ring |
---|
186 | |
---|
187 | for (i = size(f);i >= 1;i--) { |
---|
188 | l[i] = leadcoef(f[i]); |
---|
189 | } |
---|
190 | positive = varsigns(l); |
---|
191 | |
---|
192 | variable = var(nvars(basering)); // charpoly uses the last variable |
---|
193 | f = subst(f,variable,-variable); |
---|
194 | |
---|
195 | for (i = size(f);i >= 1;i--) { |
---|
196 | l[i] = leadcoef(f[i]); |
---|
197 | } |
---|
198 | |
---|
199 | negative = varsigns(l); |
---|
200 | return (positive - negative); |
---|
201 | } |
---|
202 | example |
---|
203 | { |
---|
204 | echo = 2; |
---|
205 | ring r = 0,(x,y),dp; |
---|
206 | ideal i = x4-y2x,y2-13; |
---|
207 | i = std(i); |
---|
208 | ideal b = qbase(i); |
---|
209 | |
---|
210 | matrix m = matbil(1,b,i); |
---|
211 | symsignature(m); |
---|
212 | } |
---|
213 | /////////////////////////////////////////////////////////////////////////////// |
---|
214 | |
---|
215 | proc sturmquery(poly h,ideal B,ideal I) |
---|
216 | "USAGE: sturmquery(h,b,i); h poly, b,i ideal |
---|
217 | RETURN: int: the Sturm query of h in V(i) |
---|
218 | ASSUME: i is a Groebner basis, b is an ordered monomial basis |
---|
219 | of r/i, r = basering. |
---|
220 | SEE ALSO: symsignature,matbil |
---|
221 | EXAMPLE: example sturmquery; shows an example" |
---|
222 | { |
---|
223 | if (isparam(h) || isparam(B) || isparam(I)) { |
---|
224 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
225 | } |
---|
226 | |
---|
227 | return (mysymmsig(matbil(h,B,I))); |
---|
228 | } |
---|
229 | example |
---|
230 | { |
---|
231 | echo = 2; |
---|
232 | ring r = 0,(x,y),dp; |
---|
233 | ideal i = x4-y2x,y2-13; |
---|
234 | i = std(i); |
---|
235 | ideal b = qbase(i); |
---|
236 | |
---|
237 | sturmquery(1,b,i); |
---|
238 | } |
---|
239 | /////////////////////////////////////////////////////////////////////////////// |
---|
240 | |
---|
241 | static proc mysymmsig(matrix m) |
---|
242 | // returns the signature of a square symmetric matrix m |
---|
243 | { |
---|
244 | int positive, negative, i; |
---|
245 | list l; |
---|
246 | poly variable; |
---|
247 | |
---|
248 | poly f = charpoly(m); // Uses the last variable of the ring |
---|
249 | |
---|
250 | for (i = size(f);i >= 1;i--) { |
---|
251 | l[i] = leadcoef(f[i]); |
---|
252 | } |
---|
253 | positive = varsigns(l); |
---|
254 | |
---|
255 | variable = var(nvars(basering)); // charpoly uses the last variable |
---|
256 | f = subst(f,variable,-variable); |
---|
257 | |
---|
258 | for (i = size(f);i >= 1;i--) { |
---|
259 | l[i] = leadcoef(f[i]); |
---|
260 | } |
---|
261 | |
---|
262 | negative = varsigns(l); |
---|
263 | return (positive - negative); |
---|
264 | } |
---|
265 | /////////////////////////////////////////////////////////////////////////////// |
---|
266 | |
---|
267 | proc matbil(poly h,ideal B,ideal I) |
---|
268 | "USAGE: matbil(h,b,i); h poly, b,i ideal |
---|
269 | RETURN: matrix: the matrix of the bilinear form (f,g) |-> trace(m_fhg), |
---|
270 | m_fhg = multiplication with fhg on r/i |
---|
271 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
272 | r = basering |
---|
273 | SEE ALSO: matmult,tracemult |
---|
274 | EXAMPLE: example matbil; shows an example" |
---|
275 | { |
---|
276 | matrix m[size(B)][size(B)]; |
---|
277 | poly f; |
---|
278 | int k,l; |
---|
279 | //h = reduce(h,I); |
---|
280 | |
---|
281 | for (k = 1; k <= size(B); k++) { |
---|
282 | for (l = 1; l <= k; l++) { |
---|
283 | m[k,l] = tracemult(h*B[k]*B[l],B,I)[1]; |
---|
284 | m[l,k] = m[k,l]; // The matrix we are trying to compute is symmetric |
---|
285 | } |
---|
286 | } |
---|
287 | return(m); |
---|
288 | } |
---|
289 | example |
---|
290 | { |
---|
291 | echo = 2; |
---|
292 | ring r = 0,(x,y),dp; |
---|
293 | ideal i = x4-y2x,y2-13; |
---|
294 | i = std(i); |
---|
295 | ideal b = qbase(i); |
---|
296 | poly f = x3-xy+y-13+x4-y2x; |
---|
297 | |
---|
298 | matrix m = matbil(f,b,i); |
---|
299 | print(m); |
---|
300 | |
---|
301 | } |
---|
302 | /////////////////////////////////////////////////////////////////////////////// |
---|
303 | |
---|
304 | proc tracemult(poly f,ideal B,ideal I) |
---|
305 | "USAGE: tracemult(f,b,i);f poly, b,i ideal |
---|
306 | RETURN: number: the trace of the multiplication by f (m_f) on r/i, |
---|
307 | written in the monomial basis b of r/i, r = basering |
---|
308 | (faster than matmult + trace) |
---|
309 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i |
---|
310 | SEE ALSO: matmult,trace |
---|
311 | EXAMPLE: example tracemult; shows an example" |
---|
312 | { |
---|
313 | int k; // Iterates over the basis monomials |
---|
314 | int l; // Iterates over the rows of the matrix |
---|
315 | list coordinates; |
---|
316 | number m; |
---|
317 | poly g; |
---|
318 | |
---|
319 | //f = reduce(f,I); |
---|
320 | for (k = 1; k <= size(B); k++) { |
---|
321 | l=1; |
---|
322 | g = reduce(f*B[k],I); |
---|
323 | while (l <= k) { |
---|
324 | if (leadmonom(g[l]) == B[k]) { |
---|
325 | m = m + leadcoef(g[l]); |
---|
326 | break; |
---|
327 | } |
---|
328 | l++; |
---|
329 | } |
---|
330 | } |
---|
331 | return (m); |
---|
332 | } |
---|
333 | example |
---|
334 | { |
---|
335 | echo = 2; |
---|
336 | ring r = 0,(x,y),dp; |
---|
337 | ideal i = x4-y2x,y2-13; |
---|
338 | i = std(i); |
---|
339 | ideal b = qbase(i); |
---|
340 | |
---|
341 | poly f = x3-xy+y-13+x4-y2x; |
---|
342 | matrix m = matmult(f,b,i); |
---|
343 | print(m); |
---|
344 | |
---|
345 | tracemult(f,b,i); //the trace of m |
---|
346 | } |
---|
347 | /////////////////////////////////////////////////////////////////////////////// |
---|
348 | |
---|
349 | proc matmult(poly f, ideal B, ideal I) |
---|
350 | "USAGE: matmult(f,b,i); f poly, b,i ideal |
---|
351 | RETURN: matrix: the matrix of the multiplication map by f (m_f) on r/i |
---|
352 | w.r.t. to the monomial basis b of r/i (r = basering) |
---|
353 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
354 | as given by qbase(i) |
---|
355 | SEE ALSO: coords,matbil |
---|
356 | EXAMPLE: example matmult; shows an example" |
---|
357 | { |
---|
358 | int k; // Iterates over the basis monomials |
---|
359 | int l; // Iterates over the rows of the matrix |
---|
360 | list coordinates; |
---|
361 | matrix m[size(B)][size(B)]; |
---|
362 | |
---|
363 | //f = reduce(f,I); |
---|
364 | for (k = 1;k <= size(B);k++) { |
---|
365 | coordinates = coords(f*(B[k]),B,I); // f*x_k written on the basis B |
---|
366 | for (l = 1;l <= size(B);l++) { |
---|
367 | m[l,k] = coordinates[l]; |
---|
368 | } |
---|
369 | } |
---|
370 | return (m); |
---|
371 | } |
---|
372 | example |
---|
373 | { |
---|
374 | echo = 2; |
---|
375 | ring r = 0,(x,y),dp; |
---|
376 | ideal i = x4-y2x,y2-13; |
---|
377 | i = std(i); |
---|
378 | ideal b = qbase(i); |
---|
379 | |
---|
380 | poly f = x3-xy+y-13+x4-y2x; |
---|
381 | matrix m = matmult(f,b,i); |
---|
382 | print(m); |
---|
383 | } |
---|
384 | /////////////////////////////////////////////////////////////////////////////// |
---|
385 | |
---|
386 | proc coords(poly f,ideal B,ideal I) |
---|
387 | "USAGE: coords(f,b,i), f poly, b,i ideal |
---|
388 | RETURN: list of numbers: the coordinates of the class of f (mod i) |
---|
389 | in the monomial basis b |
---|
390 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
391 | r = basering |
---|
392 | SEE ALSO: matmult,matbil |
---|
393 | KEYWORDS: coordinates |
---|
394 | EXAMPLE: example coords; shows an example" |
---|
395 | { |
---|
396 | // We assume the basis is sorted according to the ring order |
---|
397 | poly g; |
---|
398 | int k,l=1,1; |
---|
399 | list coordinates; |
---|
400 | int N = size(B); |
---|
401 | |
---|
402 | // We first compute the normal form of f wrt I |
---|
403 | g = reduce(f,I); |
---|
404 | int n = size(g); //allways n <= N |
---|
405 | |
---|
406 | while (k <= N) { |
---|
407 | if (leadmonom(g[l]) == B[k]) { |
---|
408 | coordinates[k] = leadcoef(g[l]); |
---|
409 | l++; |
---|
410 | } else { |
---|
411 | coordinates[k] = number(0); |
---|
412 | } |
---|
413 | k++; |
---|
414 | } |
---|
415 | return (coordinates); |
---|
416 | } |
---|
417 | example |
---|
418 | { |
---|
419 | echo = 2; |
---|
420 | ring r = 0,(x,y),dp; |
---|
421 | ideal i = x4-y2x,y2-13; |
---|
422 | poly f = x3-xy+y-13+x4-y2x; |
---|
423 | i = std(i); |
---|
424 | ideal b = qbase(i); |
---|
425 | b; |
---|
426 | coords(f,b,i); |
---|
427 | } |
---|
428 | /////////////////////////////////////////////////////////////////////////////// |
---|
429 | |
---|
430 | static proc isSquare(matrix m) |
---|
431 | // returns 1 iff m is a square matrix |
---|
432 | { |
---|
433 | return (nrows(m)==ncols(m)); |
---|
434 | } |
---|
435 | /////////////////////////////////////////////////////////////////////////////// |
---|
436 | |
---|
437 | proc randcharpoly(ideal B,ideal I,list #) |
---|
438 | "USAGE: randcharpoly(b,i); randcharpoly(b,i,n); b,i ideal; n int |
---|
439 | RETURN: poly: the characteristic polynomial of a pseudorandom |
---|
440 | rational univariate projection having one zero per zero of i. |
---|
441 | If n<10 is given, it is the number of digits being used for the |
---|
442 | pseudorandom coefficients (default: n=5) |
---|
443 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
444 | r = basering |
---|
445 | NOTE: shows a warning if printlevel>0 (default: printlevel=0) |
---|
446 | KEYWORDS: rational univariate projection |
---|
447 | EXAMPLE: example randcharpoly; shows an example" |
---|
448 | { |
---|
449 | int pr = printlevel - voice + 2; |
---|
450 | poly p; |
---|
451 | poly generic; |
---|
452 | list l; |
---|
453 | matrix m; |
---|
454 | poly q; |
---|
455 | |
---|
456 | if (size(#) == 1) { |
---|
457 | generic = randlinpoly(#[1]); |
---|
458 | } else { |
---|
459 | generic = randlinpoly(); |
---|
460 | } |
---|
461 | |
---|
462 | p = reduce(generic,I); |
---|
463 | m = matmult(p,B,I); |
---|
464 | q = charpoly(m); |
---|
465 | |
---|
466 | dbprint(pr,"*********************************************************************"); |
---|
467 | dbprint(pr,"* WARNING: This polynomial was obtained using pseudorandom numbers.*"); |
---|
468 | dbprint(pr,"* If you want to verify the result, please use the command *"); |
---|
469 | dbprint(pr,"* *"); |
---|
470 | dbprint(pr,"* verify(p,b,i) *"); |
---|
471 | dbprint(pr,"* *"); |
---|
472 | dbprint(pr,"* where p is the polynomial I returned, b is the monomial basis *"); |
---|
473 | dbprint(pr,"* used, and i the Groebner basis of the ideal *"); |
---|
474 | dbprint(pr,"*********************************************************************"); |
---|
475 | |
---|
476 | return(q); |
---|
477 | } |
---|
478 | example |
---|
479 | { |
---|
480 | echo = 2; |
---|
481 | ring r = 0,(x,y,z),dp; |
---|
482 | ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2; |
---|
483 | i = std(i); |
---|
484 | ideal b = qbase(i); |
---|
485 | poly p = randcharpoly(b,i); |
---|
486 | p; |
---|
487 | nrroots(p); // See nrroots in urrcount.lib |
---|
488 | |
---|
489 | int pr = printlevel; |
---|
490 | printlevel = pr+2; |
---|
491 | p = randcharpoly(b,i,5); |
---|
492 | nrroots(p); |
---|
493 | printlevel = pr; |
---|
494 | } |
---|
495 | |
---|
496 | /////////////////////////////////////////////////////////////////////////////// |
---|
497 | |
---|
498 | proc verify(poly p,ideal b,ideal i) |
---|
499 | "USAGE: verify(p,b,i);p poly, b,i,ideal |
---|
500 | RETURN: integer: 1 iff the polynomial p splits the points of V(i). |
---|
501 | It's used to check the result of randcharpoly |
---|
502 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
503 | r = basering |
---|
504 | NOTE: comments the result if printlevel>0 (default: printlevel=0) |
---|
505 | SEE ALSO: randcharpoly |
---|
506 | EXAMPLE: example verify; shows an example" |
---|
507 | { |
---|
508 | int pr = printlevel - voice + 2; |
---|
509 | poly sqr_free; |
---|
510 | int correct; |
---|
511 | poly variable; |
---|
512 | |
---|
513 | if (isparam(p) || isparam(b) || isparam(i)) { |
---|
514 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
515 | } |
---|
516 | |
---|
517 | variable = isuni(p); |
---|
518 | sqr_free = p/gcd(p,diff(p,variable)); |
---|
519 | correct = (mat_rk(matbil(1,b,i)) == deg(sqr_free)); |
---|
520 | |
---|
521 | if (correct) { |
---|
522 | dbprint(pr,"//Verification successful"); |
---|
523 | } else { |
---|
524 | dbprint(pr,"//The choice of random numbers was not useful"); |
---|
525 | dbprint(pr,"//You might want to try randcharpoly with a larger number of digits"); |
---|
526 | } |
---|
527 | return (correct); |
---|
528 | } |
---|
529 | example |
---|
530 | { |
---|
531 | echo = 2; |
---|
532 | ring r = 0,(x,y),dp; |
---|
533 | poly f = x3-xy+y-13+x4-y2x; |
---|
534 | ideal i = x4-y2x,y2-13; |
---|
535 | i = std(i); |
---|
536 | ideal b = qbase(i); |
---|
537 | poly p = randcharpoly(b,i); |
---|
538 | verify(p,b,i); |
---|
539 | } |
---|
540 | /////////////////////////////////////////////////////////////////////////////// |
---|
541 | |
---|
542 | proc randlinpoly(list #) |
---|
543 | "USAGE: randlinpoly(); randlinpoly(n); n int |
---|
544 | RETURN: poly: linear combination of the variables of the ring, with |
---|
545 | pseudorandom coefficients. If n<10 is given, it is the number of |
---|
546 | digits being used for the range of the coefficients (default: n=5) |
---|
547 | SEE ALSO: randcharpoly; |
---|
548 | EXAMPLE: example randlinpoly; shows an example" |
---|
549 | { |
---|
550 | int n,i; |
---|
551 | poly p = 0; |
---|
552 | int ndigits = 5; |
---|
553 | |
---|
554 | if (size(#) == 1) { |
---|
555 | ndigits = #[1]; |
---|
556 | } |
---|
557 | |
---|
558 | n = nvars(basering); |
---|
559 | for (i = 1;i <= n;i++) { |
---|
560 | p = p + var(i)*random(1,10^ndigits); |
---|
561 | } |
---|
562 | return (p); |
---|
563 | } |
---|
564 | example |
---|
565 | { |
---|
566 | echo = 2; |
---|
567 | ring r = 0,(x,y,z,w),dp; |
---|
568 | poly p = randlinpoly(); |
---|
569 | p; |
---|
570 | randlinpoly(5); |
---|
571 | } |
---|
572 | /////////////////////////////////////////////////////////////////////////////// |
---|
573 | |
---|
574 | proc powersums(poly f,ideal B,ideal I) |
---|
575 | "USAGE: powersums(f,b,i); f poly; b,i ideal |
---|
576 | RETURN: list: the powersums of the results of evaluating f at the zeros of I |
---|
577 | ASSUME: i is a Groebner basis and b is an ordered monomial basis of r/i, |
---|
578 | r = basering |
---|
579 | SEE ALSO: symmfunc |
---|
580 | EXAMPLE: example symmfunc; shows an example" |
---|
581 | { |
---|
582 | int N,k; |
---|
583 | list sums; |
---|
584 | |
---|
585 | N = size(B); |
---|
586 | for (k = 1;k <= N;k++) { |
---|
587 | sums = sums + list(leadcoef(trace(matmult(f^k,B,I)))); |
---|
588 | } |
---|
589 | return (sums); |
---|
590 | } |
---|
591 | example |
---|
592 | { |
---|
593 | echo = 2; |
---|
594 | ring r = 0,(x,y,z),dp; |
---|
595 | |
---|
596 | ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5) |
---|
597 | i = std(i); |
---|
598 | |
---|
599 | ideal b = qbase(i); |
---|
600 | poly f = x+y+z; |
---|
601 | list psums = list(-2-3,4+9); // f evaluated at V(I) gives {-3,-2} |
---|
602 | list l = powersums(f,b,i); |
---|
603 | psums; |
---|
604 | l; |
---|
605 | } |
---|
606 | /////////////////////////////////////////////////////////////////////////////// |
---|
607 | |
---|
608 | proc symmfunc(list S) |
---|
609 | "USAGE: symmfunc(s); s list |
---|
610 | RETURN: list: the symmetric functions of the roots of a polynomial, given |
---|
611 | the power sums of those roots. |
---|
612 | SEE ALSO: powersums |
---|
613 | EXAMPLE: example symmfunc; shows an example" |
---|
614 | { |
---|
615 | // Takes the list of power sums and returns the symmetric functions |
---|
616 | list a; |
---|
617 | int j,l,N; |
---|
618 | number sum; |
---|
619 | |
---|
620 | N = size(S); |
---|
621 | a[N+1] = 1; // We set the length of the list and initialize its last element. |
---|
622 | |
---|
623 | for (l = N - 1;l >= 0;l--) { |
---|
624 | sum = 0; |
---|
625 | for (j = l + 1;j <= N;j++) { |
---|
626 | sum = sum + ((a[j+1])*(S[j-l])); |
---|
627 | } |
---|
628 | sum = -sum; |
---|
629 | a[l+1] = sum/(N-l); |
---|
630 | } |
---|
631 | |
---|
632 | a = reverse(a); |
---|
633 | return (a); |
---|
634 | } |
---|
635 | example |
---|
636 | { |
---|
637 | echo = 2; |
---|
638 | ring r = 0,x,dp; |
---|
639 | poly p = (x-1)*(x-2)*(x-3); |
---|
640 | list psums = list(1+2+3,1+4+9,1+8+27); |
---|
641 | list l = symmfunc(psums); |
---|
642 | l; |
---|
643 | p; // Compare p with the elements of l |
---|
644 | } |
---|
645 | /////////////////////////////////////////////////////////////////////////////// |
---|
646 | |
---|
647 | proc univarpoly(list l) |
---|
648 | "USAGE: univarpoly(l); l list |
---|
649 | RETURN: poly: a polynomial p on the first variable of basering, say x, |
---|
650 | with p = l[1] + l[2]*x + l[3]*x^2 + ... |
---|
651 | EXAMPLE: example univarpoly; shows an example" |
---|
652 | { |
---|
653 | poly p; |
---|
654 | int i,n; |
---|
655 | |
---|
656 | n = size(l); |
---|
657 | for (i = 1;i <= n;i++) { |
---|
658 | p = p + l[i]*var(1)^(n-i); |
---|
659 | } |
---|
660 | return (p); |
---|
661 | } |
---|
662 | example |
---|
663 | { |
---|
664 | echo = 2; |
---|
665 | ring r = 0,x,dp; |
---|
666 | list l = list(1,2,3,4,5); |
---|
667 | poly p = univarpoly(l); |
---|
668 | p; |
---|
669 | } |
---|
670 | /////////////////////////////////////////////////////////////////////////////// |
---|
671 | |
---|
672 | proc qbase(ideal i) |
---|
673 | "USAGE: qbase(I); I zero-dimensional ideal |
---|
674 | RETURN: ideal: A monomial basis of the quotient between the basering and the |
---|
675 | ideal I, sorted according to the basering order. |
---|
676 | SEE ALSO: kbase |
---|
677 | KEYWORDS: zero-dimensional |
---|
678 | EXAMPLE: example qbase; shows an example" |
---|
679 | { |
---|
680 | ideal b; |
---|
681 | |
---|
682 | b = kbase(i); |
---|
683 | b = reverseideal(sort(b)[1]); // sort sorts in ascending order |
---|
684 | return (b); |
---|
685 | } |
---|
686 | example |
---|
687 | { |
---|
688 | echo = 2; |
---|
689 | ring r = 0,(x,y,z),dp; |
---|
690 | |
---|
691 | ideal i = 2x2,-y2,z3; |
---|
692 | i = std(i); |
---|
693 | ideal b = qbase(i); |
---|
694 | b; |
---|
695 | b = kbase(i); |
---|
696 | b; // Compare this with the result of qbase |
---|
697 | } |
---|
698 | /////////////////////////////////////////////////////////////////////////////// |
---|
699 | |
---|
700 | static proc reverseideal(ideal b) // Returns b reversed |
---|
701 | { |
---|
702 | int i; |
---|
703 | ideal result; |
---|
704 | |
---|
705 | result = b[1]; |
---|
706 | for (i = 2;i <= size(b);i++) { |
---|
707 | result = b[i], result; |
---|
708 | } |
---|
709 | return (result); |
---|
710 | } |
---|
711 | /////////////////////////////////////////////////////////////////////////////// |
---|
712 | |
---|