// $Id: toric.lib,v 1.2 2000-05-11 09:03:41 Singular Exp $ // // author : Christine Theis // //version="$Id: toric.lib,v 1.2 2000-05-11 09:03:41 Singular Exp $"; /////////////////////////////////////////////////////////////////////////////// info=" LIBRARY: toric.lib PROCEDURES FOR COMPUTING TORIC IDEALS Let A an integral (mxn)-matrix. The toric ideal I_A of A is defined as the ideal I_A:=< x^u - x^v | u,v integral and nonnegative, u-v in the kernel of A > in the ring of n variables x:=x1,...,xn. Toric ideals play an important role in polyhedral geometry and may also be used for integer programming. They are generated by binomials with relatively prime monomials. Buchberger's algorithm can be specialized to these structures in a way that considerably speeds up computation. toric_ideal(intmat A, string alg); toric_ideal(intmat A, string alg, intvec prsv); procedures for computing the toric ideal of A They return the standard basis of the toric ideal of A with respect to the term ordering in the actual basering. When calling this procedure, a ring with n variables should be active. Not all term orderings are supported: The usual global term orderings may be used, but no block orderings combining them. One may call the procedure with several different algorithms: - the algorithm of Conti/Traverso using elimination (ect), - the algorithm of Pottier (pt), - an algorithm of Bigatti/La Scala/Robbiano (blr), - the algorithm of Hosten/Sturmfels (hs), - the algorithm of DiBiase/Urbanke (du). The last two seem to be the fastest in the actual implementation. `alg' should be the abbreviation (in brackets) for an algorithm as above. If `alg' is chosen to be `blr' or `hs', the algorithm needs a vector with positive coefficcients in the row space of A. If no row of A contains only positive entries, one must use the second version of toric_ideal which takes such a vector in the third argument. toric_std(ideal I); computes the standard basis of I using the specialized Buchberger algorithm. The generating system by which I is given has to consist of binomials of the form x^u-x^v (although there are other generating systems of toric ideals). There is no real check if I is toric. If the generator list of I contains a binomial whose monomials are not relatively prime, the procedure outputs a warning. If I is generated by binomials of the above form, but not toric, toric_std computes an ideal `between' I and its saturation with respect to all variables. "; /////////////////////////////////////////////////////////////////////////////// static proc toric_ideal_1(intmat A, string alg) { ideal I; // to be returned // check suitability of actual basering if(nvars(basering)0) { break; } } if(found==0) { "ERROR: algorithm needs positive vector in the row space of the matrix"; close(MATRIX); dummy=system("sh","rm -f "+matrixfile); return(I); } write(MATRIX,"positive row space vector:"); for(j=1;j<=ncols(A);j++) { write(MATRIX,A[found,j]); } } close(MATRIX); // call external program dummy=system("sh","toric_ideal -alg "+alg+" "+matrixfile); // read toric ideal from created file link TORIC_IDEAL=":r "+matrixfile+".GB."+alg; string toric_ideal=read(TORIC_IDEAL); int generators; pos=find(toric_ideal,"size"); pos=find(toric_ideal,":",pos); pos++; while(toric_ideal[pos]==" " || toric_ideal[pos]==newline) { pos++; } number_string=""; while(toric_ideal[pos]!=" " && toric_ideal[pos]!=newline) { number_string=number_string+toric_ideal[pos]; pos++; } execute("generators="+number_string+";"); intvec v; poly head; poly tail; pos=find(toric_ideal,"basis"); pos=find(toric_ideal,":",pos); pos++; for(i=1;i<=generators;i++) { head=1; tail=1; for(j=1;j<=ncols(A);j++) { while(toric_ideal[pos]==" " || toric_ideal[pos]==newline) { pos++; } number_string=""; while(toric_ideal[pos]!=" " && toric_ideal[pos]!=newline) { number_string=number_string+toric_ideal[pos]; pos++; } execute("v[j]="+number_string+";"); if(v[j]<0) { tail=tail*var(j)^(-v[j]); } if(v[j]>0) { head=head*var(j)^v[j]; } } I[i]=head-tail; } // delete all created files dummy=system("sh","rm -f "+matrixfile); dummy=system("sh","rm -f "+matrixfile+".GB."+alg); return(I); } static proc toric_ideal_2(intmat A, string alg, intvec prsv) { ideal I; // to be returned // check arguments if(size(prsv)0) { head=head*var(j)^v[j]; } } I[i]=head-tail; } // delete all created files dummy=system("sh","rm -f "+matrixfile); dummy=system("sh","rm -f "+matrixfile+".GB."+alg); return(I); } proc toric_ideal "USAGE: toric_ideal(A,alg); A intmat, alg string toric_ideal(A,alg,prsv); A intmat, alg string, prsv intvec RETURN: toric ideal of A as explained in toric.lib return type = ideal EXAMPLE: example toric_ideal; shows an example" { if(size(#)==2) { return(toric_ideal_1(#[1],#[2])); } else { return(toric_ideal_2(#[1],#[2],#[3])); } } example { "EXAMPLE"; echo=2; ring r=0,(x,y,z),dp; // call with two arguments intmat A[2][3]=1,1,0,0,1,1; A; ideal I=toric_ideal(A,"du"); I; I=toric_ideal(A,"blr"); I; // call with three arguments intvec prsv=1,2,1; I=toric_ideal(A,"blr",prsv); I; } proc toric_std(ideal I) "USAGE: toric_std(I); I ideal RETURN: standard basis of I as explained in toric.lib return type = ideal EXAMPLE: example toric_std; shows an example" { ideal J; // to be returned // check suitability of actual term ordering // the "module ordering" c or C is ignored string singular_ord=ordstr(basering); string test_ord; string external_ord=""; int i,j; intvec weightvec; if(find(singular_ord,"lp")==1) { external_ord="W_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=0; } test_ord="lp("+string(nvars(basering))+"),"; if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c")) { "Warning: block orderings are not supported; lp used for computation"; } } if(external_ord=="" && find(singular_ord,"lp")==3) { external_ord="W_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=0; } test_ord="lp("+string(nvars(basering))+"),"; if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord)) { "Warning: block orderings are not supported; lp used for computation"; } } if(external_ord=="" && find(singular_ord,"dp")==1) { external_ord="W_REV_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=1; } test_ord="dp("+string(nvars(basering))+"),"; if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c")) { "Warning: block orderings are not supported; dp used for computation"; } } if(external_ord=="" && find(singular_ord,"dp")==3) { external_ord="W_REV_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=1; } test_ord="dp("+string(nvars(basering))+"),"; if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord)) { "Warning: block orderings are not supported; dp used for computation"; } } if(external_ord=="" && find(singular_ord,"Dp")==1) { external_ord="W_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=1; } test_ord="Dp("+string(nvars(basering))+"),"; if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c")) { "Warning: block orderings are not supported; Dp used for computation"; } } if(external_ord=="" && find(singular_ord,"Dp")==3) { external_ord="W_LEX"; for(i=1;i<=nvars(basering);i++) { weightvec[i]=1; } test_ord="Dp("+string(nvars(basering))+"),"; if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord)) { "Warning: block orderings are not supported; Dp used for computation"; } } int pos; string number_string; if(external_ord=="" && find(singular_ord,"wp")==1) { external_ord="W_REV_LEX"; pos=3; for(i=1;i<=nvars(basering);i++) { pos++; number_string=""; while(singular_ord[pos]!="," && singular_ord[pos]!=")") { number_string=number_string+singular_ord[pos]; pos++; } execute("weightvec["+string(i)+"]="+number_string+";"); } test_ord="wp("+string(weightvec)+"),"; if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c")) { "Warning: block orderings are not supported; wp("+string(weightvec)+") used for computation"; } } if(external_ord=="" && find(singular_ord,"wp")==3) { external_ord="W_REV_LEX"; pos=5; for(i=1;i<=nvars(basering);i++) { pos++; number_string=""; while(singular_ord[pos]!=",") { number_string=number_string+singular_ord[pos]; pos++; } execute("weightvec[i]="+number_string); } test_ord="wp("+string(weightvec)+"),"; if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord)) { "Warning: block orderings are not supported; wp("+string(weightvec)+") used for computation"; } } if(external_ord=="" && find(singular_ord,"Wp")==1) { external_ord="W_LEX"; pos=3; for(i=1;i<=nvars(basering);i++) { pos++; number_string=""; while(singular_ord[pos]!=",") { number_string=number_string+singular_ord[pos]; pos++; } execute("weightvec[i]="+number_string); } test_ord="Wp("+string(weightvec)+"),"; if(singular_ord!=(test_ord+"C") && singular_ord!=(test_ord+"c")) { "Warning: block orderings are not supported; Wp("+string(weightvec)+") used for computation"; } } if(external_ord=="" && find(singular_ord,"Wp")==3) { external_ord="W_LEX"; pos=5; for(i=1;i<=nvars(basering);i++) { pos++; number_string=""; while(singular_ord[pos]!=",") { number_string=number_string+singular_ord[pos]; pos++; } execute("weightvec[i]="+number_string); } test_ord="Wp("+string(weightvec)+"),"; if(singular_ord!=("C"+test_ord) && singular_ord!=("c"+test_ord)) { "Warning: block orderings are not supported; Wp("+string(weightvec)+") used for computation"; } } if(external_ord=="") { "ERROR: term ordering of actual basering not supported"; return(I); } // create first temporary file with which the external program is called int dummy; int process=system("pid"); string groebnerfile="temp_GROEBNER"+string(process); link GROEBNER=":w "+groebnerfile; open(GROEBNER); write(GROEBNER,"GROEBNER","computed with algorithm:","pt","term ordering:","elimination block",0,"weighted block",nvars(basering),external_ord); // algorithm is totally unimportant, only required by the external program for(i=1;i<=nvars(basering);i++) { write(GROEBNER,weightvec[i]); } write(GROEBNER,"size:",size(I),"Groebner basis:"); poly head; poly tail; poly rest; intvec v; for(j=1;j<=size(I);j++) { // test suitability of generator j rest=I[j]; head=lead(rest); rest=rest-head; tail=lead(rest); rest=rest-tail; if(head==0 && tail==0 && rest!=0) { "ERROR: generator "+string(j)+" of input ideal is no binomial"; close(GROEBNER); dummy=system("sh","rm -f "+groebnerfile); return(J); } if(leadcoef(tail)!=-leadcoef(head)) // generator no difference of monomials (or a constant multiple) { "ERROR: generator "+string(j)+" of input ideal is no difference of monomials"; close(GROEBNER); dummy=system("sh","rm -f "+groebnerfile); return(J); } if(gcd(head,tail)!=1) { "Warning: monomials of generator "+string(j)+" of input ideal are not relatively prime"; } // write vector representation of generator j into the file v=leadexp(head)-leadexp(tail); for(i=1;i<=nvars(basering);i++) { write(GROEBNER,v[i]); } } close(GROEBNER); // create second temporary file string newcostfile="temp_NEW_COST"+string(process); link NEW_COST=":w "+newcostfile; open(NEW_COST); write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:"); for(i=1;i<=nvars(basering);i++) { write(NEW_COST,weightvec[i]); } // call external program dummy=system("sh","change_cost "+groebnerfile+" "+newcostfile); // read toric standard basis from created file link TORIC_IDEAL=":r "+newcostfile+".GB.pt"; string toric_ideal=read(TORIC_IDEAL); int generators; pos=find(toric_ideal,"size"); pos=find(toric_ideal,":",pos); pos++; while(toric_ideal[pos]==" " || toric_ideal[pos]==newline) { pos++; } number_string=""; while(toric_ideal[pos]!=" " && toric_ideal[pos]!=newline) { number_string=number_string+toric_ideal[pos]; pos++; } execute("generators="+number_string+";"); pos=find(toric_ideal,"basis"); pos=find(toric_ideal,":",pos); pos++; for(j=1;j<=generators;j++) { head=1; tail=1; for(i=1;i<=nvars(basering);i++) { while(toric_ideal[pos]==" " || toric_ideal[pos]==newline) { pos++; } number_string=""; while(toric_ideal[pos]!=" " && toric_ideal[pos]!=newline) { number_string=number_string+toric_ideal[pos]; pos++; } execute("v[i]="+number_string+";"); if(v[i]<0) { tail=tail*var(i)^(-v[i]); } if(v[i]>0) { head=head*var(i)^v[i]; } } J[j]=head-tail; } // delete all created files dummy=system("sh","rm -f "+groebnerfile); dummy=system("sh","rm -f "+groebnerfile+".GB.pt"); dummy=system("sh","rm -f "+newcostfile); return(J); } example { "EXAMPLE"; echo=2; ring r=0,(x,y,z),wp(3,2,1); // call with toric ideal (of the matrix A=(1,1,1)) ideal I=x-y,x-z; ideal J=toric_std(I); J; // call with the same ideal, but badly chosen generators: // not only binomials I=x-y,2x-y-z; J=toric_std(I); // binomials whose monomials are not relatively prime I=x-y,xy-yz,y-z; J=toric_std(I); J; // call with a non-toric ideal that seems to be toric I=x-yz,xy-z; J=toric_std(I); J; // comparison with real standard basis and saturation ideal H=std(I); H; LIB "elim.lib"; sat(H,xyz); } //////////////////////////////////////////////////////////////////////////////