# source:git/Singular/LIB/rootsur.lib@894734

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