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