////////////////////////////////////////////////////////////////////////////// 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 ------------- setring r; 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 chosen --------------- 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 * \ 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)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; */