singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / signcond.lib

///////////////////////////////////////////////////////////////////////////
version="version signcond.lib 4.0.0.0 Jun_2013 "; // $Id: 5d7118844e1de435273225d8cd9cb0962cb86d8e $
category="Symbolic-numerical solving";
info="
LIBRARY: signcond.lib Routines for computing realizable sign conditions
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 by polynomials of P on a V(I)
  psigncnd(P,l)  Pretty prints the output of signcnd (l)
  firstoct(I)    The number of elements of V(I) with every coordinate > 0

KEYWORDS: real roots,sign conditions
";

LIB "rootsmr.lib";
LIB "linalg.lib";
///////////////////////////////////////////////////////////////////////////////

proc firstoct(ideal I)
"USAGE:    firstoct(I); I ideal
RETURN:   number: the number of points of V(I) lying in the first octant
ASSUME:   I is given by a Groebner basis.
SEE ALSO: signcnd
EXAMPLE:  example firstoct; shows an example"
{
  ideal firstoctant;
  int j;
  list result;
  int n;

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

  for (j = nvars(basering);j > 0;j--) {
    firstoctant = firstoctant + var(j);
  }

  result = signcnd(firstoctant,I);

  list fst;
  for (j = nvars(basering);j > 0;j--) {
    fst[j] = 1;
  }

  n = isIn(fst,result[1]);

  if (n != -1) {
    return (result[2][n]);
  } else {
    return (0);
  }
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = (x-2)*(x+3)*x,y*(y-1);
  firstoct(i);
}
///////////////////////////////////////////////////////////////////////////////

proc signcnd(ideal P,ideal I)
"USAGE:     signcnd(P,I); ideal P,I
RETURN:    list: the sign conditions realized by the polynomials of P on V(I).
           The output of signcnd is a list of two lists. Both lists have the
           same length. This length is the number of sign conditions realized
           by the polynomials of P on the set V(i).
           Each element of the first list indicates a sign condition of the
           polynomials of P.
           Each element of the second list indicates how many elements of V(I)
           give rise to the sign condition expressed by the same position on
           the first list.
           See the example for further explanations of the output.
ASSUME:    I is a Groebner basis.
NOTE:      The procedure psigncnd performs some pretty printing of this output.
SEE ALSO:  firstoct, psigncnd
EXAMPLE:   example signcnd; shows an example"
{
  ideal B;

  // Cumulative stuff
  matrix M;
  matrix SQs;
  matrix C;
  list Signs;
  list Exponents;

  // Used to store the precalculated SQs
  list SQvalues;
  list SQpositions;

  int i;

  // Variables for each step
  matrix Mi;
  matrix M3x3[3][3];
  matrix M3x3inv[3][3]; // Constant matrices
  matrix c[3][1];
  matrix sq[3][1];
  int j;
  list exponentsi;
  list signi;
  int numberOfNonZero;

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

  M3x3 = matrix(1,3,3);
  M3x3 = 1,1,1,0,1,-1,0,1,1; // The 3x3 matrix
  M3x3inv = inverse(M3x3);

  // First, we compute sturmquery(1,V(I))
  I = groebner(I);
  B = qbase(I);
  sq[1,1] = sturmquery(1,B,I); // Number of real roots in V(I)
  SQvalues = SQvalues + list(sq[1,1]);
  SQpositions = SQpositions + list(1);

  // We initialize the cumulative variables
  M = matrix(1,1,1);
  Exponents = list(list());
  Signs = list(list());

  i = 1;

  while (i <= size(P)) { // for each poly in P

    sq[2,1] = sturmquery(P[i],B,I);
    sq[3,1] = sturmquery(P[i]^2,B,I);


    c = M3x3inv*sq;

    // We have to eliminate the 0 elements in c
    exponentsi = list();
    signi = list();


    // We determine the list of signs which correspond to a nonzero
    // number of roots
    numberOfNonZero = 3;

    if (c[1,1] != 0) {
      signi = list(0);
    } else {
      numberOfNonZero--;
    }

    if (c[2,1] != 0) {
      signi = signi + list(1);
    } else {
      numberOfNonZero--;
    }

    if (c[3,1] != 0) {
      signi = signi + list(-1);
    } else {
      numberOfNonZero--;
    }

    // We now determine the little matrix we'll work with,
    // and the list of exponents
    if (numberOfNonZero == 3) {
      Mi = M3x3;
      exponentsi = list(0,1,2);
    } else {if (numberOfNonZero == 2) {
      Mi = matrix(1,2,2);
      Mi[1,2] = 1;
      if (c[1,1] != 0 && c[2,1] != 0) { // 0,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;
      } else { // 1,-1
        Mi[2,1] = 1;
        Mi[2,2] = -1;
      }}
      exponentsi = list(0,1);
    } else {if (numberOfNonZero == 1) {
      Mi = matrix(1,1,1);
      exponentsi = list(0);
    }}}

    // We store the Sturm Queries we'll need later
    if (numberOfNonZero == 2) {
      SQvalues = SQvalues + list(sq[2,1]);
      SQpositions = SQpositions + list(size(Exponents)+1);
    } else {if (numberOfNonZero == 3) {
      SQvalues = SQvalues + list(sq[2,1],sq[3,1]);
      SQpositions = SQpositions + list(size(Exponents)+1,size(Exponents)*2+1);
    }}

    // Now, we accumulate information
    M = tensor(Mi,M);
    Signs = expprod(Signs,signi);
    Exponents = expprod(Exponents,exponentsi);

    i++;
  }

  // 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.

  M = inverse(M);

  // We have to compute the constants vector (the Sturm Queries)

  SQs = matrix(1,size(Exponents),1);

  j = 1; // We'll iterate over the presaved SQs

  for (i = 1;i <= size(Exponents);i++) {
    if (j <= size(SQvalues)) {
      if (SQpositions[j] == i) {
        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);
    }
  }

  C = M*SQs;

  list result;
  result[2] = list();
  result[1] = list();

  // We have to filter the 0 elements of C
  for (i = 1;i <= size(Signs);i++) {
    if (C[i,1] != 0) {
      result[1] = result[1] + list(Signs[i]);
      result[2] = result[2] + list(C[i,1]);
    }
  }

  return (result);
}
example
{ echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = (x-2)*(x+3)*x,y*(y-1);
  ideal P = x,y;
  list l = signcnd(P,i);

  size(l[1]);     // = the number of sign conditions of P on V(i)

  //Each element of l[1] indicates a sign condition of the polynomials of P.
  //The following means P[1] > 0, P[2] = 0:
  l[1][2];

  //Each element of l[2] indicates how many elements of V(I) give rise to
  //the sign condition expressed by the same position on the first list.
  //The following means that exactly 1 element of V(I) gives rise to the
  //condition P[1] > 0, P[2] = 0:
  l[2][2];
}
///////////////////////////////////////////////////////////////////////////////

