// stein01.lib, v 0.01 2005/09/15
//(Arnaud Bodin, last modified 2005/09/30)
///////////////////////////////////////////////////////////////////////////////
version="critic01.lib, v 0.01 2005/09/15";
category="Singularities";   // line may be commented for some old version of Singular
info="
LIBRARY:  stein01.lib      Around Stein factorisation of polynomials
AUTHORS:  Arnaud Bodin : email : Arnaud.Bodin@math.univ-lille1.fr

PROCEDURES:
  proc clean_gaussmatriw(matrix A);   clean a Gauss reduced matrix
  proc findminors(matrix A);   find the minors of maximal size
  proc checkhypo(poly f);  check if the algorithm applies
  proc linearfrompde(poly f);  make a linear system (via a pde due to Gao-Ruppert) to find reducible factors
  proc reducible(poly f);  main procedure; give reducible fibers as the zero of a polynomial
";  

LIB "matrix.lib";


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

//////////////////////////////////////////
////////    Matrix operations    /////////
//////////////////////////////////////////




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


proc clean_gaussmatrix(matrix A)
"USAGE:   clean_gaussmatrix(A); A matrix (being gauss-reduced)
RETURN:  matrix, obtained from A by deleting all non necessary rows and columns 
EXAMPLE: example clean_gaussmatrix; shows an example
"
{
  A = simplify(A,2+4);
  A = transpose(simplify(transpose(A),2+4));
  intvec Li, Lj;
  int i,j;
  // A priori we keep all
  for (i=1; i <= nrows(A); i++) {Li = Li,i;};
  for (j=1; j<=ncols(A);  j++) {Lj = Lj,j;};

  for (j=ncols(A); j>0; j--)
  {
    for (i=nrows(A); i>0; i--)
    {
      if(A[i,j]==1) {Li = takeoff(Li,i); Lj = takeoff(Lj,j); break;};
      if(A[i,j]!=0) {break;};
    } 
  }

  if (size(Li) > 1)
  {
    Li = Li[2..size(Li)]; Lj = Lj[2..size(Lj)];
    A = submat(A, Li,Lj);
    A = simplify(A,2+4);
    A = transpose(simplify(transpose(A),2+4));
  }
  return(A);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[6][8]=
   0,  0,    1,   0, B,0,  0,0,
   0,  0,  -C2,   0, 0,0,  7,0,
   0,  0,    0,   0, 0,0,  0,0,
   0,  0,    B,  -A, 0,2A, 0,0,
   0,  0,    0,   0, 0,0,  1,0,
   0,  0,    0,   0, 0,0,  0,1;
   print(clean_gaussmatrix(m));
}


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

proc findminors(matrix A)
"USAGE:   findminors(A); A matrix 
RETURN:  ideal, of minors of maximal size 
EXAMPLE: example findminors; shows an example
"
{
  ideal I;
  if (ncols(A) > nrows(A))
  {
    I = minor(A,nrows(A));
  }
  else
  {
    I = minor(A,ncols(A));
  };
  I = std(I);
  return(I);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(A,B,C),dp;
   matrix m[6][8]=
   0,  0,    1,   0, B,0,  0,0,
   0,  1,  -C2,   0, 0,0,  7,0,
   1,  0,    0,   0, 0,0,  0,0,
   0,  0,    B,  -A, 0,2A, 0,0,
   1,  1,    0,   0, 0,0,  1,0,
   1,  0,    0,   0, 0,0,  0,1;
   print(findminors(m));
}

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


//////////////////////////////////////////
//////  Linear system from a pde     /////
//////////////////////////////////////////

proc checkhypo(poly f)
"USAGE:   checkhypo(f); f poly
RETURN:  messages, if hypotheses are not verified
EXAMPLE: example checkhypo; shows an example
"
{
  def bring = basering;

// Number of variables 2 + 1 parameters
if (nvars(bring)!= 3)
{"Warning : incorrect number of variables !";
 "Should be three variables : 2 true variables and 1 for the result ex: (x,y,c).";};

// gcd(f-c,fx)=1
number f0= number(subst(subst(f,x,0),y,0));
ideal J = jacob(f);
if ((gcd(f-f0,diff(f,var(1)))!=1) or  (dim(std(f))>1))
{"Warning : the polynomial should verified  gcd(f-c,fx)=1 for all c !";};

// My warning : should be for all c gcd (f-c, ....
// Ok for complex number
// 
// To be changed : in fact it works if some f-c is not reduced
// But it doesn't give the right multiplicity
// Hence I have to exclude the case of composite polynomials

// Characteristic 
  intvec v= 1,0,0;
  int m = deg(f,v);
  intvec w = 0,1,1;
  int n = deg(f, w);
  if ((char(bring)>0) and (char(bring)<= (2*m-1)*n))
{"Warning : the chararcteristic of the ring is not enough big !";
 "char(r) should be zero or greater than (2deg(f,x)-1)deg(f,y)."; };
 

}
example
{ "EXAMPLE:"; echo = 2;
   ring r=3, (x,y,c), dp;
   poly f = x2-y4;
   print(checkhypo(f));
   ring r=0, (x,y,c), dp;
   poly f = x2*(xy-1)-1;
   print(checkhypo(f));
}


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

