=================================================================== RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v retrieving revision 1.14 retrieving revision 1.30 diff -u -p -r1.14 -r1.30 --- OpenXM/src/k097/lib/minimal/minimal.k 2000/06/09 08:04:54 1.14 +++ OpenXM/src/k097/lib/minimal/minimal.k 2000/11/19 05:50:30 1.30 @@ -1,13 +1,28 @@ -/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.13 2000/06/08 08:37:53 takayama Exp $ */ +/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.29 2000/08/22 05:34:06 takayama Exp $ */ #define DEBUG 1 -/* #define ORDINARY 1 */ +Sordinary = false; /* If you run this program on openxm version 1.1.2 (FreeBSD), make a symbolic link by the command ln -s /usr/bin/cpp /lib/cpp */ #define OFFSET 0 -#define TOTAL_STRATEGY 1 /* #define OFFSET 20*/ +Sverbose = false; /* Be extreamly verbose */ +Sverbose2 = true; /* Don't be quiet and show minimal information */ +def Sprintln(s) { + if (Sverbose) Println(s); +} +def Sprint(s) { + if (Sverbose) Print(s); +} +def Sprintln2(s) { + if (Sverbose2) Println(s); +} +def Sprint2(s) { + if (Sverbose2) Print(s); + sm1(" [(flush)] extension "); +} + /* Test sequences. Use load["minimal.k"];; @@ -35,6 +50,7 @@ def load_tower() { if (Boundp("k0-tower.sm1.loaded")) { }else{ sm1(" [(parse) (k0-tower.sm1) pushfile ] extension "); + sm1(" [(parse) (new.sm1) pushfile ] extension "); sm1(" /k0-tower.sm1.loaded 1 def "); } sm1(" oxNoX "); @@ -50,6 +66,108 @@ def Reverse(f) { def Sgroebner(f) { sm1(" [f] groebner /FunctionValue set"); } + +def Sinvolutive(f,w) { + local g,m; + if (IsArray(f[0])) { + m = NewArray(Length(f[0])); + }else{ + m = [0]; + } + g = Sgroebner(f); + /* This is a temporary code. */ + sm1(" g 0 get { w m init_w} map /FunctionValue set "); +} + + + +def Error(s) { + sm1(" s error "); +} + +def IsNull(s) { + if (Stag(s) == 0) return(true); + else return(false); +} + +def MonomialPart(f) { + sm1(" [(lmonom) f] gbext /FunctionValue set "); +} + +def Warning(s) { + Print("Warning: "); + Println(s); +} +def RingOf(f) { + local r; + if (IsPolynomial(f)) { + if (f != Poly("0")) { + sm1(f," (ring) dc /r set "); + }else{ + sm1(" [(CurrentRingp)] system_variable /r set "); + } + }else{ + Warning("RingOf(f): the argument f must be a polynomial. Return the current ring."); + sm1(" [(CurrentRingp)] system_variable /r set "); + } + return(r); +} + +def Ord_w_m(f,w,m) { + sm1(" f w m ord_w { (universalNumber) dc } map /FunctionValue set "); +} +HelpAdd(["Ord_w_m", +["Ord_w_m(f,w,m) returns the order of f with respect to w with the shift m.", + "Note that the order of the ring and the weight w must be the same.", + "When f is zero, it returns -intInfinity = -999999999.", + "Example: Sweyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ", + " Ord_w_m([x*Dx+1,Dx^2+x^5],[\"x\",-1,\"Dx\",1],[2,0]):"]]); + +def Init_w_m(f,w,m) { + sm1(" f w m init_w /FunctionValue set "); +} +HelpAdd(["Init_w_m", +["Init_w_m(f,w,m) returns the initial of f with respect to w with the shift m.", + "Note that the order of the ring and the weight w must be the same.", + "Example: Sweyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ", + " Init_w_m([x*Dx+1,Dx^2+x^5],[\"x\",-1,\"Dx\",1],[2,0]):"]]); + +def Max(v) { + local i,t,n; + n = Length(v); + if (n == 0) return(null); + t = v[0]; + for (i=0; i t) { t = v[i];} + } + return(t); +} +HelpAdd(["Max", +["Max(v) returns the maximal element in v."]]); + +def Kernel(f) { + sm1(" [f] syz /FunctionValue set "); +} +def Syz(f) { + sm1(" [f] syz /FunctionValue set "); +} +HelpAdd(["Kernel", +["Kernel(f) returns the syzygy of f.", + "Return value [b, c]: b is a set of generators of the syzygies of f", + " : c=[gb, backward transformation, syzygy without", + " dehomogenization", + "Example: Weyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ", + " s=Kernel([x*Dx+1,Dx^2+x^5]); s[0]:"]]); +/* cf. sm1_syz in cohom.k */ +def Gb(f) { + sm1(" [f] gb /FunctionValue set "); +} +HelpAdd(["Gb", +["Gb(f) returns the Groebner basis of f.", + "cf. Kernel, Weyl."]]); + + +/* End of standard functions that should be moved to standard libraries. */ def test0() { local f; Sweyl("x,y,z"); @@ -69,7 +187,6 @@ def test1() { } - def Sweyl(v,w) { /* extern WeightOfSweyl ; */ local ww,i,n; @@ -137,11 +254,34 @@ sm1(" [(AvoidTheSameRing)] pushEnv def SresolutionFrameWithTower(g,opt) { local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof, - gbasis; + gbasis, nohomog,i,n; + /* extern Sordinary */ + nohomog = false; + count = -1; Sordinary = false; /* default value for options. */ if (Length(Arglist) >= 2) { - if (IsInteger(opt)) count = opt; - }else{ - count = -1; + if (IsArray(opt)) { + n = Length(opt); + for (i=0; i= 1 in SpairAndReduction."); p = skel[level,ii]; myindex = p[0]; i = myindex[0]; j = myindex[1]; bases = freeRes[level-1]; - Println(["p and bases ",p,bases]); + Sprintln(["p and bases ",p,bases]); if (IsNull(bases[i]) || IsNull(bases[j])) { - Println([level,i,j,bases[i],bases[j]]); + Sprintln([level,i,j,bases[i],bases[j]]); Error("level, i, j : bases[i], bases[j] must not be NULL."); } tower2 = StowerOf(tower,level-1); SsetTower(tower2); - Println(["level=",level]); - Println(["tower2=",tower2]); + Sprintln(["level=",level]); + Sprintln(["tower2=",tower2]); /** sm1(" show_ring "); */ gi = Stoes_vec(bases[i]); @@ -720,23 +872,23 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) sj = ssp[0,1]; syzHead = si*es^i; /* This will be the head term, I think. But, double check. */ - Println([si*es^i,sj*es^j]); + Sprintln([si*es^i,sj*es^j]); - Print("[gi, gj] = "); Println([gi,gj]); - sm1(" [(Homogenize)] system_variable message "); - Print("Reduce the element "); Println(si*gi+sj*gj); - Print("by "); Println(bases); + Sprint("[gi, gj] = "); Sprintln([gi,gj]); + sm1(" [(Homogenize)] system_variable "); + Sprint("Reduce the element "); Sprintln(si*gi+sj*gj); + Sprint("by "); Sprintln(bases); tmp = Sreduction(si*gi+sj*gj, bases); - Print("result is "); Println(tmp); + Sprint("result is "); Sprintln(tmp); /* This is essential part for V-minimal resolution. */ /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */ vdeg = SvDegree(si*gi,tower,level-1,ww); vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww); - Print("vdegree of the original = "); Println(vdeg); - Print("vdegree of the remainder = "); Println(vdeg_reduced); + Sprint("vdegree of the original = "); Sprintln(vdeg); + Sprint("vdegree of the remainder = "); Sprintln(vdeg_reduced); t_syz = tmp[2]; si = si*tmp[1]+t_syz[i]; @@ -752,7 +904,7 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced]; /* pos is the place to put syzygy at level. */ /* pos2 is the place to put a new GB at level-1. */ - Println(ans); + Sprintln(ans); return(ans); } @@ -785,24 +937,6 @@ def Sreduction(f,myset) { return([tmp[0],tmp[1],t_syz]); } -def Warning(s) { - Print("Warning: "); - Println(s); -} -def RingOf(f) { - local r; - if (IsPolynomial(f)) { - if (f != Poly("0")) { - sm1(f," (ring) dc /r set "); - }else{ - sm1(" [(CurrentRingp)] system_variable /r set "); - } - }else{ - Warning("RingOf(f): the argument f must be a polynomial. Return the current ring."); - sm1(" [(CurrentRingp)] system_variable /r set "); - } - return(r); -} def Sfrom_es(f,size) { local c,ans, i, d, myes, myee, j,n,r,ans2; @@ -861,7 +995,10 @@ def Sbases_to_vec(bases,size) { } HelpAdd(["Sminimal", -["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm", +["It constructs the V-minimal free resolution by LaScala's algorithm", + "option: \"homogenized\" (no automatic homogenization)", + " : \"Sordinary\" (no (u,v)-minimal resolution)", + "Options should be given as an array.", "Example: Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);", " v=[[2*x*Dx + 3*y*Dy+6, 0],", " [3*x^2*Dy + 2*y*Dx, 0],", @@ -873,15 +1010,25 @@ HelpAdd(["Sminimal", " IsExact_h(b,[x,y]):", "Note: a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]); -def Sminimal(g) { +def Sminimal(g,opt) { local r, freeRes, redundantTable, reducer, maxLevel, minRes, seq, maxSeq, level, betti, q, bases, dr, - betti_levelplus, newbases, i, j,qq, tminRes; - r = SlaScala(g); + betti_levelplus, newbases, i, j,qq, tminRes,bettiTable, ansSminimal; + if (Length(Arglist) < 2) { + opt = null; + } + /* Sordinary is set in SlaScala(g,opt) --> SresolutionFrameWithTower */ + + ScheckIfSchreyer("Sminimal:0"); + r = SlaScala(g,opt); /* Should I turn off the tower?? */ + ScheckIfSchreyer("Sminimal:1"); freeRes = r[0]; redundantTable = r[1]; reducer = r[2]; + bettiTable = SbettiTable(redundantTable); + Sprintln2("BettiTable ------"); + if (Sverbose || Sverbose2) {sm1_pmat(bettiTable);} minRes = SnewArrayOfFormat(freeRes); seq = 0; maxSeq = SgetMaxSeq(redundantTable); @@ -891,12 +1038,12 @@ def Sminimal(g) { } seq=maxSeq+1; while (seq > 1) { - seq--; + seq--; Sprint2(seq); for (level = 0; level < maxLevel; level++) { betti = Length(freeRes[level]); for (q = 0; q The betti numbers are 3, 2. + a=Sannfs2("x^3-y^2-x"); + a=Sannfs2("x*y*(x-y)"); +*/ -/* Do not forget to turn on TOTAL_STRATEGY */ -def Sannfs2_laScala(f) { - local p,pp; - p = Sannfs(f,"x,y"); - /* Do not make laplace transform. - sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set "); - p = [p]; - */ - Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); - pp = Map(p[0],"Spoly"); - return(Sminimal(pp)); -} -def Sannfs2_laScala2(f) { - local p,pp; - p = Sannfs(f,"x,y"); - sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set "); - p = [p]; - Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1], - ["x",-1,"y",-1,"Dx",1,"Dy",1]]); - pp = Map(p[0],"Spoly"); - return(Sminimal(pp)); -} - def Sannfs3(f) { local p,pp; p = Sannfs(f,"x,y,z"); sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set "); Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]); pp = Map(p,"Spoly"); - return(Sminimal_v(pp)); + return(Sminimal(pp)); } HelpAdd(["Sannfs3", ["Sannfs3(f) constructs the V-minimal free resolution for the weight (-1,1)", "of the Laplace transform of the annihilating ideal of the polynomial f in x,y,z.", - "See also Sminimal_v, Sannfs2.", + "See also Sminimal, Sannfs2.", "Example: a=Sannfs3(\"x^3-y^2*z^2\");", " b=a[0]; sm1_pmat(b);", " b[1]*b[0]: b[2]*b[1]:"]]); -/* - The betti numbers of most examples are 2,1. (0-th and 1-th). - a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2. - a=Sannfs2("x^3-y^2-x"); : it causes an error. It should be fixed. - a=Sannfs2("x*y*(x-y)"); : it causes an error. It should be fixed. - -*/ -def Sannfs3_laScala2(f) { - local p,pp; - p = Sannfs(f,"x,y,z"); - sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set "); - Sweyl("x,y,z",[["x",1,"y",1,"z",1,"Dx",1,"Dy",1,"Dz",1,"h",1], - ["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]); - pp = Map(p,"Spoly"); - return(Sminimal(pp)); -} - -/* The below does not use LaScala-Stillman's algorithm. */ -def Sschreyer(g) { - local rf, tower, reductionTable, skel, redundantTable, bases, - strategy, maxOfStrategy, height, level, n, i, - freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww, - redundantTable_ordinary, redundant_seq_ordinary, - reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase; - /* extern WeightOfSweyl; */ - ww = WeightOfSweyl; - Print("WeghtOfSweyl="); Println(WeightOfSweyl); - rf = SresolutionFrameWithTower(g); - redundant_seq = 1; redundant_seq_ordinary = 1; - tower = rf[1]; - reductionTable = SgenerateTable(tower); - skel = rf[2]; - redundantTable = SnewArrayOfFormat(rf[1]); - redundantTable_ordinary = SnewArrayOfFormat(rf[1]); - reducer = SnewArrayOfFormat(rf[1]); - freeRes = SnewArrayOfFormat(rf[1]); - bettiTable = SsetBettiTable(rf[1],g); - - height = Length(reductionTable); - for (level = 0; level < height; level++) { - n = Length(reductionTable[level]); - for (i=0; i f[5]) { - /* Zero in the gr-module */ - Print("v-degree of [org,remainder] = "); - Println([f[4],f[5]]); - Print("[level,i] = "); Println([level,i]); - redundantTable[level-1,place] = 0; - }else{ - redundantTable[level-1,place] = redundant_seq; - redundant_seq++; - } -#endif - redundantTable_ordinary[level-1,place] - =redundant_seq_ordinary; - redundant_seq_ordinary++; - bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */ - /* redundantTable[level,i] = 0; - redundantTable_ordinary[level,i] = 0; */ - /* i must be equal to f[2], I think. Double check. */ - - /* Correction Of Constant */ - /* Correction of syzygy */ - c2 = f[6]; /* or -f[6]? Double check. */ - Print("c2="); Println(c2); - nn = Length(bases); - for (ii=0; ii=0; ii--) { - if (!IsNull(reducerBase[ii])) { - for (jj=ii-1; jj>=0; jj--) { - if (!IsNull(reducerBase[jj])) { - if (!IsZero(reducerBase[jj,ii])) { - /* reducerBase[ii,ii] should be always constant. */ - reducerBase[jj] = reducerBase[ii,ii]*reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii]; - } - } - } - } - } - Println("New reducer"); - sm1_pmat(reducerBase); - reducer[level-1] = reducerBase; - } - - } /* level loop */ - n = Length(freeRes); - freeResV = SnewArrayOfFormat(freeRes); - for (i=0; i= 1 in SpairAndReduction."); - p = skel[level,ii]; - myindex = p[0]; - i = myindex[0]; j = myindex[1]; - bases = freeRes[level-1]; - Println(["p and bases ",p,bases]); - if (IsNull(bases[i]) || IsNull(bases[j])) { - Println([level,i,j,bases[i],bases[j]]); - Error("level, i, j : bases[i], bases[j] must not be NULL."); - } - - tower2 = StowerOf(tower,level-1); - SsetTower(tower2); - Println(["level=",level]); - Println(["tower2=",tower2]); - /** sm1(" show_ring "); */ - - gi = Stoes_vec(bases[i]); - gj = Stoes_vec(bases[j]); - - ssp = Sspolynomial(gi,gj); - si = ssp[0,0]; - sj = ssp[0,1]; - syzHead = si*es^i; - /* This will be the head term, I think. But, double check. */ - Println([si*es^i,sj*es^j]); - - Print("[gi, gj] = "); Println([gi,gj]); - sm1(" [(Homogenize)] system_variable message "); - Print("Reduce the element "); Println(si*gi+sj*gj); - Print("by "); Println(bases); - - tmp = Sreduction(si*gi+sj*gj, bases); - - Print("result is "); Println(tmp); - if (!IsZero(tmp[0])) { - Print("Error: base = "); - Println(Map(bases,"Stoes_vec")); - Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring."); - } - t_syz = tmp[2]; - si = si*tmp[1]+t_syz[i]; - sj = sj*tmp[1]+t_syz[j]; - t_syz[i] = si; - t_syz[j] = sj; - - c2 = null; - /* tmp[0] must be zero */ - n = Length(t_syz); - for (i=0; i=0; q--) { - if (redundantTable[level,q] > 0) { - Print("[seq,level,q]="); Println([seq,level,q]); - if (level < maxLevel-1) { - bases = freeRes[level+1]; - dr = reducer[level,q]; - /* dr[q] = -1; We do not need this in our reducer format. */ - /* dr[q] should be a non-zero constant. */ - newbases = SnewArrayOfFormat(bases); - betti_levelplus = Length(bases); - /* - bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j] - */ - for (i=0; i (homogenized Weyl algebra)", + "cf. ReParse" +]]); + def ReParse(a) { local c; if (IsArray(a)) { @@ -1730,3 +1492,125 @@ HelpAdd(["ReParse", "In this case, ReParse the outputs from these functions.", "cf. IsExaxt_h" ]]); + +def ScheckIfSchreyer(s) { + local ss; + sm1(" (report) (grade) switch_function /ss set "); + if (ss != "module1v") { + Print("ScheckIfSchreyer: from "); Println(s); + Error("grade is not module1v"); + } + /* + sm1(" (report) (mmLarger) switch_function /ss set "); + if (ss != "tower") { + Print("ScheckIfSchreyer: from "); Println(s); + Error("mmLarger is not tower"); + } + */ + sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set "); + if (ss != 1) { + Print("ScheckIfSchreyer: from "); Printl(s); + Error("Schreyer order is not set."); + } + /* More check will be necessary. */ + return(true); +} + +def SgetShift(mat,w,m) { + local omat; + sm1(" mat { w m ord_w {(universalNumber) dc}map } map /omat set"); + return(Map(omat,"Max")); +} +HelpAdd(["SgetShift", +["SgetShift(mat,w,m) returns the shift vector of mat with respect to w with the shift m.", + "Note that the order of the ring and the weight w must be the same.", + "Example: Sweyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ", + " SgetShift([[x*Dx+1,Dx^2+x^5],[Poly(\"0\"),x],[x,x]],[\"x\",-1,\"Dx\",1],[2,0]):"]]); + +def SgetShifts(resmat,w) { + local i,n,ans,m0; + n = Length(resmat); + ans = NewArray(n+1); + m0 = NewArray(Length(resmat[0,0])); + ans[0] = m0; + for (i=0; i