source:git/Singular/LIB/urrcount.lib@edfd36

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