// defpol01.lib, v 0.01 2003/04/11
//(Arnaud Bodin, last modified 2003/06/11)
///////////////////////////////////////////////////////////////////////////////
version="defpol01.lib, v 0.01 2003/04/11";
category="Singularities";   // line may be commented for some old version of Singular
info="
LIBRARY:  defpol01.lib      Critical parameters for defomations of polynomials
AUTHORS:  Arnaud Bodin : email : Arnaud.Bodin@agat.univ-lille1.fr

PROCEDURES:
  proc parCritDegree(poly F);        detects the change in the degree
  proc parCritMu (poly F);           detects the change in the affine Milnor number
  proc parCritBaff (poly F);         detects the change in the affine critical values
  proc parCritInfx (poly F, int k);  detects the change at infinity (lambda or Binf) in a chart
  proc parCritInf (poly F);          detects the change at infinity (lambda or Binf)
  proc parCritB (ideal I1, ideal I2) detects the change in B = Baff union Binf
  proc parCritIdeal (ideal I);       detects the change in a braid
";  

LIB "elim.lib";
LIB "primdec.lib";

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



//////////////////////////////////////////
///// Critical parameters for degree /////
//////////////////////////////////////////

// Critical parameters for the degree
// Input  : a polynomial
// Output : polynomial whose roots detect the change in the degree

proc parCritDegree(poly F)
"USAGE:  parCritDegree(F); F = poly
RETURN:  polynomial whose roots detect the change in the degree
EXAMPLE: example parCritDegree; shows an example
"
{
  // Change of ring
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  ring r = 0, (x(1..n),s,z,t), dp;
  poly f = fetch(bring,F);

  // Now s, t are parameters 
  // In fact there is no t, no z, only x(i) and s

  ring rr=(0,s,t), (x(1..n),z), dp;
  poly ff = imap(r,f);
  
  poly ffd = ff - jet(ff,deg(ff)-1);

  poly DD = factorize(ffd)[1][1];  

  
  setring r;
  poly D = imap(rr,DD);

  setring bring;
  poly DDD = fetch(r,D);
  return(polyred(DDD));

}
example
{ "EXAMPLE:"; echo = 2;
   ring myr   = 0, (x,y,s), dp;
   poly F  = (s^2-1)*x^3+(s-1)*x*y^2+x*y;
   parCritDegree(F);
}
///////////////////////////////////////////////////////////////////////////////




//////////////////////////////////////////
//////  Critical parameters for mu ///////
//////////////////////////////////////////

// Critical parameters for the affine Milnor number
// Input  : a polynomial
// Output : polynomial whose roots are potential
//          critical parameters for mu

proc parCritMu (poly F)
"USAGE:  parCritMu(F); F = poly
RETURN:  polynomial whose roots detect the change in mu
EXAMPLE: example parCritMu; shows an example
"
{
  // Change of ring
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  ring r = 0, (x(1..n),s,z,t), dp;
  poly f = fetch(bring,F);

  // Now s is a parameter, there is no t
  ring rr=(0,s), (x(1..n),t,z), dp;
  poly ff = imap(r,f);

  // Jacobian ideal
  ideal J = jacob(ff);

  // Homogenization
  ideal JJH = homog(J,z);

  // Now s is a variable
  setring r;

  // Ideal JJH in the new ring
  ideal JH = imap(rr,JJH);

  // At infinity
  ideal Jinf = JH, z;

  // JH\Jinf
  ideal Jaff = sat(JH,Jinf)[1];

  // At infinity
  ideal Jaffinf = Jaff,z;

  // Projectivization
  ideal Pr = x(1..n),z;

  // Elimination of the solution x(1)=x(2)=...=x(n)=z=0
  ideal JJaffinf = sat(Jaffinf,Pr)[1];

  // Projection to s 
  poly V = 1;
  for (int i = 1; i<=n; i=i+1) {V = V * var(i);};
  ideal ICrit = eliminate(JJaffinf,V*z);
  poly PCritMu = ICrit[1];

  // Answer in user's ring
  setring bring;
  poly PPPCritMu = fetch(r,PCritMu);
  return(polyred(PPPCritMu));
}
example
{ "EXAMPLE:"; echo = 2;
   ring myr   = 0, (x,y,s),dp;
   poly P  =  x*(x-s)*y;
   parCritMu(P);
}
///////////////////////////////////////////////////////////////////////////////



