source: git/Singular/LIB/rootsmr.lib @ a57b65

/////////////////////////////////////////////////////////////////////////////
version="version rootsmr.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Teaching";
info="
LIBRARY: rootsmr.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. Two methods are implemented: deterministic
           computation of the number of roots, via the signature of a certain
           bilinear form (nrRootsDeterm); and a rational univariate projection,
           using a pseudorandom polynomial (nrRootsProbab). It also includes a
           command to verify the correctness of the pseudorandom answer.
           References: Basu, Pollack, Roy, \"Algorithms in Real Algebraic
           Geometry\", Springer, 2003.

PROCEDURES:
 nrRootsProbab(I)    Number of real roots of 0-dim ideal (probabilistic)
 nrRootsDeterm(I)    Number of real roots of 0-dim ideal (deterministic)
 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
";
///////////////////////////////////////////////////////////////////
LIB "linalg.lib";   // We use charpoly
LIB "rootsur.lib"; // We use varsigns

proc nrRootsProbab(ideal I, list #)
"USAGE:     nrRootsProbab(I,[n]); ideal I, int n
RETURN:    int: the number of real roots of the ideal I by a probabilistic
           algorithm
ASSUME:    If I is not a Groebner basis, then a Groebner basis will be computed
           by using std. If I is already a Groebner basis (i.e. if
           attrib(I,"isSB"); returns 1) then this Groebner basis will be
           used, hence it must be one w.r.t. (any) global ordering. This may
           be useful if the ideal is known to be a Groebner basis or if it
           can be computed faster by a different method.
NOTE:      If n<10 is given, n is the number of digits being used for
           constructing a random characteristic polynomial, a bigger n is
           more safe but slower (default: n=5).
           If printlevel>0 the number of complex solutions is displayed
           (default: printlevel=0).
SEE ALSO:  nrroots, nrRootsDeterm, randcharpoly, solve
EXAMPLE:   example nrRootsProbab; shows an example"
{
  //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)).
  //Then the algorithm needs:
  //1 std(I) and ~n NF computations (of randcharpoly w.r.t. I)

  if (isparam(I)) {
    ERROR("This procedure cannot operate with parametric arguments");
  }
  int pr = printlevel-voice+2;
  int v;
  int n=5;
  if (size(#) == 1) {
    n=#[1];
  }
  if (attrib(I,"isSB")!=1) {
    I = std(I);
  }

  ideal b = qbase(I);
  v = size(b);
  if (v == 0) {
    ERROR("ideal is not 0-dimensional");
  }
  dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity");

  poly p = randcharpoly(b,I,n);

  return (nrroots(p));
}

example
{
  echo = 2;
  ring r = 0,(x,y,z),lp;
  ideal i = (x-1)*(x-2),(y-1)^3*(x-y),(z-1)*(z-2)*(z-3)^2;
  nrRootsProbab(i);       //no of real roots (using internally std)

  i = groebner(i);        //using the hilbert driven GB computation
  int pr = printlevel;
  printlevel = 2;
  nrRootsProbab(i);
  printlevel = pr;
}
///////////////////////////////////////////////////////////////////////////////

proc nrRootsDeterm(ideal I)
"USAGE:     nrRootsDeterm(I); ideal I
RETURN:    int: the number of real roots of the ideal I by a deterministic
           algorithm
ASSUME:    If I is not a Groebner basis, then a Groebner basis will be computed
           by using std. If I is already a Groebner basis (i.e. if
           attrib(I,"isSB"); returns 1) then this Groebner basis will be
           used, hence it must be one w.r.t. (any) global ordering. This may
           be useful if the ideal is known to be a Groebner basis or if it
           can be computed faster by a different method.
NOTE:      If printlevel>0 the number of complex solutions is displayed
           (default: printlevel=0). The procedure nrRootsProbab is usually faster.
SEE ALSO:  nrroots, nrRootsProbab, sturmquery, solve
EXAMPLE:   example nrRootsDeterm; shows an example"
{
  //Note on complexity: Let n = no of complex roots of I (= vdim(std(I)).
  //Then the algotithm needs:
  //1 std(I) and (1/2)n*(n+1)^2 ~ 1/2n^3 NF computations (of monomials w.r.t. I)

  if (isparam(I)) {
    ERROR("This procedure cannot operate with parametric arguments");
  }
  int pr = printlevel-voice+2;
  int v;

  if (attrib(I,"isSB")!=1) {
    I = std(I);
  }

  ideal b = qbase(I);
  v = size(b);
  if (v == 0) {
    ERROR("ideal is not 0-dimensional");
  }
  dbprint(pr,"//ideal has " +string(v)+ " complex solutions, counted with multiplicity");

  return (sturmquery(1,b,I));
}

example
{
  echo = 2;
  ring r = 0,(x,y,z),lp;
  ideal I = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2;
  nrRootsDeterm(I);       //no of real roots (using internally std)

  I = groebner(I);        //using the hilbert driven GB computation
  int pr = printlevel;
  printlevel = 2;
  nrRootsDeterm(I);
  printlevel = pr;
}
///////////////////////////////////////////////////////////////////////////////

proc symsignature(matrix m)
"USAGE:     symsignature(m); m matrix. m must be symmetric.
RETURN:    int: the signature of m
SEE ALSO:  matbil,sturmquery
EXAMPLE:   example symsignature; shows an example"
{
  int positive, negative, i, j;
  list l;
  poly variable;

  if (isparam(m)) {
    ERROR("This procedure cannot operate with parametric arguments");
  }

  if (!isSquare(m)) {
    ERROR ("m must be a square matrix");
  }

  // We check whether m is symmetric
  for (i = 1;i <= nrows(m);i++) {
    for (j = i;j <= nrows(m);j++) {
      if (m[i,j] != m[j,i]) {
        ERROR ("m must be a symmetric matrix");
      }
    }
  }

  poly f = charpoly(m); // Uses the last variable of the ring

  for (i = size(f);i >= 1;i--) {
    l[i] = leadcoef(f[i]);
  }
  positive = varsigns(l);

  variable = var(nvars(basering)); // charpoly uses the last variable
  f = subst(f,variable,-variable);

  for (i = size(f);i >= 1;i--) {
    l[i] = leadcoef(f[i]);
  }

  negative = varsigns(l);
  return (positive - negative);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);

  matrix m = matbil(1,b,i);
  symsignature(m);
}
///////////////////////////////////////////////////////////////////////////////

proc sturmquery(poly h,ideal B,ideal I)
"USAGE:     sturmquery(h,b,i); h poly, b,i ideal
RETURN:    int: the Sturm query of h in V(i)
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"
{
  if (isparam(h) || isparam(B) || isparam(I)) {
    ERROR("This procedure cannot operate with parametric arguments");
  }

  return (mysymmsig(matbil(h,B,I)));
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);

  sturmquery(1,b,i);
}
///////////////////////////////////////////////////////////////////////////////

