source: git/Singular/LIB/elim.lib @ 7f3ad4

// $Id: elim.lib,v 1.26 2009-01-14 16:07:04 Singular Exp $
///////////////////////////////////////////////////////////////////////////////
version="$Id: elim.lib,v 1.26 2009-01-14 16:07:04 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
 elimRing(p);          create ring with block ordering for elimating vars in p
 elim(id,..);          variables .. eliminated from id (ideal/module)
 elim1(id,p);          p=product of vars to be eliminated from id
 elim2(id,..);         variables .. eliminated from id (ideal/module)
 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 elimRing ( poly vars, list #)
"USAGE:   elimRing(vars [,w]); vars = product of variables to be eliminated
         (type poly), w = intvec (specifying weights for all variables)
RETURN:  a ring, say R, s.t. the monomial ordering of R has 2 blocks.
         The first block corresponds to the (given) variables to be eliminated
         and has ordering dp if these variables have weight all 1; if w is
         given or if not all variables in vars have weight 1 the ordering is
         wp(w1) where w1 is the intvec of weights of the variables to be
         eliminated.
         The second block corresponds to variables not to be eliminated.
@*       If the first variable not to be eliminated is global (i.e. > 1),
         resp. local (i.e. < 1), the second block has ordering dp, resp. ds,
         (or wp(w2), resp. ws(w2), where w2 is the intvec of weights of the
         variables not to be eliminated).
@*       If the basering is a quotient ring P/Q, then R a quotient ring
         with Q replaced by a standard basis of Q w.r.t. the new ordering
         (parameters are not touched).
NOTE:    The ordering in R is an elimination ordering for the variables
         appearing in vars.
PURPOSE: Prepare a ring for eliminating vars from an ideal/moduel by
         computing a standard basis in R with a fast monomial ordering.
         This procedure is used by the procedure elim.
EXAMPLE: example elimRing; shows an example
"
{
  def BR = basering;
  int nvarBR = nvars(BR);
  list BRlist = ringlist(BR);
  intvec @w;                   //to store weights of all variables
  @w[nvarBR] = 0;
  @w = @w + 1;                 //initialize @w as 1..1
  if (size(#) == 1)
  {
    if ( typeof(#[1]) == "intvec" )
    {
       @w = #[1];              //take the given weights
    }
  }
  else
  {
     @w = ringweights(BR);     //compute the ring weights (proc from ring.lib)
  }

  //--- get variables to be eliminated and ringweights:
  intvec w1,w2;  //for ringweights of first (w1) and second (w2) block
  list v1,v2;    //for variables of first (to be liminated) and second block

  int ii;
  for( ii=1; ii<=nvarBR; ii++ )
  {
     if( vars/var(ii)==0 )    //treat variables not to be eliminated
     {
        w2 = w2,@w[ii];
        v2 = v2+list(string(var(ii)));
        if ( defined(local) == 0 )
        {
           int local = (var(ii) < 1);
         }
     }
     else
     {
        w1 = w1,@w[ii];
        v1 = v1+list(string(var(ii)));
     }
  }

  int l1, l2 = size(w1), size(w2);
  if ( l1 <= 1 )
  {
    ERROR("no elimination ?");
    //return(BR);
  }
  if ( l2 <= 1 )
  {
    ERROR("## elimination of all variables is not possible");
  }

  w1 = w1[2..size(w1)];
  w2 = w2[2..size(w2)];

  //--- put variables to be eliminated in front:
  BRlist[2] = v1 + v2;

  //--- create a block ordering with two blocks and weights:
  int nblock = size(BRlist[3]);      //number of blocks
  list BR3 =  BRlist[3];             //save ordering
  BRlist[3] = list();
  list B3;

  if( w1==1 )
  {
     B3[1] = list("dp", w1);
  }
  else
  {
     B3[1] = list("wp", w1);
  }

  if( w2==1 )
  {
     if ( local==1 )
     {
        B3[2] = list("ds", w2);
     }
     else
     {
        B3[2] = list("dp", w2);
     }
  }
  else
  {
     if ( local==1 )
     {
        B3[2] = list("ws", w2);
     }
     else
     {
        B3[2] = list("wp", w2);
     }
  }

  BRlist[3] = B3;

  //Module ordering stays in front resp. at the end:
  if( BR3[nblock][1] =="c" || BR3[nblock][1] =="C" )
  {
    BRlist[3] = insert(BRlist[3],BR3[nblock],size(B3));
  }
  else
  {
    BRlist[3] = insert(BRlist[3],BR3[1]);
  }

  def eRing = ring(quotientList(BRlist));
  return (eRing);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R = 0,(x,y,z,u,v),(c,lp);
   def P = elimRing(yu);  P;
   intvec w = 1,1,3,4,5;
   elimRing(yu,w);

   ring S =  (0,a),(x,y,z,u,v),ws(1,2,3,4,5);
   minpoly = a2+1;
   qring T = std(ideal(x+y2+v3,(x+v)^2));
   def Q = elimRing(yv);
   setring Q; Q;
}
///////////////////////////////////////////////////////////////////////////////

proc elim (id, list #)
"USAGE:   elim(id,arg[,\"withWeights\"]);  id ideal/module, arg can be either
        an intvec vor a product p of variables (type poly)
RETURN: ideal/module obtained from id by eliminating either the variables
        with indices appearing in v or the variables appearing in p.
        Works also in a qring.
METHOD: elim uses elimRing to create a ring with block ordering with two
        blocks where the first block contains the variables to be eliminated
        and then uses groebner. If the variables in the basering have weights
        these weights are used in elimRing.
@*      If a string \"withWeigts\" as second, optional argument is given,
        Singular computes weights for the variables to make the input as
        homogeneous as possible.
@*      The method is different from that used by eliminate and elim1;
        in some examples elim can be significantly faster.
NOTE:   No special monomial ordering is required, i.e. the ordering can be
        local or mixed. The result is a SB with respect to the ordering of
        the second block used by elimRing. E.g. if the first var not to be
        eliminated is global, resp. local, this ordering is dp, resp. ds
        (or wp, resp. ws, with the given weights for these variables).
        If printlevel > 0 the ring for which the output is a SB is shown.
SEE ALSO: eliminate, elim1
EXAMPLE: example elim; shows an example
"
{
  if (size(#) == 0)
  {
    ERROR("## specify variables to be eliminated");
  }
  int pr = printlevel - voice + 2;   //for ring display if printlevel > 0
  def BR = basering;
//-------------------------------- check input -------------------------------
  poly vars;
  int ii;
  if (size(#) > 0)
  {
    if ( typeof(#[1]) == "poly" )
    {
      vars = #[1];
    }
    if ( typeof(#[1]) == "intvec")
    {
      vars=1;
      for( ii=1; ii<=size(#[1]); ii++ )
      {
        vars=vars*var(#[1][ii]);
      }
    }
  }
  if (size(#) == 2)
  {
    if ( typeof(#[2]) == "string" )
    {
       if ( #[2] == "withWeights" )
       {
         intvec @w = weight(id);
       }
    }
  }

//-------------- create new ring and map objects to new ring ------------------
  if ( defined(@w) )
  {
     def ER = elimRing(vars,@w);
  }
  else
  {
     def ER = elimRing(vars);
  }
  setring ER;
  def id = imap(BR,id);
  poly vars = imap(BR,vars);
//---------- now eliminate in new ring and map back to old ring ---------------
  id = groebner(id);
  id = nselect(id,1..size(ringlist(ER)[3][1][2]));
  if ( pr > 0 )
  {
     "// result is a SB in the following ring:";
     ER;
  }
  setring BR;
  return(imap(ER,id));
}
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);
   elim(i,uv);
   int p = printlevel;
   printlevel = 2;
   elim(i,uv,"withWeights");
   printlevel = p;

   ring S =  (0,a),(x,y,z,u,v),ws(1,2,3,4,5);
   minpoly = a2+1;
   qring T = std(ideal(ax+y2+v3,(x+v)^2));
   ideal i=x-u,y-u2,az-u3,v-x+ay3;
   module m=i*gen(1)+i*gen(2);
   m=elim(m,xy);
   show(m);
}
///////////////////////////////////////////////////////////////////////////////

proc elim2 (id, intvec va)
"USAGE:   elim2(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, elim
EXAMPLE: example elim2; 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;
   elim2(i,3..4);
   module m=i*gen(1)+i*gen(2);
   m=elim2(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
METHOD:  elim1 calls eliminate but in a ring with ordering dp (resp. ls)
         if the first var not to be eliminated belongs to a -p (resp. -s)
         ordering.
NOTE:    no special monomial ordering is required.
         This proc uses 'execute' or calls a procedure using 'execute'.
SEE ALSO: elim, eliminate
EXAMPLE: example elim1; shows examples
"
{
   def br = basering;
   if ( size(ideal(br)) != 0 )
   {
      ERROR ("cannot eliminate in a qring");
   }
//---- 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 ------------------
   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, intvec v)
"USAGE:   nselect(id,v); id = ideal, module or matrix, v = intvec
RETURN:  generators (or columns) of id not containing the variables with index
         an entry of v
SEE ALSO: select, select1
EXAMPLE: example nselect; shows examples
"{
   if (typeof(id)!="ideal")
   {
     if (typeof(id)=="module" || typeof(id)=="matrix")
     {
       module id1 = module(id);
     }
     else
     {
       ERROR("// *** input must be of type ideal or module or matrix");
     }
   }
   else
   {
      ideal id1 = id;
   }
   int j,k;
   int n,m = size(v), ncols(id1);
   for( k=1; k<=m; k++ )
   {
      for( j=1; j<=n; j++ )
      {
        if( size(id1[k]/var(v[j]))!=0 )
        {
           id1[k]=0; break;
        }
      }
   }
   id1=simplify(id1,2);
   if(typeof(id)=="matrix")
   {
      return(matrix(id1));
   }
   return(id1);
}
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));
   m;
   nselect(m,3..4);
   nselect(matrix(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/matrix, v = intvec
RETURN:  generators/columns of id containing all variables with index
         an entry of v
NOTE:    use 'select1' for selecting generators/columns containing at least
         one of the variables with index an entry of v
SEE ALSO: select1, nselect
EXAMPLE: example select; shows examples
"{
   if (typeof(id)!="ideal")
   {
     if (typeof(id)=="module" || typeof(id)=="matrix")
     {
       module id1 = module(id);
     }
     else
     {
       ERROR("// *** input must be of type ideal or module or matrix");
     }
   }
   else
   {
     ideal id1 = id;
   }
   int j,k;
   int n,m = size(v), ncols(id1);
   for( k=1; k<=m; k++ )
   {
      for( j=1; j<=n; j++ )
      {
         if( size(id1[k]/var(v[j]))==0)
         {
            id1[k]=0; break;
         }
      }
   }
   if(typeof(id)=="matrix")
   {
      return(matrix(simplify(id1,2)));
   }
   return(simplify(id1,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);
   select(matrix(m),1..2);
}
///////////////////////////////////////////////////////////////////////////////

proc select1 (id, intvec v)
"USAGE:   select1(id,v); id = ideal/module/matrix, v = intvec
RETURN:  generators/columns of id containing at least one of the variables
         with index an entry of v
NOTE:    use 'select' for selecting generators/columns containing all variables
         with index an entry of v
SEE ALSO: select, nselect
EXAMPLE: example select1; shows examples
"{
   if (typeof(id)!="ideal")
   {
     if (typeof(id)=="module" || typeof(id)=="matrix")
     {
       module id1 = module(id);
       module I;
     }
     else
     {
       ERROR("// *** input must be of type ideal or module or matrix");
     }
   }
   else
   {
     ideal id1 = id;
     ideal I;
   }
   int j,k;
   int n,m = size(v), ncols(id1);
   for( k=1; k<=m; k++ )
   {  for( j=1; j<=n; j++ )
      {
         if( size(subst(id1[k],var(v[j]),0)) != size(id1[k]) )
         {
            I = I,id1[k]; break;
         }
      }
   }
   I=simplify(I,2);
   if(typeof(id)=="matrix")
   {
      return(matrix(I));
   }
   return(I);
}
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);
   select1(matrix(m),1..2);
}
///////////////////////////////////////////////////////////////////////////////