//////////////////////////////////////////
////// Critical parameters for Baff //////
//////////////////////////////////////////

// Critical parameters for the affine critical values
// Input  : a polynomial
// Output : [1] polynomial whose roots are potential
//              critical parameters for Baff
//          [2] the braid of affine critical values

proc parCritBaff (poly F)
"USAGE:  parCritBaff(F); F = poly
RETURN:  [1] polynomial whose roots detect the change in Baff, 
         [2] the braid of affine critical values (in the 2 first variables)
EXAMPLE: example parCritBaff; shows an example
"
{
  // Change of ring
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  ring r = 0, (x(1..n),s,z,t), dp;
  poly f = fetch(bring,F);

  // Now s is a parameter
  ring rr=(0,s), (x(1..n),z,t), dp;
  poly ff = imap(r,f);

  // Jacobian ideal + graph
  ideal J = jacob(ff), ff-t;

  // Now s is a variable
  setring r;
  ideal JJ = imap(rr,J);

  // Projection
  poly V = 1;
  for (int i = 1; i<=n; i=i+1) {V = V * var(i);};
  ideal K = eliminate(JJ,V);

  // Now consider the bifuraction set
  ring rpar = 0, (t,s), dp;
  ideal Kpar =  imap(r,K);

  poly Ppar = parCritIdeal(Kpar);

  // Back to the user's ring

  setring r;
  poly PCritBaff = fetch(rpar,Ppar);

  setring bring;
  poly PPCritBaff = imap(r,PCritBaff);
  ideal KK = fetch(rpar,Kpar);

  poly PPPCritBaff = polyred(PPCritBaff);
  return(list(PPPCritBaff,KK));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r   = 0, (x,y,s),dp;
   poly P  =  x*(x-s)*y;
   parCritBaff(P)[1];
}
///////////////////////////////////////////////////////////////////////////////



//////////////////////////////////////////
///// Critical parameters at infinity ////
//////////////////////////////////////////

// Critical parameters  at infinity
// Input  : a polynomial, a chart k : x(k)=1
// Output : [1] polynomial whose roots are potential
//              critical parameters for Binf or lambda
//          [2] the braid of critical values at infinity

proc parCritInfx (poly F, int k)
"USAGE:  parCritInfx(F,k); F = poly, k integer
RETURN:  [1] polynomial whose roots detect the change at infinity in the chart xk=1, 
         [2] the braid of critical values at infinity (in the 2 first variables)
EXAMPLE: example parCritInfx; shows an example
"
{
  // Change of ring
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  ring r = 0, (x(1..n),s,z,t), dp;
  poly f = fetch(bring,F);

 
  // Now s is a parameter (for homogeneization)
  ring rr=(0,s), (x(1..n),z,t), dp;
  poly ff = imap(r,f);
  
  // X = (fH=0)
  poly ffH = homog(ff,z)-t*z^deg(ff);   // t is a parameter not a variable


  // Polar Curve with respect to (xk=0)
  list L = list();
  for (int i = 1; i<k; i=i+1) {L[i] = diff(ff,x(i));}
  for (i = k+1; i<=n; i=i+1) {L[i-1] = diff(ff,x(i));}
  ideal P = L[1..n-1];   // Affine polar curve
  ideal PPH = homog(P,z);


  // Now s is a variable
  setring r;

  // X = (fH=0)
  poly fH = imap(rr,ffH);
  ideal X = fH;
  ideal Xinf = z,fH;

  // Polar Curve with respect to (xk=0)
  ideal PH = imap(rr,PPH);

  // Cloture of the intersection of polar curve and X\Xinf and (xk<>0)
  // We localise outside xk=0 projectively, i.e. we set xk=1 
  ideal Cbar= x(k)-1,X,PH;

  // We exclude the part of the polar curve that comes from Xinf
  ideal C = sat(Cbar,Xinf)[1];  
  
  // We calculate the intersection of the "real" polar curve C
  // with the hyperplane at infinity
  ideal Cinf = z,x(k)-1,C;

  // We go in x1=0, x2=0, x{k-1}=0, xk=1 to avoid redundancy
  list LL = list();
  for (i = 1; i<k; i=i+1) {LL[i] = x(i);}
  if (k==1) {ideal J = 0;} else {ideal J = LL[1..k-1];};
 
  ideal Cred = Cinf, J; 

  // We get the critical values at infinity by elimination
  poly V = 1;
  for (i = 1; i<=n; i=i+1) {V = V * var(i);};
  ideal K = eliminate(Cred,V*z);

  // Critical parameters for lambda
  ideal CK = C, K; 
  ideal CinfK = Cinf, K;
  ideal D = sat(CK,CinfK)[1];
  ideal Dinf = z, x(k)-1, J, D;
  ideal L =  eliminate(Dinf,V*z*t);
  poly PCritLambda = L[1];

  // Critical parameters for Binf
  ring rpar = 0, (t,s), dp;
  ideal Kpar =  imap(r,K);
  poly Ppar = parCritIdeal(Kpar);

  // Back to the user's ring
  setring r;
  poly PCritBinf = imap(rpar,Ppar);

  setring bring;
  poly PPCritBinf = fetch(r,PCritBinf);
  poly PPCritLambda = fetch(r,PCritLambda);
  poly PCritInf = PPCritLambda*PPCritBinf;
 
  ideal KK = fetch(rpar,Kpar);

  return(list(polyred(PCritInf),KK));
}
example
{ "EXAMPLE:"; echo = 2;
   ring myr   = 0, (x,y,s),dp;
   poly P = (x-s)*(xy+1);
   parCritInfx(P,1)[1];
   parCritInfx(P,2)[1];
}
///////////////////////////////////////////////////////////////////////////////




