// $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)
///////////////////////////////////////////////////////////////////////////////
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 the 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 "mrrcount.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 orthant
ASSUME:   i is 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. i must be a grobner basis
RETURN:    list: the sign conditions realized by the polynomials of P on V(I).
           See the example for an explanation of the output.
SEE ALSO:  firstoct
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);
  echo = 0;

  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
is");
  print("print(size(l[1]));");
  print(size(l[1]));
  print("Each element of the first list indicates a sign condition of the");
  print("polynomials of P. For example,");
  print("print(l[1][2]);");
  print(l[1][2]);
  print("means P[1] > 0,P[2] = 0");
  print("Each element of the second list indicates how many elements of V(I)");
  print("give rise to the sign condition expressed by the same position on the");
  print("first list. For example");
  print("print(l[2][2]);");
  print(l[2][2]);
  print("indicates that exactly 1 elemnt of V(I) gives rise to the condition");
  print("P[1] > 0,P[2] = 0.");
  print("The procedure psigncnd performs some pretty printing on this output.");
}
///////////////////////////////////////////////////////////////////////////////

proc psigncnd(ideal P,list l)
"USAGE:     psigncnd(P,I); 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);
}
///////////////////////////////////////////////////////////////////////////////