Index: Singular/LIB/mrrcount.lib
===================================================================
--- Singular/LIB/mrrcount.lib	(revision 58f0e9648d6829b34ac7378d6842b1a441057ae6)
+++ Singular/LIB/mrrcount.lib	(revision d3b3ab061482ee56c8453f942b9806ce6f1ef50d)
@@ -1,28 +1,34 @@
-// E. Tobis  12.Nov.2004
-// last change 16. Apr. 2005 (G.-M. Greuel)
+// $Id: mrrcount.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
+// E. Tobis  12.Nov.2004, April 2004
+// last change 1. May 2005 (G.-M. Greuel)
 ///////////////////////////////////////////////////////////////////////////////
 category="Symbolic-numerical solving"
 info="
-LIBRARY: mrrcount.lib Counting number of real roots of polynomial systems
-AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar 
+LIBRARY: mrrcount.lib Counting the number of real roots of polynomial systems
+AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
 
 OVERVIEW:  Routines for counting the number of real roots of a multivariate
-           polynomial system.
-           References:
+           polynomial system. Two methods are implemented: deterministic
+           computation of the number of roots, via the signature of a certain
+           bilinear form; and a rational univariate projection, using a
+           pseudorandom polynomial. Also includes a command to verify the
+           correctness of the pseudorandom answer.
+           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
+           Geometry\", Springer, 2003.
 
 PROCEDURES:
- symsignature(m)   Signature of the symmetric matrix m
- sturmquery(h,B,I) Sturm query of h on V(I)
- matbil(h,B,I)     Matrix of the bilinear form on R/I associated to h
- matmult(f,B,I)    Matrix of multiplication by f (m_f) on R/I in the basis B
- tracemult(f,B,I)  Trace of m_f (B is an ordered basis of R/I)
- coords(f,B,I)     Coordinates of f in the ordered basis B
- randcharpoly(B,I) Pseudorandom charpoly of univ. projection
- verify(p,B,i)     Verifies the result of randcharpoly
- randlinpoly()     Pseudorandom linear polynomial 
- powersums(f,B,I)  Powersums of the roots of a char polynomial
- symmfunc(S)       Symmetric functions from the powersums S
- univarpoly(l)     Polynomial with coefficients from l
- qbase(i)          Like kbase, but the monomials are ordered
+ symsignature(m)     Signature of the symmetric matrix m
+ sturmquery(h,B,I)   Sturm query of h on V(I)
+ matbil(h,B,I)       Matrix of the bilinear form on R/I associated to h
+ matmult(f,B,I)      Matrix of multiplication by f (m_f) on R/I in the basis B
+ tracemult(f,B,I)    Trace of m_f (B is an ordered basis of R/I)
+ coords(f,B,I)       Coordinates of f in the ordered basis B
+ randcharpoly(B,I,n) Pseudorandom charpoly of univ. projection, n optional
+ verify(p,B,i)       Verifies the result of randcharpoly
+ randlinpoly(n)      Pseudorandom linear polynomial, n optional
+ powersums(f,B,I)    Powersums of the roots of a char polynomial
+ symmfunc(S)         Symmetric functions from the powersums S
+ univarpoly(l)       Polynomial with coefficients from l
+ qbase(i)            Like kbase, but the monomials are ordered
 
 KEYWORDS: real roots, univariate projection
@@ -54,5 +60,5 @@
     for (j = i;j <= nrows(m);j++) {
       if (m[i,j] != m[j,i]) {
-	ERROR ("m must be a symmetric matrix");
+        ERROR ("m must be a symmetric matrix");
       }
     }
@@ -81,5 +87,5 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
 
@@ -92,5 +98,6 @@
 "USAGE:     sturmquery(h,b,i); h poly, b,i ideal
 RETURN:    number: the Sturm query of h in V(i)
-ASSUME:    b is an ordered monomial basis of r/i, r = basering
+ASSUME:    i is a Groebner basis, b is an ordered monomial basis
+           of r/i, r = basering.
 SEE ALSO:  symsignature,matbil
 EXAMPLE:   example sturmquery; shows an example"
@@ -107,5 +114,5 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
 
@@ -113,5 +120,5 @@
 }
 
