Changeset 16fba9 in git

Timestamp:

Oct 8, 2010, 2:57:36 PM (14 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

10c87563326684cb871f508482250e5d4481798e

Parents:

db357053d1d7c7b0f68829f8597192533f97171b

Message:

rewritten sagbi.lib (was sagbi2)

git-svn-id: file:///usr/local/Singular/svn/trunk@13423 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• Singular/LIB/sagbi.lib (modified) (4 diffs)

Singular/LIB/sagbi.lib

@@ -3 +3 @@
 category="Commutative Algebra";
 info="
-LIBRARY:  sagbi.lib  Compute subalgebra bases analogous to Groebner bases for ideals
-AUTHORS: Gerhard Pfister,        pfister@mathematik.uni-kl.de,
-@*       Anen Lakhal,            alakhal@mathematik.uni-kl.de
+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
 
 PROCEDURES:
- sagbiRreduction(p,I);  perform one step subalgebra reducton (for short S-reduction) of p w.r.t I
- sagbiSPoly(I);   compute the S-polynomials of the Subalgebra defined by the genartors of I
- sagbiNF(id,I);   perform iterated S-reductions in order to compute Subalgebras normal forms
- sagbi(I);        construct SAGBI basis for the Subalgebra defined by I
- sagbiPart(I);    construct partial SAGBI basis for the Subalgebra defined by I
+ sagbiSPoly(A [,r,m]); SAGBI S-polynomials of A
+ sagbiReduce(I,A);     subalgebra reduction of I by A
+ sagbi(A [,m,t]);      SAGBI basis for A
+ sagbiPart(A,k[,m]);   partial SAGBI basis for A
+ algDep(I,it);         iteration of SAGBI for algebraic dependencies of I
+
+SEE ALSO: algebra_lib
 ";
 
+
+//AUXILIARY PROCEDURES:
+//uniqueVariableName(s)                              adds character "@" at the beginning of string s until this is a unique new variable name
+//extendRing(r, ...)                                         creates a new ring, which is an extension of r
+//stdKernPhi(kernNew, kernOld,...) computes Groebner basis of kernNew+kernOld
+//spolynomialsGB(A,...)                                  computes the SAGBI S-polynomials of the subalgebra defined by the generators in A using Groebner bases
+//spolynomialsToric(A)                              computes the SAGBI S-polynomials of the subalgebra defined by the generators in A using toric.lib
+//reductionGB(ideal F, ideal A,....)         performs subalgebra reduction of F by A
+//reduceByMonomials(A)                                        performs subalgebra reduction of all polynomials in A by the subset of monomials
+
+LIB "elim.lib"; //for nselect
+LIB "toric.lib"; //for toric_ideal
 LIB "algebra.lib";
-LIB "elim.lib";
-
-///////////////////////////////////////////////////////////////////////////////
-proc sagbiSPoly(id ,list #)
-"USAGE: sagbiSPoly(id [,n]); id ideal, n positive integer.
-RETURN: an ideal
-@format
-      - If (n=0 or default) an ideal, whose generators are the S-polynomials.
-      - If (n=1) a list of size 2:
-        the first element of this list is the ideal of S-polynomials.
-        the second element of this list is the ring in which the ideal of algebraic
-        relations is defined.
+//////////////////////////////////////////////////////////////////////////////
+
+
+
+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 algDep 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 Groebner basis
+   * 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
+   */
+  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)];
+  }
+  //-------------------------- calulate kernel of Phi depending on method choosen -----------------------
+
+  attrib(kernOld,"isSB",1);
+  ideal kern=stdKernPhi(kernNew,kernOld,leadTermsAlgebra,method);
+  //-------------------------- calulate algebraic relations -----------------------
+  ideal algebraicRelations=nselect(kern,1..nvars(br));
+  attrib(algebraicRelationsOld,"isSB",1);
+  ideal algebraicRelationsNew=reduce(algebraicRelations,algebraicRelationsOld);
+  /*        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])
+   */
+  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]);
+  }
+  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);
+  //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
+  ideal algebraicRelations=toric_ideal(A,"ect");
+  //Suitable substitution
+  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));
+    }
+  }
+  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 Groebner basis 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.
+   */
+  def br=basering;
+  int numVarsBasering=nvars(br);
+  ideal varsBasering=maxideal(1);
+  int i;
+
+  if (numVarsBasering==nvars(r))
+  {
+    /* 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 done (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);
+  }
+  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] ---------------------------------------
+  for (i=1; i<=ncols(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) { 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;
+    }
+  }
+  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.
+ * Note that parRed is used for testing a special modification and can be ignored (assume parRed==0).
+ */
+{
+  def br=basering;
+  int i=1;
+
+  ideal reducedParameters;
+  int numReducedParameters=1;
+  int j;
+  if (parRed==0)
+  {
+    algebra=reduceByMonomials(algebra);
+    algebra=simplify(simplify(algebra,3),4);
+  }
+  int step=1;
+  if (iterations==-1) //case: infintitly many iterations
+  {
+    step=0;
+    iterations=1;
+  }
+  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 && i<=iterations)
+  {
+    def rNew=spolynomialsGB(algebra,r,method);
+    kill r;
+    def r=rNew;
+    kill rNew;
+    phi=r,varsBasering,algebra;
+    spolynomialsNew=simplify(phi(algebraicRelationsNew),6);
+    //By construction spolynomialsNew only contains the spolynomials, that have not already
+    //been calculated in the steps before.
+    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);
+    }
+    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+step;
+  }
+  if (step==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);
+  return(algebra);
+}
+
+
+proc sagbiSPoly(ideal algebra,list #)
+"USAGE:   sagbiSPoly(A[, returnRing, meth]);  A is an ideal, returnRing and meth are integers.
+RETURN:   Returns SAGBI S-polynomials of the leading terms of 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; show an example "
-{
-  if(size(#)==0)
-  {
-    #[1]=0;
-  }
-  degBound=0;
-  def bsr=basering;
-  ideal vars=maxideal(1);
-  ideal B=ideal(bsr);//when the basering is quotient ring this "type casting"
-                    //gives th quotient ideal.
-  int b=size(B);
-
-  //In quotient rings,SINGULAR does not reduce polynomials w.r.t the
-  //quotient ideal,therefore we should first 'reduce';if it is necessary for
-  //computations to have a uniquely determined representant for each equivalent
-  //class,which is the case of this procedure.
-
-  if(b!=0)
-  {
-    id =reduce(id,groebner(0));
-  }
-  int n,m=nvars(bsr),ncols(id);
-  int z;
-  string mp=string(minpoly);
-  ideal P;
-  list L;
-
-  if(id==0)
-  {
-    if(#[1]==0)
-    {
-      return(P);
-    }
-    else
-    {
-      return(L);
-    }
+EXAMPLE:  example sagbiSPoly; shows an example"
+{
+  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
   {
-    //==================create anew ring with extra variables================
-
-    execute("ring R1=("+charstr(bsr)+"),("+varstr(bsr)+",@y(1..m)),(dp(n),dp(m));");
-    execute("minpoly=number("+mp+");");
-    ideal id=imap(bsr,id);
-    ideal A;
-
-    for(z=1;z<=m;z++)
-    {
-      A[z]=lead(id[z])-@y(z);
-    }
-
-    A=groebner(A);
-    ideal kern=nselect(A,1..n);// "kern" is the kernel of the ring map:
-                        // R1----->bsr ;y(z)----> lead(id[z]).
-                        //"kern" is the ideal of algebraic relations between
-                        // lead(id[z]).
-    export kern,A;// we export:
-                  // * the ideal A  to avoid useless computations
-                  //   between 2 steps in sagbi procedure.
-                  // * the ideal kern : some times we can get intresting
-                  //   informations from this ideal.
-    setring bsr;
-    map phi=R1,vars,id;
-
-    // the sagbiSPolynomials are the image by phi of the generators of kern
-
-    P=simplify(phi(kern),1);
-    if(#[1]==0)
-    {
-      return(P);
-    }
-    else
-    {
-      L=P,R1;
-      kill phi,vars;
-
-      dbprint(printlevel-voice+3,"
-// 'sagbiSPoly' created a ring as 2nd element of the list.
-// The ring contains the ideal 'kern'  of algebraic relations between the
-//leading terms of the generators of I.
-// To access to this ring and see 'kern' you should give the ring a name,
-// e.g.:
-               def S = L[2]; setring S; kern;
-      ");
-      return(L);
-    }
+    setring rNew;
+    ideal spolynomials=fetch(br,spolynomials);
+    export spolynomials;
+    setring br;
+    return(rNew);
   }
 }
 example
 { "EXAMPLE:"; echo = 2;
- ring r=0, (x,y),dp;
- poly f1,f2,f3,f4=x2,y2,xy+y,2xy2;
- ideal I=f1,f2,f3,f4;
- sagbiSPoly(I);
- list L=sagbiSPoly(I,1);
- L[1];
- def S= L[2]; setring S; kern;
-}
-
-///////////////////////////////////////////////////////////////////////////////
-static proc std1(ideal J,ideal I,list #)
-  // I is contained in J, and it is assumed to be a standard bases!
-  //This procedure computes a Standard basis for J from I one's
-  //This procedure is essential for Spoly1 procedure.
+  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:   Returns the remainder(s) of I after SAGBI reduction by A
+@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"
+{
+  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: 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"
+{
+  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.");
+    }
+  }
+  algebra=sagbiConstruction(algebra,-1,tailreduction,method,0);
+  algebra=simplify(algebra,7);
+  algebra=interreduced(algebra);
+  return(algebra);
+}
+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: Performs k iterations of the SAGBI construction algorithm for the subalgebra given by the generators in 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"
+{
+  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.");
+  }
+  algebra=sagbiConstruction(algebra,iterations,tailreduction,method,0);
+  algebra=simplify(algebra,7);
+  algebra=interreduced(algebra);
+  return(algebra);
+}
+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);
+}
+
+
+//DONE? 1nsen als werden noch mit ausgegeben, aus algebraischen AbhÀngigkeiten entfernen
+// VL: finished the documentation
+// TO comapre with algDependent from algebra_lib
+proc algDep(ideal I,int iterations)
+"USAGE: algDep(I,it); I an an ideal, it is an integer
+RETURN: In 'it' iterations, compute algebraic dependencies between elements of I
+EXAMPLE: example algDep; shows an example"
 {
   def br=basering;
-  int tt;
-  ideal Res,@result;
-
-  if(size(#)>0) {tt=#[1];}
-
-  if(size(I)==0) {@result=groebner(J);}
-
-  if((size(I)!=0) && (size(J)-size(I)>=1))
-  {
-    qring Q=I;
-    ideal J=fetch(br,J);
-    J=groebner(J);
-    setring br;
-    Res=fetch(Q,J);// Res contains the generators that we add to I
-                   // to get the generators of std(J);
-    @result=Res+I;
-  }
-
-  if(tt==0) { return(@result);}
-  else      { return(Res);}
-}
-
-///////////////////////////////////////////////////////////////////////////////
-
-static proc Spoly1(list l,ideal I,ideal J,int a)
-  //an implementation of SAGBI construction Algorithm using Spoly
-  //procedure leads to useless computations and affect the efficiency
-  //of SAGBI bases computations. This procedure is a variant of Spoly
-  //in order to avoid these useless compuations.
-{
-  degBound=0;
-  def br=basering;
-  ideal vars=maxideal(1);
-  ideal B=ideal(br);
-  int b=size(B);
-
-  if(b!=0)
-  {
-    I=reduce(I,groebner(0));
-    J=reduce(J,groebner(0));
-  }
-  int n,ii,jj=nvars(br),ncols(I),ncols(J);
-  int z;
-  list @L;
-  string mp =string(minpoly);
-
-  if(size(J)==0)
-  {
-    @L =sagbiSPoly(I,1);
-  }
-  else
-  {
-    ideal @sum=I+J;
-    ideal P1;
-    ideal P=l[1];//P is the ideal of spolynomials of I;
-    def R=l[2];setring R;int kk=nvars(R);
-    ideal J=fetch(br,J);
-
-    //================create a new ring with extra variables==============
-    execute("ring R1=("+charstr(R)+"),("+varstr(R)+",@y((ii+1)..(ii+jj))),(dp(n),dp(kk+jj-n));");
-    // *levandov: would it not be easier and better to use
-    // ring @Y = char(R),(@y((ii+1)..(ii+jj))),dp;
-    // def R1 = R + @Y;
-    // setring R1;
-    // -> thus
-    ideal kern1;
-    ideal A=fetch(R,A);
-    attrib(A,"isSB",1);
-    ideal J=fetch(R,J);
-    ideal kern=fetch(R,kern);
-    ideal A1;
-    for(z=1;z<=jj;z++)
-    {
-      A1[z]=lead(J[z])-var(z+kk);
-    }
-    A1=A+A1;
-    ideal @Res=std1(A1,A,1);// the generators of @Res are whose we have to add
-                            // to A to get std(A1).
-    A=A+@Res;
-    kern1=nselect(@Res,1..n);
-    kern=kern+kern1;
-    export kern,kern1,A;
-    setring br;
-    map phi=R1,vars,@sum;
-    P1=simplify(phi(kern1),1);//P1 is th ideal we add to P to get the ideal
-                              //of Spolynomials of @sum.
-    P=P+P1;
-
-    if (a==1)
-    {
-      @L=P,R1;
-      kill phi,vars;
-      dbprint(printlevel-voice+3,"
-// 'Spoly1' created a ring as 2nd element of the list.
-// The ring contains the ideal 'kern'  of algebraic relations between the
-//generators of I+J.
-// To access to this ring and see 'kern' you should give the ring a name,
-// e.g.:
-               def @ring = L[2]; setring @ring ; kern;
-      ");
-    }
-    if(a==2)
-    {
-      @L=P1,R1;
-      kill phi,vars;
-    }
-  }
-  return(@L);
+  int 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))));
+  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,")");
+  }
+  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;
+  //         }
+
+
+  ideal I=fetch(br,I);
+  ideal algebra;
+  for (i=1; i<=ncols(I);i++)
+  {
+    algebra[i]=I[i]-par(i);
+  }
+  algebra=sagbiConstruction(algebra, iterations,0,0,1);
+  int j=1;
+  ideal algDep;
+  for (i=1; i<= ncols(algebra); i++)
+  {
+    if (leadmonom(algebra[i])==1)
+    {
+      algDep[j]=algebra[i];
+      j++;
+    }
+  }
+  setring extendVarRing;
+  ideal algDep=imap(parameterRing,algDep);
+  ideal algebra=imap(parameterRing,algebra);
+  export algDep,algebra;
+  //print(algDep);
+  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 = algDep(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 = algDep(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)
+{
+  ideal J,B;
+  int i,j,k;
+  poly f;
+  for(k=1;k<=ncols(I);k++)
+  {
+    f=I[k];
+    I[k]=0;
+    f=sagbiReduce(f,I,1);
+    I[k]=f;
+  }
+  I=simplify(I,2);
+  return(I);
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -386 +1036 @@
 }
 
-///////////////////////////////////////////////////////////////////////////////
-static proc completeReduction(poly p,ideal dom,list#)//reduction
-{
-  poly p1=p;
-  poly p2=sagbiReduction(p,dom,#);
-  while (p1!=p2)
-  {
-    p1=p2;
-    p2=sagbiReduction(p1,dom,#);
-  }
-  return(p2);
-}
-///////////////////////////////////////////////////////////////////////////////
-
-static proc completeReduction1(poly p,ideal dom,list #) //tail reduction
-{
-  poly p1,p2,re;
-  p1=p;
-  while(p1!=0)
-  {
-    p2=sagbiReduction(p1,dom,#);
-    if(p2!=p1)
-    {
-      p1=p2;
-    }
-    else
-    {
-      re=re+lead(p2);
-      p1=p2-lead(p2);
-    }
-  }
-  return(re);
-}
-
-
-
-///////////////////////////////////////////////////////////////////////////////
-
 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.
@@ -439 +1051 @@
 EXAMPLE: example sagbiNF; show example "
 {
-  int z;
-  ideal Red;
-  poly re;
-  if(typeof(id)=="ideal")
-  {
-    int i=ncols(id);
-    for(z=1;z<=i;z++)
-    {
-      if(k==0)
-      {
-        id[z]=completeReduction(id[z],dom,#);
-      }
-      else
-      {
-        id[z]=completeReduction1(id[z],dom,#);//tail reduction.
-      }
-    }
-    Red=simplify(id,7);
-    return(Red);
-  }
-  if(typeof(id)=="poly")
-  {
-    if(k==0)
-    {
-      re=completeReduction(id,dom,#);
-    }
-    else
-    {
-      re=completeReduction1(id,dom,#);
-    }
-    return(re);
-  }
+   return(sagbiReduce(id,dom,k));
 }
 example
@@ -483 +1064 @@
 }
 
-
-///////////////////////////////////////////////////////////////////////////////
-
-static proc intRed(id,int k, list #)
-{
-  int i,z;
-  ideal Rest,intRed;
-  z=ncols(id);
-  for(i=1;i<=z;i++)
-  {
-    Rest=id;
-    Rest[i]=0;
-    Rest=simplify(Rest,2);
-    if(k==0)
-    {
-      intRed[i]=completeReduction(id[i],Rest,#);
-    }
-    else
-    {
-      intRed[i]=completeReduction1(id[i],Rest,#);
-    }
-  }
-  intRed=simplify(intRed,7);//1+2+4 in simplify command
-  return(intRed);
-}
-//////////////////////////////////////////////////////////////////////////////
-
-proc sagbi(id,int k,list#)
-"USAGE: sagbi(id,k[,n]); id ideal, k and n positive integers.
-RETURN: A SAGBI basis for the subalgebra defined by the generators of id.
-@format
-    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, and S-interreduced SAGBI basis
-       is returned.
-    Three algorithm variants are used to perform subalgebra reduction.
-    The positive interger n determine which variant should be used.
-    n may take the values (0 or default),1 or 2.
-@end format
-NOTE: SAGBI bases computations may be performed in polynomial rings or quotients
-      thereof.
-EXAMPLE: example sagbi; show example "
-{
-  degBound=0;
-  ideal S,oldS,Red;
-  list L;
-  S=intRed(id,k,#);
-  while(size(S)!=size(oldS))
-  {
-    L=Spoly1(L,S,Red,2);
-    Red=L[1];
-    Red=sagbiNF(Red,S,k,#);
-    oldS=S;
-    S=S+Red;
-  }
-  return(interreduced(S));
-}
-example
-{ "EXAMPLE:"; echo = 2;
- ring r= 0,(x,y),dp;
- ideal I=x2,y2,xy+y;
- sagbi(I,1,1);
-}
-
-proc interreduced(ideal I)
-{
-  ideal J,B;
-  int i,j,k;
-  poly f;
-  for(k=1;k<=ncols(I);k++)
-  {
-    f=I[k];
-    I[k]=0;
-    f=sagbiNF(f,I,1);
-    I[k]=f;
-  }
-  I=simplify(I,2);
-  return(I);
-}
-
-///////////////////////////////////////////////////////////////////////////////
-proc sagbiPart(id,int k,int c,list #)
-"USAGE: sagbiPart(id,k,c[,n]); id ideal, k, c and n positive integers
-RETURN: A partial SAGBI basis for the subalgebra defined by the generators of id.
-@format
-     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, and S-intereduced SAGBI basis
-      is returned.
-     c determines, after how many loops the Sagbi basis computation should stop.
-     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:- SAGBI bases computations may be performed in polynomial rings or quotients
-       thereof.
-     - This version of sagbi is interesting in the case of subalgebras with infinte
-       SAGBI basis. In this case, it may be used to check, if the elements of this
-       basis have a particular form.
-EXAMPLE: example sagbiPart; show example "
-{
-  degBound=0;
-  ideal S,oldS,Red;
-  int counter;
-  list L;
-  S=intRed(id,k,#);
-  while((size(S)!=size(oldS))&&(counter<=c))
-  {
-    L=Spoly1(L,S,Red,2);
-    Red=L[1];
-    Red=sagbiNF(Red,S,k,#);
-    oldS=S;
-    S=S+Red;
-    counter=counter+1;
-  }
-  return(S);
-}
-example
-{ "EXAMPLE:"; echo = 2;
- ring r= 0,(x,y),dp;
- ideal I=x,xy-y2,xy2;//the corresponding Subalgebra has an infinte SAGBI basis
- sagbiPart(I,1,3);// computations should stop after 3 turns.
-}
-//////////////////////////////////////////////////////////////////////////////
+/*
+  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 = algDep(J,2);
+  setring DI; algDep;
+*/