proc linearfrompde(poly f)
"USAGE:   linearfrompde(f); f poly (ring being r= 0, (x,y,c), dp)
RETURN:  matrix, whose rank give reducible fibers of f, based on Ruppert-Gao pde
EXAMPLE: example linearfrompde; shows an example
"
{

  def bring = basering;

  // Partial degrees
  intvec v= 1,0,0;
  int m = deg(f,v);
  intvec w = 0,1,1;
  int n = deg(f, w);

  // New ring with unknowns
  def bring = basering;      // user's ring, of type r=0, (x,y,c), dp;
  poly theminpoly = minpoly;
  string thechar = charstr(bring);
  string thevar = varstr(bring);
  string rstr = "ring rr = (" +  thechar + "), (aa(1..m*(n+1)),bb(1..(m+1)*n),"+ thevar  +"), dp;";  
  execute(rstr);
  minpoly = number(fetch(bring,theminpoly));

  poly ff = imap(bring,f);       // f the polynomial in (x,y)

  poly varx = var(2*m*n+m+n+1);  // x variable
  poly vary = var(2*m*n+m+n+2);  // y variable
  poly varc = var(2*m*n+m+n+3);  // c parameter


  /// g;
  poly g =0; 
  int i,j;
  for (i=1; i <= m; i++)
  { for(j=1; j<= n+1; j++)
    { g = g + aa(i+(j-1)*m)*varx^(i-1)*vary^(j-1);    }
  };

  /// h;
  poly h =0; 
  for (i=1; i <= m+1; i++)
  { for(j=1; j<= n; j++)
    { h = h + bb(i+(j-1)*(m+1))*varx^(i-1)*vary^(j-1); }
  };

  /// PDE Equation;
  poly E = (ff-varc)*(diff(g,vary)-diff(h,varx))+h*diff(ff,varx)-g*diff(ff,vary);

  // Make a linear system in aa(1),..., bb(m*n)
  // Note : not linear in c !!
  matrix C = coeffs(E,vary);
  int degy = size(C);
  ideal EE;
  for (i=1; i<= degy; i++)
  { EE = EE + ideal(coeffs(C[i,1],varx)); }


  // Make a matrix
  int p = size(EE);
  int q = 2*m*n+m+n;
  matrix M[p][q];

  for (i=1; i<=p; i++)
  {
    for (j=1; j<= m*(n+1); j++)
    {
      C = coef(EE[i],aa(j));
      if (C[1,1]!=1) {M[i,j] = C[2,1];};
    };
    for (j=m*(n+1)+1; j<= q; j++)
    {
      C = coef(EE[i],bb(j-m*(n+1)));
      if (C[1,1]!=1) {M[i,j] = C[2,1]; };
    };
  }

  // Back to user's ring
  setring bring;
  matrix MM;
  MM  = imap(rr,M);
  return(MM);

}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0, (x,y,c), dp;
   poly f = x2-y4;
   print(linearfrompde(f));
}



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

proc reducible(poly f);  main procedure; give reducible fibers as the zero of a polynomial
"USAGE:   reducible(f); f poly (ring being r= 0, (x,y,c), dp)
RETURN:  poly, whose zeroes c are corresponds to reducible fibers (f=c)
EXAMPLE: example reducible; shows an example
"
{

  checkhypo(f);
  matrix M;
  M  = linearfrompde(f);
  M = gauss_col(M);
  M = clean_gaussmatrix(M);

// Should be check (by me): M should now really be simplified
// otherwise I think it could false

  ideal I;
  I = findminors(M);
  poly g = I[1];
  return(g);

 
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0, (x,y,c), dp;
   poly f = x+y*(x-1)*(x-2)*(x-3)*(x-4);
  poly g = reducible(f);
  "Reducible fiber are given by the zeroes of :", g;
  "After factorization, reducible fibers are :"; 
  factorize(g,1);
}

//////////////////////////////////////////
//////    Internal procedures        /////
//////////////////////////////////////////

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

