Changeset d3b3ab in git


Ignore:
Timestamp:
May 2, 2005, 2:24:16 PM (18 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
ebce66fe9c5de71af2446287b25216fc9f8634ac
Parents:
58f0e9648d6829b34ac7378d6842b1a441057ae6
Message:
*gmg: new libs
*hannes: examples fixed


git-svn-id: file:///usr/local/Singular/svn/trunk@7991 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular/LIB
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/mrrcount.lib

    r58f0e9 rd3b3ab  
    1 // E. Tobis  12.Nov.2004
    2 // last change 16. Apr. 2005 (G.-M. Greuel)
     1// $Id: mrrcount.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
     2// E. Tobis  12.Nov.2004, April 2004
     3// last change 1. May 2005 (G.-M. Greuel)
    34///////////////////////////////////////////////////////////////////////////////
    45category="Symbolic-numerical solving"
    56info="
    6 LIBRARY: mrrcount.lib Counting number of real roots of polynomial systems
    7 AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar 
     7LIBRARY: mrrcount.lib Counting the number of real roots of polynomial systems
     8AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
    89
    910OVERVIEW:  Routines for counting the number of real roots of a multivariate
    10            polynomial system.
    11            References:
     11           polynomial system. Two methods are implemented: deterministic
     12           computation of the number of roots, via the signature of a certain
     13           bilinear form; and a rational univariate projection, using a
     14           pseudorandom polynomial. Also includes a command to verify the
     15           correctness of the pseudorandom answer.
     16           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
     17           Geometry\", Springer, 2003.
    1218
    1319PROCEDURES:
    14  symsignature(m)   Signature of the symmetric matrix m
    15  sturmquery(h,B,I) Sturm query of h on V(I)
    16  matbil(h,B,I)     Matrix of the bilinear form on R/I associated to h
    17  matmult(f,B,I)    Matrix of multiplication by f (m_f) on R/I in the basis B
    18  tracemult(f,B,I)  Trace of m_f (B is an ordered basis of R/I)
    19  coords(f,B,I)     Coordinates of f in the ordered basis B
    20  randcharpoly(B,I) Pseudorandom charpoly of univ. projection
    21  verify(p,B,i)     Verifies the result of randcharpoly
    22  randlinpoly()     Pseudorandom linear polynomial
    23  powersums(f,B,I)  Powersums of the roots of a char polynomial
    24  symmfunc(S)       Symmetric functions from the powersums S
    25  univarpoly(l)     Polynomial with coefficients from l
    26  qbase(i)          Like kbase, but the monomials are ordered
     20 symsignature(m)     Signature of the symmetric matrix m
     21 sturmquery(h,B,I)   Sturm query of h on V(I)
     22 matbil(h,B,I)       Matrix of the bilinear form on R/I associated to h
     23 matmult(f,B,I)      Matrix of multiplication by f (m_f) on R/I in the basis B
     24 tracemult(f,B,I)    Trace of m_f (B is an ordered basis of R/I)
     25 coords(f,B,I)       Coordinates of f in the ordered basis B
     26 randcharpoly(B,I,n) Pseudorandom charpoly of univ. projection, n optional
     27 verify(p,B,i)       Verifies the result of randcharpoly
     28 randlinpoly(n)      Pseudorandom linear polynomial, n optional
     29 powersums(f,B,I)    Powersums of the roots of a char polynomial
     30 symmfunc(S)         Symmetric functions from the powersums S
     31 univarpoly(l)       Polynomial with coefficients from l
     32 qbase(i)            Like kbase, but the monomials are ordered
    2733
    2834KEYWORDS: real roots, univariate projection
    … …  
    5460    for (j = i;j <= nrows(m);j++) {
    5561      if (m[i,j] != m[j,i]) {
    56         ERROR ("m must be a symmetric matrix");
     62        ERROR ("m must be a symmetric matrix");
    5763      }
    5864    }
    … …  
    8187  ring r = 0,(x,y),dp;
    8288  ideal i = x4-y2x,y2-13;
    83   i = groebner(i);
     89  i = std(i);
    8490  ideal b = qbase(i);
    8591
    … …  
    9298"USAGE:     sturmquery(h,b,i); h poly, b,i ideal
    9399RETURN:    number: the Sturm query of h in V(i)
    94 ASSUME:    b is an ordered monomial basis of r/i, r = basering
     100ASSUME:    i is a Groebner basis, b is an ordered monomial basis
     101           of r/i, r = basering.
    95102SEE ALSO:  symsignature,matbil
    96103EXAMPLE:   example sturmquery; shows an example"
    … …  
    107114  ring r = 0,(x,y),dp;
    108115  ideal i = x4-y2x,y2-13;
    109   i = groebner(i);
     116  i = std(i);
    110117  ideal b = qbase(i);
    111118
    … …  
    113120}
    114121
    115 static proc mysymmsig(matrix m) 
     122static proc mysymmsig(matrix m)
    116123// returns the signature of a square symmetric matrix m
    117124{
    … …  
    166173  ring r = 0,(x,y),dp;
    167174  ideal i = x4-y2x,y2-13;
    168   i = groebner(i);
     175  i = std(i);
    169176  ideal b = qbase(i);
    170177  poly f = x3-xy+y-13+x4-y2x;
    … …  
    178185proc tracemult(poly f,ideal B,ideal I)
    179186"USAGE:     tracemult(f,b,i);f poly, b,i ideal
    180 RETURN:    number: the trace of the multiplication by f (m_f) on r/i, 
    181            written in the monomial basis b of r/i, r = basering 
     187RETURN:    number: the trace of the multiplication by f (m_f) on r/i,
     188           written in the monomial basis b of r/i, r = basering
    182189           (faster than matmult + trace)
    183190ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i
    … …  
    192199
    193200  g = reduce(f,I);
    194  
     201
    195202  m = 0;
    196203  for (k = 1;k <= size(B);k++) {
    … …  
    205212  ring r = 0,(x,y),dp;
    206213  ideal i = x4-y2x,y2-13;
    207   i = groebner(i);
     214  i = std(i);
    208215  ideal b = qbase(i);
    209216
    … …  
    231238
    232239  g = reduce(f,I);
    233  
     240
    234241  for (k = 1;k <= size(B);k++) {
    235242    coordinates = coords(g*(B[k]),B,I); // f*x_k written on the basis B
    … …  
    245252  ring r = 0,(x,y),dp;
    246253  ideal i = x4-y2x,y2-13;
    247   i = groebner(i);
     254  i = std(i);
    248255  ideal b = qbase(i);
    249256
    … …  
    290297  ring r = 0,(x,y),dp;
    291298  ideal i = x4-y2x,y2-13;
    292   i = groebner(i);
    293   ideal b = qbase(i);
    294  
     299  i = std(i);
     300  ideal b = qbase(i);
     301
    295302  coords(x3-xy+y-13+x4-y2x,b,i);
    296303  b;
    … …  
    298305///////////////////////////////////////////////////////////////////////////////
    299306
    300 static proc isSquare(matrix m) 
     307static proc isSquare(matrix m)
    301308// returns 1 iff m is a square matrix
    302309{
    … …  
    305312///////////////////////////////////////////////////////////////////////////////
    306313
    307 proc randcharpoly(ideal B,ideal I)
    308 "USAGE:     randcharpoly(b,i); b,i ideal
     314proc randcharpoly(ideal B,ideal I,list #)
     315"USAGE:     randcharpoly(b,i); randcharpoly(b,i,n); b,i ideal; n int
    309316RETURN:    poly: the characteristic polynomial of a pseudorandom
    310            rational univariate projection having one zero per zero of i
     317           rational univariate projection having one zero per zero of i.
     318           If n is given, it is the number of digits being used for the
     319           pseudorandom coefficients (default: n=2)
    311320ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
    312321           r = basering
    … …  
    320329  poly q;
    321330
    322   generic = randlinpoly();
     331  if (size(#) == 1) {
     332    generic = randlinpoly(#[1]);
     333  } else {
     334    generic = randlinpoly();
     335  }
    323336
    324337  p = reduce(generic,I);
    … …  
    332345  print("* verify(p,b,i)                                                     *");
    333346  print("*                                                                   *");
    334   print("* where p is the polynomial, b is the monomial basis used, and i    *");
    335   print("* the Groebner basis of the ideal                                   *");
     347  print("* where p is the polynomial I returned, b is the monomial basis     *");
     348  print("* used, and i the Groebner basis of the ideal                       *");
    336349  print("*********************************************************************");
    337350
    … …  
    343356  ring r = 0,(x,y,z),dp;
    344357  ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2;
    345   i = groebner(i);
     358  i = std(i);
    346359  ideal b = qbase(i);
    347360  poly p = randcharpoly(b,i);
    348361  p;
    349362  nrroots(p); // See nrroots in urrcount.lib
    350 }
     363  p = randcharpoly(b,i,5);
     364  p;
     365  nrroots(p);
     366}
     367
    351368///////////////////////////////////////////////////////////////////////////////
    352369
    353370proc verify(poly p,ideal b,ideal i)
    354371"USAGE:     verify(p,b,i);p poly, b,i,ideal
    355 RETURN:    integer: 1 iff the polynomial p splits the points of V(i). 
     372RETURN:    integer: 1 iff the polynomial p splits the points of V(i).
    356373           It's used to check the result of randcharpoly
    357374ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
    … …  
    376393  } else {
    377394    print("The choice of random numbers was not useful");
     395    print("You might want to try randcharpoly with a larger number of digits");
    378396  }
    379397  return (correct);
    … …  
    385403  poly f = x3-xy+y-13+x4-y2x;
    386404  ideal i = x4-y2x,y2-13;
    387   i = groebner(i);
     405  i = std(i);
    388406  ideal b = qbase(i);
    389407  poly p = randcharpoly(b,i);
    … …  
    392410///////////////////////////////////////////////////////////////////////////////
    393411
    394 proc randlinpoly()
    395 "USAGE:     randlinpoly();
     412proc randlinpoly(list #)
     413"USAGE:     randlinpoly(); randlinpoly(n); n int
    396414RETURN:    poly: a polynomial linear in each variable of the ring, with
    397            pseudorandom coefficients
     415           pseudorandom coefficients. If n is given, it is the number of
     416           digits being used for the range of the coefficients (default: n=2)
    398417SEE ALSO:  randcharpoly;
    399418EXAMPLE:   example randlinpoly; shows an example"
    … …  
    401420  int n,i;
    402421  poly p = 0;
    403  
     422  int ndigits = 2;
     423
     424  if (size(#) == 1) {
     425    ndigits = #[1];
     426  }
     427
    404428  n = nvars(basering);
    405429  for (i = 1;i <= n;i++) {
    406     p = p + var(i)*random(1,1000000);
     430    p = p + var(i)*random(1,10^ndigits);
    407431  }
    408432  return (p);
    … …  
    413437  ring r = 0,(x,y,z,w),dp;
    414438  poly p = randlinpoly();
     439  p;
     440  p = randlinpoly(5);
    415441  p;
    416442}
    … …  
    424450EXAMPLE:   example symmfunc; shows an example"
    425451{
    426   int N,k; 
     452  int N,k;
    427453  list sums;
    428454
    … …  
    439465
    440466  ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5)
    441   i = groebner(i);
     467  i = std(i);
    442468
    443469  ideal b = qbase(i);
    … …  
    450476///////////////////////////////////////////////////////////////////////////////
    451477
    452 proc symmfunc(list S) 
     478proc symmfunc(list S)
    453479// Takes the list of power sums and returns the symmetric functions
    454480"USAGE:     symmfunc(s); s list
    … …  
    491517proc univarpoly(list l)
    492518"USAGE:     univarpoly(l); l list
    493 RETURN:    poly: a polynomial p on the first variable of basering, say x, 
     519RETURN:    poly: a polynomial p on the first variable of basering, say x,
    494520           with p = l[1] + l[2]*x + l[3]*x^2 + ...
    495521EXAMPLE:  example univarpoly; shows an example"
    … …  
    497523  poly p;
    498524  int i,n;
    499  
     525
    500526  n = size(l);
    501527  for (i = 1;i <= n;i++) {
    … …  
    534560
    535561  ideal i = 2x2,-y2,z3;
    536   i = groebner(i);
     562  i = std(i);
    537563  ideal b = qbase(i);
    538564  b;

  • Singular/LIB/signcond.lib

    r58f0e9 rd3b3ab  
     1// $Id: signcond.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
    12// E. Tobis  12.Nov.2004
    23// last change 16. Apr. 2005 (G.-M. Greuel)
    … …  
    56info="
    67LIBRARY: signcond.lib Routines for computing realizable sign conditions
    7 AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
     8AUTHOR:               Enrique A. Tobis, etobis@dc.uba.ar
     9
     10OVERVIEW:  Routines to determine the number of solutions of a multivariate
     11           polynomial system which satisfy a given sign configuration.
     12           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
     13           Geometry\", Springer, 2003.
    814
    915PROCEDURES:
    1016  signcnd(P,I)   The sign conditions realized the polynomials of P on a V(I)
    1117  psigncnd(P,l)  Pretty prints the output of signcnd (l)
    12   firstoct(P)    The number of elements of a V(I) with every coordinate > 0
     18  firstoct(I)    The number of elements of V(I) with every coordinate > 0
    1319
    1420KEYWORDS: real roots,sign conditions
    … …  
    2228"USAGE:    firstoct(i); i ideal
    2329RETURN:   number: the number of points of V(i) lying in the first orthant
    24 ASSUME:   i is a Groebner basis 
     30ASSUME:   i is a Groebner basis
    2531SEE ALSO: signcnd
    2632EXAMPLE:  example firstoct; shows an example"
    … …  
    112118
    113119  // We initialize the cumulative variables
    114   M = matrix(1,1,1); 
     120  M = matrix(1,1,1);
    115121  Exponents = list(list());
    116122  Signs = list(list());
    … …  
    131137
    132138
    133     // We determine the list of signs which correspond to a nonzero 
     139    // We determine the list of signs which correspond to a nonzero
    134140    // number of roots
    135141    numberOfNonZero = 3;
    … …  
    146152      numberOfNonZero--;
    147153    }
    148    
     154
    149155    if (c[3,1] != 0) {
    150156      signi = signi + list(-1);
    … …  
    152158      numberOfNonZero--;
    153159    }
    154    
    155     // We now determine the little matrix we'll work with, 
     160
     161    // We now determine the little matrix we'll work with,
    156162    // and the list of exponents
    157163    if (numberOfNonZero == 3) {
    … …  
    162168      Mi[1,2] = 1;
    163169      if (c[1,1] != 0 && c[2,1] != 0) { // 0,1
    164         Mi[2,1] = 0;
    165         Mi[2,2] = 1;
     170        Mi[2,1] = 0;
     171        Mi[2,2] = 1;
    166172      } else {if (c[1,1] != 0 && c[3,1] != 0) { // 0,-1
    167         Mi[2,1] = 0;
    168         Mi[2,2] = -1;
     173        Mi[2,1] = 0;
     174        Mi[2,2] = -1;
    169175      } else { // 1,-1
    170         Mi[2,1] = 1;
    171         Mi[2,2] = -1;
     176        Mi[2,1] = 1;
     177        Mi[2,2] = -1;
    172178      }}
    173179      exponentsi = list(0,1);
    … …  
    176182      exponentsi = list(0);
    177183    }}}
    178    
     184
    179185    // We store the Sturm Queries we'll need later
    180186    if (numberOfNonZero == 2) {
    … …  
    194200  }
    195201
    196   // At this point, we have the cumulative matrix, 
    197   // the vector of exponents and the matching sign conditions. 
     202  // At this point, we have the cumulative matrix,
     203  // the vector of exponents and the matching sign conditions.
    198204  // We have to solve the big linear system to finish.
    199205
    … …  
    209215    if (j <= size(SQvalues)) {
    210216      if (SQpositions[j] == i) {
    211         SQs[i,1] = SQvalues[j];
    212         j++;
     217        SQs[i,1] = SQvalues[j];
     218        j++;
    213219      } else {
    214220      SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
    215       } 
     221      }
    216222    } else {
    217         SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
     223        SQs[i,1] = sturmquery(evalp(Exponents[i],P),B,I);
    218224    }
    219225  }
    … …  
    244250  echo = 0;
    245251
    246   print("The output of signcnd is a list of two lists. Both lists have the 
     252  print("The output of signcnd is a list of two lists. Both lists have the
    247253same");
    248254  print("length. That length is the number of sign conditions realized by the");
    249   print ("polynomials of P on the set V(i). In this example, that number 
     255  print ("polynomials of P on the set V(i). In this example, that number
    250256is");
    251257  print("print(size(l[1]));");
    … …  
    335341  int i,j;
    336342  int found;
    337  
     343
    338344  found = 0;
    339345  i = 1;
    … …  
    345351    } else {
    346352      while(j <= size(a)) {
    347         found = found && a[j] == b[i][j];
    348         j++;
     353        found = found && a[j] == b[i][j];
     354        j++;
    349355      }
    350356    }
    351357    i++;
    352358  }
    353  
     359
    354360  if (found) {
    355361    return (i-1);
    … …  
    375381  result[la*lb] = 0;
    376382
    377  
     383
    378384  for (i = 0;i < lb;i++) {
    379385    for (j = 0;j < la;j++) {
    … …  
    398404///////////////////////////////////////////////////////////////////////////////
    399405
    400 static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate 
     406static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate
    401407{
    402408  int i;

  • Singular/LIB/urrcount.lib

    r58f0e9 rd3b3ab  
    1 // E. Tobis  12.Nov.2004
    2 // last change 16. Apr. 2005 (G.-M. Greuel)
     1// $Id: urrcount.lib,v 1.2 2005-05-02 12:24:16 Singular Exp $
     2// E. Tobis  12.Nov.2004, April 2004
     3// last change 1. May 2005 (G.-M. Greuel)
    34///////////////////////////////////////////////////////////////////////////////
    45category="Symbolic-numerical solving"
    56info="
    67LIBRARY: urrcount.lib   Counting number of real roots of univariate polynomial
    7 AUTHOR:                 Enrique A. Tobis, etobis@dc.uba.ar 
     8AUTHOR:                 Enrique A. Tobis, etobis@dc.uba.ar
    89
    910OVERVIEW:  Routines for bounding and counting the number of real roots of a
    10            univariate polynomial by using Strum and Sturm-Habicht sequences.
    11            References:
     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
    1221
    1322PROCEDURES:
    … …  
    2635  sturmhaseq(p)    A Sturm-Habicht Sequence of a polynomial
    2736  reverse(l)       Reverses a list
    28   nrroots(p)        The number of real roots of p
     37  nrroots(p)       The number of real roots of p
    2938  isparam(p)       Returns 0 iff the polynomial has non-parametric coefficients
    3039
    … …  
    7685      i = size(ar);
    7786      while (!ispar && (i >= 1)) {
    78         ispar = ispar || (isparametric(ar[i]));
    79         i--;
     87        ispar = ispar || (isparametric(ar[i]));
     88        i--;
    8089      }
    8190    } else {
    … …  
    8392      int j;
    8493      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--;
     94      while (!ispar && (i >= 1)) {
     95        j = nrows(m);
     96        while (!ispar && (j >= 1)) {
     97          ispar = ispar || (isparametric(m[i,j]));
     98          j--;
     99        }
     100        i--;
    92101      }
    93102    }
    … …  
    97106    while (!ispar && (i >= 1)) {
    98107      if ((typeof(#[i]) != "poly") && (typeof(#[i]) != "number") &&
    99           typeof(#[i]) != "int") {
    100         ERROR("This procedure only works with lists of polynomials");
     108          typeof(#[i]) != "int") {
     109              ERROR("This procedure only works with lists of polynomials");
    101110      }
    102111      ispar = ispar || (isparametric(#[i]));
    … …  
    161170SEE ALSO:  isuni
    162171EXAMPLE:   example whichvariable; shows an example"
    163 { 
     172{
    164173  int i;
    165174  int found;
    … …  
    219228    ERROR("This procedure cannot operate with parametric arguments");
    220229  }
    221  
     230
    222231  lastsign = sign(l[1]);
    223232
    … …  
    240249}
    241250///////////////////////////////////////////////////////////////////////////////
    242 proc boundBuFou(poly p,number a,number b) 
     251proc boundBuFou(poly p,number a,number b)
    243252"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 
     253RETURN:    number: an upper bound for the number of real roots of p in (a,b],
     254           with the same parity as the actual number of roots (using the
    246255           Budan-Fourier Theorem)
    247256ASSUME:    p is a univarite polynomials with rational coefficients
    … …  
    273282
    274283  d = deg(p);
    275  
     284
    276285  // We calculate the list of derivatives
    277286
    … …  
    425434proc allrealst(poly p)
    426435"USAGE:     allrealst(p); poly p
    427 RETURN:    int: 1 iff all the roots of p are real, 0 otherwise 
     436RETURN:    int: 1 iff all the roots of p are real, 0 otherwise
    428437           Checks whether all the roots of p are real by using Sturm's Theorem
    429438ASSUME:    p is a univarite polynomials with rational coefficients
    … …  
    495504}
    496505///////////////////////////////////////////////////////////////////////////////
    497 proc sturm(poly p,number a,number b) 
     506proc sturm(poly p,number a,number b)
    498507"USAGE:     sturm(p,a,b); poly p, number a,b
    499508RETURN:    number: the number of real roots of p in (a,b]
    … …  
    562571RETURN:    list: a Sturm sequence of p
    563572ASSUME:    p is a univarite polynomials with rational coefficients
     573THEORY:    The Sturm sequence of p (also called remainder sequence) is the
     574           sequence begininng with p, p' and goes on with minus the remainder
     575           of the two previous polynomials, until the remainder is zero.
     576           See: Basu, Pollack, Roy, "Algorithms in Real Algebraic Geometry",
     577           Springer, 2003.
    564578SEE ALSO:  sturm,sturmhaseq
    565579EXAMPLE:   example sturmseq; shows an example"
    … …  
    568582  poly variable;
    569583  int i;
    570  
     584
    571585  variable = isuni(p);
    572  
     586
    573587  if (isparam(p)) {
    574588    ERROR("This procedure cannot operate with parametric arguments");
    … …  
    578592    ERROR ("p must be a univariate polynomial");
    579593  }
    580  
     594
    581595  // The two first polynomials in Sturm's sequence
    582596  stseq = list();
    … …  
    595609
    596610  // Right now, we have gcd(P,P') in stseq[size(stseq)];
    597  
     611
    598612  for (i = size(stseq)-1;i >= 1;i--) {
    599613    stseq[i] = stseq[i]/(sign(leadcoef(stseq[size(stseq)]))*stseq[size(stseq)]);
    600614    stseq[i] = stseq[i]/abs(leadcoef(stseq[i]));
    601615  }
    602  
     616
    603617  // We divide the gcd by itself
    604618  stseq[size(stseq)] = sign(leadcoef(stseq[size(stseq)]));
    605  
     619
    606620  return (stseq);
    607621}
    … …  
    650664proc sturmha(poly P,number a,number b)
    651665"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 
     666RETURN:    number: the number of real roots of p in (a,b) (using a
    653667           Sturm-Habicht sequence)
    654668SEE ALSO:  sturm,allreal
    … …  
    717731///////////////////////////////////////////////////////////////////////////////
    718732proc sturmhaseq(poly P)
    719 "USAGE:     sturmhaseq(P); P poly. P must be univariate
     733"USAGE:     sturmhaseq(P); P poly.
    720734RETURN:    list: the nonzero polynomials of the Sturm-Habicht sequence of P
     735ASSUME:    P is a univariate polynomial.
     736THEORY:    The Sturm-Habicht sequence (also subresultant sequence) is closely
     737           related to the Sturm sequence, but behaves better with respect to
     738           the size of the coefficients. It is defined via subresultants.
     739           See: Basu, Pollack, Roy, "Algorithms in Real Algebraic Geometry",
     740           Springer, 2003.
    721741SEE ALSO:  sturm,sturmseq,sturmha
    722742EXAMPLE:   example sturmhaseq; shows an example"
    … …  
    767787//      T[k-1+2] = SR[k-1+2];
    768788      srbar[k-1+2] = leadcoef(SR[k-1+2]);
    769     } 
     789    }
    770790    if (k < j-1) {
    771791      // Computation of sr[k+2]
    772792      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];
     793        srbar[j-l-1+2] = ((-1)^l)*(srbar[j-1+2]*srbar[j-l+2])/sr[j+2];
    774794    }
    775795    sr[k+2] = srbar[k+2];
    … …  
    795815    if (typeof(SR[i]) != "none") {
    796816      if (SR[i] != 0) {
    797         filtered = insert(filtered,SR[i]);
     817        filtered = insert(filtered,SR[i]);
    798818      }
    799819    }
    … …  
    847867    av = -x;
    848868  }
    849  
     869
    850870  return (av);
    851871}
    … …  
    9931013    } else {
    9941014      temp = lastsign * sign(l[i]);
    995      
     1015
    9961016      if (temp < 0) {
    997         sc++;
     1017        sc++;
    9981018      } else {
    999         if (nofzeros == 2) {
    1000           sc = sc + 2;
    1001         }
     1019        if (nofzeros == 2) {
     1020          sc = sc + 2;
     1021        }
    10021022      }
    10031023      nofzeros = 0;
    … …  
    10091029}
    10101030///////////////////////////////////////////////////////////////////////////////
    1011 ///////////////////////////////////////////////////////////////////////////////
    1012 
     1031
     1032
Note: See TracChangeset for help on using the changeset viewer.