1 | // $Id: urrcount.lib,v 1.4 2005-05-06 09:42:01 Singular Exp $ |
---|
2 | // E. Tobis 12.Nov.2004, April 2004 |
---|
3 | // last change 1. May 2005 (G.-M. Greuel) |
---|
4 | /////////////////////////////////////////////////////////////////////////////// |
---|
5 | category="Symbolic-numerical solving" |
---|
6 | info=" |
---|
7 | LIBRARY: urrcount.lib Counting number of real roots of univariate polynomial |
---|
8 | AUTHOR: Enrique A. Tobis, etobis@dc.uba.ar |
---|
9 | |
---|
10 | OVERVIEW: Routines for bounding and counting the number of real roots of a |
---|
11 | univariate polynomial, by means of several different methods, namely |
---|
12 | Descartes' rule of signs, the Budan-Fourier theorem, Sturm sequences |
---|
13 | and Sturm-Habicht sequences. The first two give bounds on the number |
---|
14 | of roots. The other two compute the actual number of roots of the |
---|
15 | polynomial. There are several wrapper functions, to simplify the |
---|
16 | application of the aforesaid theorems and some functions |
---|
17 | to determine whether a given polynomial is univariate. |
---|
18 | References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic |
---|
19 | Geometry\", Springer, 2003. |
---|
20 | |
---|
21 | |
---|
22 | PROCEDURES: |
---|
23 | isuni(p) Checks whether a polynomials is univariate |
---|
24 | whichvariable(p) The only variable of a univariate monomial (or 0) |
---|
25 | varsigns(p) Number of sign changes in a list |
---|
26 | boundBuFou(p,a,b) Bound for number of real roots of poly p in interval (a,b) |
---|
27 | boundposDes(p) Bound for the number of positive real roots of poly p |
---|
28 | boundDes(p) Bound for the number of real roots of poly p |
---|
29 | allrealst(p) Checks whether all the roots of a poly are real (via Sturm) |
---|
30 | maxabs(p) A bound for the maximum absolute value of a root of a poly |
---|
31 | allreal(p) Checks whether all the roots of a poly are real (via St-Ha) |
---|
32 | sturm(p,a,b) Number of real roots of a poly on an interval (via Sturm) |
---|
33 | sturmseq(p) Sturm sequence of a polynomial |
---|
34 | sturmha(p,a,b) Number of real roots of a poly in (a,b) (via Sturm-Habicht) |
---|
35 | sturmhaseq(p) A Sturm-Habicht Sequence of a polynomial |
---|
36 | reverse(l) Reverses a list |
---|
37 | nrroots(p) The number of real roots of p |
---|
38 | isparam(p) Returns 0 iff the polynomial has non-parametric coefficients |
---|
39 | |
---|
40 | KEYWORDS: real roots, univariate polynomial |
---|
41 | "; |
---|
42 | /////////////////////////////////////////////////////////////////////////////// |
---|
43 | |
---|
44 | static proc isparametric(poly p) |
---|
45 | { |
---|
46 | int ispar; |
---|
47 | def ba = basering; |
---|
48 | |
---|
49 | // If the basering has parameters declared |
---|
50 | if (npars(basering) != 0) { |
---|
51 | // If we were given just a polynomial |
---|
52 | list lba = ringlist(ba); |
---|
53 | lba[1]=0; |
---|
54 | def rba = ring(lba); setring rba; |
---|
55 | poly p1 = imap(ba,p); |
---|
56 | setring ba; |
---|
57 | poly p1 = imap(rba,p1); |
---|
58 | ispar = ispar || (size(p-p1)!=0); |
---|
59 | } |
---|
60 | return (ispar); |
---|
61 | } |
---|
62 | /////////////////////////////////////////////////////////////////////////////// |
---|
63 | proc isparam(list #) |
---|
64 | "USAGE: isparam(ideal/module/poly/list); |
---|
65 | RETURN: int: 0 if the argument is not parametric and 1 if it |
---|
66 | EXAMPLE: example isparam; shows an example" |
---|
67 | { |
---|
68 | int i; |
---|
69 | int ispar; |
---|
70 | |
---|
71 | def ar = #[1]; |
---|
72 | |
---|
73 | ispar = 0; |
---|
74 | |
---|
75 | // It we were given only one argument (not a list) |
---|
76 | if (size(#) == 1) { |
---|
77 | if (typeof(ar) == "number") { |
---|
78 | ispar = (pardeg(ar) > 0); |
---|
79 | } else { |
---|
80 | if (typeof(ar) == "poly") { |
---|
81 | ispar = isparametric(ar); |
---|
82 | } else { |
---|
83 | if (typeof(ar) == "ideal" || typeof(ar) == "module") { |
---|
84 | // Ciclo que revisa cada polinomio |
---|
85 | i = size(ar); |
---|
86 | while (!ispar && (i >= 1)) { |
---|
87 | ispar = ispar || (isparametric(ar[i])); |
---|
88 | i--; |
---|
89 | } |
---|
90 | } else { |
---|
91 | if (typeof(ar) == "matrix" || typeof(ar) == "intmat") { |
---|
92 | int j; |
---|
93 | i = nrows(m); |
---|
94 | while (!ispar && (i >= 1)) { |
---|
95 | j = nrows(m); |
---|
96 | while (!ispar && (j >= 1)) { |
---|
97 | ispar = ispar || (isparametric(m[i,j])); |
---|
98 | j--; |
---|
99 | } |
---|
100 | i--; |
---|
101 | } |
---|
102 | } |
---|
103 | }}}} else { |
---|
104 | if (size(#) > 1) { |
---|
105 | i = size(#); |
---|
106 | while (!ispar && (i >= 1)) { |
---|
107 | if ((typeof(#[i]) != "poly") && (typeof(#[i]) != "number") && |
---|
108 | typeof(#[i]) != "int") { |
---|
109 | ERROR("This procedure only works with lists of polynomials"); |
---|
110 | } |
---|
111 | ispar = ispar || (isparametric(#[i])); |
---|
112 | i--; |
---|
113 | } |
---|
114 | }} |
---|
115 | return (ispar); |
---|
116 | } |
---|
117 | example |
---|
118 | { |
---|
119 | echo = 2; |
---|
120 | ring r = 0,x,dp; |
---|
121 | isparam(2x3-56x+2); |
---|
122 | ring s = (0,a,b,c),x,dp; |
---|
123 | isparam(2x3-56x+2); |
---|
124 | isparam(2x3-56x+abc); |
---|
125 | } |
---|
126 | /////////////////////////////////////////////////////////////////////////////// |
---|
127 | proc isuni(poly p) |
---|
128 | "USAGE: isuni(p); poly p; |
---|
129 | RETURN: poly: if p is a univariate polynomial, it returns the variable. If |
---|
130 | not, zero. |
---|
131 | SEE ALSO: whichvariable |
---|
132 | EXAMPLE: example isuni; shows an example" |
---|
133 | { |
---|
134 | poly variable; |
---|
135 | poly g; |
---|
136 | int i; |
---|
137 | int sofarsogood; |
---|
138 | |
---|
139 | variable = whichvariable(leadmonom(p)); |
---|
140 | |
---|
141 | if (variable != 0) { |
---|
142 | i = 1; |
---|
143 | sofarsogood = 1; |
---|
144 | while (sofarsogood && i <= size(p)) { |
---|
145 | g = p[i]/(variable^(deg(p[i]))); |
---|
146 | sofarsogood = sofarsogood && deg(g) == 0; |
---|
147 | i++; |
---|
148 | } |
---|
149 | } |
---|
150 | |
---|
151 | if (sofarsogood) { |
---|
152 | return (variable); |
---|
153 | } else { |
---|
154 | return (0); |
---|
155 | } |
---|
156 | } |
---|
157 | example |
---|
158 | { |
---|
159 | echo = 2; |
---|
160 | ring r = 0,(x,y),dp; |
---|
161 | poly p = 6x7-3x2+2x-15/7; |
---|
162 | isuni(p); |
---|
163 | isuni(p*y); |
---|
164 | } |
---|
165 | /////////////////////////////////////////////////////////////////////////////// |
---|
166 | proc whichvariable(poly p) |
---|
167 | "USAGE: whichvariable(p); poly p |
---|
168 | RETURN: poly: if p is a univariate monomial, the variable. Otherwise 0. |
---|
169 | ASSUME: p is a monomial |
---|
170 | SEE ALSO: isuni |
---|
171 | EXAMPLE: example whichvariable; shows an example" |
---|
172 | { |
---|
173 | int i; |
---|
174 | int found; |
---|
175 | poly g; |
---|
176 | poly variable; |
---|
177 | |
---|
178 | if (size(p) != 1) { |
---|
179 | ERROR("p must be a monomial"); |
---|
180 | } |
---|
181 | |
---|
182 | found = 0; |
---|
183 | |
---|
184 | i = 1; |
---|
185 | while (!found && i <= nvars(basering)) { |
---|
186 | g = diff(p,var(i)); |
---|
187 | if (deg(g) == deg(p) - 1) { |
---|
188 | // We suspect p is a polynomial in var(i) |
---|
189 | variable = var(i); |
---|
190 | found = 1; |
---|
191 | } |
---|
192 | i++; |
---|
193 | } |
---|
194 | |
---|
195 | // We now have to confirm our suspicion |
---|
196 | if (found) { |
---|
197 | if (leadmonom(p/(variable^deg(p))) == 1) { //If p is a monomial on variable |
---|
198 | return (variable); |
---|
199 | } else { if (deg(p) <= 0) { // If we received a constant polynomial |
---|
200 | return (var(1)); |
---|
201 | } else { |
---|
202 | return (0); |
---|
203 | }} |
---|
204 | } else { if (deg(p) <= 0) { // If we received a constant polynomial |
---|
205 | return (var(1)); |
---|
206 | } else { |
---|
207 | return (0); |
---|
208 | }} |
---|
209 | } |
---|
210 | example |
---|
211 | { |
---|
212 | echo = 2; |
---|
213 | ring r = 0,(x,y),dp; |
---|
214 | whichvariable(x5); |
---|
215 | whichvariable(x3y); |
---|
216 | } |
---|
217 | /////////////////////////////////////////////////////////////////////////////// |
---|
218 | proc varsigns(list l) |
---|
219 | "USAGE: varsigns(l); list l. |
---|
220 | RETURN: number: the number of sign changes in the list l |
---|
221 | SEE ALSO: boundposDes |
---|
222 | EXAMPLE: example varsigns; shows an example" |
---|
223 | { |
---|
224 | int lastsign; |
---|
225 | int numberofchanges = 0; |
---|
226 | |
---|
227 | if (isparam(l)) { |
---|
228 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
229 | } |
---|
230 | |
---|
231 | lastsign = sign(l[1]); |
---|
232 | |
---|
233 | for (int i = 1; i <= size(l); i = i + 1) { |
---|
234 | if (sign(l[i]) != lastsign && sign(l[i]) != 0) { |
---|
235 | numberofchanges++; |
---|
236 | lastsign = sign(l[i]); |
---|
237 | } |
---|
238 | } |
---|
239 | return (numberofchanges); |
---|
240 | } |
---|
241 | example |
---|
242 | { |
---|
243 | echo = 2; |
---|
244 | ring r = 0,x,dp; |
---|
245 | list l = 1,2,3; |
---|
246 | varsigns(l); |
---|
247 | l = 1,-1,2,-2,3,-3; |
---|
248 | varsigns(l); |
---|
249 | } |
---|
250 | /////////////////////////////////////////////////////////////////////////////// |
---|
251 | proc boundBuFou(poly p,number a,number b) |
---|
252 | "USAGE: boundBuFou(p,a,b); p poly, a,b number |
---|
253 | RETURN: number: an upper bound for the number of real roots of p in (a,b], |
---|
254 | with the same parity as the actual number of roots (using the |
---|
255 | Budan-Fourier Theorem) |
---|
256 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
257 | a, b are rational numbers with a < b |
---|
258 | SEE ALSO: boundposDes,varsigns |
---|
259 | EXAMPLE: example boundBuFou; shows an example" |
---|
260 | { |
---|
261 | int i; |
---|
262 | poly variable; |
---|
263 | list Der; |
---|
264 | list Dera,Derb; |
---|
265 | int d; |
---|
266 | number bound; |
---|
267 | |
---|
268 | variable = isuni(p); |
---|
269 | |
---|
270 | if (isparam(p) || isparam(a) || isparam(b)) { |
---|
271 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
272 | } |
---|
273 | |
---|
274 | // p must be a univariate polynomial |
---|
275 | if (variable == 0) { |
---|
276 | ERROR("p must be a univariate polynomial"); |
---|
277 | } |
---|
278 | |
---|
279 | if (a >= b) { |
---|
280 | ERROR("a must be smaller than b"); |
---|
281 | } |
---|
282 | |
---|
283 | d = deg(p); |
---|
284 | |
---|
285 | // We calculate the list of derivatives |
---|
286 | |
---|
287 | Der[d+1] = p; |
---|
288 | |
---|
289 | for (i = 0;i < d;i++) { |
---|
290 | Der[d-i] = diff(Der[d-i+1],variable); |
---|
291 | } |
---|
292 | |
---|
293 | // Then evaluate that list |
---|
294 | |
---|
295 | for (i = d+1;i >= 1;i--) { |
---|
296 | Dera [i] = leadcoef(subst(Der[i],variable,a)); |
---|
297 | Derb [i] = leadcoef(subst(Der[i],variable,b)); |
---|
298 | } |
---|
299 | |
---|
300 | // Finally we calculate the sign variations |
---|
301 | |
---|
302 | bound = varsigns(Dera) - varsigns(Derb); |
---|
303 | |
---|
304 | return(bound); |
---|
305 | } |
---|
306 | example |
---|
307 | { |
---|
308 | echo = 2; |
---|
309 | ring r = 0,x,dp; |
---|
310 | poly p = (x+2)*(x-1)*(x-5); |
---|
311 | boundBuFou(p,-3,5); |
---|
312 | boundBuFou(p,-2,5); |
---|
313 | } |
---|
314 | /////////////////////////////////////////////////////////////////////////////// |
---|
315 | proc boundposDes(poly p) |
---|
316 | "USAGE: boundposDes(p); poly p |
---|
317 | RETURN: number: an upper bound for the number of positive roots of p, with |
---|
318 | the same parity as the actual number of positive roots of p. |
---|
319 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
320 | SEE ALSO: boundBuFou |
---|
321 | EXAMPLE: example boundposDes; shows an example" |
---|
322 | { |
---|
323 | poly g; |
---|
324 | number nroots; |
---|
325 | poly variable; |
---|
326 | list coefficients; |
---|
327 | int i; |
---|
328 | |
---|
329 | variable = isuni(p); |
---|
330 | |
---|
331 | if (isparam(p)) { |
---|
332 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
333 | } |
---|
334 | |
---|
335 | // p must be a univariate polynomial |
---|
336 | if (variable == 0) { |
---|
337 | ERROR("p must be a univariate polynomial"); |
---|
338 | } |
---|
339 | |
---|
340 | g = p; // We will work with g |
---|
341 | |
---|
342 | // We check whether 0 is a root of g, and if so, remove it |
---|
343 | if (subst(g,variable,0) == 0) { |
---|
344 | g = g/variable^(deg(g[size[g]])); |
---|
345 | } |
---|
346 | |
---|
347 | // We count the number of positive roots |
---|
348 | i = size(g); |
---|
349 | while (i >= 1) { |
---|
350 | coefficients[i] = leadcoef(g[i]); |
---|
351 | i--; |
---|
352 | } |
---|
353 | nroots = varsigns(coefficients); |
---|
354 | |
---|
355 | return(nroots); |
---|
356 | } |
---|
357 | example |
---|
358 | { |
---|
359 | echo = 2; |
---|
360 | ring r = 0,x,dp; |
---|
361 | poly p = (x+2)*(x-1)*(x-5); |
---|
362 | boundposDes(p); |
---|
363 | |
---|
364 | p = p*(x2+1); |
---|
365 | |
---|
366 | boundposDes(p); |
---|
367 | } |
---|
368 | /////////////////////////////////////////////////////////////////////////////// |
---|
369 | proc boundDes(poly p) |
---|
370 | "USAGE: boundDes(p); poly p |
---|
371 | RETURN: number: an upper bound for the number of real roots of p, with |
---|
372 | the same parity as the actual number of real roots of p. |
---|
373 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
374 | SEE ALSO: boundBuFou |
---|
375 | EXAMPLE: example boundDes; shows an example" |
---|
376 | { |
---|
377 | poly g; |
---|
378 | number nroots; |
---|
379 | poly variable; |
---|
380 | list coefficients; |
---|
381 | int i; |
---|
382 | |
---|
383 | variable = isuni(p); |
---|
384 | |
---|
385 | if (isparam(p)) { |
---|
386 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
387 | } |
---|
388 | |
---|
389 | // p must be a univariate polynomial |
---|
390 | if (variable == 0) { |
---|
391 | ERROR("p must be a univariate polynomial"); |
---|
392 | } |
---|
393 | |
---|
394 | g = p; // We will work with g |
---|
395 | |
---|
396 | nroots = 0; |
---|
397 | // We check whether 0 is a root of g, and if so, remove it |
---|
398 | if (subst(g,variable,0) == 0) { |
---|
399 | g = g/variable^(deg(g[size[g]])); |
---|
400 | nroots++; |
---|
401 | } |
---|
402 | |
---|
403 | // We count the number of positive roots |
---|
404 | i = size(g); |
---|
405 | while (i >= 1) { |
---|
406 | coefficients[i] = leadcoef(g[i]); |
---|
407 | i--; |
---|
408 | } |
---|
409 | nroots = nroots + varsigns(coefficients); |
---|
410 | |
---|
411 | // We count the number of negative roots |
---|
412 | g = subst(g,variable,-variable); |
---|
413 | i = size(g); |
---|
414 | while (i >= 1) { |
---|
415 | coefficients[i] = leadcoef(g[i]); |
---|
416 | i--; |
---|
417 | } |
---|
418 | nroots = nroots + varsigns(coefficients); |
---|
419 | |
---|
420 | return(nroots); |
---|
421 | } |
---|
422 | example |
---|
423 | { |
---|
424 | echo = 2; |
---|
425 | ring r = 0,x,dp; |
---|
426 | poly p = (x+2)*(x-1)*(x-5); |
---|
427 | boundDes(p); |
---|
428 | |
---|
429 | p = p*(x2+1); |
---|
430 | |
---|
431 | boundDes(p); |
---|
432 | } |
---|
433 | /////////////////////////////////////////////////////////////////////////////// |
---|
434 | proc allrealst(poly p) |
---|
435 | "USAGE: allrealst(p); poly p |
---|
436 | RETURN: int: 1 iff all the roots of p are real, 0 otherwise |
---|
437 | Checks whether all the roots of p are real by using Sturm's Theorem |
---|
438 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
439 | SEE ALSO: allreal,sturm,sturmha |
---|
440 | EXAMPLE: example allrealst; shows an example" |
---|
441 | { |
---|
442 | number upper,lower; |
---|
443 | poly sqfp; // The square-free part of p |
---|
444 | poly variable; |
---|
445 | |
---|
446 | variable = isuni(p); |
---|
447 | |
---|
448 | if (isparam(p)) { |
---|
449 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
450 | } |
---|
451 | if (variable == 0) { |
---|
452 | ERROR ("p must be a univariate polynomial"); |
---|
453 | } |
---|
454 | |
---|
455 | sqfp = p/gcd(p,diff(p,variable)); |
---|
456 | |
---|
457 | upper = maxabs(sqfp); // By adding one we ensure that sqfp(upper) != 0 |
---|
458 | lower = -upper; |
---|
459 | |
---|
460 | return (sturm(sqfp,lower,upper) == deg(sqfp)); |
---|
461 | } |
---|
462 | example |
---|
463 | { |
---|
464 | echo = 2; |
---|
465 | ring r = 0,x,dp; |
---|
466 | poly p = (x+2)*(x-1)*(x-5); |
---|
467 | allrealst(p); |
---|
468 | p = p*(x2+1); |
---|
469 | allrealst(p); |
---|
470 | } |
---|
471 | /////////////////////////////////////////////////////////////////////////////// |
---|
472 | proc maxabs(poly p) |
---|
473 | "USAGE: maxabs(p); poly p |
---|
474 | RETURN: number: an upper bound for the largest absolute value of a root of p |
---|
475 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
476 | SEE ALSO: sturm |
---|
477 | EXAMPLE: example maxabs; shows an example" |
---|
478 | { |
---|
479 | number maximum; |
---|
480 | poly monic; |
---|
481 | int i; |
---|
482 | |
---|
483 | if (isparam(p)) { |
---|
484 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
485 | } |
---|
486 | |
---|
487 | monic = simplify(p,1); |
---|
488 | |
---|
489 | maximum = 0; |
---|
490 | |
---|
491 | for (i = 1; i <= size(monic); i = i + 1) { |
---|
492 | maximum = max(abs(leadcoef(p[i])),maximum); |
---|
493 | } |
---|
494 | |
---|
495 | return (maximum + 1); |
---|
496 | } |
---|
497 | example |
---|
498 | { |
---|
499 | echo = 2; |
---|
500 | echo = 2; |
---|
501 | ring r = 0,x,dp; |
---|
502 | poly p = (x+2)*(x-1)*(x-5); |
---|
503 | maxabs(p); |
---|
504 | } |
---|
505 | /////////////////////////////////////////////////////////////////////////////// |
---|
506 | proc sturm(poly p,number a,number b) |
---|
507 | "USAGE: sturm(p,a,b); poly p, number a,b |
---|
508 | RETURN: number: the number of real roots of p in (a,b] |
---|
509 | ASSUME: p is a univarite polynomials with rational coefficients, |
---|
510 | a, b are rational numbers with a < b |
---|
511 | SEE ALSO: sturmha,allrealst,allreal |
---|
512 | EXAMPLE: example sturm; shows an example" |
---|
513 | { |
---|
514 | list l; |
---|
515 | list pa; |
---|
516 | list pb; |
---|
517 | int signsA,signsB; |
---|
518 | int i; |
---|
519 | int nroots; |
---|
520 | poly variable; |
---|
521 | |
---|
522 | if (isparam(p)) { |
---|
523 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
524 | } |
---|
525 | |
---|
526 | variable = isuni(p); |
---|
527 | |
---|
528 | if (variable == 0) { |
---|
529 | ERROR ("p must be a univariate polynomial"); |
---|
530 | } |
---|
531 | |
---|
532 | if (a >= b) { |
---|
533 | ERROR("a must be lower than b"); |
---|
534 | } |
---|
535 | |
---|
536 | if (subst(p,variable,a) == 0 || subst(p,variable,b) == 0) { |
---|
537 | ERROR ("Neither a nor b can be roots of P"); |
---|
538 | } |
---|
539 | |
---|
540 | l = sturmseq(p); |
---|
541 | |
---|
542 | i = size(l); |
---|
543 | |
---|
544 | while (i >= 1) { // We build the sequences |
---|
545 | pa[i] = leadcoef(subst(l[i],variable,a)); |
---|
546 | pb[i] = leadcoef(subst(l[i],variable,b)); |
---|
547 | i--; |
---|
548 | } |
---|
549 | |
---|
550 | signsA = varsigns(pa); |
---|
551 | signsB = varsigns(pb); |
---|
552 | |
---|
553 | nroots = signsA - signsB + nroots; |
---|
554 | |
---|
555 | return (nroots); |
---|
556 | } |
---|
557 | example |
---|
558 | { |
---|
559 | echo = 2; |
---|
560 | ring r = 0,x,dp; |
---|
561 | poly p = (x+2)*(x-1)*(x-5); |
---|
562 | sturm(p,-3,6); |
---|
563 | p = p*(x2+1); |
---|
564 | sturm(p,-3,6); |
---|
565 | p = p*(x+2); |
---|
566 | sturm(p,-3,6); |
---|
567 | } |
---|
568 | /////////////////////////////////////////////////////////////////////////////// |
---|
569 | proc sturmseq(poly p) |
---|
570 | "USAGE: sturmseq(p); p poly |
---|
571 | RETURN: list: a Sturm sequence of p |
---|
572 | ASSUME: p is a univarite polynomials with rational coefficients |
---|
573 | THEORY: The Sturm sequence of p (also called remainder sequence) is the |
---|
574 | sequence begininng with p, p' and goes on with minus the remainder |
---|
575 | of the two previous polynomials, until the remainder is zero. |
---|
576 | See: Basu, Pollack, Roy, Algorithms in Real Algebraic Geometry, |
---|
577 | Springer, 2003. |
---|
578 | SEE ALSO: sturm,sturmhaseq |
---|
579 | EXAMPLE: example sturmseq; shows an example" |
---|
580 | { |
---|
581 | list stseq; |
---|
582 | poly variable; |
---|
583 | int i; |
---|
584 | |
---|
585 | variable = isuni(p); |
---|
586 | |
---|
587 | if (isparam(p)) { |
---|
588 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
589 | } |
---|
590 | |
---|
591 | if (variable == 0) { |
---|
592 | ERROR ("p must be a univariate polynomial"); |
---|
593 | } |
---|
594 | |
---|
595 | // The two first polynomials in Sturm's sequence |
---|
596 | stseq = list(); |
---|
597 | stseq[1] = p; |
---|
598 | stseq[2] = diff(p,variable); |
---|
599 | poly q = stseq[2]; |
---|
600 | i = 3; |
---|
601 | |
---|
602 | while (q <>0) { |
---|
603 | q=reduce(stseq[i-2],std(stseq[i-1])); |
---|
604 | if (q<>0) { |
---|
605 | stseq[i] = q; |
---|
606 | } |
---|
607 | i++; |
---|
608 | } |
---|
609 | |
---|
610 | // Right now, we have gcd(P,P') in stseq[size(stseq)]; |
---|
611 | |
---|
612 | for (i = size(stseq)-1;i >= 1;i--) { |
---|
613 | stseq[i] = stseq[i]/(sign(leadcoef(stseq[size(stseq)]))*stseq[size(stseq)]); |
---|
614 | stseq[i] = stseq[i]/abs(leadcoef(stseq[i])); |
---|
615 | } |
---|
616 | |
---|
617 | // We divide the gcd by itself |
---|
618 | stseq[size(stseq)] = sign(leadcoef(stseq[size(stseq)])); |
---|
619 | |
---|
620 | return (stseq); |
---|
621 | } |
---|
622 | example |
---|
623 | { |
---|
624 | echo = 2; |
---|
625 | ring r = 0,(z,x),dp; |
---|
626 | poly p = x5-3x4+12x3+7x-153; |
---|
627 | sturmseq(p); |
---|
628 | } |
---|
629 | /////////////////////////////////////////////////////////////////////////////// |
---|
630 | proc allreal(poly p) |
---|
631 | "USAGE: allreal(p); |
---|
632 | RETURN: bool: 1 iff all the roots of p are real |
---|
633 | SEE ALSO: allrealst |
---|
634 | EXAMPLE: example allreal; shows an example" |
---|
635 | { |
---|
636 | number upper,lower; |
---|
637 | poly sqfp; // The square-free part of p |
---|
638 | poly variable; |
---|
639 | |
---|
640 | if (isparam(p)) { |
---|
641 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
642 | } |
---|
643 | |
---|
644 | variable = isuni(p); |
---|
645 | |
---|
646 | if (variable == 0) { |
---|
647 | ERROR ("p must be a univariate polynomial"); |
---|
648 | } |
---|
649 | |
---|
650 | sqfp = p/gcd(p,diff(p,variable)); |
---|
651 | |
---|
652 | return (sturmha(sqfp,-maxabs(p),maxabs(p)) == deg(sqfp)); |
---|
653 | } |
---|
654 | example |
---|
655 | { |
---|
656 | echo = 2; |
---|
657 | ring r = 0,x,dp; |
---|
658 | poly p = (x+2)*(x-1)*(x-5); |
---|
659 | allreal(p); |
---|
660 | p = p*(x2+1); |
---|
661 | allreal(p); |
---|
662 | } |
---|
663 | /////////////////////////////////////////////////////////////////////////////// |
---|
664 | proc sturmha(poly P,number a,number b) |
---|
665 | "USAGE: sturmha(p,a,b); poly p, number a,b |
---|
666 | RETURN: number: the number of real roots of p in (a,b) (using a |
---|
667 | Sturm-Habicht sequence) |
---|
668 | SEE ALSO: sturm,allreal |
---|
669 | EXAMPLE: example sturmha; shows an example" |
---|
670 | { |
---|
671 | list seq; |
---|
672 | int i; |
---|
673 | list seqa,seqb; |
---|
674 | poly variable; |
---|
675 | number bound; |
---|
676 | number result; |
---|
677 | |
---|
678 | if (isparam(P) || isparam(a) || isparam(b)) { |
---|
679 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
680 | } |
---|
681 | |
---|
682 | variable = isuni(P); |
---|
683 | |
---|
684 | if (variable == 0) { |
---|
685 | ERROR ("P must be a univariate polynomial"); |
---|
686 | } |
---|
687 | |
---|
688 | if (a >= b) { |
---|
689 | ERROR("a must be lower than b"); |
---|
690 | } |
---|
691 | |
---|
692 | if (subst(P,variable,a) == 0 || subst(P,variable,b) == 0) { |
---|
693 | ERROR ("Neither a nor b can be roots of P"); |
---|
694 | } |
---|
695 | |
---|
696 | seq = sturmhaseq(P); |
---|
697 | |
---|
698 | bound = maxabs(P); |
---|
699 | |
---|
700 | if (a < -bound) { |
---|
701 | a = -bound; |
---|
702 | } |
---|
703 | |
---|
704 | if (b > bound) { |
---|
705 | b = bound; |
---|
706 | } |
---|
707 | |
---|
708 | // if (a == -bound && b == bound) { |
---|
709 | // for (i = size(seq);i >= 1;i--) { |
---|
710 | // seq[i] = leadcoef(seq[i]); |
---|
711 | // } |
---|
712 | // result = D(seq); |
---|
713 | // } else { |
---|
714 | for (i = size(seq);i >= 1;i--) { |
---|
715 | seqa[i] = leadcoef(subst(seq[i],variable,a)); |
---|
716 | seqb[i] = leadcoef(subst(seq[i],variable,b)); |
---|
717 | } |
---|
718 | result = (W(seqa) - W(seqb)); |
---|
719 | // } |
---|
720 | return (result); |
---|
721 | } |
---|
722 | example |
---|
723 | { |
---|
724 | echo = 2; |
---|
725 | ring r = 0,x,dp; |
---|
726 | poly p = (x+2)*(x-1)*(x-5); |
---|
727 | sturmha(p,-3,6); |
---|
728 | p = p*(x2+1); |
---|
729 | sturmha(p,-3,6); |
---|
730 | } |
---|
731 | /////////////////////////////////////////////////////////////////////////////// |
---|
732 | proc sturmhaseq(poly P) |
---|
733 | "USAGE: sturmhaseq(P); P poly. |
---|
734 | RETURN: list: the nonzero polynomials of the Sturm-Habicht sequence of P |
---|
735 | ASSUME: P is a univariate polynomial. |
---|
736 | THEORY: The Sturm-Habicht sequence (also subresultant sequence) is closely |
---|
737 | related to the Sturm sequence, but behaves better with respect to |
---|
738 | the size of the coefficients. It is defined via subresultants. |
---|
739 | See: Basu, Pollack, Roy, Algorithms in Real Algebraic Geometry, |
---|
740 | Springer, 2003. |
---|
741 | SEE ALSO: sturm,sturmseq,sturmha |
---|
742 | EXAMPLE: example sturmhaseq; shows an example" |
---|
743 | { |
---|
744 | poly Q; |
---|
745 | poly variable; |
---|
746 | int p,q,i,j,k,l; |
---|
747 | list SR; |
---|
748 | list sr; |
---|
749 | list srbar; |
---|
750 | list T; |
---|
751 | |
---|
752 | if (isparam(P)) { |
---|
753 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
754 | } |
---|
755 | |
---|
756 | variable = isuni(P); |
---|
757 | |
---|
758 | if (variable == 0) { |
---|
759 | ERROR ("P must be a univariate polynomial"); |
---|
760 | } |
---|
761 | |
---|
762 | p = deg(P); |
---|
763 | Q = diff(P,variable); |
---|
764 | q = deg(Q); |
---|
765 | |
---|
766 | // Initialization |
---|
767 | SR[p+2] = sign(leadcoef(P)^(p-q-1))*P; |
---|
768 | // T[p+2] = SR[p+2]; |
---|
769 | |
---|
770 | srbar[p+2] = sign(leadcoef(P)^(p-q)); |
---|
771 | sr[p+2] = srbar[p+2]; |
---|
772 | |
---|
773 | SR[p-1+2] = sign(leadcoef(P)^(p-q+1))*Q; |
---|
774 | // T[p-1+2] = SR[p-1+2]; |
---|
775 | srbar[p-1+2] = sign(leadcoef(P)^(p-q+1))*leadcoef(Q); |
---|
776 | |
---|
777 | i = p+1; |
---|
778 | j = p; |
---|
779 | |
---|
780 | while (SR[j-1+2] != 0) { |
---|
781 | k = deg(SR[j-1+2]); |
---|
782 | if (k == j-1) { |
---|
783 | sr[j-1+2] = srbar[j-1+2]; |
---|
784 | SR[k-1+2] = -(reduce(sr[j-1+2]^2*SR[i-1+2], |
---|
785 | std(SR[j-1+2])))/(sr[j+2]*srbar[i-1+2]); |
---|
786 | |
---|
787 | // T[k-1+2] = SR[k-1+2]; |
---|
788 | srbar[k-1+2] = leadcoef(SR[k-1+2]); |
---|
789 | } |
---|
790 | if (k < j-1) { |
---|
791 | // Computation of sr[k+2] |
---|
792 | for (l = 1;l <= j-k-1;l++) { |
---|
793 | srbar[j-l-1+2] = ((-1)^l)*(srbar[j-1+2]*srbar[j-l+2])/sr[j+2]; |
---|
794 | } |
---|
795 | sr[k+2] = srbar[k+2]; |
---|
796 | |
---|
797 | // Computation of SR[k-1+2] |
---|
798 | SR[k-1+2] = -reduce(srbar[j-1+2]*sr[k+2]*SR[i-1+2], |
---|
799 | std(SR[j-1+2]))/(sr[j+2]*srbar[i-1+2]); |
---|
800 | |
---|
801 | srbar[k-1+2] = leadcoef(SR[k-1+2]); |
---|
802 | |
---|
803 | SR[k+2] = SR[j-1+2] * ( sr[k+2] / leadcoef(SR[j-1+2])); |
---|
804 | } |
---|
805 | i = j; |
---|
806 | j = k; |
---|
807 | } |
---|
808 | |
---|
809 | // We build a new list, discarding the undefined and zero elements |
---|
810 | // Plus, we reverse the elements |
---|
811 | |
---|
812 | list filtered; |
---|
813 | i = size(SR); |
---|
814 | while (i >= 1) { |
---|
815 | if (typeof(SR[i]) != "none") { |
---|
816 | if (SR[i] != 0) { |
---|
817 | filtered = insert(filtered,SR[i]); |
---|
818 | } |
---|
819 | } |
---|
820 | i--; |
---|
821 | } |
---|
822 | |
---|
823 | return (filtered); |
---|
824 | } |
---|
825 | example |
---|
826 | { |
---|
827 | echo = 2; |
---|
828 | ring r = 0,x,dp; |
---|
829 | poly p = x5-x4+x-3/2; |
---|
830 | list l = sturmhaseq(p); |
---|
831 | l; |
---|
832 | } |
---|
833 | /////////////////////////////////////////////////////////////////////////////// |
---|
834 | proc nrroots(poly p) |
---|
835 | "USAGE: nrroots(p); poly p |
---|
836 | RETURN: number: the number of real roots of p |
---|
837 | SEE ALSO: boundposDes, sturm, sturmha |
---|
838 | EXAMPLE: example nrroots; shows an example" |
---|
839 | { |
---|
840 | if (isparam(p)) { |
---|
841 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
842 | } |
---|
843 | |
---|
844 | number a = maxabs(p); |
---|
845 | |
---|
846 | return (sturmha(p,-a,a)); |
---|
847 | |
---|
848 | } |
---|
849 | example |
---|
850 | { |
---|
851 | echo = 2; |
---|
852 | ring r = 0,x,dp; |
---|
853 | poly p = (x+2)*(x-1)*(x-5); |
---|
854 | nrroots(p); |
---|
855 | p = p*(x2+1); |
---|
856 | nrroots(p) |
---|
857 | } |
---|
858 | /////////////////////////////////////////////////////////////////////////////// |
---|
859 | static proc abs(number x) |
---|
860 | // Returns the absolute value of x |
---|
861 | { |
---|
862 | number av; |
---|
863 | |
---|
864 | if (x >= 0) { |
---|
865 | av = x; |
---|
866 | } else { |
---|
867 | av = -x; |
---|
868 | } |
---|
869 | |
---|
870 | return (av); |
---|
871 | } |
---|
872 | /////////////////////////////////////////////////////////////////////////////// |
---|
873 | proc sign(number x) |
---|
874 | { |
---|
875 | int sgn; |
---|
876 | |
---|
877 | if (isparam(x)) { |
---|
878 | print(x); |
---|
879 | ERROR("This procedure cannot operate with parameters"); |
---|
880 | } |
---|
881 | |
---|
882 | if (x > 0) { |
---|
883 | sgn = 1; |
---|
884 | } else { if (x < 0) { |
---|
885 | sgn = -1; |
---|
886 | } else { |
---|
887 | sgn = 0; |
---|
888 | }} |
---|
889 | |
---|
890 | return (sgn); |
---|
891 | } |
---|
892 | /////////////////////////////////////////////////////////////////////////////// |
---|
893 | static proc max(number a,number b) |
---|
894 | { |
---|
895 | number maximum; |
---|
896 | |
---|
897 | if (isparam(a) || isparam(b)) { |
---|
898 | ERROR("This procedure cannot operate with parameters"); |
---|
899 | } |
---|
900 | |
---|
901 | if (a > b) { |
---|
902 | maximum = a; |
---|
903 | } else { |
---|
904 | maximum = b; |
---|
905 | } |
---|
906 | |
---|
907 | return (maximum); |
---|
908 | } |
---|
909 | /////////////////////////////////////////////////////////////////////////////// |
---|
910 | proc reverse(list l) |
---|
911 | "USAGE: reverse(l); l list |
---|
912 | RETURN: list: l reversed. |
---|
913 | EXAMPLE: example reverse; shows an example" |
---|
914 | { |
---|
915 | int i; |
---|
916 | list result; |
---|
917 | |
---|
918 | for (i = 1;i <= size(l);i++) { |
---|
919 | result = list(l[i]) + result; |
---|
920 | } |
---|
921 | return (result); |
---|
922 | } |
---|
923 | example |
---|
924 | { |
---|
925 | echo = 2; |
---|
926 | ring r = 0,x,dp; |
---|
927 | list l = 1,2,3,4,5; |
---|
928 | list rev = reverse(l); |
---|
929 | l; |
---|
930 | rev; |
---|
931 | } |
---|
932 | /////////////////////////////////////////////////////////////////////////////// |
---|
933 | static proc D(list l) |
---|
934 | { |
---|
935 | int p; |
---|
936 | int q; |
---|
937 | int i; |
---|
938 | int sc; // The modified number of sign changes |
---|
939 | |
---|
940 | if (l[size(l)] == 0) { |
---|
941 | ERROR("l[size(l)] cannot be 0"); |
---|
942 | } |
---|
943 | |
---|
944 | sc = 0; |
---|
945 | |
---|
946 | // We know that l[size(l)]] != 0 |
---|
947 | p = size(l); |
---|
948 | q = p - 1; |
---|
949 | |
---|
950 | while (searchnot(l,q,-1,0)) { |
---|
951 | q = searchnot(l,q,-1,0); |
---|
952 | if ((p - q) % 2 == 1) { // if p-q is odd |
---|
953 | sc = sc + ((-1)^(((p-q)*(p-q-1)) / 2))*sign(l[p]*l[q]); |
---|
954 | } |
---|
955 | p = q; |
---|
956 | q = p - 1; |
---|
957 | } |
---|
958 | |
---|
959 | return (sc); |
---|
960 | } |
---|
961 | /////////////////////////////////////////////////////////////////////////////// |
---|
962 | static proc search(list l,int from,int dir,number element) |
---|
963 | { |
---|
964 | int i; |
---|
965 | int result; |
---|
966 | i = from; |
---|
967 | |
---|
968 | result = 0; |
---|
969 | |
---|
970 | while (i + dir >= 0 && i + dir <= size(l) + 1 && !result) { |
---|
971 | if (l[i] == element) { |
---|
972 | result = i; |
---|
973 | } |
---|
974 | i = i + dir; |
---|
975 | } |
---|
976 | |
---|
977 | return (result); |
---|
978 | } |
---|
979 | /////////////////////////////////////////////////////////////////////////////// |
---|
980 | static proc searchnot(list l,int from,int dir,number element) |
---|
981 | { |
---|
982 | int i; |
---|
983 | int result; |
---|
984 | i = from; |
---|
985 | |
---|
986 | result = 0; |
---|
987 | |
---|
988 | while (i + dir >= 0 && i + dir <= size(l) + 1 && !result) { |
---|
989 | if (l[i] != element) { |
---|
990 | result = i; |
---|
991 | } |
---|
992 | i = i + dir; |
---|
993 | } |
---|
994 | |
---|
995 | return (result); |
---|
996 | } |
---|
997 | /////////////////////////////////////////////////////////////////////////////// |
---|
998 | static proc W(list l) |
---|
999 | { |
---|
1000 | int i,temp,sc,lastsign,nofzeros,n; |
---|
1001 | |
---|
1002 | n = size(l); |
---|
1003 | sc = 0; |
---|
1004 | nofzeros = 0; |
---|
1005 | i = 1; |
---|
1006 | lastsign = sign(l[i]); |
---|
1007 | |
---|
1008 | i++; |
---|
1009 | |
---|
1010 | while (i <= n) { |
---|
1011 | if (l[i] == 0) { |
---|
1012 | nofzeros++; |
---|
1013 | } else { |
---|
1014 | temp = lastsign * sign(l[i]); |
---|
1015 | |
---|
1016 | if (temp < 0) { |
---|
1017 | sc++; |
---|
1018 | } else { |
---|
1019 | if (nofzeros == 2) { |
---|
1020 | sc = sc + 2; |
---|
1021 | } |
---|
1022 | } |
---|
1023 | nofzeros = 0; |
---|
1024 | lastsign = temp/lastsign; |
---|
1025 | } |
---|
1026 | i++; |
---|
1027 | } |
---|
1028 | return (sc); |
---|
1029 | } |
---|
1030 | /////////////////////////////////////////////////////////////////////////////// |
---|
1031 | |
---|
1032 | |
---|