// $Id: standard.lib,v 1.38 1999-07-07 14:39:14 obachman Exp $ ////////////////////////////////////////////////////////////////////////////// version="$Id: standard.lib,v 1.38 1999-07-07 14:39:14 obachman Exp $"; info=" LIBRARY: standard.lib PROCEDURES WHICH ARE ALWAYS LOADED AT START-UP PROCEDURES: stdfglm(ideal[,ord]) standard basis of ideal via fglm [and ordering ord] stdhilb(ideal[,h]) standard basis of ideal using the Hilbert function groebner(ideal/module) standard basis using a heuristically choosen method quot(any,any[,n]) quotient using heuristically choosen method sprintf(fmt,...) returns fomatted string fprintf(link,fmt,..) writes formatted string to link printf(fmt,...) displays formatted string "; ////////////////////////////////////////////////////////////////////////////// proc stdfglm (ideal i, list #) "USAGE: stdfglm(i[,s]); i ideal, s string (any allowed ordstr of a ring) RETURN: stdfglm(i): standard basis of i in the basering, calculated via fglm from ordering \"dp\" to the ordering of the basering. stdfglm(i,s): standard basis of i in the basering, calculated via fglm from ordering s to the ordering of the basering. SEE ALSO: stdhilb, std, groebner EXAMPLE: example stdfglm; shows an example" { string os; def dr= basering; if( (size(#)==0) or (typeof(#[1]) != "string") ) { os = "dp(" + string( nvars(dr) ) + ")"; if ( (find( ordstr(dr), os ) != 0) and (find( ordstr(dr), "a") == 0) ) { os= "Dp"; } else { os= "dp"; } } else { os = #[1]; } execute "ring sr=("+charstr(dr)+"),("+varstr(dr)+"),"+os+";"; ideal i= fetch(dr,i); intvec opt= option(get); option(redSB); i=std(i); option(set,opt); setring dr; return (fglm(sr,i)); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; ideal i = y3+x2, x2y+x2, x3-x2, z4-x2-y; ideal i1= stdfglm(i); //uses fglm from "dp" to "lp" i1; ideal i2= stdfglm(i,"Dp"); //uses fglm from "Dp" to "lp" i2; } ///////////////////////////////////////////////////////////////////////////// proc stdhilb(ideal i,list #) "USAGE: stdhilb(i); i ideal stdhilb(i,v); i homogeneous ideal, v intvec (the Hilbert function) RETURN: stdhilb(i): a standard basis of i (computing v internally) stdhilb(i,v): standard basis of i, using the given Hilbert function SEE ALSO: stdfglm, std, groebner EXAMPLE: example stdhilb; shows an example" { def R=basering; if((homog(i)==1)||(ordstr(basering)[1]=="d")) { if ((size(#)!=0)&&(homog(i)==1)) { return(std(i,#[1])); } return(std(i)); } execute "ring S = ("+charstr(R)+"),("+varstr(R)+",@t),dp;"; ideal i=homog(imap(R,i),@t); intvec v=hilb(std(i),1); execute "ring T = ("+charstr(R)+"),("+varstr(R)+",@t),("+ordstr(R)+");"; ideal i=fetch(S,i); ideal a=std(i,v); setring R; map phi=T,maxideal(1),1; ideal a=phi(a); int k,j; poly m; int c=size(i); for(j=1;j0) { m=lead(a[j]); for(k=j+1;k<=c;k++) { if(size(lead(a[k])/m)>0) { a[k]=0; } } } } a=simplify(a,2); attrib(a,"isSB",1); return(a); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),lp; ideal i = y3+x2, x2y+x2, x3-x2, z4-x2-y; ideal i1= stdhilb(i); i1; // is in this case equivalent to: intvec v=1,0,0,-3,0,1,0,3,-1,-1; ideal i2=stdhilb(i,v); } ////////////////////////////////////////////////////////////////////////// proc groebner(def i, list #) "USAGE: groebner(i[, wait]) i -- ideal/module; wait -- int RETURNS: Standard basis of ideal or module which is computed using a heuristically choosen method: If the ordering of the current ring is a local ordering, or if it is a non-block ordering and the current ring has no parameters, then std(i) is returned. Otherwise, i is mapped into a ring with no parameters and ordering dp, where its Hilbert series is computed. This is followed by a Hilbert-series based std computation in the original ring. NOTE: If a 2nd argument 'wait' is given, then the computation proceeds at most 'wait' seconds. That is, if no result could be computed in 'wait' seconds, then the computation is interrupted, 0 is returned, a warning message is displayed, and the global variable 'groebner_error' is defined. SEE ALSO: stdhilb, stdfglm, std EXAMPLE: example groebner; shows an example" { def P=basering; // we have two arguments -- try to use MPfork links if (size(#) > 0) { if (system("with", "MP")) { if (typeof(#[1]) == "int") { int wait = #[1]; int j = 10; string bs = nameof(basering); link l_fork = "MPtcp:fork"; open(l_fork); write(l_fork, quote(system("pid"))); int pid = read(l_fork); write(l_fork, quote(groebner(eval(i)))); // sleep in small intervalls for appr. one second if (wait > 0) { while(j < 1000000) { if (status(l_fork, "read", "ready", j)) {break;} j = j + j; } } // sleep in intervalls of one second from now on j = 1; while (j < wait) { if (status(l_fork, "read", "ready", 1000000)) {break;} j = j + 1; } if (status(l_fork, "read", "ready")) { def result = read(l_fork); if (bs != nameof(basering)) { def PP = basering; setring P; def result = imap(PP, result); kill PP; } if (defined(groebner_error)) { kill(groebner_error); } kill (l_fork); } else { ideal result; if (! defined(groebner_error)) { int groebner_error = 1; export groebner_error; } "// ** groebner did not finish"; j = system("sh", "kill " + string(pid)); } return (result); } else { "// ** groebner needs int as 2nd arg"; } } else { "// ** groebner with two args is not supported in this configuration"; } } // we are still here -- do the actual computation string ordstr_P = ordstr(P); if (find(ordstr_P,"s") > 0) { //spaeter den lokalen fall ueber lp oder aehnlich behandeln return(std(i)); } int IsSimple_P; if (system("nblocks") <= 2) { if (find(ordstr_P, "M") <= 0) { IsSimple_P = 1; } } int npars_P = npars(P); // return std if no parameters and (dp or wp) if ((npars_P <= 1) && IsSimple_P) { if (find(ordstr_P, "d") > 0) { return (std(i)); } if (find(ordstr_P,"w") > 0) { return (std(i)); } } // reset options intvec opt=option(get); int p_opt; string s_opt = option(); option(none); // turn on option(prot) and/or option(mem), if previously set if (find(s_opt, "prot")) { option(prot); p_opt = 1; } if (find(s_opt, "mem")) { option(mem); } // construct ring in which first std computation is done string varstr_P = varstr(P); string parstr_P = parstr(P); int is_homog = (homog(i) && (npars_P <= 1)); int add_vars = 0; string ri = "ring Phelp ="; // more than one parameters are converted to ring variables if (npars_P > 1) { ri = ri + string(char(P)) + ",(" + varstr_P + "," + parstr_P; add_vars = npars_P; } else { ri = ri + "(" + charstr(P) + "),(" + varstr_P; } // a homogenizing variable is added, if necessary if (! is_homog) { ri = ri + ",@t"; add_vars = add_vars + 1; } // ordering is set to (dp, C) ri = ri + "),(dp,C);"; // change the ring execute(ri); // get ideal from previous ring if (is_homog) { ideal qh = imap(P, i); } else { // and homogenize ideal qh=homog(imap(P,i),@t); } // compute std and hilbert series if (p_opt) { "std in " + ri[13, size(ri) - 13]; } ideal qh1=std(qh); intvec hi=hilb(qh1,1); if (add_vars == 0) { // no additional variables were introduced setring P; // can immediately change to original ring // simply compute std with hilbert series in original ring if (p_opt) { "std with hilb in basering"; } i = std(i, hi); } else { // additional variables were introduced // need another intermediate ring ri = "ring Phelp1 = (" + charstr(Phelp) + "),(" + varstr(Phelp) + "),(" + ordstr_P; // for lp wit at most one parameter, we do not need a block ordering if ( ! (IsSimple_P && (add_vars <2) && find(ordstr_P, "l"))) { // need block ordering ri = ri + ", dp(" + string(add_vars) + ")"; } ri = ri + ");"; // change to intermediate ring execute(ri); ideal qh = imap(Phelp, qh); kill Phelp; if (p_opt) { "std with hilb in " + ri[14,size(ri)-14]; } // compute std with Hilbert series qh = std(qh, hi); // subst 1 for homogenizing var if (!is_homog) { if (p_opt) { "dehomogenization"; } qh = subst(qh, @t, 1); } // go back to original ring setring P; // get ideal, delete zeros and clean SB if (p_opt) { "imap to original ring"; } i = imap(Phelp1,qh); if (p_opt) { "simplification"; } i = simplify(i, 34); kill Phelp1; } // clean-up time option(set, opt); if (find(s_opt, "redSB") > 0) { if (p_opt) { "interreduction"; } i=interred(i); } attrib(i, "isSB", 1); return (i); } example { "EXAMPLE: "; echo = 2; ring r = 0, (a,b,c,d), lp; option(prot); ideal i = a+b+c+d, ab+ad+bc+cd, abc+abd+acd+bcd, abcd-1; // cyclic 4 groebner(i); ring rp = (0, a, b), (c,d), lp; ideal i = imap(r, i); ideal j = groebner(i); option(noprot); j; simplify(j, 1); std(i); if (system("with", "MP")) {groebner(i, 0);} defined(groebner_error); } ////////////////////////////////////////////////////////////////////////// proc res(list #) { def P=basering; def m=#[1]; //the ideal or module int i=#[2]; //the length of the resolution //if size(#)>2 a minimal resolution is computed string varstr_P = varstr(P); if(size(ideal(basering)) > 0) { // the quick hack for qrings - seems to fit most needs // (lres is not implemented for qrings, sres is not so efficient) return(nres(m,i)); } //LaScala for the homogeneous case if(homog(m)==1) { resolution re; if ((i==0) or (i>=nvars(basering))) { re=lres(m,i); if(size(#)>2) { re=minres(re); } } else { if(size(#)>2) { re=mres(m,i); } else { re=sres(std(m),i); } } return(re); } //mres for the global non homogeneous case if(find(ordstr(P),"s")==0) { string ri= "ring Phelp =" +string(char(P))+",("+varstr_P+"),(dp,C);"; execute(ri); def m=imap(P,m); list re=mres(m,i); setring P; resolution result=imap(Phelp,re); return(result); } //sres for the local case and not minimal resolution if(size(#)<=2) { string ri= "ring Phelp =" +string(char(P))+",("+varstr_P+"),(ls,c);"; execute(ri); def m=imap(P,m); m=std(m); list re=sres(m,i); setring P; resolution result=imap(Phelp,re); return(result); } //mres for the local case and minimal resolution string ri= "ring Phelp =" +string(char(P))+",("+varstr_P+"),(ls,C);"; execute(ri); def m=imap(P,m); list re=mres(m,i); setring P; resolution result=imap(Phelp,re); return(result); } proc quot (m1,m2,list #) "USAGE: quot(m1, m2[, n]); m1, m2 two submodules of k^s, n (optional) integer (1<= n <=5) RETURN: the quotient of m1 and m2 SEE ALSO: quotient EXAMPLE: example quot; shows an example" { if (((typeof(m1)!="ideal") and (typeof(m1)!="module")) or ((typeof(m2)!="ideal") and (typeof(m2)!="module"))) { "USAGE: quot(m1, m2[, n]); m1, m2 two submodules of k^s,"; " n (optional) integer (1<= n <=5)"; "RETURN: the quotient of m1 and m2"; "EXAMPLE: example quot; shows an example"; return(); } if (typeof(m1)!=typeof(m2)) { return(quotient(m1,m2)); } if (size(#)>0) { if (typeof(#[1])=="int" ) { return(quot1(m1,m2,#[1])); } } else { return(quot1(m1,m2,2)); } } example { "EXAMPLE:"; echo = 2; ring r=181,(x,y,z),(c,ls); ideal id1=maxideal(4); ideal id2=x2+xyz,y2-z3y,z3+y5xz; option(prot); ideal id6=quotient(id1,id2); id6; ideal id7=quot(id1,id2,1); id7; ideal id8=quot(id1,id2,2); id8; } static proc quot1 (module m1, module m2,int n) "USAGE: quot1(m1, m2, n); m1, m2 two submodules of k^s, n integer (1<= n <=5) RETURN: the quotient of m1 and m2 EXAMPLE: example quot1; shows an example" { if (n==1) { return(quotient1(m1,m2)); } else { if (n==2) { return(quotient2(m1,m2)); } else { if (n==3) { return(quotient3(m1,m2)); } else { if (n==4) { return(quotient4(m1,m2)); } else { if (n==5) { return(quotient5(m1,m2)); } else { return(quotient(m1,m2)); } } } } } } example { "EXAMPLE:"; echo = 2; ring r=181,(x,y,z),(c,ls); ideal id1=maxideal(4); ideal id2=x2+xyz,y2-z3y,z3+y5xz; option(prot); ideal id6=quotient(id1,id2); id6; ideal id7=quot1(id1,id2,1); id7; ideal id8=quot1(id1,id2,2); id8; } static proc quotient0(module a,module b) { module mm=b+a; resolution rs=lres(mm,0); list I=list(rs); matrix M=I[2]; matrix A[1][nrows(M)]=M[1..nrows(M),1]; ideal i=A; return (i); } proc quotient1(module a,module b) //17sec "USAGE: quotient1(m1, m2); m1, m2 two submodules of k^s, RETURN: the quotient of m1 and m2" { int i; a=std(a); module dummy; module B=NF(b,a)+dummy; ideal re=quotient(a,module(B[1])); for(i=2;i<=size(B);i++) { re=intersect1(re,quotient(a,module(B[i]))); } return(re); } proc quotient2(module a,module b) //13sec "USAGE: quotient2(m1, m2); m1, m2 two submodules of k^s, RETURN: the quotient of m1 and m2" { a=std(a); module dummy; module bb=NF(b,a)+dummy; int i=size(bb); ideal re=quotient(a,module(bb[i])); bb[i]=0; module temp; module temp1; module bbb; int mx; i=i-1; while (1) { if (i==0) break; temp = a+bb*re; temp1 = lead(interred(temp)); mx=ncols(a); if (ncols(temp1)>ncols(a)) { mx=ncols(temp1); } temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx); temp1 = dummy+temp1; if (deg(temp1[1])<0) break; re=intersect1(re,quotient(a,module(bb[i]))); bb[i]=0; i = i-1; } return(re); } proc quotient3(module a,module b) //89sec "USAGE: quotient3(m1, m2); m1, m2 two submodules of k^s, only for global rings RETURN: the quotient of m1 and m2" { string s="ring @newr=("+charstr(basering)+ "),("+varstr(basering)+",@t,@w),dp;"; def @newP=basering; execute s; module b=imap(@newP,b); module a=imap(@newP,a); int i; int j=size(b); vector @b; for(i=1;i<=j;i++) { @b=@b+@t^(i-1)*@w^(j-i+1)*b[i]; } ideal re=quotient(a,module(@b)); setring @newP; ideal re=imap(@newr,re); return(re); } proc quotient5(module a,module b) //89sec "USAGE: quotient5(m1, m2); m1, m2 two submodules of k^s, only for global rings RETURN: the quotient of m1 and m2" { string s="ring @newr=("+charstr(basering)+ "),("+varstr(basering)+",@t),dp;"; def @newP=basering; execute s; module b=imap(@newP,b); module a=imap(@newP,a); int i; int j=size(b); vector @b; for(i=1;i<=j;i++) { @b=@b+@t^(i-1)*b[i]; } @b=homog(@b,@w); ideal re=quotient(a,module(@b)); setring @newP; ideal re=imap(@newr,re); return(re); } proc quotient4(module a,module b) //95sec "USAGE: quotient4(m1, m2); m1, m2 two submodules of k^s, only for global rings RETURN: the quotient of m1 and m2" { string s="ring @newr=("+charstr(basering)+ "),("+varstr(basering)+",@t),dp;"; def @newP=basering; execute s; module b=imap(@newP,b); module a=imap(@newP,a); int i; vector @b=b[1]; for(i=2;i<=size(b);i++) { @b=@b+@t^(i-1)*b[i]; } matrix sy=modulo(@b,a); ideal re=sy; setring @newP; ideal re=imap(@newr,re); return(re); } static proc intersect1(ideal i,ideal j) { def R=basering; execute "ring gnir = ("+charstr(basering)+"), ("+varstr(basering)+",@t),(C,dp);"; ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j); ideal j=eliminate(i,var(nvars(basering))); setring R; map phi=gnir,maxideal(1); return(phi(j)); } ////////////////////////////////////////////////////////////////// /// /// sprintf, fprintf printf /// proc sprintf(string fmt, list #) "USAGE: sprintf(fmt, ...) fmt string RETURN: string PURPOSE: sprintf performs output formatting. The first argument is a format control string. Additional arguments may be required, depending on the contents of the control string. A series of output characters is generated as directed by the control string; these characters are returned as a string. The control string is simply text to be copied, except that the string may contain conversion specifications. Do 'help print:' for a listing of valid conversion specifications. As an addition to the conversions of 'print', the '%n' and '%2' conversion specification does not consume an additional argument, but simply generates a newline character. NOTE: If one of the additional arguments is a list, then it should be enclosed once more into a list() command, since passing a list as an argument flattens the list by one level. SEE ALSO: fprintf, printf, print, string EXAMPLE : example sprintf; shows an example " { int sfmt = size(fmt); if (sfmt <= 1) { return (fmt); } int next, l, nnext; string ret; list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2"; while (1) { if (size(#) <= 0) { return (ret + fmt); } nnext = 0; while (nnext < sfmt) { nnext = find(fmt, "%", nnext + 1); if (nnext == 0) { next = 0; break; } l = 1; while (l <= size(formats)) { next = find(fmt, formats[l], nnext); if (next == nnext) break; l++; } if (next == nnext) break; } if (next == 0) { return (ret + fmt); } if (formats[l] != "%2" && formats[l] != "%n") { ret = ret + fmt[1, next - 1] + print(#[1], formats[l]); # = delete(#, 1); } else { ret = ret + fmt[1, next - 1] + print("", "%2s"); } if (size(fmt) <= (next + size(formats[l]) - 1)) { return (ret); } fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1]; } } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; module m=[1,y],[0,x+z]; intmat M=betti(mres(m,0)); list l = r, m, M; string s = sprintf("s:%s,%n l:%l", 1, 2); s; s = sprintf("s:%n%s", l); s; s = sprintf("s:%2%s", list(l)); s; s = sprintf("2l:%n%2l", list(l)); s; s = sprintf("%p", list(l)); s; s = sprintf("%;", list(l)); s; s = sprintf("%b", M); s; } proc printf(string fmt, list #) "USAGE: printf(fmt, ...) fmt string RETURN: none PURPOSE: printf performs output formatting. The first argument is a format control string. Additional arguments may be required, depending on the contents of the control string. A series of output characters is generated as directed by the control string; these characters are displayed (i.e. printed to standard out). The control string is simply text to be copied, except that the string may contain conversion specifications. Do 'help print:' for a listing of valid conversion specifications. As an addition to the conversions of 'print', the '%n' and '%2' conversion specification does not consume an additional argument, but simply generates a newline character. NOTE: If one of the additional arguments is a list, then it should be enclosed once more into a list() command, since passing a list as an argument flattens the list by one level. SEE ALSO: sprintf, fprintf, print, string EXAMPLE : example printf; shows an example " { write("", sprintf(fmt, #)); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; module m=[1,y],[0,x+z]; intmat M=betti(mres(m,0)); list l = r, m, M; printf("s:%s, l:%l", 1, 2); printf("s:%s", l); printf("s:%s", list(l)); printf("2l:%2l", list(l)); printf("%p", list(l)); printf("%;", list(l)); printf("%b", M); } proc fprintf(link l, string fmt, list #) "USAGE: fprintf(l, fmt, ...) l link; fmt string RETURN: none PURPOSE: fprintf performs output formatting. The second argument is a format control string. Additional arguments may be required, depending on the contents of the control string. A series of output characters is generated as directed by the control string; these characters are written to the link l. The control string is simply text to be copied, except that the string may contain conversion specifications. Do 'help print:' for a listing of valid conversion specifications. As an addition to the conversions of 'print', the '%n' and '%2' conversion specification does not consume an additional argument, but simply generates a newline character. NOTE: If one of the additional arguments is a list, then it should be enclosed once more into a list() command, since passing a list as an argument flattens the list by one level. SEE ALSO: sprintf, printf, print, string EXAMPLE : example fprintf; shows an example " { write(l, sprintf(fmt, #)); } example { "EXAMPLE:"; echo=2; ring r=0,(x,y,z),dp; module m=[1,y],[0,x+z]; intmat M=betti(mres(m,0)); list l = r, m, M; link li = ""; // link to stdout fprintf(li, "s:%s, l:%l", 1, 2); fprintf(li, "s:%s", l); fprintf(li, "s:%s", list(l)); fprintf(li, "2l:%2l", list(l)); fprintf(li, "%p", list(l)); fprintf(li, "%;", list(l)); fprintf(li, "%b", M); } /* proc minres(list #) { if (size(#) == 2) { if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module") { if (typeof(#[2] == "int")) { return (res(#[1],#[2],1)); } } } if (typeof(#[1]) == "resolution") { return minimizeres(#[1]); } else { return minimizeres(#); } } */