// $Id: elim.lib,v 1.23 2008-10-06 17:04:26 Singular Exp $
// (GMG, last modified 22.06.96)
///////////////////////////////////////////////////////////////////////////////
version="$Id: elim.lib,v 1.23 2008-10-06 17:04:26 Singular Exp $";
category="Commutative Algebra";
info="
LIBRARY:  elim.lib      Elimination, Saturation and Blowing up

PROCEDURES:
 blowup0(j[,s1,s2]);    create presentation of blownup ring of ideal j
 elim(id,v);            eliminate variables in v from id (ideal/module)
 elim1(id,p);           p=product of vars to be eliminated from id
 nselect(id,v);         select generators not containing variables given by v
 sat(id,j);             saturated quotient of ideal/module id by ideal j
 select(id,v);          select generators containing all variables given by v
 select1(id,v);         select generators containing one variable given by v

           (parameters in square brackets [] are optional)
";

LIB "inout.lib";
LIB "general.lib";
LIB "poly.lib";

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

proc blowup0 (ideal J,ideal C, list #)
"USAGE:   blowup0(J,C [,W]); J,C,W ideals
@*       C = ideal of center of blowup, J = ideal to be blown up,
         W = ideal of ambient space
ASSUME:  inclusion of ideals : W in J, J in C.
         If not, the procedure replaces J by J+W and C by C+J+W
RETURN:  a ring, say B, containing the ideals C,J,W and the ideals
@*         - bR (ideal defining the blown up basering)
@*         - aS (ideal of blown up ambient space)
@*         - eD (ideal of exceptional divisor)
@*         - tT (ideal of total transform)
@*         - sT (ideal of strict transform)
@*         - bM (ideal of the blowup map from basering to B)
@*       such that B/bR is isomorphic to the blowup ring BC.
PURPOSE: compute the projective blowup of the basering in the center C, the
         exceptional locus, the total and strict tranform of J,
         and the blowup map.
         The projective blowup is a presentation of the blowup ring
         BC = R[C] = R + t*C + t^2*C^2 + ... (also called Rees ring) of the
         ideal C in the ring basering R.
THEORY:  If basering = K[x1,...,xn] and C = <f1,...,fk> then let
         B = K[x1,...,xn,y1,...,yk] and aS the preimage in B of W
         under the map B -> K[x1,...,xn,t], xi -> xi, yi -> t*fi.
         aS is homogeneous in the variables yi and defines a variety
         Z=V(aS) in  A^n x P^(k-1), the ambient space of the blowup of V(W).
         The projection Z -> A^n is an isomorphism outside the preimage
         of the center V(C) in A^n and is called the blowup of the center.
         The preimage of V(C) is called the exceptional set, the preimage of
         V(J) is called the total transform of V(J). The strict transform
         is the closure of (total transform - exceptional set).
@*       If C = <x1,...,xn> then aS = <yi*xj - yj*xi | i,j=1,...,n>
         and Z is the blowup of A^n in 0, the exceptional set is P^(k-1).
NOTE:    The procedure creates a new ring with variables y(1..k) and x(1..n)
         where n=nvars(basering) and k=ncols(C). The ordering is a block
         ordering where the x-block has the ordering of the basering and
         the y-block has ordering dp if C is not homogeneous
         resp. the weighted ordering wp(b1,...bk) if C is homogeneous
         with deg(C[i])=bi.
SEE ALSO:blowUp
EXAMPLE: example blowup0; shows examples
"{
   def br = basering;
   list l = ringlist(br);
   int n,k,i = nvars(br),ncols(C),0;
   ideal W;
   if (size(#) !=0)
   { W = #[1];}
   J = J,W;
   //J = interred(J+W);
//------------------------- create rings for blowup ------------------------
//Create rings tr = K[x(1),...,x(n),t] and nr = K[x(1),...,x(n),y(1),...,y(k)]
//and map Bl: nr --> tr, x(i)->x(i), y(i)->t*fi.
//Let ord be the ordering of the basering.
//We change the ringlist l by changing l[2] and l[3]
//For K[t,x(1),...,x(n),t]
// - l[2]: the variables to x(1),...,x(n),t
// - l[3]: the ordering to a block ordering (ord,dp(1))
//For K[x(1),...,x(n),y(1),...,y(k)]
// - l[2]: the variables to x(1),...,x(n),y(1),...,y(k),
// - l[3]: the ordering to a block ordering (ord,dp) if C is
//         not homogeneous or to (ord,wp(b1,...bk),ord) if C is
//         homogeneous with deg(C[i])=bi;

//--------------- create tr = K[x(1),...,x(n),t] ---------------------------
   int s = size(l[3]);
   for ( i=1; i<=n; i++)
   {
      l[2][i]="x("+string(i)+")";
   }
   l[2]=insert(l[2],"t",n);
   l[3]=insert(l[3],list("dp",1),s-1);
   def tr = ring(l);

//--------------- create nr = K[x(1),...,x(n),y(1),...,y(k)] ---------------
   l[2]=delete(l[2],n+1);
   l[3]=delete(l[3],s);
   for ( i=1; i<=k; i++)
   {
      l[2][n+i]="y("+string(i)+")";
   }

   //---- change l[3]:
   l[3][s+1] = l[3][s];         // save the module ordering of the basering
   intvec w;
   w[k]=0; w=w+1;
   intvec v;                    // containing the weights for the varibale
   if( homog(C) )
   {
      for( i=1; i<=k; i++)
      {
         v[i]=deg(C[i]);
      }
      if (v != w)
      {
         l[3][s]=list("wp",v);
      }
      else
      {
         l[3][s]=list("dp",v);
      }
   }
   else
   {
      for( i=1; i<=k; i++)
      {
         v[i]=1;
      }
      l[3][s]=list("dp",v);
   }
   def nr = ring(l);

//-------- create blowup map Bl: nr --> tr, x(i)->x(i), y(i)->t*fi ---------
   setring tr;
   ideal C = fetch(br,C);
   ideal bl = x(1..n);
   for( i=1; i<=k; i++) { bl = bl,t*C[i]; }
   map Bl = nr,bl;
   ideal Z;
//------------------ compute blown up objects and return  -------------------
   setring nr;
   ideal bR = preimage(tr,Bl,Z);   //ideal of blown up affine space A^n
   ideal C = fetch(br,C);
   ideal J = fetch(br,J);
   ideal W = fetch(br,W);
   ideal aS = interred(bR+W);                //ideal of ambient space
   ideal tT = interred(J+bR+W);              //ideal of total transform
   ideal eD = interred(C+J+bR+W);            //ideal of exceptional divisor
   ideal sT = sat(tT,C)[1];       //ideal of strict transform
   ideal bM = x(1..n);            //ideal of blowup map br --> nr

   export(bR,C,J,W,aS,tT,eD,sT,bM);
   return(nr);
}
example
{ "EXAMPLE:"; echo = 2;
   ring r  = 0,(x,y),dp;
   poly f  = x2+y3;
   ideal C = x,y;           //center of blowup
   def B1 = blowup0(f,C);
   setring B1;
   aS;                      //ideal of blown up ambient space
   tT;                      //ideal of total transform of f
   sT;                      //ideal of strict transform of f
   eD;                      //ideal of exceptional divisor
   bM;                      //ideal of blowup map r --> B1

   ring R  = 0,(x,y,z),ds;
   poly f  = y2+x3+z5;
   ideal C = y2,x,z;
   ideal W = z-x;
   def B2 = blowup0(f,C,W);
   setring B2;
   B2;                       //weighted ordering
   bR;                       //ideal of blown up R
   aS;                       //ideal of blown up R/W
   sT;                       //strict transform of f
   eD;                       //ideal of exceptional divisor
   //Note that the different affine charts are {y(i)=1}
 }
///////////////////////////////////////////////////////////////////////////////


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

proc elim (id, intvec va)
"USAGE:   elim(id,v);  id ideal/module, v intvec
RETURNS: ideal/module obtained from id by eliminating variables in v
NOTE:    no special monomial ordering is required, result is a SB with
         respect to ordering dp (resp. ls) if the first var not to be
         eliminated belongs to a -p (resp. -s) blockordering
         This proc uses 'execute' or calls a procedure using 'execute'.
SEE ALSO: elim1, eliminate
EXAMPLE: example elim; shows examples
"
{
//---- get variables to be eliminated and create string for new ordering ------
   int ii; poly vars=1;
   for( ii=1; ii<=size(va); ii++ ) { vars=vars*var(va[ii]); }
   if(  attrib(basering,"global")) { string ordering = "),dp;"; }
   else { string ordering = "),ls;"; }
   string mpoly=string(minpoly);
//-------------- create new ring and map objects to new ring ------------------
   def br = basering;
   string str = "ring @newr = ("+charstr(br)+"),("+varstr(br)+ordering;
   execute(str);
   if (mpoly!="0") { execute("minpoly="+mpoly+";"); }
   def i = imap(br,id);
   poly vars = imap(br,vars);
//---------- now eliminate in new ring and map back to old ring ---------------
   i = eliminate(i,vars);
   setring br;
   return(imap(@newr,i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,u,v,w),dp;
   ideal i=x-u,y-u2,w-u3,v-x+y3;
   elim(i,3..4);
   module m=i*gen(1)+i*gen(2);
   m=elim(m,3..4);show(m);
}
///////////////////////////////////////////////////////////////////////////////

proc elim1 (id, poly vars)
"USAGE:   elim1(id,p); id ideal/module, p product of vars to be eliminated
RETURN:  ideal/module obtained from id by eliminating vars occuring in poly
NOTE:    no special monomial ordering is required, result is a SB with
         respect to ordering dp (resp. ls) if the first var not to be
         eliminated belongs to a -p (resp. -s) blockordering
         This proc uses 'execute' or calls a procedure using 'execute'.
SEE ALSO: elim, eliminate
EXAMPLE: example elim1; shows examples
"
{
//---- get variables to be eliminated and create string for new ordering ------
   int ii;
   for( ii=1; ii<=nvars(basering); ii++ )
   {
      if( vars/var(ii)==0 ) { poly p = 1+var(ii); break;}
   }
   if( ord(p)==0 ) { string ordering = "),ls;"; }
   if( ord(p)>0 ) { string ordering = "),dp;"; }
//-------------- create new ring and map objects to new ring ------------------
   def br = basering;
   string str = "ring @newr = ("+charstr(br)+"),("+varstr(br)+ordering;
   execute(str);
   def id = fetch(br,id);
   poly vars = fetch(br,vars);
//---------- now eliminate in new ring and map back to old ring ---------------
   id = eliminate(id,vars);
   setring br;
   return(imap(@newr,id));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,t,s,z),dp;
   ideal i=x-t,y-t2,z-t3,s-x+y3;
   elim1(i,ts);
   module m=i*gen(1)+i*gen(2);
   m=elim1(m,st); show(m);
}
///////////////////////////////////////////////////////////////////////////////

proc nselect (id_inp, intvec v)
"USAGE:   nselect(id,v); id = module or ideal, v = intvec
RETURN:  generators of id not containing the variables with index an entry of v
SEE ALSO: select, select1
EXAMPLE: example nselect; shows examples
"
{
  if (typeof(id_inp)!="ideal") { module id=module(id_inp); }
  else                         { ideal id=id_inp;         }
  int j,k;
  int n,m = size(v), ncols(id);
  //listvar(id);
  for( k=1; k<=m; k++ )
  {  
    for( j=1; j<=n; j++ )
    {  
      if( size(id[k]/var(v[j]))!=0 ) 
      { 
         id[k]=0; break; 
      }
    }
  }
  //if(typeof(id)=="ideal") { return(simplify(id,2)); }
  //return(simplify(module(id),2));
  id=(simplify(module(id),2));
  //listvar(id);
  return(id);
  //return(simplify(id,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,t,s,z),(c,dp);
   ideal i=x-y,y-z2,z-t3,s-x+y3;
   nselect(i,3);
   module m=i*(gen(1)+gen(2));
   show(m);
   show(nselect(m,3..4));
}
///////////////////////////////////////////////////////////////////////////////

proc sat (id, ideal j)
"USAGE:   sat(id,j);  id=ideal/module, j=ideal
RETURN:  list of an ideal/module [1] and an integer [2]:
         [1] = saturation of id with respect to j (= union_(k=1...) of id:j^k)
         [2] = saturation exponent (= min( k | id:j^k = id:j^(k+1) ))
NOTE:    [1] is a standard basis in the basering
DISPLAY: saturation exponent during computation if printlevel >=1
EXAMPLE: example sat; shows an example
"{
   int ii,kk;
   def i=id;
   id=std(id);
   int p = printlevel-voice+3;  // p=printlevel+1 (default: p=1)
   while( ii<=size(i) )
   {
      dbprint(p-1,"// compute quotient "+string(kk+1));
      i=quotient(id,j);
      for( ii=1; ii<=size(i); ii++ )
      {
         if( reduce(i[ii],id,1)!=0 ) break;
      }
      id=std(i); kk++;
   }
   dbprint(p-1,"// saturation becomes stable after "+string(kk-1)+" iteration(s)","");
   list L = id,kk-1;
   return (L);
}
example
{ "EXAMPLE:"; echo = 2;
   int p      = printlevel;
   ring r     = 2,(x,y,z),dp;
   poly F     = x5+y5+(x-y)^2*xyz;
   ideal j    = jacob(F);
   sat(j,maxideal(1));
   printlevel = 2;
   sat(j,maxideal(2));
   printlevel = p;
}
///////////////////////////////////////////////////////////////////////////////

proc select (id, intvec v)
"USAGE:   select(id,n[,m]); id = ideal/module, v= intvec
RETURN:  generators of id containing all variables with index an entry of v
NOTE:    use 'select1' for selecting generators containing at least one of the
         variables with index an entry of v
SEE ALSO: select1, nselect
EXAMPLE: example select; shows examples
"{
   int j,k;
   int n,m = size(v), ncols(id);
   for( k=1; k<=m; k++ )
   {  
      for( j=1; j<=n; j++ )
      {   
         if( size(id[k]/var(v[j]))==0) 
         { 
            id[k]=0; break; 
         }
      }
   }
   return(simplify(id,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,t,s,z),(c,dp);
   ideal i=x-y,y-z2,z-t3,s-x+y3;
   ideal j=select(i,1);
   j;
   module m=i*(gen(1)+gen(2));
   m;
   select(m,1..2);
}
///////////////////////////////////////////////////////////////////////////////

proc select1 (id, intvec v)
"USAGE:   select1(id,v); id = ideal/module, v = intvec
RETURN:  generators of id containing at least one of the variables with
         index an entry of v
NOTE:    use 'select' for selecting generators containing all the
         variables with index an entry of v
SEE ALSO: select, nselect
EXAMPLE: example select1; shows examples
"{
   int j,k;
   int n,m = size(v), ncols(id);
   execute (typeof(id)+" I;");
   for( k=1; k<=m; k++ )
   {  for( j=1; j<=n; j++ )
      {
         if( size(subst(id[k],var(v[j]),0)) != size(id[k]) )
         {  
            I = I,id[k]; break; 
         }
      }
   }
   return(simplify(I,2));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y,t,s,z),(c,dp);
   ideal i=x-y,y-z2,z-t3,s-x+y3;
   ideal j=select1(i,1);
   j;
   module m=i*(gen(1)+gen(2));
   m;
   select1(m,1..2);
}
///////////////////////////////////////////////////////////////////////////////