-static proc mysymmsig(matrix m) 
+static proc mysymmsig(matrix m)
 // returns the signature of a square symmetric matrix m
 {
@@ -166,5 +173,5 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
   poly f = x3-xy+y-13+x4-y2x;
@@ -178,6 +185,6 @@
 proc tracemult(poly f,ideal B,ideal I)
 "USAGE:     tracemult(f,b,i);f poly, b,i ideal
-RETURN:    number: the trace of the multiplication by f (m_f) on r/i, 
-           written in the monomial basis b of r/i, r = basering 
+RETURN:    number: the trace of the multiplication by f (m_f) on r/i,
+           written in the monomial basis b of r/i, r = basering
            (faster than matmult + trace)
 ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i
@@ -192,5 +199,5 @@
 
   g = reduce(f,I);
-  
+
   m = 0;
   for (k = 1;k <= size(B);k++) {
@@ -205,5 +212,5 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
 
@@ -231,5 +238,5 @@
 
   g = reduce(f,I);
-  
+
   for (k = 1;k <= size(B);k++) {
     coordinates = coords(g*(B[k]),B,I); // f*x_k written on the basis B
@@ -245,5 +252,5 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
 
@@ -290,7 +297,7 @@
   ring r = 0,(x,y),dp;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
-  ideal b = qbase(i);
-  
+  i = std(i);
+  ideal b = qbase(i);
+
   coords(x3-xy+y-13+x4-y2x,b,i);
   b;
@@ -298,5 +305,5 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc isSquare(matrix m) 
+static proc isSquare(matrix m)
 // returns 1 iff m is a square matrix
 {
@@ -305,8 +312,10 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-proc randcharpoly(ideal B,ideal I)
-"USAGE:     randcharpoly(b,i); b,i ideal
+proc randcharpoly(ideal B,ideal I,list #)
+"USAGE:     randcharpoly(b,i); randcharpoly(b,i,n); b,i ideal; n int
 RETURN:    poly: the characteristic polynomial of a pseudorandom
-           rational univariate projection having one zero per zero of i
+           rational univariate projection having one zero per zero of i.
+           If n is given, it is the number of digits being used for the
+           pseudorandom coefficients (default: n=2)
 ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
            r = basering
@@ -320,5 +329,9 @@
   poly q;
 
-  generic = randlinpoly();
+  if (size(#) == 1) {
+    generic = randlinpoly(#[1]);
+  } else {
+    generic = randlinpoly();
+  }
 
   p = reduce(generic,I);
@@ -332,6 +345,6 @@
   print("* verify(p,b,i)                                                     *");
   print("*                                                                   *");
-  print("* where p is the polynomial, b is the monomial basis used, and i    *");
-  print("* the Groebner basis of the ideal                                   *");
+  print("* where p is the polynomial I returned, b is the monomial basis     *");
+  print("* used, and i the Groebner basis of the ideal                       *");
   print("*********************************************************************");
 
@@ -343,15 +356,19 @@
   ring r = 0,(x,y,z),dp;
   ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
   poly p = randcharpoly(b,i);
   p;
   nrroots(p); // See nrroots in urrcount.lib
-}
+  p = randcharpoly(b,i,5);
+  p;
+  nrroots(p);
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 
 proc verify(poly p,ideal b,ideal i)
 "USAGE:     verify(p,b,i);p poly, b,i,ideal
-RETURN:    integer: 1 iff the polynomial p splits the points of V(i). 
+RETURN:    integer: 1 iff the polynomial p splits the points of V(i).
            It's used to check the result of randcharpoly
 ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
@@ -376,4 +393,5 @@
   } else {
     print("The choice of random numbers was not useful");
+    print("You might want to try randcharpoly with a larger number of digits");
   }
   return (correct);
@@ -385,5 +403,5 @@
   poly f = x3-xy+y-13+x4-y2x;
   ideal i = x4-y2x,y2-13;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
   poly p = randcharpoly(b,i);
@@ -392,8 +410,9 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-proc randlinpoly() 
-"USAGE:     randlinpoly();
+proc randlinpoly(list #)
+"USAGE:     randlinpoly(); randlinpoly(n); n int
 RETURN:    poly: a polynomial linear in each variable of the ring, with
-           pseudorandom coefficients
+           pseudorandom coefficients. If n is given, it is the number of
+           digits being used for the range of the coefficients (default: n=2)
 SEE ALSO:  randcharpoly;
 EXAMPLE:   example randlinpoly; shows an example"
@@ -401,8 +420,13 @@
   int n,i;
   poly p = 0;
-  
+  int ndigits = 2;
+
+  if (size(#) == 1) {
+    ndigits = #[1];
+  }
+
   n = nvars(basering);
   for (i = 1;i <= n;i++) {
-    p = p + var(i)*random(1,1000000);
+    p = p + var(i)*random(1,10^ndigits);
   }
   return (p);
@@ -413,4 +437,6 @@
   ring r = 0,(x,y,z,w),dp;
   poly p = randlinpoly();
+  p;
+  p = randlinpoly(5);
   p;
 }
@@ -424,5 +450,5 @@
 EXAMPLE:   example symmfunc; shows an example"
 {
-  int N,k; 
+  int N,k;
   list sums;
 
@@ -439,5 +465,5 @@
 
   ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5)
-  i = groebner(i);
+  i = std(i);
 
   ideal b = qbase(i);
@@ -450,5 +476,5 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-proc symmfunc(list S) 
+proc symmfunc(list S)
 // Takes the list of power sums and returns the symmetric functions
 "USAGE:     symmfunc(s); s list
@@ -491,5 +517,5 @@
 proc univarpoly(list l)
 "USAGE:     univarpoly(l); l list
-RETURN:    poly: a polynomial p on the first variable of basering, say x, 
+RETURN:    poly: a polynomial p on the first variable of basering, say x,
            with p = l[1] + l[2]*x + l[3]*x^2 + ...
 EXAMPLE:  example univarpoly; shows an example"
@@ -497,5 +523,5 @@
   poly p;
   int i,n;
-  
+
   n = size(l);
   for (i = 1;i <= n;i++) {
@@ -534,5 +560,5 @@
 
   ideal i = 2x2,-y2,z3;
-  i = groebner(i);
+  i = std(i);
   ideal b = qbase(i);
   b;
Index: Singular/LIB/signcond.lib
===================================================================
--- Singular/LIB/signcond.lib	(revision 58f0e9648d6829b34ac7378d6842b1a441057ae6)
+++ Singular/LIB/signcond.lib	(revision d3b3ab061482ee56c8453f942b9806ce6f1ef50d)
@@ -1,2 +1,3 @@
+// $Id: signcond.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
 // E. Tobis  12.Nov.2004
 // last change 16. Apr. 2005 (G.-M. Greuel)
@@ -5,10 +6,15 @@
 info="
 LIBRARY: signcond.lib Routines for computing realizable sign conditions
-AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar 
+AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
+
+OVERVIEW:  Routines to determine the number of solutions of a multivariate
+           polynomial system which satisfy a given sign configuration.
+           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
+           Geometry\", Springer, 2003.
 
 PROCEDURES:
   signcnd(P,I)   The sign conditions realized the polynomials of P on a V(I)
   psigncnd(P,l)  Pretty prints the output of signcnd (l)
-  firstoct(P)    The number of elements of a V(I) with every coordinate > 0
+  firstoct(I)    The number of elements of V(I) with every coordinate > 0
 
 KEYWORDS: real roots,sign conditions
@@ -22,5 +28,5 @@
 "USAGE:    firstoct(i); i ideal
 RETURN:   number: the number of points of V(i) lying in the first orthant
-ASSUME:   i is a Groebner basis 
+ASSUME:   i is a Groebner basis
 SEE ALSO: signcnd
 EXAMPLE:  example firstoct; shows an example"
@@ -112,5 +118,5 @@
 
   // We initialize the cumulative variables
-  M = matrix(1,1,1); 
+  M = matrix(1,1,1);
   Exponents = list(list());
   Signs = list(list());
@@ -131,5 +137,5 @@
 
 
-    // We determine the list of signs which correspond to a nonzero 
+    // We determine the list of signs which correspond to a nonzero
     // number of roots
     numberOfNonZero = 3;
@@ -146,5 +152,5 @@
       numberOfNonZero--;
     }
-    
+
     if (c[3,1] != 0) {
       signi = signi + list(-1);
@@ -152,6 +158,6 @@
       numberOfNonZero--;
     }
-    
-    // We now determine the little matrix we'll work with, 
+
+    // We now determine the little matrix we'll work with,
     // and the list of exponents
     if (numberOfNonZero == 3) {
@@ -162,12 +168,12 @@
       Mi[1,2] = 1;
       if (c[1,1] != 0 && c[2,1] != 0) { // 0,1
-	Mi[2,1] = 0;
-	Mi[2,2] = 1;
+        Mi[2,1] = 0;
+        Mi[2,2] = 1;
       } else {if (c[1,1] != 0 && c[3,1] != 0) { // 0,-1
-	Mi[2,1] = 0;
-	Mi[2,2] = -1;
+        Mi[2,1] = 0;
+        Mi[2,2] = -1;
       } else { // 1,-1
-	Mi[2,1] = 1;
-	Mi[2,2] = -1;
+        Mi[2,1] = 1;
+        Mi[2,2] = -1;
       }}
       exponentsi = list(0,1);
@@ -176,5 +182,5 @@
       exponentsi = list(0);
     }}}
-    
+
     // We store the Sturm Queries we'll need later
     if (numberOfNonZero == 2) {
@@ -194,6 +200,6 @@
   }
 
-  // At this point, we have the cumulative matrix, 
-  // the vector of exponents and the matching sign conditions. 
+  // At this point, we have the cumulative matrix,
+  // the vector of exponents and the matching sign conditions.
   // We have to solve the big linear system to finish.
 
@@ -209,11 +215,11 @@
     if (j <= size(SQvalues)) {
       if (SQpositions[j] == i) {
-	SQs[i,1] = SQvalues[j];
-	j++;
+        SQs[i,1] = SQvalues[j];
+        j++;
       } else {
       SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
-      } 
+      }
     } else {
-	SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
+        SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
     }
   }
@@ -244,8 +250,8 @@
   echo = 0;
 
-  print("The output of signcnd is a list of two lists. Both lists have the 
+  print("The output of signcnd is a list of two lists. Both lists have the
 same");
   print("length. That length is the number of sign conditions realized by the");
-  print ("polynomials of P on the set V(i). In this example, that number 
+  print ("polynomials of P on the set V(i). In this example, that number
 is");
   print("print(size(l[1]));");
@@ -335,5 +341,5 @@
   int i,j;
   int found;
-  
+
   found = 0;
   i = 1;
@@ -345,11 +351,11 @@
     } else {
       while(j <= size(a)) {
-	found = found && a[j] == b[i][j];
-	j++;
+        found = found && a[j] == b[i][j];
+        j++;
       }
     }
     i++;
   }
-  
+
   if (found) {
     return (i-1);
@@ -375,5 +381,5 @@
   result[la*lb] = 0;
 
-  
+
   for (i = 0;i < lb;i++) {
     for (j = 0;j < la;j++) {
@@ -398,5 +404,5 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate 
+static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate
 {
   int i;
Index: Singular/LIB/urrcount.lib
===================================================================
--- Singular/LIB/urrcount.lib	(revision 58f0e9648d6829b34ac7378d6842b1a441057ae6)
+++ Singular/LIB/urrcount.lib	(revision d3b3ab061482ee56c8453f942b9806ce6f1ef50d)
@@ -1,13 +1,22 @@
-// E. Tobis  12.Nov.2004
-// last change 16. Apr. 2005 (G.-M. Greuel)
+// $Id: urrcount.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
+// E. Tobis  12.Nov.2004, April 2004
+// last change 1. May 2005 (G.-M. Greuel)
 ///////////////////////////////////////////////////////////////////////////////
 category="Symbolic-numerical solving"
 info="
 LIBRARY: urrcount.lib   Counting number of real roots of univariate polynomial
-AUTHOR:                 Enrique A. Tobis, etobis@dc.uba.ar 
+AUTHOR:                 Enrique A. Tobis, etobis@dc.uba.ar
 
 OVERVIEW:  Routines for bounding and counting the number of real roots of a
-           univariate polynomial by using Strum and Sturm-Habicht sequences.
-           References:
+           univariate polynomial, by means of several different methods, namely
+           Descartes' rule of signs, the Budan-Fourier theorem, Sturm sequences
+           and Sturm-Habicht sequences. The first two give bounds on the number
+           of roots. The other two compute the actual number of roots of the
+           polynomial. There are several wrapper functions, to simplify the
+           application of the aforesaid theorems and some functions
+           to determine whether a given polynomial is univariate.
+           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
+           Geometry\", Springer, 2003.
+
 
 PROCEDURES:
@@ -26,5 +35,5 @@
   sturmhaseq(p)    A Sturm-Habicht Sequence of a polynomial
   reverse(l)       Reverses a list
-  nrroots(p)        The number of real roots of p
+  nrroots(p)       The number of real roots of p
   isparam(p)       Returns 0 iff the polynomial has non-parametric coefficients
 
@@ -76,6 +85,6 @@
       i = size(ar);
       while (!ispar && (i >= 1)) {
-	ispar = ispar || (isparametric(ar[i]));
-	i--;
+        ispar = ispar || (isparametric(ar[i]));
+        i--;
       }
     } else {
@@ -83,11 +92,11 @@
       int j;
       i = nrows(m);
-      while (!isparam && (i >= 1)) {
-	j = nrows(m);
-	while (!isparam && (j >= 1)) {
-	  isparam = isparam || (isparametric(m[i,j]));
-	  j--;
-	}
-	i--;
+      while (!ispar && (i >= 1)) {
+        j = nrows(m);
+        while (!ispar && (j >= 1)) {
+          ispar = ispar || (isparametric(m[i,j]));
+          j--;
+        }
+        i--;
       }
     }
@@ -97,6 +106,6 @@
     while (!ispar && (i >= 1)) {
       if ((typeof(#[i]) != "poly") && (typeof(#[i]) != "number") &&
-	  typeof(#[i]) != "int") {
-      	ERROR("This procedure only works with lists of polynomials");
+          typeof(#[i]) != "int") {
+              ERROR("This procedure only works with lists of polynomials");
       }
       ispar = ispar || (isparametric(#[i]));
@@ -161,5 +170,5 @@
 SEE ALSO:  isuni
 EXAMPLE:   example whichvariable; shows an example"
-{ 
+{
   int i;
   int found;
@@ -219,5 +228,5 @@
     ERROR("This procedure cannot operate with parametric arguments");
   }
-  
+
   lastsign = sign(l[1]);
 
@@ -240,8 +249,8 @@
 }
 ///////////////////////////////////////////////////////////////////////////////
-proc boundBuFou(poly p,number a,number b) 
+proc boundBuFou(poly p,number a,number b)
 "USAGE:     boundBuFou(p,a,b); p poly, a,b number
-RETURN:    number: an upper bound for the number of real roots of p in (a,b], 
-           with the same parity as the actual number of roots (using the 
+RETURN:    number: an upper bound for the number of real roots of p in (a,b],
+           with the same parity as the actual number of roots (using the
            Budan-Fourier Theorem)
 ASSUME:    p is a univarite polynomials with rational coefficients
@@ -273,5 +282,5 @@
 
   d = deg(p);
-  
+
   // We calculate the list of derivatives
 
@@ -425,5 +434,5 @@
 proc allrealst(poly p)
 "USAGE:     allrealst(p); poly p
-RETURN:    int: 1 iff all the roots of p are real, 0 otherwise 
+RETURN:    int: 1 iff all the roots of p are real, 0 otherwise
            Checks whether all the roots of p are real by using Sturm's Theorem
 ASSUME:    p is a univarite polynomials with rational coefficients
@@ -495,5 +504,5 @@
 }
 ///////////////////////////////////////////////////////////////////////////////
-proc sturm(poly p,number a,number b) 
+proc sturm(poly p,number a,number b)
 "USAGE:     sturm(p,a,b); poly p, number a,b
 RETURN:    number: the number of real roots of p in (a,b]
@@ -562,4 +571,9 @@
 RETURN:    list: a Sturm sequence of p
 ASSUME:    p is a univarite polynomials with rational coefficients
+THEORY:    The Sturm sequence of p (also called remainder sequence) is the
+           sequence begininng with p, p' and goes on with minus the remainder
+           of the two previous polynomials, until the remainder is zero.
+           See: Basu, Pollack, Roy, "Algorithms in Real Algebraic Geometry",
+           Springer, 2003.
 SEE ALSO:  sturm,sturmhaseq
 EXAMPLE:   example sturmseq; shows an example"
@@ -568,7 +582,7 @@
   poly variable;
   int i;
-  
+
   variable = isuni(p);
-  
+
   if (isparam(p)) {
     ERROR("This procedure cannot operate with parametric arguments");
@@ -578,5 +592,5 @@
     ERROR ("p must be a univariate polynomial");
   }
-  
+
   // The two first polynomials in Sturm's sequence
   stseq = list();
@@ -595,13 +609,13 @@
 
   // Right now, we have gcd(P,P') in stseq[size(stseq)];
- 
+
   for (i = size(stseq)-1;i >= 1;i--) {
     stseq[i] = stseq[i]/(sign(leadcoef(stseq[size(stseq)]))*stseq[size(stseq)]);
     stseq[i] = stseq[i]/abs(leadcoef(stseq[i]));
   }
-  
+
   // We divide the gcd by itself
   stseq[size(stseq)] = sign(leadcoef(stseq[size(stseq)]));
-  
+
   return (stseq);
 }
@@ -650,5 +664,5 @@
 proc sturmha(poly P,number a,number b)
 "USAGE:     sturmha(p,a,b); poly p, number a,b
-RETURN:    number: the number of real roots of p in (a,b) (using a 
+RETURN:    number: the number of real roots of p in (a,b) (using a
            Sturm-Habicht sequence)
 SEE ALSO:  sturm,allreal
@@ -717,6 +731,12 @@
 ///////////////////////////////////////////////////////////////////////////////
 proc sturmhaseq(poly P)
-"USAGE:     sturmhaseq(P); P poly. P must be univariate
+"USAGE:     sturmhaseq(P); P poly.
 RETURN:    list: the nonzero polynomials of the Sturm-Habicht sequence of P
+ASSUME:    P is a univariate polynomial.
+THEORY:    The Sturm-Habicht sequence (also subresultant sequence) is closely
+           related to the Sturm sequence, but behaves better with respect to
+           the size of the coefficients. It is defined via subresultants.
+           See: Basu, Pollack, Roy, "Algorithms in Real Algebraic Geometry",
+           Springer, 2003.
 SEE ALSO:  sturm,sturmseq,sturmha
 EXAMPLE:   example sturmhaseq; shows an example"
@@ -767,9 +787,9 @@
 //      T[k-1+2] = SR[k-1+2];
       srbar[k-1+2] = leadcoef(SR[k-1+2]);
-    } 
+    }
     if (k < j-1) {
       // Computation of sr[k+2]
       for (l = 1;l <= j-k-1;l++) {
-	srbar[j-l-1+2] = ((-1)^l)*(srbar[j-1+2]*srbar[j-l+2])/sr[j+2];
+        srbar[j-l-1+2] = ((-1)^l)*(srbar[j-1+2]*srbar[j-l+2])/sr[j+2];
     }
     sr[k+2] = srbar[k+2];
@@ -795,5 +815,5 @@
     if (typeof(SR[i]) != "none") {
       if (SR[i] != 0) {
-	filtered = insert(filtered,SR[i]);
+        filtered = insert(filtered,SR[i]);
       }
     }
@@ -847,5 +867,5 @@
     av = -x;
   }
- 
+
   return (av);
 }
@@ -993,11 +1013,11 @@
     } else {
       temp = lastsign * sign(l[i]);
-      
+
       if (temp < 0) {
-	sc++;
+        sc++;
       } else {
-	if (nofzeros == 2) {
-	  sc = sc + 2;
-	}
+        if (nofzeros == 2) {
+          sc = sc + 2;
+        }
       }
       nofzeros = 0;
@@ -1009,4 +1029,4 @@
 }
 ///////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////
-
+
+