proc takeoff(intvec L, int n); take off n from L, size(L) >1
{
  int i; intvec LL;
  for (i=1; i<= size(L); i++)
  {      
    if (L[i]!= n)  {LL = LL,L[i];};
  }
  if (size(LL)>1) {LL = LL[2..size(LL)];};
  return(LL); 
}


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



//////////////////////////////////////////
//////     Rational fraction         /////
//////////////////////////////////////////


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

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

proc linearfrompdeF(poly p, poly q)
"USAGE:   linearfrompdeF(polyp, poly q); p/q fraction (ring being r= 0, (x,y,c), dp)
RETURN:  matrix, whose rank give reducible fibers of p-cq, based on Ruppert-Gao pde
EXAMPLE: example linearfrompdeR; shows an example
"
{

  def bring = basering;

  // Partial degrees
  intvec v= 1,0,0;
  if (deg(p,v)>deg(q,v)) {int m = deg(p,v);} else {int m = deg(q,v);};
  intvec w = 0,1,1;
  if (deg(p,w)>deg(q,w)) {int n = deg(p,w);} else {int n = deg(q,w);};


  // New ring with unknowns
  def bring = basering;      // user's ring, of type r=0, (x,y,c), dp;
  poly theminpoly = minpoly;
  string thechar = charstr(bring);
  string thevar = varstr(bring);
  string rstr = "ring rr = (" +  thechar + "), (aa(1..m*(n+1)),bb(1..(m+1)*n),"+ thevar  +"), dp;";  
  execute(rstr);
  minpoly = number(fetch(bring,theminpoly));

  poly pp = imap(bring,p);       // p the polynomial in (x,y)
  poly qq = imap(bring,q);       // q the polynomial in (x,y)

  poly varx = var(2*m*n+m+n+1);  // x variable
  poly vary = var(2*m*n+m+n+2);  // y variable
  poly varc = var(2*m*n+m+n+3);  // c parameter


  /// g;
  poly g =0; 
  int i,j;
  for (i=1; i <= m; i++)
  { for(j=1; j<= n+1; j++)
    { g = g + aa(i+(j-1)*m)*varx^(i-1)*vary^(j-1);    }
  };

  /// h;
  poly h =0; 
  for (i=1; i <= m+1; i++)
  { for(j=1; j<= n; j++)
    { h = h + bb(i+(j-1)*(m+1))*varx^(i-1)*vary^(j-1); }
  };

  /// PDE Equation;
  poly E = (pp-varc*qq)*(diff(g,vary)-diff(h,varx))+h*diff(pp-varc*qq,varx)-g*diff(pp-varc*qq,vary);

  // Make a linear system in aa(1),..., bb(m*n)
  // Note : not linear in c !!
  matrix C = coeffs(E,vary);
  int degy = size(C);
  ideal EE;
  for (i=1; i<= degy; i++)
  { EE = EE + ideal(coeffs(C[i,1],varx)); }


  // Make a matrix
  int ip = size(EE);
  int iq = 2*m*n+m+n;
  matrix M[ip][iq];

  for (i=1; i<=ip; i++)
  {
    for (j=1; j<= m*(n+1); j++)
    {
      C = coef(EE[i],aa(j));
      if (C[1,1]!=1) {M[i,j] = C[2,1];};
    };
    for (j=m*(n+1)+1; j<= iq; j++)
    {
      C = coef(EE[i],bb(j-m*(n+1)));
      if (C[1,1]!=1) {M[i,j] = C[2,1]; };
    };
  }

  // Back to user's ring
  setring bring;
  matrix MM;
  MM  = imap(rr,M);
  return(MM);

}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0, (x,y,c), dp;
   poly p = x3-y3; poly q = 1-y^3;
   print(linearfrompdeF(f));
}



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

proc reducibleF(poly p, poly q);  main procedure; give reducible fibers as the zero of a polynomial
"USAGE:   reducibleF(p,q); p,q poly (ring being r= 0, (x,y,c), dp)
RETURN:  poly, whose zeroes c are corresponds to reducible fibers (p-cq=0)
EXAMPLE: example reducibleF; shows an example
"
{

//  checkhypo(f);
  matrix M;
  M  = linearfrompdeF(p,q);
  M = gauss_col(M);
  M = clean_gaussmatrix(M);
  ideal I;
print(M);
  I = findminors(M);
  poly g = I[1];
  return(g); 
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0, (x,y,c), dp;
  poly p = x3-y3;
  poly q = 1-y3;
  poly g = reducible(p,q);
  "Reducible fiber are given by the zeroes of :", g;
  "After factorization, reducible fibers are :"; 
  factorize(g,1);
}