//////////////////////////////////////////
///// Critical parameters at infinity ////
//////////////////////////////////////////

// Critical parameters  at infinity
// Input  : a polynomial
// Output : [1] polynomial whose roots are potential
//              critical parameters for Binf or lambda
//          [2] the braid of critical values at infinity 
proc parCritInf (poly F)
"USAGE:  parCritInf(F); F = poly
RETURN:  [1] polynomial whose roots detect the change at infinity,  
         [2] the braid of critical values at infinity (in the 2 first variables)
EXAMPLE: example parCritInf; shows an example
"
{
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  list Linf;
  poly Pinf = 1;
  ideal Iinf = 1;

  for (int i=1; i<=n; i++) 
  {
    Linf = parCritInfx(F,i);
    Pinf = Pinf * Linf[1];
    Iinf = Iinf * Linf[2];
  };

  return(list(polyred(Pinf),Iinf));
}
example
{ "EXAMPLE:"; echo = 2;
   ring myr = 0, (x,y,s),dp;
   poly P = (x-s)*(xy+1);
   parCritInf(P)[1];
}
///////////////////////////////////////////////////////////////////////////////



/////////////////////////////////////////////////////
///// Critical parameters of the critical values ////
/////////////////////////////////////////////////////

// Critical parameters  at infinity
// Input  : the critical values in affine space and at infinity
// Output :  polynomial whose roots are potential
//              critical parameters for B = Baff union Binf
proc parCritB (ideal I1, ideal I2)
"USAGE:  parCritB(I1,I2); I1,I2 = ideals
RETURN:  polynomial whose roots detect the change in B
EXAMPLE: example parCritB; shows an example
"
{
  // Here the braids are in the two first variables(x(1), x(2)) (yes, it is quiet strange !)

  // Change of ring
  def bring = basering;      // user's ring
  int n = nvars(bring)-1;
  ring r = 0, (x(1..n),s,z,t), dp;
  ideal II1 = fetch(bring,I1);
  ideal II2 = fetch(bring,I2);
  // Now the braids are in (x(1), x(2)) (yes, it is quiet strange !)

  ring rpar = 0, (t,s), dp;
  ideal K1 =  fetch(r,II1);
  ideal K2 =  fetch(r,II2);
  //Now the braids are in (t,s) as it should be

  poly ParB = parCritIdeal(K1*K2);

  // Back to the user's ring
  setring r;
  poly PrB = imap(rpar,ParB);

  setring bring;
  poly PB = fetch(r,PrB);

  return(polyred(PB));
}
example
{ "EXAMPLE:"; echo = 2;
   ring rpar = 0, (x,y,s), dp;
   ideal I1 = x-y;
   ideal I2 = x+y;
   parCritB(I1,I2);
}
///////////////////////////////////////////////////////////////////////////////




//////////////////////////////////////////
////////    Critical parameters //////////
//////////////////////////////////////////

