source: git/Singular/LIB/sagbi.lib @ e2fd5d

//////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Commutative Algebra";
info="
LIBRARY: sagbi.lib  Compute SAGBI basis (subalgebra bases analogous to Groebner bases for ideals) of a subalgebra
AUTHORS: Jan Hackfeld,     Jan.Hackfeld@rwth-aachen.de
         Gerhard Pfister,  pfister@mathematik.uni-kl.de
         Viktor Levandovskyy,     levandov@math.rwth-aachen.de

OVERVIEW:
SAGBI stands for 'subalgebra bases analogous to Groebner bases for ideals'.
SAGBI bases provide important tools for working with finitely presented
subalgebras of a polynomial ring. Note, that in contrast to Groebner
bases, SAGBI bases may be infinite.

REFERENCES:
Ana Bravo: Some Facts About Canonical Subalgebra Bases, 
MSRI Publications  51, p. 247-254

PROCEDURES:
 sagbiSPoly(A [,r,m]); computes SAGBI S-polynomials of A
 sagbiReduce(I,A [,t,mt]);  performs subalgebra reduction of I by A
 sagbi(A [,m,t]);      computes SAGBI basis for A
 sagbiPart(A,k[,m]);   computes partial SAGBI basis for A
 algebraicDependence(I,it); performs iterations of SAGBI for algebraic dependencies of I

SEE ALSO: algebra_lib
";

LIB "elim.lib";
LIB "toric.lib";
LIB "algebra.lib";
//////////////////////////////////////////////////////////////////////////////

static proc assumeQring()
{
  if (ideal(basering) != 0)
  {
    ERROR("This function has not yet been implemented over qrings.");
  }
}


static proc uniqueVariableName (string variableName)
{
  //Adds character "@" at the beginning of variableName until this name ist unique
  //(not contained in the names of the ring variables or description of the coefficient field)
  string ringVars = charstr(basering) + "," + varstr(basering);
  while (find(ringVars,variableName) <> 0)
  {
    variableName="@"+variableName;
  }
  return(variableName);
}

static proc extendRing(r, ideal leadTermsAlgebra, int method) {
  /* Extends ring r with additional variables. If k=ncols(leadTermsAlgebra) and
   * r contains already m additional variables @y, the procedure adds k-m variables
   * @y(m+1)...@y(k) to the ring.
   * The monomial ordering of the extended ring depends on method.
   * Important: When calling this function, the basering (where algebra is defined) has to be active
   */
  def br=basering;
  int i;
  ideal varsBasering=maxideal(1);
  int numTotalAdditionalVars=ncols(leadTermsAlgebra);
  string variableName=uniqueVariableName("@y");
  //get a variable name different from existing variables

  //-------- extend current baserring r with new variables @y,
  // one for each new element in ideal algebra  -------------
  list l = ringlist(r);
  for (i=nvars(r)-nvars(br)+1; i<=numTotalAdditionalVars;i++)
  {
    l[2][i+nvars(br)]=string(variableName,"(",i,")");
  }
  if (method>=0 && method<=1)
  {
    if (nvars(r)==nvars(br))
    {        //first run of spolynomialGB in sagbi construction algorithms
      l[3][size(l[3])+1]=l[3][size(l[3])]; //save module ordering
      l[3][size(l[3])-1]=list("dp",intvec(1:numTotalAdditionalVars));
    }
    else
    {        //overwrite existing order for @y(i) to only get one block for the @y
      l[3][size(l[3])-1]=list("dp",intvec(1:numTotalAdditionalVars));
    }
  }
  // VL : todo noncomm case: correctly use l[5] and l[6]
  // that is update matrices
  // at the moment this is troublesome, so use nc_algebra call
  // see how it done in algebraicDependence proc // VL
  def rNew=ring(l);
  setring br;
  return(rNew);
}


static proc stdKernPhi(ideal kernNew, ideal kernOld, ideal leadTermsAlgebra,int method)
{
  /* Computes Groebner basis of kernNew+kernOld, where kernOld already is a GB
   * and kernNew contains elements of the form @y(i)-leadTermsAlgebra[i] added to it.
   * The techniques chosen is specified by the integer method
   */
  ideal kern;
  attrib(kernOld,"isSB",1);
  if (method==0)
  {
    kernNew=reduce(kernNew,kernOld);
    kern=kernOld+kernNew;
    kern=std(kern);
    //kern=std(kernOld,kernNew); //Found bug using this method.
    // TODO Change if bug is removed
    //this call of std return Groebner Basis of ideal kernNew+kernOld
    // given that kernOld is a Groebner basis
  }
  if (method==1)
  {
    kernNew=reduce(kernNew,kernOld);
    kern=slimgb(kernNew+kernOld);
  }
  return(kern);
}


