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 = i + 1) { |
---|
182 | if (sign(l[i]) != lastsign && sign(l[i]) != 0) { |
---|
183 | numberofchanges++; |
---|
184 | lastsign = sign(l[i]); |
---|
185 | } |
---|
186 | } |
---|
187 | return (numberofchanges); |
---|
188 | } |
---|
189 | example |
---|
190 | { |
---|
191 | echo = 2; |
---|
192 | ring r = 0,x,dp; |
---|
193 | list l = 1,2,3; |
---|
194 | varsigns(l); |
---|
195 | l = 1,-1,2,-2,3,-3; |
---|
196 | varsigns(l); |
---|
197 | } |
---|
198 | /////////////////////////////////////////////////////////////////////////////// |
---|
199 | proc boundBuFou(poly p,number a,number b) |
---|
200 | "USAGE: boundBuFou(p,a,b); p poly, a,b number |
---|
201 | RETURN: int: an upper bound for the number of real roots of p in (a,b], |
---|
202 | with the same parity as the actual number of roots (using the |
---|
203 | Budan-Fourier Theorem) |
---|
204 | ASSUME: - p is a univariate polynomial with rational coefficients@* |
---|
205 | - a, b are rational numbers with a < b |
---|
206 | SEE ALSO: boundposDes,varsigns |
---|
207 | EXAMPLE: example boundBuFou; shows an example" |
---|
208 | { |
---|
209 | int i; |
---|
210 | poly variable; |
---|
211 | list Der; |
---|
212 | list Dera,Derb; |
---|
213 | int d; |
---|
214 | number bound; |
---|
215 | |
---|
216 | variable = isuni(p); |
---|
217 | |
---|
218 | if (isparam(p) || isparam(a) || isparam(b)) { |
---|
219 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
220 | } |
---|
221 | |
---|
222 | // p must be a univariate polynomial |
---|
223 | if (variable == 0) { |
---|
224 | ERROR("p must be a univariate polynomial"); |
---|
225 | } |
---|
226 | |
---|
227 | if (a >= b) { |
---|
228 | ERROR("a must be smaller than b"); |
---|
229 | } |
---|
230 | |
---|
231 | d = deg(p); |
---|
232 | |
---|
233 | // We calculate the list of derivatives |
---|
234 | |
---|
235 | Der[d+1] = p; |
---|
236 | |
---|
237 | for (i = 0;i < d;i++) { |
---|
238 | Der[d-i] = diff(Der[d-i+1],variable); |
---|
239 | } |
---|
240 | |
---|
241 | // Then evaluate that list |
---|
242 | |
---|
243 | for (i = d+1;i >= 1;i--) { |
---|
244 | Dera [i] = leadcoef(subst(Der[i],variable,a)); |
---|
245 | Derb [i] = leadcoef(subst(Der[i],variable,b)); |
---|
246 | } |
---|
247 | |
---|
248 | // Finally we calculate the sign variations |
---|
249 | |
---|
250 | bound = varsigns(Dera) - varsigns(Derb); |
---|
251 | |
---|
252 | return(bound); |
---|
253 | } |
---|
254 | example |
---|
255 | { |
---|
256 | echo = 2; |
---|
257 | ring r = 0,x,dp; |
---|
258 | poly p = (x+2)*(x-1)*(x-5); |
---|
259 | boundBuFou(p,-3,5); |
---|
260 | boundBuFou(p,-2,5); |
---|
261 | } |
---|
262 | /////////////////////////////////////////////////////////////////////////////// |
---|
263 | proc boundposDes(poly p) |
---|
264 | "USAGE: boundposDes(p); poly p |
---|
265 | RETURN: int: an upper bound for the number of positive roots of p, with |
---|
266 | the same parity as the actual number of positive roots of p. |
---|
267 | ASSUME: p is a univariate polynomial with rational coefficients |
---|
268 | SEE ALSO: boundBuFou |
---|
269 | EXAMPLE: example boundposDes; shows an example" |
---|
270 | { |
---|
271 | poly g; |
---|
272 | number nroots; |
---|
273 | poly variable; |
---|
274 | list coefficients; |
---|
275 | int i; |
---|
276 | |
---|
277 | variable = isuni(p); |
---|
278 | |
---|
279 | if (isparam(p)) { |
---|
280 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
281 | } |
---|
282 | |
---|
283 | // p must be a univariate polynomial |
---|
284 | if (variable == 0) { |
---|
285 | ERROR("p must be a univariate polynomial"); |
---|
286 | } |
---|
287 | |
---|
288 | g = p; // We will work with g |
---|
289 | |
---|
290 | // We check whether 0 is a root of g, and if so, remove it |
---|
291 | if (subst(g,variable,0) == 0) { |
---|
292 | g = g/variable^(deg(g[size[g]])); |
---|
293 | } |
---|
294 | |
---|
295 | // We count the number of positive roots |
---|
296 | i = size(g); |
---|
297 | while (i >= 1) { |
---|
298 | coefficients[i] = leadcoef(g[i]); |
---|
299 | i--; |
---|
300 | } |
---|
301 | nroots = varsigns(coefficients); |
---|
302 | |
---|
303 | return(nroots); |
---|
304 | } |
---|
305 | example |
---|
306 | { |
---|
307 | echo = 2; |
---|
308 | ring r = 0,x,dp; |
---|
309 | poly p = (x+2)*(x-1)*(x-5); |
---|
310 | boundposDes(p); |
---|
311 | |
---|
312 | p = p*(x2+1); |
---|
313 | |
---|
314 | boundposDes(p); |
---|
315 | } |
---|
316 | /////////////////////////////////////////////////////////////////////////////// |
---|
317 | proc boundDes(poly p) |
---|
318 | "USAGE: boundDes(p); poly p |
---|
319 | RETURN: int: an upper bound for the number of real roots of p, with |
---|
320 | the same parity as the actual number of real roots of p. |
---|
321 | ASSUME: p is a univariate polynomial with rational coefficients |
---|
322 | SEE ALSO: boundBuFou |
---|
323 | EXAMPLE: example boundDes; shows an example" |
---|
324 | { |
---|
325 | poly g; |
---|
326 | number nroots; |
---|
327 | poly variable; |
---|
328 | list coefficients; |
---|
329 | int i; |
---|
330 | |
---|
331 | variable = isuni(p); |
---|
332 | |
---|
333 | if (isparam(p)) { |
---|
334 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
335 | } |
---|
336 | |
---|
337 | // p must be a univariate polynomial |
---|
338 | if (variable == 0) { |
---|
339 | ERROR("p must be a univariate polynomial"); |
---|
340 | } |
---|
341 | |
---|
342 | g = p; // We will work with g |
---|
343 | |
---|
344 | nroots = 0; |
---|
345 | // We check whether 0 is a root of g, and if so, remove it |
---|
346 | if (subst(g,variable,0) == 0) { |
---|
347 | g = g/variable^(deg(g[size[g]])); |
---|
348 | nroots++; |
---|
349 | } |
---|
350 | |
---|
351 | // We count the number of positive roots |
---|
352 | i = size(g); |
---|
353 | while (i >= 1) { |
---|
354 | coefficients[i] = leadcoef(g[i]); |
---|
355 | i--; |
---|
356 | } |
---|
357 | nroots = nroots + varsigns(coefficients); |
---|
358 | |
---|
359 | // We count the number of negative roots |
---|
360 | g = subst(g,variable,-variable); |
---|
361 | i = size(g); |
---|
362 | while (i >= 1) { |
---|
363 | coefficients[i] = leadcoef(g[i]); |
---|
364 | i--; |
---|
365 | } |
---|
366 | nroots = nroots + varsigns(coefficients); |
---|
367 | |
---|
368 | return(nroots); |
---|
369 | } |
---|
370 | example |
---|
371 | { |
---|
372 | echo = 2; |
---|
373 | ring r = 0,x,dp; |
---|
374 | poly p = (x+2)*(x-1)*(x-5); |
---|
375 | boundDes(p); |
---|
376 | |
---|
377 | p = p*(x2+1); |
---|
378 | |
---|
379 | boundDes(p); |
---|
380 | } |
---|
381 | /////////////////////////////////////////////////////////////////////////////// |
---|
382 | proc allrealst(poly p) |
---|
383 | "USAGE: allrealst(p); poly p |
---|
384 | RETURN: int: 1 if and only if all the roots of p are real, 0 otherwise. |
---|
385 | Checks by using Sturm's Theorem whether all the roots of p are real |
---|
386 | ASSUME: p is a univariate polynomial with rational coefficients |
---|
387 | SEE ALSO: allreal,sturm,sturmha |
---|
388 | EXAMPLE: example allrealst; shows an example" |
---|
389 | { |
---|
390 | number upper,lower; |
---|
391 | poly sqfp; // The square-free part of p |
---|
392 | poly variable; |
---|
393 | |
---|
394 | variable = isuni(p); |
---|
395 | |
---|
396 | if (isparam(p)) { |
---|
397 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
398 | } |
---|
399 | if (variable == 0) { |
---|
400 | ERROR ("p must be a univariate polynomial"); |
---|
401 | } |
---|
402 | |
---|
403 | sqfp = p/gcd(p,diff(p,variable)); |
---|
404 | |
---|
405 | upper = maxabs(sqfp); // By adding one we ensure that sqfp(upper) != 0 |
---|
406 | lower = -upper; |
---|
407 | |
---|
408 | return (sturm(sqfp,lower,upper) == deg(sqfp)); |
---|
409 | } |
---|
410 | example |
---|
411 | { |
---|
412 | echo = 2; |
---|
413 | ring r = 0,x,dp; |
---|
414 | poly p = (x+2)*(x-1)*(x-5); |
---|
415 | allrealst(p); |
---|
416 | p = p*(x2+1); |
---|
417 | allrealst(p); |
---|
418 | } |
---|
419 | /////////////////////////////////////////////////////////////////////////////// |
---|
420 | proc maxabs(poly p) |
---|
421 | "USAGE: maxabs(p); poly p |
---|
422 | RETURN: number: an upper bound for the largest absolute value of a root of p |
---|
423 | ASSUME: p is a univariate polynomial with rational coefficients |
---|
424 | SEE ALSO: sturm |
---|
425 | EXAMPLE: example maxabs; shows an example" |
---|
426 | { |
---|
427 | number maximum; |
---|
428 | poly monic; |
---|
429 | int i; |
---|
430 | |
---|
431 | if (isparam(p)) { |
---|
432 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
433 | } |
---|
434 | |
---|
435 | monic = simplify(p,1); |
---|
436 | |
---|
437 | maximum = 0; |
---|
438 | |
---|
439 | for (i = 1; i <= size(monic); i = i + 1) { |
---|
440 | maximum = max(abs(leadcoef(p[i])),maximum); |
---|
441 | } |
---|
442 | |
---|
443 | return (maximum + 1); |
---|
444 | } |
---|
445 | example |
---|
446 | { |
---|
447 | echo = 2; |
---|
448 | echo = 2; |
---|
449 | ring r = 0,x,dp; |
---|
450 | poly p = (x+2)*(x-1)*(x-5); |
---|
451 | maxabs(p); |
---|
452 | } |
---|
453 | /////////////////////////////////////////////////////////////////////////////// |
---|
454 | proc sturm(poly p,number a,number b) |
---|
455 | "USAGE: sturm(p,a,b); poly p, number a,b |
---|
456 | RETURN: int: the number of real roots of p in (a,b] |
---|
457 | ASSUME: p is a univariate polynomial with rational coefficients,@* |
---|
458 | a, b are rational numbers with a < b |
---|
459 | SEE ALSO: sturmha,allrealst,allreal |
---|
460 | EXAMPLE: example sturm; shows an example" |
---|
461 | { |
---|
462 | list l; |
---|
463 | list pa; |
---|
464 | list pb; |
---|
465 | int signsA,signsB; |
---|
466 | int i; |
---|
467 | int nroots; |
---|
468 | poly variable; |
---|
469 | |
---|
470 | if (isparam(p)) { |
---|
471 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
472 | } |
---|
473 | |
---|
474 | variable = isuni(p); |
---|
475 | |
---|
476 | if (variable == 0) { |
---|
477 | ERROR ("p must be a univariate polynomial"); |
---|
478 | } |
---|
479 | |
---|
480 | if (a >= b) { |
---|
481 | ERROR("a must be lower than b"); |
---|
482 | } |
---|
483 | |
---|
484 | if (subst(p,variable,a) == 0 || subst(p,variable,b) == 0) { |
---|
485 | ERROR ("Neither a nor b can be roots of P"); |
---|
486 | } |
---|
487 | |
---|
488 | l = sturmseq(p); |
---|
489 | |
---|
490 | i = size(l); |
---|
491 | |
---|
492 | while (i >= 1) { // We build the sequences |
---|
493 | pa[i] = leadcoef(subst(l[i],variable,a)); |
---|
494 | pb[i] = leadcoef(subst(l[i],variable,b)); |
---|
495 | i--; |
---|
496 | } |
---|
497 | |
---|
498 | signsA = varsigns(pa); |
---|
499 | signsB = varsigns(pb); |
---|
500 | |
---|
501 | nroots = signsA - signsB + nroots; |
---|
502 | |
---|
503 | return (nroots); |
---|
504 | } |
---|
505 | example |
---|
506 | { |
---|
507 | echo = 2; |
---|
508 | ring r = 0,x,dp; |
---|
509 | poly p = (x+2)*(x-1)*(x-5); |
---|
510 | sturm(p,-3,6); |
---|
511 | p = p*(x2+1); |
---|
512 | sturm(p,-3,6); |
---|
513 | p = p*(x+2); |
---|
514 | sturm(p,-3,6); |
---|
515 | } |
---|
516 | /////////////////////////////////////////////////////////////////////////////// |
---|
517 | proc sturmseq(poly p) |
---|
518 | "USAGE: sturmseq(p); p poly |
---|
519 | RETURN: list: a Sturm sequence of p |
---|
520 | ASSUME: p is a univariate polynomial with rational coefficients |
---|
521 | THEORY: The Sturm sequence of p (also called remainder sequence) is the |
---|
522 | sequence beginning with p, p' and goes on with the negative part of |
---|
523 | the remainder of the two previous polynomials, until the remainder |
---|
524 | is zero. |
---|
525 | See: Basu, Pollack, Roy, Algorithms in Real Algebraic Geometry, |
---|
526 | Springer, 2003. |
---|
527 | SEE ALSO: sturm,sturmhaseq |
---|
528 | EXAMPLE: example sturmseq; shows an example" |
---|
529 | { |
---|
530 | list stseq; |
---|
531 | poly variable; |
---|
532 | int i; |
---|
533 | |
---|
534 | variable = isuni(p); |
---|
535 | |
---|
536 | if (isparam(p)) { |
---|
537 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
538 | } |
---|
539 | |
---|
540 | if (variable == 0) { |
---|
541 | ERROR ("p must be a univariate polynomial"); |
---|
542 | } |
---|
543 | |
---|
544 | // The two first polynomials in Sturm's sequence |
---|
545 | stseq = list(); |
---|
546 | stseq[1] = p; |
---|
547 | stseq[2] = diff(p,variable); |
---|
548 | |
---|
549 | poly q = -reduce(stseq[1],std(stseq[2])); |
---|
550 | i = 3; |
---|
551 | |
---|
552 | while (q <> 0) { |
---|
553 | stseq[i] = q; |
---|
554 | q = -reduce(stseq[i-1],std(stseq[i])); |
---|
555 | i++; |
---|
556 | } |
---|
557 | |
---|
558 | // Right now, we have gcd(P,P') in stseq[size(stseq)]; |
---|
559 | |
---|
560 | for (i = size(stseq)-1;i >= 1;i--) { |
---|
561 | stseq[i] = stseq[i]/(sign(leadcoef(stseq[size(stseq)]))*stseq[size(stseq)]); |
---|
562 | stseq[i] = stseq[i]/abs(leadcoef(stseq[i])); |
---|
563 | } |
---|
564 | |
---|
565 | // We divide the gcd by itself |
---|
566 | stseq[size(stseq)] = sign(leadcoef(stseq[size(stseq)])); |
---|
567 | |
---|
568 | return (stseq); |
---|
569 | } |
---|
570 | example |
---|
571 | { |
---|
572 | echo = 2; |
---|
573 | ring r = 0,(z,x),dp; |
---|
574 | poly p = x5-3x4+12x3+7x-153; |
---|
575 | sturmseq(p); |
---|
576 | } |
---|
577 | /////////////////////////////////////////////////////////////////////////////// |
---|
578 | proc allreal(poly p) |
---|
579 | "USAGE: allreal(p); |
---|
580 | RETURN: int: 1 if and only if all the roots of p are real, 0 otherwise |
---|
581 | SEE ALSO: allrealst |
---|
582 | EXAMPLE: example allreal; shows an example" |
---|
583 | { |
---|
584 | number upper,lower; |
---|
585 | poly sqfp; // The square-free part of p |
---|
586 | poly variable; |
---|
587 | |
---|
588 | if (isparam(p)) { |
---|
589 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
590 | } |
---|
591 | |
---|
592 | variable = isuni(p); |
---|
593 | |
---|
594 | if (variable == 0) { |
---|
595 | ERROR ("p must be a univariate polynomial"); |
---|
596 | } |
---|
597 | |
---|
598 | sqfp = p/gcd(p,diff(p,variable)); |
---|
599 | |
---|
600 | return (sturmha(sqfp,-maxabs(p),maxabs(p)) == deg(sqfp)); |
---|
601 | } |
---|
602 | example |
---|
603 | { |
---|
604 | echo = 2; |
---|
605 | ring r = 0,x,dp; |
---|
606 | poly p = (x+2)*(x-1)*(x-5); |
---|
607 | allreal(p); |
---|
608 | p = p*(x2+1); |
---|
609 | allreal(p); |
---|
610 | } |
---|
611 | /////////////////////////////////////////////////////////////////////////////// |
---|
612 | proc sturmha(poly P,number a,number b) |
---|
613 | "USAGE: sturmha(p,a,b); poly p, number a,b |
---|
614 | RETURN: int: the number of real roots of p in (a,b) (using a Sturm-Habicht sequence) |
---|
615 | SEE ALSO: sturm,allreal |
---|
616 | EXAMPLE: example sturmha; shows an example" |
---|
617 | { |
---|
618 | list seq; |
---|
619 | int i; |
---|
620 | list seqa,seqb; |
---|
621 | poly variable; |
---|
622 | number bound; |
---|
623 | //number result; |
---|
624 | int result; |
---|
625 | |
---|
626 | if (isparam(P) || isparam(a) || isparam(b)) { |
---|
627 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
628 | } |
---|
629 | |
---|
630 | variable = isuni(P); |
---|
631 | |
---|
632 | if (variable == 0) { |
---|
633 | ERROR ("P must be a univariate polynomial"); |
---|
634 | } |
---|
635 | |
---|
636 | if (a >= b) { |
---|
637 | ERROR("a must be lower than b"); |
---|
638 | } |
---|
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) { |
---|
649 | a = -bound; |
---|
650 | } |
---|
651 | |
---|
652 | if (b > bound) { |
---|
653 | b = bound; |
---|
654 | } |
---|
655 | |
---|
656 | // if (a == -bound && b == bound) { |
---|
657 | // for (i = size(seq);i >= 1;i--) { |
---|
658 | // seq[i] = leadcoef(seq[i]); |
---|
659 | // } |
---|
660 | // result = D(seq); |
---|
661 | // } else { |
---|
662 | for (i = size(seq);i >= 1;i--) { |
---|
663 | seqa[i] = leadcoef(subst(seq[i],variable,a)); |
---|
664 | seqb[i] = leadcoef(subst(seq[i],variable,b)); |
---|
665 | } |
---|
666 | result = (W(seqa) - W(seqb)); |
---|
667 | // } |
---|
668 | return (result); |
---|
669 | } |
---|
670 | example |
---|
671 | { |
---|
672 | echo = 2; |
---|
673 | ring r = 0,x,dp; |
---|
674 | poly p = (x+2)*(x-1)*(x-5); |
---|
675 | sturmha(p,-3,6); |
---|
676 | p = p*(x2+1); |
---|
677 | sturmha(p,-3,6); |
---|
678 | } |
---|
679 | /////////////////////////////////////////////////////////////////////////////// |
---|
680 | proc sturmhaseq(poly P) |
---|
681 | "USAGE: sturmhaseq(P); P poly. |
---|
682 | RETURN: list: the non-zero polynomials of the Sturm-Habicht sequence of P |
---|
683 | ASSUME: P is a univariate polynomial. |
---|
684 | THEORY: The Sturm-Habicht sequence (also subresultant sequence) is closely |
---|
685 | related to the Sturm sequence, but behaves better with respect to |
---|
686 | the size of the coefficients. It is defined via subresultants. |
---|
687 | See: Basu, Pollack, Roy, Algorithms in Real Algebraic Geometry, |
---|
688 | Springer, 2003. |
---|
689 | SEE ALSO: sturm,sturmseq,sturmha |
---|
690 | EXAMPLE: example sturmhaseq; shows an example" |
---|
691 | { |
---|
692 | poly Q; |
---|
693 | poly variable; |
---|
694 | int p,q,i,j,k,l; |
---|
695 | list SR; |
---|
696 | list sr; |
---|
697 | list srbar; |
---|
698 | list T; |
---|
699 | |
---|
700 | if (isparam(P)) { |
---|
701 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
702 | } |
---|
703 | |
---|
704 | variable = isuni(P); |
---|
705 | |
---|
706 | if (variable == 0) { |
---|
707 | ERROR ("P must be a univariate polynomial"); |
---|
708 | } |
---|
709 | |
---|
710 | p = deg(P); |
---|
711 | Q = diff(P,variable); |
---|
712 | q = deg(Q); |
---|
713 | |
---|
714 | // Initialization |
---|
715 | SR[p+2] = sign(leadcoef(P)^(p-q-1))*P; |
---|
716 | // T[p+2] = SR[p+2]; |
---|
717 | |
---|
718 | srbar[p+2] = sign(leadcoef(P)^(p-q)); |
---|
719 | sr[p+2] = srbar[p+2]; |
---|
720 | |
---|
721 | SR[p-1+2] = sign(leadcoef(P)^(p-q+1))*Q; |
---|
722 | // T[p-1+2] = SR[p-1+2]; |
---|
723 | srbar[p-1+2] = sign(leadcoef(P)^(p-q+1))*leadcoef(Q); |
---|
724 | |
---|
725 | i = p+1; |
---|
726 | j = p; |
---|
727 | |
---|
728 | while (SR[j-1+2] != 0) { |
---|
729 | k = deg(SR[j-1+2]); |
---|
730 | if (k == j-1) { |
---|
731 | sr[j-1+2] = srbar[j-1+2]; |
---|
732 | SR[k-1+2] = -(reduce(sr[j-1+2]^2*SR[i-1+2], |
---|
733 | std(SR[j-1+2])))/(sr[j+2]*srbar[i-1+2]); |
---|
734 | |
---|
735 | // T[k-1+2] = SR[k-1+2]; |
---|
736 | srbar[k-1+2] = leadcoef(SR[k-1+2]); |
---|
737 | } |
---|
738 | if (k < j-1) { |
---|
739 | // Computation of sr[k+2] |
---|
740 | for (l = 1;l <= j-k-1;l++) { |
---|
741 | srbar[j-l-1+2] = ((-1)^l)*(srbar[j-1+2]*srbar[j-l+2])/sr[j+2]; |
---|
742 | } |
---|
743 | sr[k+2] = srbar[k+2]; |
---|
744 | |
---|
745 | // Computation of SR[k-1+2] |
---|
746 | SR[k-1+2] = -reduce(srbar[j-1+2]*sr[k+2]*SR[i-1+2], |
---|
747 | std(SR[j-1+2]))/(sr[j+2]*srbar[i-1+2]); |
---|
748 | |
---|
749 | srbar[k-1+2] = leadcoef(SR[k-1+2]); |
---|
750 | |
---|
751 | SR[k+2] = SR[j-1+2] * ( sr[k+2] / leadcoef(SR[j-1+2])); |
---|
752 | } |
---|
753 | i = j; |
---|
754 | j = k; |
---|
755 | } |
---|
756 | |
---|
757 | // We build a new list, discarding the undefined and zero elements |
---|
758 | // Plus, we reverse the elements |
---|
759 | |
---|
760 | list filtered; |
---|
761 | i = size(SR); |
---|
762 | while (i >= 1) { |
---|
763 | if (typeof(SR[i]) != "none") { |
---|
764 | if (SR[i] != 0) { |
---|
765 | filtered = insert(filtered,SR[i]); |
---|
766 | } |
---|
767 | } |
---|
768 | i--; |
---|
769 | } |
---|
770 | |
---|
771 | return (filtered); |
---|
772 | } |
---|
773 | example |
---|
774 | { |
---|
775 | echo = 2; |
---|
776 | ring r = 0,x,dp; |
---|
777 | poly p = x5-x4+x-3/2; |
---|
778 | list l = sturmhaseq(p); |
---|
779 | l; |
---|
780 | } |
---|
781 | /////////////////////////////////////////////////////////////////////////////// |
---|
782 | proc nrroots(poly p) |
---|
783 | "USAGE: nrroots(p); poly p |
---|
784 | RETURN: int: the number of real roots of p |
---|
785 | SEE ALSO: boundposDes, sturm, sturmha |
---|
786 | EXAMPLE: example nrroots; shows an example" |
---|
787 | { |
---|
788 | if (isparam(p)) { |
---|
789 | ERROR("This procedure cannot operate with parametric arguments"); |
---|
790 | } |
---|
791 | |
---|
792 | number a = maxabs(p); |
---|
793 | |
---|
794 | return (sturmha(p,-a,a)); |
---|
795 | |
---|
796 | } |
---|
797 | example |
---|
798 | { |
---|
799 | echo = 2; |
---|
800 | ring r = 0,x,dp; |
---|
801 | poly p = (x+2)*(x-1)*(x-5); |
---|
802 | nrroots(p); |
---|
803 | p = p*(x2+1); |
---|
804 | nrroots(p); |
---|
805 | } |
---|
806 | /////////////////////////////////////////////////////////////////////////////// |
---|
807 | static proc abs(number x) |
---|
808 | // Returns the absolute value of x |
---|
809 | { |
---|
810 | number av; |
---|
811 | |
---|
812 | if (x >= 0) { |
---|
813 | av = x; |
---|
814 | } else { |
---|
815 | av = -x; |
---|
816 | } |
---|
817 | |
---|
818 | return (av); |
---|
819 | } |
---|
820 | /////////////////////////////////////////////////////////////////////////////// |
---|
821 | proc sign(number x) |
---|
822 | { |
---|
823 | int sgn; |
---|
824 | |
---|
825 | if (isparam(x)) { |
---|
826 | print(x); |
---|
827 | ERROR("This procedure cannot operate with parameters"); |
---|
828 | } |
---|
829 | |
---|
830 | if (x > 0) { |
---|
831 | sgn = 1; |
---|
832 | } else { if (x < 0) { |
---|
833 | sgn = -1; |
---|
834 | } else { |
---|
835 | sgn = 0; |
---|
836 | }} |
---|
837 | |
---|
838 | return (sgn); |
---|
839 | } |
---|
840 | /////////////////////////////////////////////////////////////////////////////// |
---|
841 | static proc max(number a,number b) |
---|
842 | { |
---|
843 | number maximum; |
---|
844 | |
---|
845 | if (isparam(a) || isparam(b)) { |
---|
846 | ERROR("This procedure cannot operate with parameters"); |
---|
847 | } |
---|
848 | |
---|
849 | if (a > b) { |
---|
850 | maximum = a; |
---|
851 | } else { |
---|
852 | maximum = b; |
---|
853 | } |
---|
854 | |
---|
855 | return (maximum); |
---|
856 | } |
---|
857 | /////////////////////////////////////////////////////////////////////////////// |
---|
858 | proc reverse(list l) |
---|
859 | "USAGE: reverse(l); l list |
---|
860 | RETURN: list: l reversed. |
---|
861 | EXAMPLE: example reverse; shows an example" |
---|
862 | { |
---|
863 | int i; |
---|
864 | list result; |
---|
865 | |
---|
866 | for (i = 1;i <= size(l);i++) { |
---|
867 | result = list(l[i]) + result; |
---|
868 | } |
---|
869 | return (result); |
---|
870 | } |
---|
871 | example |
---|
872 | { |
---|
873 | echo = 2; |
---|
874 | ring r = 0,x,dp; |
---|
875 | list l = 1,2,3,4,5; |
---|
876 | list rev = reverse(l); |
---|
877 | l; |
---|
878 | rev; |
---|
879 | } |
---|
880 | /////////////////////////////////////////////////////////////////////////////// |
---|
881 | static proc D(list l) |
---|
882 | { |
---|
883 | int p; |
---|
884 | int q; |
---|
885 | int i; |
---|
886 | int sc; // The modified number of sign changes |
---|
887 | |
---|
888 | if (l[size(l)] == 0) { |
---|
889 | ERROR("l[size(l)] cannot be 0"); |
---|
890 | } |
---|
891 | |
---|
892 | sc = 0; |
---|
893 | |
---|
894 | // We know that l[size(l)]] != 0 |
---|
895 | p = size(l); |
---|
896 | q = p - 1; |
---|
897 | |
---|
898 | while (searchnot(l,q,-1,0)) { |
---|
899 | q = searchnot(l,q,-1,0); |
---|
900 | if ((p - q) % 2 == 1) { // if p-q is odd |
---|
901 | sc = sc + ((-1)^(((p-q)*(p-q-1)) / 2))*sign(l[p]*l[q]); |
---|
902 | } |
---|
903 | p = q; |
---|
904 | q = p - 1; |
---|
905 | } |
---|
906 | |
---|
907 | return (sc); |
---|
908 | } |
---|
909 | /////////////////////////////////////////////////////////////////////////////// |
---|
910 | static proc search(list l,int from,int dir,number element) |
---|
911 | { |
---|
912 | int i; |
---|
913 | int result; |
---|
914 | i = from; |
---|
915 | |
---|
916 | result = 0; |
---|
917 | |
---|
918 | while (i + dir >= 0 && i + dir <= size(l) + 1 && !result) { |
---|
919 | if (l[i] == element) { |
---|
920 | result = i; |
---|
921 | } |
---|
922 | i = i + dir; |
---|
923 | } |
---|
924 | |
---|
925 | return (result); |
---|
926 | } |
---|
927 | /////////////////////////////////////////////////////////////////////////////// |
---|
928 | static proc searchnot(list l,int from,int dir,number element) |
---|
929 | { |
---|
930 | int i; |
---|
931 | int result; |
---|
932 | i = from; |
---|
933 | |
---|
934 | result = 0; |
---|
935 | |
---|
936 | while (i + dir >= 0 && i + dir <= size(l) + 1 && !result) { |
---|
937 | if (l[i] != element) { |
---|
938 | result = i; |
---|
939 | } |
---|
940 | i = i + dir; |
---|
941 | } |
---|
942 | |
---|
943 | return (result); |
---|
944 | } |
---|
945 | /////////////////////////////////////////////////////////////////////////////// |
---|
946 | static proc W(list l) |
---|
947 | { |
---|
948 | int i,temp,sc,lastsign,nofzeros,n; |
---|
949 | |
---|
950 | n = size(l); |
---|
951 | sc = 0; |
---|
952 | nofzeros = 0; |
---|
953 | i = 1; |
---|
954 | lastsign = sign(l[i]); |
---|
955 | |
---|
956 | i++; |
---|
957 | |
---|
958 | while (i <= n) { |
---|
959 | if (l[i] == 0) { |
---|
960 | nofzeros++; |
---|
961 | } else { |
---|
962 | temp = lastsign * sign(l[i]); |
---|
963 | |
---|
964 | if (temp < 0) { |
---|
965 | sc++; |
---|
966 | } else { |
---|
967 | if (nofzeros == 2) { |
---|
968 | sc = sc + 2; |
---|
969 | } |
---|
970 | } |
---|
971 | nofzeros = 0; |
---|
972 | lastsign = temp div lastsign; |
---|
973 | } |
---|
974 | i++; |
---|
975 | } |
---|
976 | return (sc); |
---|
977 | } |
---|
978 | /////////////////////////////////////////////////////////////////////////////// |
---|
979 | |
---|
980 | |
---|
981 | |
---|
982 | |
---|