// Critical parameters
// Input  : a polynomial
// Output : polynomial whose roots detect the change in the topology
// Warning some parameter may not be critical !

proc parCrit(poly F)
"USAGE:  parCrit(F); F = poly
RETURN:  polynomial whose roots detect the change in the topology
EXAMPLE: example parCrit; shows an example
"
{  
  def bring = basering;
  int m = nvars(bring);

  "Family  : " + string(F);
  "Parameter : " + string(var(m));

  poly PDeg = parCritDegree(F);
  poly PMu = parCritMu(F);
  list LBaff = parCritBaff(F);
  poly PBaff = LBaff[1];
  list LInf = parCritInf(F);
  poly PInf =  LInf[1];
  poly PB = parCritB(LBaff[2],LInf[2]);
  poly PCrit = PDeg*PMu*PBaff*PInf*PB;

  PCrit = polyred(PCrit);

  "Critical parameters are included in the roots of " + string(PCrit);

  return(PCrit);
}
example
{ "EXAMPLE:"; echo = 2;
   ring myr   = 0, (x,y,s), dp;
   poly F  = (s^2-1)*x^3+(s-1)*x*y^2+(s^3+2)*x*y;
   poly C = parCrit(F);
   factorize (C,1);
}
///////////////////////////////////////////////////////////////////////////////









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

//////////////////////////////////////////
/////////////     Tools      /////////////
//////////////////////////////////////////


//////////////////////////////////////////
//// Critical parameters for an ideal ////
//////////////////////////////////////////

// Critical parameters for an ideal
// Input  : a ideal I of K[t,s]
// Output : critical parameters s where
//          I(s) is not a braid

proc parCritIdeal (ideal I)
"USAGE:  parCritIdeal(I); I = ideal
RETURN:  polynomial whose roots detect the change in I(s)
EXAMPLE: example parCritIdeal; shows an example
"
{
  def userr = basering;      // user's ring

// Starting ring is like
// ring rpar=0, (t,s), dp;
// t = var(1); s=var(2);

  poly PCrit = 1;
  ideal K;

 if (dim(std(I)) > -1)
 {

  // Decomposition into prime ideals
  list J = minAssChar(I);

  // Remember the dimensions
  list thedim;
  for (int i=1; i<= size(J); i++)
  {
     thedim[i] = dim(std(J[i]));
  }
  
  // 1 // Isolated points
  for (i=1; i<= size(J); i++)
  {
     if (thedim[i]==0)
     {
        K = eliminate(J[i],var(1));
        PCrit = PCrit * K[1];
     } 
  }
  // One dimensionnal components
  // Now s is a parameter
  ring rrpar = (0,var(2)), var(1), dp;
  list JJ = imap(userr,J);
  poly PPCrit = 1;
  poly Q;
  poly totQ=1;
  for (i=1; i<= size(JJ); i++)
  {
     if (thedim[i]==1) 
     // so JJ[i] has only one equation
     {
       Q = JJ[i][1];  // the irreducible component
       totQ = totQ * Q;

       // 2 // Vertical lines (s=cst)
       if (deg(Q)==0)
       {
         PPCrit = PPCrit * Q;
       }
       if (deg(Q)>0)
       {  
       // 3 // Asymptotic parameters
       PPCrit = PPCrit * leadcoef(Q);
       }
     }
  }

  // 4 // Singular parameters
  poly R = resultant(totQ,diff(totQ,var(1)),var(1));
  if (deg(R)>-1)
  {
    PPCrit = PPCrit * leadcoef(R);
  }

  setring userr;
  PCrit = PCrit * imap(rrpar,PPCrit);
 }  

  return(polyred(PCrit));
}
example
{ "EXAMPLE:"; echo = 2; 
   ring rpar = 0, (t,s),dp;
   ideal I = t,s;
   ideal J= t-1, s-2;
   ideal K = s2+1;
   ideal L = (t-s-1)*(t+s+1);
   factorize(parCritIdeal(I*J*K*L),1);
}
///////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////
//// Reduce a polynomial ////
//////////////////////////////////////////

static proc polyred(poly B)
{ ideal L = factorize(B,1);
  poly Bred = 1;
  for (int i=1; i<= size(L); i=i+1) 
   {Bred = Bred * L[i];};
  return(Bred);
}


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