static proc spolynomialsGB(ideal algebra,r,int method)
{
  /* This procedure does the actual S-polynomial calculation using Groebner basis methods and is
   * called by the procedures sagbiSPoly,sagbi and sagbiPart. As this procedure is called
   * at each step of the SAGBI construction algorithm, we can reuse the information already calculated
   * which is contained in the ring r. This is done in the following order
   * 1. If r already contain m additional variables and m'=number of elements in algebra, extend r with variables @y(m+1),...,@y(m')
   * 2. Transfer all objects to this ring, kernOld=kern is the Groebnerbasis already computed
   * 3. Define ideal kernNew containing elements of the form leadTermsAlgebra(m+1)-@y(m+1),...,leadTermsAlgebra(m')-@y(m')
   * 4. Compute Groebnerbasis of kernOld+kernNew
   * 5. Compute the new algebraic relations
   */
  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"//Spoly-1- initialisation and precomputation");  
  def br=basering;
  ideal varsBasering=maxideal(1);
  ideal leadTermsAlgebra=lead(algebra);
  //save leading terms as ordering in ring extension
  //may not be compatible with ordering in basering
  int numGenerators=ncols(algebra);

  def rNew=extendRing(r,leadTermsAlgebra,method);
  // important: br has to be active here
  setring r;
  if (!defined(kern))
  //only true for first run of spolynomialGB in sagbi construction algorithms
  {
    ideal kern=0;
    ideal algebraicRelations=0;
  }
  setring rNew;
  //-------------------------- transfer object to new ring rNew ----------------------
  ideal varsBasering=fetch(br,varsBasering);
  ideal kernOld,algebraicRelationsOld;
  kernOld=fetch(r,kern); //kern is Groebner basis of the kernel of the map Phi:r->K[x_1,...,x_n], x(i)->x(i), @y(i)->leadTermsAlgebra(i)
  algebraicRelationsOld=fetch(r,algebraicRelations);
  ideal leadTermsAlgebra=fetch(br,leadTermsAlgebra);
  ideal listOfVariables=maxideal(1);
  //---------define kernNew containing elements to be added to the ideal kern --------
  ideal kernNew;
  for (int i=nvars(r)-nvars(br)+1; i<=numGenerators; i++)
  {
    kernNew[i-nvars(r)+nvars(br)]=leadTermsAlgebra[i]-listOfVariables[i+nvars(br)];
  }
  //--------------- calculate kernel of Phi depending on method choosen ---------------
  dbprint(ppl,"//Spoly-2- Groebner basis computation"); 
  attrib(kernOld,"isSB",1);
  ideal kern=stdKernPhi(kernNew,kernOld,leadTermsAlgebra,method);
  dbprint(ppl-2,"//Spoly-2-1- ideal kern",kern);        
  //-------------------------- calulate algebraic relations -----------------------
  dbprint(ppl,"//Spoly-3- computing new algebraic relations");  
  ideal algebraicRelations=nselect(kern,1..nvars(br));
  attrib(algebraicRelationsOld,"isSB",1);
  ideal algebraicRelationsNew=reduce(algebraicRelations,algebraicRelationsOld);
  /* canonicalizing: */
  algebraicRelationsNew=canonicalform(algebraicRelationsNew);
  dbprint(ppl-2,"//Spoly-3-1- ideal of new algebraic relations",algebraicRelationsNew);
  /*        algebraicRelationsOld is a groebner basis by construction (as variable
   *        ordering is
   *        block ordering we have an elemination ordering for the varsBasering)
   *        Therefore, to only get the new algebraic relations, calculate
   *        <algebraicRelations>\<algebraicRelationsOld> using groebner reduction
   */
  kill kernOld,kernNew,algebraicRelationsOld,listOfVariables;
  export algebraicRelationsNew,algebraicRelations,kern;
  setring br;
  return(rNew);
}

static proc spolynomialsToric(ideal algebra) {
  /* This procedure does the actual S-polynomial calculation using toric.lib for
   * computation of a Groebner basis for the toric ideal kern(phi), where
   * phi:K[y_1,...,y_m]->K[x_1,...,x_n], y_i->leadmonom(algebra[i])
   * By suitable substitutions we obtain the kernel of the map
   * K[y_1,...,y_m]->K[x_1,...,x_n], x(i)->x(i), @y(i)->leadterm(algebra[i])
   */
  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"//Spoly-1- initialisation and precomputation");  
  def br=basering;
  int m=ncols(algebra);
  int n=nvars(basering);
  intvec tempVec;
  int i,j;
  ideal leadCoefficients;
  for (i=1;i<=m; i++)
  {
    leadCoefficients[i]=leadcoef(algebra[i]);
  }
  dbprint(ppl-2,"//Spoly-1-1- Vector of leading coefficients",leadCoefficients);
  int k=1;
  for (i=1;i<=n;i++)
  {
    for (j=1; j<=m; j++)
    {
      tempVec[k]=leadexp(algebra[j])[i];
      k++;
    }
  }
  //The columns of the matrix A are now the exponent vectors
  //of the leadings monomials in algebra.
  intmat A[n][m]=intmat(tempVec,n,m);
  dbprint(ppl-2,"//Spoly-1-2- Matrix A",A);
  //Create the preimage ring K[@y(1),...,@y(m)], where m=ncols(algebra).
  string variableName=uniqueVariableName("@y");
  list l = ringlist(basering);
  for (i=1; i<=m;i++)
  {
    l[2][i]=string(variableName,"(",i,")");
  }
  l[3][2]=l[3][size(l[3])];
  l[3][1]=list("dp",intvec(1:m));
  def rNew=ring(l);
  setring rNew;
  //Use toric_ideal to compute the kernel
   dbprint(ppl,"//Spoly-2- call of toric_ideal");
  ideal algebraicRelations=toric_ideal(A,"ect");
  //Suitable substitution
        dbprint(ppl,"//Spoly-3- substitutions");        
  ideal leadCoefficients=fetch(br,leadCoefficients);
  for (i=1; i<=m; i++)
  {
    if (leadCoefficients[i]!=0)
    {
      algebraicRelations=subst(algebraicRelations,var(i),1/leadCoefficients[i]*var(i));
    }
  }
        dbprint(ppl-2,"//Spoly-3-1- algebraic relations",algebraicRelations);   
  export algebraicRelations;
  return(rNew);
}


