source: git/Singular/LIB/elim.lib @ b0ee97

//////////////////////////////////////////////////////////////////////////////
version="version elim.lib 4.0.0.0 Jun_2013 "; // $Id$
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)          variables .. eliminated from id (different algorithm)
 elim2(id,..)         variables .. eliminated from id (different algorithm)
 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";
LIB "ring.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 minus the 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, blowUp2
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=n; i>=1; 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=k; i>=1; 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=1:k;
   intvec v;                    // containing the weights for the varibale
   if( homog(C) )
   {
      for( i=k; i>=1; i--)
      {
         v[i]=deg(C[i]);
      }
      if (v != w)
      {
         l[3][s]=list("wp",v);
      }
      else
      {
         l[3][s]=list("dp",v);
      }
   }
   else
   {
      v=1:k;
      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,str]); vars = product of variables to be eliminated
         (type poly), w = intvec (specifying weights for all variables),
         str = string either \"a\" or \"b\" (default: w=ringweights, str=\"a\")
RETURN:  a list, say L, with R:=L[1] a ring and L[2] an intvec.
         The ordering in R is an elimination ordering for the variables
         appearing in vars depending on \"a\" resp. \"b\". Let w1 (resp. w2)
         be the intvec of weights of the variables to be eliminated (resp. not
         to be eliminated).
         The monomial ordering of R has always 2 blocks, the first
         block corresponds to the (given) variables to be eliminated.
@*       If str = \"a\" the first block is a(w1,0..0) and the second block is
         wp(w) resp. ws(w) if the first variable not to be eliminated is local.
@*       If str = \"b\" the 1st block has ordering wp(w1) and the 2nd block
         is wp(w2) resp. ws(w2) if the first variable not to be eliminated is
         local.
@*       If the basering is a quotient ring P/Q, then R is also a quotient ring
         with Q replaced by a standard basis of Q w.r.t. the new ordering
         (parameters are not touched).