proc psigncnd(ideal P,list l)
"USAGE:     psigncnd(P,l); ideal P, list l
RETURN:    list: a formatted version of l
SEE ALSO:  signcnd
EXAMPLE:   example psigncnd; shows an example"
{
  string s;
  int n = size(l[1]);
  int i;

  for (i = 1;i <= n;i++) {
    s = s + string(l[2][i]) + " elements of V(I) satisfy " + psign(P,l[1][i])
        + sprintf("%n",12);
  }
  return(s);
}
example
{
  echo = 2;
  ring r = 0,(x,y),dp;
  ideal i = (x-2)*(x+3)*x,(y-1)*(y+2)*(y+4);
  ideal P = x,y;
  list l = signcnd(P,i);
  psigncnd(P,l);
}
///////////////////////////////////////////////////////////////////////////////

static proc psign(ideal P,list s)
{
  int i;
  int n = size(P);
  string output;

  output = "{P[1]";

  if (s[1] == -1) {
    output = output + " < 0";
  }
  if (s[1] == 0) {
    output = output + " = 0";
  }
  if (s[1] == 1) {
    output = output + " > 0";
  }

  for (i = 2;i <= n;i++) {
    output = output + ",";
    output = output + "P[" + string(i) + "]";
    if (s[i] == -1) {
      output = output + " < 0";
    }
    if (s[i] == 0) {
      output = output + " = 0";
    }
    if (s[i] == 1) {
      output = output + " > 0";
    }

  }
  output = output + "}";
  return (output);
}
///////////////////////////////////////////////////////////////////////////////

static proc isIn(list a,list b) //a is a list. b is a list of lists
{
  int i,j;
  int found;

  found = 0;
  i = 1;
  while (i <= size(b) && !found) {
    j = 1;
    found = 1;
    if (size(a) != size(b[i])) {
      found = 0;
    } else {
      while(j <= size(a)) {
        found = found && a[j] == b[i][j];
        j++;
      }
    }
    i++;
  }

  if (found) {
    return (i-1);
  } else {
    return (-1);
  }
}
///////////////////////////////////////////////////////////////////////////////

static proc expprod(list A,list B) // Computes the product of the list of lists A and the list B.
{
  int i,j;
  list result;
  int la,lb;

  if (size(A) == 0) {
    A = list(list());
  }

  la = size(A);
  lb = size(B);

  result[la*lb] = 0;


  for (i = 0;i < lb;i++) {
    for (j = 0;j < la;j++) {
      result[i*la+j+1] = A[j+1] + list(B[i+1]);
    }
  }

  return (result);
}
///////////////////////////////////////////////////////////////////////////////

static proc initlist(int n) // Returns an n-element list of 0s.
{
  list l;
  int i;
  l[n] = 0;
  for (i = 1;i < n;i++) {
    l[i] = 0;
  }
  return(l);
}
///////////////////////////////////////////////////////////////////////////////

static proc evalp(list exp,ideal P) // Elevates each polynomial in P to the appropriate
{
  int i;
  int n;
  poly result;

  n = size(exp);
  result = 1;

  for (i = 1;i <= n; i++) {
    result = result * (P[i]^exp[i]);
  }
  return (result);
}
///////////////////////////////////////////////////////////////////////////////

static proc incexp(list exp)
{
  int k;

  k = 1;

  while (exp[k] == 2) { // We assume exp is not the last exponent (i.e. 2,...,2)
    exp[k] = 0;
    k++;
  }

  // exp[k] < 2
  exp[k] = exp[k] + 1;

  return (exp);
}
///////////////////////////////////////////////////////////////////////////////