static proc mysymmsig(matrix m)
// returns the signature of a square symmetric matrix m
{
  int positive, negative, i;
  list l;
  poly variable;

  poly f = charpoly(m); // Uses the last variable of the ring

  for (i = size(f);i >= 1;i--) {
    l[i] = leadcoef(f[i]);
  }
  positive = varsigns(l);

  variable = var(nvars(basering)); // charpoly uses the last variable
  f = subst(f,variable,-variable);

  for (i = size(f);i >= 1;i--) {
    l[i] = leadcoef(f[i]);
  }

  negative = varsigns(l);
  return (positive - negative);
}
///////////////////////////////////////////////////////////////////////////////

proc matbil(poly h,ideal B,ideal I)
"USAGE:    matbil(h,b,i); h poly, b,i ideal
RETURN:    matrix: the matrix of the bilinear form (f,g) |-> trace(m_fhg),
           m_fhg = multiplication with fhg on r/i
ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
           r = basering
SEE ALSO:  matmult,tracemult
EXAMPLE:   example matbil; shows an example"
{
  matrix m[size(B)][size(B)];
  poly f;
  int k,l;
  //h = reduce(h,I);

  for (k = 1; k <= size(B); k++) {
    for (l = 1; l <= k; l++) {
      m[k,l] = tracemult(h*B[k]*B[l],B,I)[1];
      m[l,k] = m[k,l]; // The matrix we are trying to compute is symmetric
    }
   }
  return(m);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);
  poly f = x3-xy+y-13+x4-y2x;

  matrix m = matbil(f,b,i);
  print(m);

}
///////////////////////////////////////////////////////////////////////////////

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 (faster than matmult + trace)
ASSUME:    I is given by a Groebner basis and B is an ordered monomial basis of r/I
SEE ALSO:  matmult,trace
EXAMPLE:   example tracemult; shows an example"
{
  int k; // Iterates over the basis monomials
  int l; // Iterates over the rows of the matrix
  list coordinates;
  number m;
  poly g;

  //f = reduce(f,I);
  for (k = 1; k <= size(B); k++) {
    l=1;
    g = reduce(f*B[k],I);
    while (l <= k) {
      if (leadmonom(g[l]) == B[k]) {
        m = m + leadcoef(g[l]);
        break;
      }
      l++;
    }
  }
  return (m);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);

  poly f = x3-xy+y-13+x4-y2x;
  matrix m = matmult(f,b,i);
  print(m);

  tracemult(f,b,i);            //the trace of m
}
///////////////////////////////////////////////////////////////////////////////