@*       The intvec L[2] is the intvec of variable weights (or the given w)
         with weights <= 0 replaced by 1.
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);
  //------------------ set resp. compute ring weights ----------------------
  int ii;
  intvec @w;              //to store weights of all variables
  @w[nvarBR] = 0;
  @w = @w + 1;            //initialize @w as 1..1
  string str = "a";       //default for specifying elimination ordering
  if (size(#) == 0)       //default values
  {
     @w = ringweights(BR);     //compute the ring weights (proc from ring.lib)
  }

  if (size(#) == 1)
  {
    if ( typeof(#[1]) == "intvec" )
    {
       @w = #[1];              //take the given weights
    }
    if ( typeof(#[1]) == "string" )
    {
       str = #[1];             //string for specifying elimination ordering
    }
  }

  if (size(#) >= 2)
  {
    if ( typeof(#[1]) == "intvec" and typeof(#[2]) == "string" )
    {
       @w = #[1];              //take the given weights
       str = #[2];             //string for specifying elimination ordering

    }
    if ( typeof(#[1]) == "string" and typeof(#[2]) == "intvec" )
    {
       str = #[1];             //string for specifying elimination ordering
       @w = #[2];              //take the given weights
    }
  }

  for ( ii=1; ii<=size(@w); ii++ )
  {
    if ( @w[ii] <= 0 )
    {
       @w[ii] = 1;             //replace non-positive weights by 1
    }
  }

  //------ get variables to be eliminated together with their weights -------
  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

  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) )
        {
           int local = (var(ii) < 1);
         }
     }
     else
     {
        w1 = w1,@w[ii];
        v1 = v1+list(string(var(ii)));
     }
  }

  if ( size(w1) <= 1 )
  {
    return(BR);
  }
  if ( size(w2) <= 1 )
  {
    ERROR("## elimination of all variables is not possible");
  }

  w1 = w1[2..size(w1)];
  w2 = w2[2..size(w2)];
  BRlist[2] = v1 + v2;         //put variables to be eliminated in front

  //-------- create elimination ordering with two blocks and weights ---------
  //Assume that the first r of the n variables are to be eliminated.
  //Then, in case of an a-ordering (default), the new ring ordering will be
  //of the form (a(1..1,0..0),dp) with r 1's and n-r 0's or (a(w1,0..0),wp(@w))
  //if there are varaible weights which are not 1.
  //In the case of a b-ordering the ordering will be a block ordering with two
  //blocks of the form (dp(r),dp(n-r))  resp. (wp(w1),dp(w2))

  list B3;                     //this will become the list for new ordering

  //----- b-ordering case:
  if ( str == "b" )
  {
    if( w1==1 )              //weights for vars to be eliminated are all 1
    {
       B3[1] = list("dp", w1);
    }
    else
    {
       B3[1] = list("wp", w1);
    }

    if( w2==1 )              //weights for vars not to be eliminated are all 1
    {
       if ( local )
       {
          B3[2] = list("ds", w2);
       }
       else
       {
          B3[2] = list("dp", w2);
       }
    }
    else
    {
       if ( local )
       {
          B3[2] = list("ws", w2);
       }
       else
       {
          B3[2] = list("wp", w2);
       }
    }
  }

  //----- a-ordering case:
  else
  {
    //define first the second block
    if( @w==1 )              //weights for all vars are 1
    {
       if ( local )
       {
          B3[2] = list("ls", @w);
       }
       else
       {
          B3[2] = list("dp", @w);
       }
    }
    else
    {
       if ( local )
       {
          B3[2] = list("ws", @w);
       }
       else
       {
          B3[2] = list("wp", @w);
       }
    }

    //define now the first a-block of the form a(w1,0..0)
    intvec @v;
    @v[nvarBR] = 0;
    @v = @v+w1;
    B3[1] = list("a", @v);
  }
  BRlist[3] = B3;

  //----------- put module ordering always at the end and return -------------

  BRlist[3] = insert(BRlist[3],list("C",intvec(0)),size(B3));

  def eRing = ring(quotientList(BRlist));
  list result = eRing, @w;
  return (result);
}
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)[1];
   setring Q; Q;
}
///////////////////////////////////////////////////////////////////////////////

proc elim (id, list #)
"USAGE:   elim(id,arg[,s]);  id ideal/module, arg can be either an intvec v or
         a product p of variables (type poly), s a string determining the
         method which can be \"slimgb\" or \"std\" or, additionally,
         \"withWeigts\".
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 an elimination ordering for
        the variables to be eliminated and then applies std if \"std\"
        is given, or slimgb if \"slimgb\" is given, or a heuristically chosen
        method.
@*      If the variables in the basering have weights these weights are used
        in elimRing. If a string \"withWeigts\" as (optional) argument is given
        @sc{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;
        depending on the example, any of these commands can be 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;
  intvec save_opt=option(get);
  list lER;                          //for list returned by elimRing
//-------------------------------- check input -------------------------------
  poly vars;
  int ne;                           //for number of vars to be eliminated
  int ii;
  if (size(#) > 0)
  {
    if ( typeof(#[1]) == "poly" )
    {
      vars = #[1];
      for( ii=1; ii<=nvars(BR); ii++ )
      {
        if ( vars/var(ii) != 0)
        { ne++; }
      }
    }
    if ( typeof(#[1]) == "intvec" or typeof(#[1]) == "int")
    {
      ne = size(#[1]);
      vars=1;
      for( ii=1; ii<=ne; ii++ )
      {
        vars=vars*var(#[1][ii]);
      }
    }
  }

  string method;                    //for "std" or "slimgb" or "withWeights"
  if (size(#) >= 2)
  {
    if ( typeof(#[2]) == "string" )
    {
       if ( #[2] == "withWeights" )
       {
         intvec @w = weight(id);        //computation of weights
       }
       if ( #[2] == "std" ) { method = "std"; }
       if ( #[2] == "slimgb" ) { method = "slimgb"; }

    }
    else
    {
      ERROR("expected `elim(ideal,intvec[,string])`");
    }

    if (size(#) == 3)
    {
       if ( typeof(#[3]) == "string" )
       {
          if ( #[3] == "withWeights" )
          {
            intvec @w = weight(id);        //computation of weights
          }
          if ( #[3] == "std" ) { method = "std"; }
          if ( #[3] == "slimgb" ) { method = "slimgb"; }
       }
    }

    if ( method == "" )
    {
      ERROR("expected \"std\" or \"slimgb\" or \"withWeights\" as the optional string parameters");
    }
  }

//-------------- create new ring and map objects to new ring ------------------
   if ( defined(@w) )
   {
      lER = elimRing(vars,@w);  //in this case lER[2] = @w
   }
   else
   {
      lER = elimRing(vars);
      intvec @w = lER[2];      //in this case w is the intvec of
                               //variable weights as computed in elimRing
   }
   def ER = lER[1];
   setring ER;
   def id = imap(BR,id);
   poly vars = imap(BR,vars);

//---------- now eliminate in new ring and map back to old ring ---------------
  //if possible apply std(id,hi,w) where hi is the first hilbert function
  //of id with respect to the weights w. If w is not defined (i.e. good weights
  //@w are computed then id is only approximately @w-homogeneous and
  //the hilbert driven std cannot be used directly; however, stdhilb
  //homogenizes first and applies the hilbert driven std to the homogenization

   option(redThrough);
   if (typeof(id)=="matrix")
   {
     id = matrix(stdhilb(module(id),method,@w));
   }
   else
   {
     id = stdhilb(id,method,@w);
   }

   //### Todo: hier sollte id = groebner(id, "hilb"); verwendet werden.
   //da z.Zt. (Jan 09) groebener bei extra Gewichtsvektor a(...) aber stets std
   //aufruft und ausserdem "withWeigts" nicht kennt, ist groebner(id, "hilb")
   //zunaechst nicht aktiviert. Ev. nach Ueberarbeitung von groebner aktivieren

   id = nselect(id,1..ne);
   if ( pr > 0 )
   {
     "// result is a SB in the following ring:";
      ER;
   }
   option(set,save_opt);
   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","slimgb");
   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, list #)
"USAGE:   elim1(id,arg); id ideal/module, arg can be either an intvec v or 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
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 ("elim1 cannot eliminate in a qring");
   }
//------------- create product vars of variables to be eliminated -------------
  poly vars;
  int ii;
  if (size(#) > 0)
  {
    if ( typeof(#[1]) == "poly" ) {  vars = #[1]; }
    if ( typeof(#[1]) == "intvec" or typeof(#[1]) == "int")
    {
      vars=1;
      for( ii=1; ii<=size(#[1]); ii++ )
      {
        vars=vars*var(#[1][ii]);
      }
    }
  }
//---- get variables to be eliminated and create string for new ordering ------
   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,3..4); 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;
        }
      }
   }
   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;
   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;
         }
      }
   }
   if(typeof(id)=="matrix")
   {
      return(matrix(simplify(I,2)));
   }
   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);
   select1(matrix(m),1..2);
}
/*
///////////////////////////////////////////////////////////////////////////////
//                      EXAMPLEs
///////////////////////////////////////////////////////////////////////////////
// Siehe auch file 'tst-elim' mit grossem Beispiel;
example blowup0;
example elimRing;
example elim;
example elim1;
example nselect;
example sat;
example select;
example select1;
//===========================================================================
//         Rationale Normalkurve vom Grad d im P^d bzw. im A^d:
//homogen s:t -> (t^d:t^(d-1)s: ...: s^d), inhomogen t ->(t^d:t^(d-1): ...:t)

//------------------- 1. Homogen:
//Varianten der Methode
int d = 5;
ring R = 0,(s,t,x(0..d)),dp;
ideal I;
for( int ii=0; ii<=d; ii++) {I = I,ideal(x(ii)-t^(d-ii)*s^ii); }

int tt = timer;
ideal eI = elim(I,1..2,"std");
ideal eI = elim(I,1..2,"slimgb");
ideal eI = elim(I,st,"withWeights");
ideal eI = elim(I,st,"std","withWeights");
//komplizierter
int d = 50;
ring R = 0,(s,t,x(0..d)),dp;
ideal I;
for( int ii=0; ii<=d; ii++) {I = I,ideal(x(ii)-t^(d-ii)*s^ii); }
int tt = timer;
ideal eI = elim(I,1..2);               //56(44)sec (slimgb 22(17),hilb 33(26))
timer-tt; tt = timer;
ideal eI = elim1(I,1..2);              //71(53)sec
timer-tt; tt = timer;
ideal eI = eliminate(I,st);            //70(51)sec (wie elim1)
timer-tt;
timer-tt; tt = timer;
ideal eI = elim(I,1..2,"withWeights"); //190(138)sec
                                //(weights73(49), slimgb43(33), hilb71(53)
timer-tt;
//------------------- 2. Inhomogen
int d = 50;
ring r = 0,(t,x(0..d)),dp;
ideal I;
for( int ii=0; ii<=d; ii++) {I = I+ideal(x(ii)-t^(d-ii)); }
int tt = timer;
ideal eI = elim(I,1,);                //20(15)sec (slimgb13(10), hilb6(5))
ideal eI = elim(I,1,"std");           //17sec (std 11, hilb 6)
timer-tt; tt = timer;
ideal eI = elim1(I,t);               //8(6)sec
timer-tt; tt = timer;
ideal eI = eliminate(I,t);           //7(6)sec
timer-tt;
timer-tt; tt = timer;
ideal eI = elim(I,1..1,"withWeights"); //189(47)sec
//(weights73(42), slimgb43(1), hilb70(2)
timer-tt;

//===========================================================================
//        Zufaellige Beispiele, homogen
system("random",37);
ring R = 0,x(1..6),lp;
ideal I = sparseid(4,3);

int tt = timer;
ideal eI = elim(I,1);                //108(85)sec (slimgb 29(23), hilb79(61)
timer-tt; tt = timer;
ideal eI = elim(I,1,"std");          //(139)sec (std 77, hilb 61)
timer-tt; tt = timer;
ideal eI = elim1(I,1);               //(nach 45 min abgebrochen)
timer-tt; tt = timer;
ideal eI = eliminate(I,x(1));        //(nach 45 min abgebrochen)
timer-tt; tt = timer;

//        Zufaellige Beispiele, inhomogen
system("random",37);
ring R = 32003,x(1..5),dp;
ideal I = sparseid(4,2,3);
option(prot,redThrough);

intvec w = 1,1,1,1,1,1;
int tt = timer;
ideal eI = elim(I,1,w);            //(nach 5min abgebr.) hilb schlaegt nicht zu
timer-tt; tt = timer;              //BUG!!!!!!

int tt = timer;
ideal eI = elim(I,1);              //(nach 5min abgebr.) hilb schlaegt nicht zu
timer-tt; tt = timer;              //BUG!!!!!!
ideal eI = elim1(I,1);             //8(7.8)sec
timer-tt; tt = timer;
ideal eI = eliminate(I,x(1));      //8(7.8)sec
timer-tt; tt = timer;

                              BUG!!!!
//        Zufaellige Beispiele, inhomogen, lokal
system("random",37);
ring R = 32003,x(1..6),ds;
ideal I = sparseid(4,1,2);
option(prot,redThrough);
int tt = timer;
ideal eI = elim(I,1);                //(haengt sich auf)
timer-tt; tt = timer;
ideal eI = elim1(I,1);               //(0)sec !!!!!!
timer-tt; tt = timer;
ideal eI = eliminate(I,x(1));        //(ewig mit ...., abgebrochen)
timer-tt; tt = timer;

ring R1 =(32003),(x(1),x(2),x(3),x(4),x(5),x(6)),(a(1,0,0,0,0,0),ds,C);
ideal I = imap(R,I);
I = std(I);                           //(haengt sich auf) !!!!!!!

ideal eI = elim(I,1..1,"withWeights"); //(47)sec (weights42, slimgb1, hilb2)
timer-tt;

ring R1 =(32003),(x(1),x(2),x(3),x(4),x(5),x(6)),(a(1,0,0,0,0,0),ds,C);
ideal I = imap(R,I);
I = std(I);                           //(haengt sich auf) !!!!!!!

ideal eI = elim(I,1..1,"withWeights"); //(47)sec (weights42, slimgb1, hilb2)
timer-tt;
*/