static proc reductionGB(ideal F, ideal algebra,r, int tailreduction,int method,int parRed)
{
  /* This procedure does the actual SAGBI/subalgebra reduction using GB methods and is
   * called by the procedures sagbiReduce,sagbi and sagbiPart
   * If r already is an extension of the basering
   * and contains the ideal kern needed for the subalgebra reduction,
   * the reduction can be started directly, at each reduction step using the fact that
   * p=reduce(leadF,kern) in K[@y(1),...,@y(m)] <=> leadF in K[lead(algebra)]
   * Otherwise some precomputation has to be done, outlined below.
   * When using sagbiReduce,sagbi and sagbiPart the integer parRed will always be zero. Only the procedure
   * algebraicDependence causes this procedure to be called with parRed<>0. The only difference when parRed<>0
   * is that the reduction algorithms returns the non-zero constants it attains (instead of just returning zero as the
         * correct remainder), as they will be expressions in parameters for an algebraic dependence.
   */
  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"//Red-1- initialisation and precomputation");    
  def br=basering;
  int numVarsBasering=nvars(br);
  ideal varsBasering=maxideal(1);
  int i;

  if (numVarsBasering==nvars(r))
  {
                dbprint(ppl-1,"//Red-1-1- Groebner basis computation");
    /* Case that ring r is the same ring as the basering. Using proc extendRing,
     * stdKernPhi
     * one construct the extension of the current baserring with new variables @y, one for each element
     * in ideal algebra and calculates the kernel of Phi, where
     * Phi: r---->br, x_i-->x_i, y_i-->f_i,
     * algebra={f_1,...f_m}, br=K[x1,...,x_n] und r=K[x1,...x_n,@y1,...@y_m]
     * This is similarly dones
     * (however step by step for each run of the SAGBI construction algorithm)
     * in the procedure spolynomialsGB
     */
    ideal leadTermsAlgebra=lead(algebra);
    kill r;
    def r=extendRing(br,leadTermsAlgebra,method);
    setring r;
    ideal listOfVariables=maxideal(1);
    ideal leadTermsAlgebra=fetch(br,leadTermsAlgebra);
    ideal kern;
    for (i=1; i<=ncols(leadTermsAlgebra); i++)
    {
      kern[i]=leadTermsAlgebra[i]-listOfVariables[numVarsBasering+i];
    }
    kern=stdKernPhi(kern,0,leadTermsAlgebra,method);
    dbprint(ppl-2,"//Red-1-1-1- Ideal kern",kern);
  }
  setring r;
  poly p,leadF;
  ideal varsBasering=fetch(br,varsBasering);
  setring br;
  map phi=r,varsBasering,algebra;
  poly p,normalform,leadF;
  intvec tempExp;
  //-------------algebraic reduction for each polynomial F[i] ------------------------
  dbprint(ppl,"//Red-2- reduction, polynomial by polynomial");
  for (i=1; i<=ncols(F);i++)
  {
                dbprint(ppl-1,"//Red-2-"+string(i)+"- starting with new polynomial");
                dbprint(ppl-2,"//Red-2-"+string(i)+"-1- Polynomial before reduction",F[i]);
    normalform=0;
    while (F[i]!=0)
    {
      leadF=lead(F[i]);
      if(leadmonom(leadF)==1)
      {
      //K is always contained in the subalgebra,
      //thus the remainder is zero in this case
        if (parRed) 
        { 
                                //If parRed<>0 save non-zero constants the reduction algorithms attains.
                                        break; 
                                }
        else 
        { 
                                        F[i]=0; 
                                        break; 
                                }
      }
      //note: as the ordering in br and r might not be compatible
      //it can be that lead(F[i]) in r is
      //different from lead(F[i]) in br.
      //To take the "correct" leading term therefore take lead(F[i])
      //in br and transfer it to the extension r
      setring r;
      leadF=fetch(br,leadF);
      p=reduce(leadF,kern);
      if (leadmonom(p)<varsBasering[numVarsBasering])
      {
        //as choosen ordering is a block ordering,
        //lm(p) in K[y_1...y_m] is equivalent to lm(p)<x_n
        //Needs to be changed, if no block ordering is used!
        setring br;
        F[i]=F[i]-phi(p);
      }
      else
      {
        if (tailreduction)
        {
          setring br;
          normalform=normalform+lead(F[i]);
          F[i]=F[i]-lead(F[i]);
        }
        else
        {
          setring br;
          break;
        }
      }
    }
    if (tailreduction)
    {
      F[i] = normalform;
    }
    dbprint(ppl-2,"//Red-2-"+string(i)+"-2- Polynomial after reduction",F[i]);
  }
  return(F);
}

static proc reduceByMonomials(ideal algebra)
/*This procedures uses the sagbiReduce procedure to reduce all polynomials in algebra,
 * which are not monomials, by the subset of all monomials.
 */
{
  ideal monomials;
  int i;
  for (i=1; i<=ncols(algebra);i++)
  {
    if(size(algebra[i])==1)
    {
      monomials[i]=algebra[i];
      algebra[i]=0;
    }
    else
    {
      monomials[i]=0;
    }
  }
  //Monomials now contains the subset of all monomials in algebra,
  //algebra contains the non-monomials.
  if(size(monomials)>0)
  {
    algebra=sagbiReduce(algebra,monomials,1);
    for (i=1; i<=ncols(algebra);i++)
    {
      if(size(monomials[i])==1)
      {
        //Put back monomials into algebra.
        algebra[i]=monomials[i];
      }
    }
  }
  return(algebra);
}