proc matmult(poly f, ideal B, ideal I)
"USAGE:     matmult(f,b,i); f poly, b,i ideal
RETURN:    matrix: the matrix of the multiplication map by f (m_f) on r/i
           w.r.t. to the monomial basis b of r/i (r = basering)
ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
           as given by qbase(i)
SEE ALSO:  coords,matbil
EXAMPLE:   example matmult; shows an example"
{
  int k; // Iterates over the basis monomials
  int l; // Iterates over the rows of the matrix
  list coordinates;
  matrix m[size(B)][size(B)];

  //f = reduce(f,I);
  for (k = 1;k <= size(B);k++) {
    coordinates = coords(f*(B[k]),B,I); // f*x_k written on the basis B
    for (l = 1;l <= size(B);l++) {
      m[l,k] = coordinates[l];
    }
  }
  return (m);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);

  poly f = x3-xy+y-13+x4-y2x;
  matrix m = matmult(f,b,i);
  print(m);
}
///////////////////////////////////////////////////////////////////////////////

proc coords(poly f,ideal B,ideal I)
"USAGE:     coords(f,b,i), f poly, b,i ideal
RETURN:    list of numbers: the coordinates of the class of f (mod i)
           in the monomial basis b
ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
           r = basering
SEE ALSO:  matmult,matbil
KEYWORDS:  coordinates
EXAMPLE:   example coords; shows an example"
{
  // We assume the basis is sorted according to the ring order
  poly g;
  int k,l=1,1;
  list coordinates;
  int N = size(B);

  // We first compute the normal form of f w.r.t. I
  g = reduce(f,I);
  int n = size(g);    //allways n <= N

  while (k <= N) {
    if (leadmonom(g[l]) == B[k]) {
      coordinates[k] = leadcoef(g[l]);
      l++;
    } else {
      coordinates[k] = number(0);
    }
    k++;
  }
  return (coordinates);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = x4-y2x,y2-13;
  poly f = x3-xy+y-13+x4-y2x;
  i = std(i);
  ideal b = qbase(i);
  b;
  coords(f,b,i);
}
///////////////////////////////////////////////////////////////////////////////

static proc isSquare(matrix m)
// returns 1 if and only if m is a square matrix
{
  return (nrows(m)==ncols(m));
}
///////////////////////////////////////////////////////////////////////////////

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.
           If n<10 is given, it is the number of digits being used for the
           pseudorandom coefficients (default: n=5)
ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
           r = basering
