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