static proc sagbiConstruction(ideal algebra, int iterations, int tailreduction, int method,int parRed)
/* This procedure is the SAGBI construction algorithm and does the actual computation
 * both for the procedure sagbi and sagbiPart.
 * - If the sagbi procedure calls this procedure, iterations==-1
 *   and this procedure only stops
 *   if all S-Polynomials reduce to zero
 *   (criterion for termination of SAGBI construction algorithm).
 * - If the sagbiPart procedure calls this procedure, iterations>=0
 *   and iterations specifies the
 *   number of iterations. A degree boundary is not used here.
 * When this method is called via the procedures sagbi and sagbiPart the integer parRed 
 * will always be zero. Only the procedure algebraicDependence calls this procedure with 
 * parRed<>0. The only difference when parRed<>0 is that the reduction algorithms returns 
 * the non-zero constants it attains (instead of just returning zero as the correct 
 * remainder), as they will be expressions in parameters for an algebraic dependence. 
 * These constants are saved in the ideal reducedParameters.
 */
{
  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"// -0- initialisation and precomputation");      
  def br=basering;
  int i=1;

  ideal reducedParameters;
  int numReducedParameters=1; //number of elements plus one in reducedParameters
  int j;
  if (parRed==0) //if parRed<>0 the algebra does not contain monomials and normalisation should be avoided
  {
    algebra=reduceByMonomials(algebra);
    algebra=simplify(simplify(algebra,3),4);
  }
  // canonicalizing the gen's:
  algebra = canonicalform(algebra);
  ideal P=1;
  //note: P is initialized this way, so that the while loop is entered.
  //P gets overriden there, anyhow.
  ideal varsBasering=maxideal(1);
  map phi;
  ideal spolynomialsNew;
  def r=br;
  while (size(P)>0)
  {
    dbprint(ppl,"// -"+string(i)+"- interation of SAGBI construction algorithm");
    dbprint(ppl-1,"// -"+string(i)+"-1- Computing algebraic relations");
    def rNew=spolynomialsGB(algebra,r,method); /* canonicalizing inside! */
    kill r;
    def r=rNew;
    kill rNew;
    phi=r,varsBasering,algebra;
    dbprint(ppl-1,"// -"+string(i)+"-2- Substituting into algebraic relations");
    spolynomialsNew=simplify(phi(algebraicRelationsNew),6);
    //By construction spolynomialsNew only contains the spolynomials,
    //that have not already
    //been calculated in the steps before.
    dbprint(ppl-1,"// -"+string(i)+"-3- SAGBI reduction");
    dbprint(ppl-2,"// -"+string(i)+"-3-1- new S-polynomials before reduction",spolynomialsNew);
    P=reductionGB(spolynomialsNew,algebra,r,tailreduction,method,parRed);
    if (parRed)
    {
      for(j=1; j<=ncols(P); j++)
      {
        if (leadmonom(P[j])==1)
        {
          reducedParameters[numReducedParameters]=P[j];
          P[j]=0;
          numReducedParameters++;
        }
      }
    }
    if (parRed==0)
    {
      P=reduceByMonomials(P);
      //Reducing with monomials is cheap and can only result in less terms
      P=simplify(simplify(P,3),4);
      //Avoid that zeros are added to the bases or one element in P more than once
    }
    else
    {
      P=simplify(P,6);
    }
    /* canonicalize ! */
    P = canonicalform(P);
    dbprint(ppl-2,"// -"+string(i)+"-3-1- new S-polynomials after reduction",P);
    algebra=algebra,P;
    //Note that elements and order of elements must in algebra must not be changed,
    //otherwise the already calculated
    //ideal in r will give wrong results. Thus it is important to use a komma here.
    i=i+1;
    if (iterations!=-1 && i>iterations) //When iterations==-1 the number of iterations is unlimited
    {
      break;
    }
  }
  if (iterations!=-1)
  { //case that sagbiPart called this procedure
    if (size(P)==0)
    {
      dbprint(4-voice,
              "//SAGBI construction algorithm terminated after "+string(i-1)
              +" iterations, as all SAGBI S-polynomials reduced to 0.
//Returned generators therefore are a SAGBI basis.");
    }
    else
    {
      dbprint(4-voice,
              "//SAGBI construction algorithm stopped as it reached the limit of "
              +string(iterations)+" iterations.
//In general the returned generators are no SAGBI basis for the given algebra.");
    }
  }
  kill r;
  if (parRed)
  {
    algebra=algebra,reducedParameters;
  }
  algebra = simplify(algebra,6);
  algebra = canonicalform(algebra);
  return(algebra);
}


proc sagbiSPoly(ideal algebra,list #)
"USAGE:   sagbiSPoly(A[, returnRing, meth]);  A is an ideal, returnRing and meth are integers.
RETURN:   ideal or ring
ASSUME: basering is not a qring
PURPOSE: Returns SAGBI S-polynomials of the leading terms of a given ideal A if returnRing=0.
@*       Otherwise returns a new ring containing the ideals algebraicRelations
@*       and spolynomials, where these objects are explained by their name.
@*       See the example on how to access these objects.
@format     The other optional argument meth determines which method is
            used for computing the algebraic relations.
            - If meth=0 (default), the procedure std is used.
            - If meth=1, the procedure slimgb is used.
            - If meth=2, the prodecure uses toric_ideal.
@end format
EXAMPLE:  example sagbiSPoly; shows an example"
{
  assumeQring();
  int returnRing;
  int method=0;
  def br=basering;
  ideal spolynomials;
  if (size(#)>=1)
  {
    if (typeof(#[1])=="int")
    {
      returnRing=#[1];
    }
    else
    {
      ERROR("Type of first optional argument needs to be int.");
    }
  }
  if (size(#)==2)
  {
    if (typeof(#[2])=="int")
    {
      if (#[2]<0 || #[2]>2)
      {
        ERROR("Type of second optional argument needs to be 0,1 or 2.");
      }
      else
      {
        method=#[2];
      }
    }
    else
    {
      ERROR("Type of second optional argument needs to be int.");
    }
  }
  if (method>=0 and method<=1)
  {
    ideal varsBasering=maxideal(1);
    def rNew=spolynomialsGB(algebra,br,method);
    map phi=rNew,varsBasering,algebra;
    spolynomials=simplify(phi(algebraicRelationsNew),7);
  }
  if(method==2)
  {
    def r2=spolynomialsToric(algebra);
    map phi=r2,algebra;
    spolynomials=simplify(phi(algebraicRelations),7);
    def rNew=extendRing(br,lead(algebra),0);
    setring rNew;
    ideal algebraicRelations=imap(r2,algebraicRelations);
    export algebraicRelations;
    setring br;
  }

  if (returnRing==0)
  {
    return(spolynomials);
  }
  else
  {
    setring rNew;
    ideal spolynomials=fetch(br,spolynomials);
    export spolynomials;
    setring br;
    return(rNew);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r= 0,(x,y),dp;
  ideal A=x*y+x,x*y^2,y^2+y,x^2+x;
  //------------------ Compute the SAGBI S-polynomials only
  sagbiSPoly(A);
  //------------------ Extended ring is to be returned, which contains
  // the ideal of algebraic relations and the ideal of the S-polynomials
  def rNew=sagbiSPoly(A,1);  setring rNew;
  spolynomials;
  algebraicRelations;
  //----------------- Now we verify that the substitution of A[i] into @y(i)
  // results in the spolynomials listed above
  ideal A=fetch(r,A);
  map phi=rNew,x,y,A;
  ideal spolynomials2=simplify(phi(algebraicRelations),1);
  spolynomials2;
}


proc sagbiReduce(idealORpoly, ideal algebra, list #)
"USAGE: sagbiReduce(I, A[, tr, mt]); I, A ideals, tr, mt optional integers
RETURN: ideal of remainders of I after SAGBI reduction by A
ASSUME: basering is not a qring
PURPOSE:
@format
    The optional argument tr=tailred determines whether tail reduction will be performed.
     - If (tailred=0), no tail reduction is done.
     - If (tailred<>0), tail reduction is done.
     The other optional argument meth determines which method is
         used for Groebner basis computations.
         - If mt=0 (default), the procedure std is used.
         - If mt=1, the procedure slimgb is used.
@end format
EXAMPLE:  example sagbiReduce; shows an example"
{
  assumeQring();
  int tailreduction=0; //Default
  int method=0; //Default
  ideal I;
  if(typeof(idealORpoly)=="ideal")
  {
    I=idealORpoly;
  }
  else
  {
    if(typeof(idealORpoly)=="poly")
    {
      I[1]=idealORpoly;
    }
    else
    {
      ERROR("Type of first argument needs to be an ideal or polynomial.");
    }
  }
  if (size(#)>=1)
  {
    if (typeof(#[1])=="int")
    {
      tailreduction=#[1];
    }
    else
    {
      ERROR("Type of optional argument needs to be int.");
    }
  }
  if (size(#)>=2 )
  {
    if (typeof(#[2])=="int")
    {
      if (#[2]<0 || #[2]>1)
      {
        ERROR("Type of second optional argument needs to be 0 or 1.");
      }
      else
      {
        method=#[2];
      }
    }
    else
    {
      ERROR("Type of optional arguments needs to be int.");
    }
  }

  def r=basering;
  I=simplify(reductionGB(I,algebra,r,tailreduction,method,0),1);

  if(typeof(idealORpoly)=="ideal")
  {
    return(I);
  }
  else
  {
    if(typeof(idealORpoly)=="poly")
    {
      return(I[1]);
    }
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=0,(x,y,z),dp;
  ideal A=x2,2*x2y+y,x3y2;
  poly p1=x^5+x2y+y;
  poly p2=x^16+x^12*y^5+6*x^8*y^4+x^6+y^4+3;
  ideal P=p1,p2;
  //---------------------------------------------
  //SAGBI reduction of polynomial p1 by algebra A.
  //Default call, that is, no tail-reduction is done.
  sagbiReduce(p1,A);
  //---------------------------------------------
  //SAGBI reduction of set of polynomials P by algebra A,
  //now tail-reduction is done.
  sagbiReduce(P,A,1);
}

proc sagbi(ideal algebra, list #)
"USAGE:   sagbi(A[, tr, mt]); A ideal, tr, mt optional integers
RETURN: ideal, a SAGBI basis for A
ASSUME: basering is not a qring
PURPOSE: Computes a SAGBI basis for the subalgebra given by the generators in A.
@format
    The optional argument tr=tailred determines whether tail reduction will be performed.
     - If (tailred=0), no tail reduction is performed,
     - If (tailred<>0), tail reduction is performed.
     The other optional argument meth determines which method is
         used for Groebner basis computations.
         - If mt=0 (default), the procedure std is used.
         - If mt=1, the procedure slimgb is used.
@end format
EXAMPLE:  example sagbi; shows an example"
{
  assumeQring();
  int tailreduction=0; //default value
  int method=0; //default value
  if (size(#)>=1)
  {
    if (typeof(#[1])=="int")
    {
      tailreduction=#[1];
    }
    else
    {
      ERROR("Type of optional argument needs to be int.");
    }
  }
  if (size(#)>=2 )
  {
    if (typeof(#[2])=="int")
    {
      if (#[2]<0 || #[2]>1)
      {
        ERROR("Type of second optional argument needs to be 0 or 1.");
      }
      else
      {
        method=#[2];
      }
    }
    else
    {
      ERROR("Type of optional arguments needs to be int.");
    }
  }
  ideal a;
  a=sagbiConstruction(algebra,-1,tailreduction,method,0);
  a=simplify(a,7);
  //  a=interreduced(a);
  return(a);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r= 0,(x,y,z),dp;
  ideal A=x2,y2,xy+y;
  //Default call, no tail-reduction is done.
  sagbi(A);
  //---------------------------------------------
  //Call with tail-reduction and method specified.
  sagbi(A,1,0);
}

proc sagbiPart(ideal algebra, int iterations, list #)
"USAGE: sagbiPart(A, k,[tr, mt]); A is an ideal, k, tr and mt are integers
RETURN: ideal
ASSUME: basering is not a qring
PURPOSE: Performs k iterations of the SAGBI construction algorithm for the subalgebra given by the generators given by A.
@format
     The optional argument tr=tailred determines if tail reduction will be performed.
     - If (tailred=0), no tail reduction is performed,
     - If (tailred<>0), tail reduction is performed.
     The other optional argument meth determines which method is
         used for Groebner basis computations.
         - If mt=0 (default), the procedure std is used.
         - If mt=1, the procedure slimgb is used.
@end format
EXAMPLE:  example sagbiPart; shows an example"
{
  assumeQring();
  int tailreduction=0; //default value
  int method=0; //default value
  if (size(#)>=1)
  {
    if (typeof(#[1])=="int")
    {
      tailreduction=#[1];
    }
    else
    {
      ERROR("Type of optional argument needs to be int.");
    }
  }
  if (size(#)>=2 )
  {
    if (typeof(#[2])=="int")
    {
      if (#[2]<0 || #[2]>3)
      {
        ERROR("Type of second optional argument needs to be 0 or 1.");
      }
      else
      {
        method=#[2];
      }
    }
    else
    {
      ERROR("Type of optional arguments needs to be int.");
    }
  }
  if (iterations<0)
  {
    ERROR("Number of iterations needs to be non-negative.");
  }
  ideal a;
  a=sagbiConstruction(algebra,iterations,tailreduction,method,0);
  a=simplify(a,6);
  //  a=interreduced(a);
  return(a);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r= 0,(x,y,z),dp;
  //The following algebra does not have a finite SAGBI basis.
  ideal A=x,xy-y2,xy2;
  //---------------------------------------------------
  //Call with two iterations, no tail-reduction is done.
  sagbiPart(A,2);
  //---------------------------------------------------
  //Call with three iterations, tail-reduction and method 0.
  sagbiPart(A,3,1,0);
}


proc algebraicDependence(ideal I,int iterations)
"USAGE: algebraicDependence(I,it); I an an ideal, it is an integer
RETURN: ring
ASSUME: basering is not a qring
PURPOSE: Returns a ring containing the ideal @code{algDep}, which contains possibly 
@*              some algebraic dependencies of the elements of I obtained through @code{it}
@*              iterations of the SAGBI construction algorithms. See the example on how 
@*              to access these objects. 
EXAMPLE: example algebraicDependence; shows an example"
{
        assumeQring();
        int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"//AlgDep-1- initialisation and precomputation"); 
  def br=basering;
  int i;
  I=simplify(I,2); //avoid that I contains zeros

        //Create two polynomial rings, which both are extensions of the current basering.
  //The first ring will contain the additional paramteres @c(1),...,@c(m), the second one
  //will contain the additional variables @c(1),...,@c(m), where m=ncols(I).
  string parameterName=uniqueVariableName("@c");
  list l = ringlist(basering);
  list parList;
  for (i=1; i<=ncols(I);i++)
  {
    parList[i]=string(parameterName,"(",i,")");
  }
  l[1]=list(l[1],parList,list(list("dp",1:ncols(I)))); //add @c(i) to the ring as paramteres
  ideal temp=0;
  l[1][4]=temp;
  // addition VL: noncomm case
  int isNCcase = 0; // default for comm
  //         if (size(l)>4)
  //         {
  //           // that is we're in the noncomm algebra
  //           isNCcase = 1; // noncomm
  //           matrix @C@ = l[5];
  //           matrix @D@ = l[6];
  //           l = l[1],l[2],l[3],l[4];
  //         }
  def parameterRing=ring(l);

  string extendVarName=uniqueVariableName("@c");
  list l2 = ringlist(basering);
  for (i=1; i<=ncols(I);i++)
  {
    l2[2][i+nvars(br)]=string(extendVarName,"(",i,")"); //add @c(i) to the rings as variables
  }
  l2[3][size(l2[3])+1]=l2[3][size(l2[3])];
  l2[3][size(l2[3])-1]=list("dp",intvec(1:ncols(I)));
  //         if (isNCcase)
  //         {
  //           // that is we're in the noncomm algebra
  //           matrix @C@2 = l2[5];
  //           matrix @D@2 = l2[6];
  //           l2 = l2[1],l2[2],l2[3],l2[4];
  //         }

  def extendVarRing=ring(l2);
  setring extendVarRing;
  // VL : this requires extended matrices
  // let's forget it for the moment
  // since this holds only for showing the answer
  //         if (isNCcase)
  //         {
  //           matrix C2=imap(br,@C@2);
  //           matrix D2=imap(br,@D@2);
  //           def er2 = nc_algebra(C2,D2);
  //           setring er2;
  //           def extendVarRing=er2;
  //         }

  setring parameterRing;

  //         if (isNCcase)
  //         {
  //           matrix C=imap(br,@C@);
  //           matrix D=imap(br,@D@);
  //           def pr = nc_algebra(C,D);
  //           setring pr;
  //           def parameterRing=pr;
  //         }

        //Compute a partial SAGBI basis of the algebra generated by I[1]-@c(1),...,I[m]-@c(m),
  //where the @c(n) are parameters
  ideal I=fetch(br,I);
  ideal algebra;
  for (i=1; i<=ncols(I);i++)
  {
    algebra[i]=I[i]-par(i);
  }
  dbprint(ppl,"//AlgDep-2- call of SAGBI construction algorithm");      
  algebra=sagbiConstruction(algebra, iterations,0,0,1);
  dbprint(ppl,"//AlgDep-3- postprocessing of results"); 
  int j=1;
  //If K[x_1,...,x_n] was the basering, then algebra is in K(@c(1),...,@c(m))[x_1,...x_n]. We intersect
  //elements in algebra with K(@c(1),..,@c(n)) to get algDep. Note that @c(i) can only appear in the numerator,
  //as the SAGBI construction algorithms just multiplies and substracts polynomials. So actually we have
  //algDep=algebra intersect K[@c(1),...,@c(m)]
  ideal algDep;
  for (i=1; i<= ncols(algebra); i++)
  {
    if (leadmonom(algebra[i])==1) //leadmonom(algebra[i])==1 iff algebra[i] in K[@c(1),...,@c(m)]
    {
      algDep[j]=algebra[i];
      j++;
    }
  }
  //Transfer algebraic dependencies to ring where @c(i) are not parameters, but now variables.
  setring extendVarRing;
  ideal algDep=imap(parameterRing,algDep);
  ideal algebra=imap(parameterRing,algebra);
  //Now get rid of constants in K that may have been added to algDep.
  for (i=1; i<=ncols(algDep); i++)
  {
                if(leadmonom(algDep[i])==1) 
                {
                                algDep[i]=0;
                }
        }
        algDep=simplify(algDep,2);
  export algDep,algebra;
  setring br;
  return(extendVarRing);
}
example
{ "EXAMPLE:"; echo = 2;
  ring r= 0,(x,y),dp;
  //The following algebra does not have a finite SAGBI basis.
  ideal I=x^2, xy-y2, xy2;
  //---------------------------------------------------
  //Call with two iterations
  def DI = algebraicDependence(I,2);
  setring DI; algDep;
  // we see that no dependency has been seen so far
  //---------------------------------------------------
  //Call with two iterations
  setring r; kill DI;
  def DI = algebraicDependence(I,3);
  setring DI; algDep;
  map F = DI,x,y,x^2, xy-y2, xy2;
  F(algDep); // we see that it is a dependence indeed
}

static proc interreduced(ideal I)
{
  /* performs subalgebra interreduction of a set of subalgebra generators */
  int ppl = printlevel-voice+3; //variable for additional printlevel-dependend information
  dbprint(ppl,"//Interred-1- starting interreduction"); 
  ideal J,B;
  int i,j,k;
  poly f;
  for(k=1;k<=ncols(I);k++)
  {
    dbprint(ppl-1,"//Interred-1-"+string(k)+"- reducing next poly");
    f=I[k];
    I[k]=0;
    f=sagbiReduce(f,I,1);
    I[k]=f;
  }
  I=simplify(I,2);
  dbprint(ppl,"//Interred-2- interreduction completed");        
  return(I);
}
///////////////////////////////////////////////////////////////////////////////

proc sagbiReduction(poly p,ideal dom,list #)
"USAGE: sagbiReduction(p,dom[,n]); p poly , dom  ideal
RETURN: polynomial, after one step of subalgebra reduction
PURPOSE:
@format
    Three algorithm variants are used to perform subalgebra reduction.
    The positive interger n determines which variant should be used.
    n may take the values 0 (default), 1 or 2.
@end format
NOTE: works over both polynomial rings and their quotients
EXAMPLE: example sagbiReduction; shows an example"
{
  def bsr=basering;
  ideal B=ideal(bsr);//When the basering is quotient ring  this type casting
                     // gives the quotient ideal.
  int b=size(B);
  int n=nvars(bsr);

  //In quotient rings, SINGULAR, usually does not reduce polynomials w.r.t the
  //quotient ideal,therefore we should first  reduce,
  //when it is necessary for computations,
  // to have a uniquely determined representant for each equivalent
  //class,which is the case of this algorithm.

  if(b !=0) //means that the basering is a quotient ring
  {
    p=reduce(p,std(0));
    dom=reduce(dom,std(0));
  }

  int i,choose;
  int z=ncols(dom);

  if((size(#)>0) && (typeof(#[1])=="int"))
  {
    choose = #[1];
  }
  if (size(#)>1)
  {
    choose =#[2];
  }

  //=======================first algorithm(default)=========================
  if ( choose == 0 )
  {
    list L = algebra_containment(lead(p),lead(dom),1);
    if( L[1]==1 )
    {
      // the ring L[2] = char(bsr),(x(1..nvars(bsr)),y(1..z)),(dp(n),dp(m)),
      // contains poly check s.t. LT(p) is of the form check(LT(f1),...,LT(fr))
      def s1 = L[2];
      map psi = s1,maxideal(1),dom;
      poly re = p - psi(check);
      // divide by the maximal power of #[1]
      if ( (size(#)>0) && (typeof(#[1])=="poly") )
      {
        while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
        {
          re=re/#[1];
        }
      }
      return(re);
    }
    return(p);
  }
  //======================2end variant of algorithm=========================
  //It uses two different commands for elimaination.
  //if(choose==1):"elimainate"command.
  //if (choose==2):"nselect" command.
  else
  {
    poly v=product(maxideal(1));

    //------------- change the basering bsr to bsr[@(0),...,@(z)] ----------
    execute("ring s=("+charstr(basering)+"),("+varstr(basering)+",@(0..z)),dp;");
    // Ev hier die Reihenfolge der Vars aendern. Dazu muss unten aber entsprechend
    // geaendert werden:
    //  execute("ring s="+charstr(basering)+",(@(0..z),"+varstr(basering)+"),dp;");

    //constructs the leading ideal of dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z))
    ideal dom=imap(bsr,dom);
    for (i=1;i<=z;i++)
    {
      dom[i]=lead(dom[i])-var(nvars(bsr)+i+1);
    }
    dom=lead(imap(bsr,p))-@(0),dom;

    //---------- eliminate the variables of the basering bsr --------------
    //i.e. computes dom intersected with K[@(0),...,@(z)].

    if(choose==1)
    {
      ideal kern=eliminate(dom,imap(bsr,v));//eliminate does not need a
                                            //standard basis as input.
    }
    if(choose==2)
    {
      ideal kern= nselect(groebner(dom),1..n);//"nselect" is combinatorial command
                                         //which uses the internal command
                                         // "simplify"
    }

    //---------  test wether @(0)-h(@(1),...,@(z)) is in ker ---------------
    // for some poly h and divide by maximal power of q=#[1]
    poly h;
    z=size(kern);
    for (i=1;i<=z;i++)
    {
      h=kern[i]/@(0);
      if (deg(h)==0)
      {
        h=(1/h)*kern[i];
        // define the map psi : s ---> bsr defined by @(i) ---> p,dom[i]
        setring bsr;
        map psi=s,maxideal(1),p,dom;
        poly re=psi(h);
        // divide by the maximal power of #[1]
        if ((size(#)>0) && (typeof(#[1])== "poly") )
        {
          while ((re!=0) && (re!=#[1]) &&(subst(re,#[1],0)==0))
          {
            re=re/#[1];
          }
        }
        return(re);
      }
    }
    setring bsr;
    return(p);
  }
}
example
{"EXAMPLE:"; echo = 2;
  ring r= 0,(x,y),dp;
  ideal dom =x2,y2,xy-y;
  poly p=x4+x3y+xy2-y2;
  sagbiReduction(p,dom);
  sagbiReduction(p,dom,2);
  // now let us see the action over quotient ring
  ideal I = xy;
  qring Q = std(I);
  ideal dom = imap(r,dom); poly p = imap(r,p);
  sagbiReduction(p,dom,1);
}

proc sagbiNF(id,ideal dom,int k,list#)
"USAGE: sagbiNF(id,dom,k[,n]); id either poly or ideal,dom ideal, k and n positive intergers.
RETURN: same as type of id; ideal or polynomial.
PURPOSE:
@format
    The integer k determines what kind of s-reduction is performed:
    - if (k=0) no tail s-reduction is performed.
    - if (k=1) tail s-reduction is performed.
    Three Algorithm variants are used to perform subalgebra reduction.
    The positive integer n determines which variant should be used.
    n may take the values (0 or default),1 or 2.
@end format
NOTE: sagbiNF works over both rings and quotient rings
EXAMPLE: example sagbiNF; show example "
{
  ideal rs;
  if (ideal(basering) == 0)
  {
    rs = sagbiReduce(id,dom,k) ;
  }
  else
  {
    rs = sagbiReduction(id,dom,k) ;
  }
  return(rs);
}
example
{"EXAMPLE:"; echo = 2;
 ring r=0,(x,y),dp;
 poly p=x4+x2y+y;
 ideal dom =x2,x2y+y,x3y2;
 sagbiNF(p,dom,1);
 ideal I= x2-xy;
 qring Q=std(I); // we go to the quotient ring
 poly p=imap(r,p);
 NF(p,std(0)); // the representative of p has changed
 ideal dom = imap(r,dom);
 print(matrix(NF(dom,std(0)))); // dom has changed as well
 sagbiNF(p,dom,0); // no tail reduction
 sagbiNF(p,dom,1);// tail subalgebra reduction is performed
}

static proc canonicalform(ideal I)
{
  /* placeholder for the canonical form of a set of gen's */
  /* for the time being we agree on content(p)=1; that is coeffs with no fractions */
  int i; ideal J=I;
  for(i=ncols(I); i>=1; i--)
  {
    J[i] = canonicalform_poly(I[i]);
  }
  return(J);
}

static proc canonicalform_poly(poly p)
{
  /* placeholder for the canonical form of a poly */
  /* for the time being we agree on content(p)=1; that is coeffs with no fractions */
  number n = content(p);
  return( p/content(p) );
}

/*
  ring r= 0,(x,y),dp;
  //The following algebra does not have a finite SAGBI basis.
  ideal J=x^2, xy-y2, xy2, x^2*(x*y-y^2)^2 - (x*y^2)^2*x^4 + 11;
  //---------------------------------------------------
  //Call with two iterations
  def DI = algebraicDependence(J,2);
  setring DI; algDep;
*/