NOTE:      shows a warning if printlevel>0 (default: printlevel=0)
KEYWORDS:  rational univariate projection
EXAMPLE:   example randcharpoly; shows an example"
{
  int pr = printlevel - voice + 2;
  poly p;
  poly generic;
  list l;
  matrix m;
  poly q;

  if (size(#) == 1) {
    generic = randlinpoly(#[1]);
  } else {
    generic = randlinpoly();
  }

  p = reduce(generic,I);
  m = matmult(p,B,I);
  q = charpoly(m);

  dbprint(pr,"*********************************************************************");
  dbprint(pr,"* WARNING: This polynomial was obtained using  pseudorandom numbers.*");
  dbprint(pr,"* If you want to verify the result, please use the command          *");
  dbprint(pr,"*                                                                   *");
  dbprint(pr,"* verify(p,b,i)                                                     *");
  dbprint(pr,"*                                                                   *");
  dbprint(pr,"* where p is the polynomial I returned, b is the monomial basis     *");
  dbprint(pr,"* used, and i the Groebner basis of the ideal                       *");
  dbprint(pr,"*********************************************************************");

  return(q);
}
example
{
  echo = 2;
  ring r = 0,(x,y,z),dp;
  ideal i = (x-1)*(x-2),(y-1),(z-1)*(z-2)*(z-3)^2;
  i = std(i);
  ideal b = qbase(i);
  poly p = randcharpoly(b,i);
  p;
  nrroots(p); // See nrroots in urrcount.lib

  int pr = printlevel;
  printlevel = pr+2;
  p = randcharpoly(b,i,5);
  nrroots(p);
  printlevel = pr;
}

///////////////////////////////////////////////////////////////////////////////

proc verify(poly p,ideal B,ideal I)
"USAGE:     verify(p,B,I); p poly, B,I,ideal
RETURN:    integer: 1 if and only if the polynomial p splits the points of V(I).
           It's used to check the result of randcharpoly
ASSUME:    I is given by a Groebner basis and B is an ordered monomial basis of r/I,
           r = basering
NOTE:      comments the result if printlevel>0 (default: printlevel=0)
SEE ALSO:  randcharpoly
EXAMPLE:   example verify; shows an example"
{
  int pr = printlevel - voice + 2;
  poly sqr_free;
  int correct;
  poly variable;

  if (isparam(p) || isparam(B) || isparam(I)) {
    ERROR("This procedure cannot operate with parametric arguments");
  }

  variable = isuni(p);
  sqr_free = p/gcd(p,diff(p,variable));
  correct = (mat_rk(matbil(1,B,I)) == deg(sqr_free));

  if (correct) {
    dbprint(pr,"//Verification successful");
  } else {
    dbprint(pr,"//The choice of random numbers was not useful");
    dbprint(pr,"//You might want to try randcharpoly with a larger number of digits");
  }
  return (correct);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  poly f = x3-xy+y-13+x4-y2x;
  ideal i = x4-y2x,y2-13;
  i = std(i);
  ideal b = qbase(i);
  poly p = randcharpoly(b,i);
  verify(p,b,i);
}
///////////////////////////////////////////////////////////////////////////////

proc randlinpoly(list #)
"USAGE:     randlinpoly(); randlinpoly(n); n int
RETURN:    poly: linear combination  of the variables of the ring, with
           pseudorandom coefficients. If n<10 is given, it is the number of
           digits being used for the range of the coefficients (default: n=5)
SEE ALSO:  randcharpoly;
EXAMPLE:   example randlinpoly; shows an example"
{
  int n,i;
  poly p = 0;
  int ndigits = 5;

  if (size(#) == 1) {
    ndigits = #[1];
  }

  n = nvars(basering);
  for (i = 1;i <= n;i++) {
    p = p + var(i)*random(1,10^ndigits);
  }
  return (p);
}
example
{
  echo = 2;
  ring r = 0,(x,y,z,w),dp;
  poly p = randlinpoly();
  p;
  randlinpoly(5);
}
///////////////////////////////////////////////////////////////////////////////

proc powersums(poly f,ideal B,ideal I)
"USAGE:     powersums(f,b,i); f poly; b,i ideal
RETURN:    list: the powersums of the results of evaluating f at the zeros of I
ASSUME:    i is a Groebner basis and b is an ordered monomial basis of r/i,
           r = basering
SEE ALSO:  symmfunc
EXAMPLE:   example symmfunc; shows an example"
{
  int N,k;
  list sums;

  N = size(B);
  for (k = 1;k <= N;k++) {
    sums = sums + list(leadcoef(trace(matmult(f^k,B,I))));
  }
  return (sums);
}
example
{
  echo = 2;
  ring r = 0,(x,y,z),dp;

  ideal i = (x-1)*(x-2),(y-1),(z+5); // V(I) = {(1,1,-5),(2,1,-5)}
  i = std(i);

  ideal b = qbase(i);
  poly f = x+y+z;
  list psums = list(-2-3,4+9); // f evaluated at V(I) gives {-3,-2}
  list l = powersums(f,b,i);
  psums;
  l;
}
///////////////////////////////////////////////////////////////////////////////

proc symmfunc(list S)
"USAGE:     symmfunc(s); s list
RETURN:    list: the symmetric functions of the roots of a polynomial, given
                 the power sums of those roots.
SEE ALSO:  powersums
EXAMPLE:   example symmfunc; shows an example"
{
  // Takes the list of power sums and returns the symmetric functions
  list a;
  int j,l,N;
  number sum;

  N = size(S);
  a[N+1] = 1; // We set the length of the list and initialize its last element.

  for (l = N - 1;l >= 0;l--) {
    sum = 0;
    for (j = l + 1;j <= N;j++) {
      sum = sum + ((a[j+1])*(S[j-l]));
    }
    sum = -sum;
    a[l+1] = sum/(N-l);
  }

  a = reverse(a);
  return (a);
}
example
{
  echo = 2;
  ring r = 0,x,dp;
  poly p = (x-1)*(x-2)*(x-3);
  list psums = list(1+2+3,1+4+9,1+8+27);
  list l = symmfunc(psums);
  l;
  p; // Compare p with the elements of l
}
///////////////////////////////////////////////////////////////////////////////

proc univarpoly(list l)
"USAGE:     univarpoly(l); l list
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"
{
  poly p;
  int i,n;

  n = size(l);
  for (i = 1;i <= n;i++) {
    p = p + l[i]*var(1)^(n-i);
  }
  return (p);
}
example
{
  echo = 2;
  ring r = 0,x,dp;
  list l = list(1,2,3,4,5);
  poly p = univarpoly(l);
  p;
}
///////////////////////////////////////////////////////////////////////////////

proc qbase(ideal i)
"USAGE:    qbase(I); I zero-dimensional ideal
RETURN:   ideal: A monomial basis of the quotient between the basering and the
          ideal I, sorted according to the basering order.
SEE ALSO: kbase
KEYWORDS: zero-dimensional
EXAMPLE:  example qbase; shows an example"
{
  ideal b;

  b = kbase(i);
  b = reverseideal(sort(b)[1]); // sort sorts in ascending order
  return (b);
}
example
{
  echo = 2;
  ring r = 0,(x,y,z),dp;

  ideal i = 2x2,-y2,z3;
  i = std(i);
  ideal b = qbase(i);
  b;
  b = kbase(i);
  b; // Compare this with the result of qbase
}
///////////////////////////////////////////////////////////////////////////////

static proc reverseideal(ideal b) // Returns b reversed
{
  int i;
  ideal result;

  result = b[1];
  for (i = 2;i <= size(b);i++) {
    result = b[i], result;
  }
  return (result);
}
///////////////////////////////////////////////////////////////////////////////