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