=================================================================== RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v retrieving revision 1.5 retrieving revision 1.11 diff -u -p -r1.5 -r1.11 --- OpenXM/src/k097/lib/minimal/minimal.k 2000/05/05 08:13:49 1.5 +++ OpenXM/src/k097/lib/minimal/minimal.k 2000/05/19 11:16:51 1.11 @@ -1,10 +1,13 @@ -/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.4 2000/05/04 11:05:20 takayama Exp $ */ +/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.10 2000/05/07 02:10:44 takayama Exp $ */ #define DEBUG 1 /* #define ORDINARY 1 */ /* 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 */ +/* #define OFFSET 20*/ /* Test sequences. Use load["minimal.k"];; @@ -34,6 +37,7 @@ def load_tower() { sm1(" [(parse) (k0-tower.sm1) pushfile ] extension "); sm1(" /k0-tower.sm1.loaded 1 def "); } + sm1(" oxNoX "); } load_tower(); SonAutoReduce = true; @@ -336,11 +340,21 @@ def test_SinitOfArray() { /* f is assumed to be a monomial with toes. */ def Sdegree(f,tower,level) { - local i; + local i,ww, wd; + /* extern WeightOfSweyl; */ + ww = WeightOfSweyl; f = Init(f); if (level <= 1) return(StotalDegree(f)); i = Degree(f,es); - return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1)); +#ifdef TOTAL_STRATEGY + return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1)); +#endif + /* Strategy must be compatible with ordering. */ + /* Weight vector must be non-negative, too. */ + /* See Sdegree, SgenerateTable, reductionTable. */ + wd = Sord_w(f,ww); + return(wd+Sdegree(tower[level-2,i],tower,level-1)); + } def SgenerateTable(tower) { @@ -351,7 +365,8 @@ def SgenerateTable(tower) { n = Length(tower[i]); ans_at_each_floor=NewArray(n); for (j=0; j<n; j++) { - ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1); + ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1) + + OFFSET; /* Println([i,j,ans_at_each_floor[j]]); */ } ans[i] = ans_at_each_floor; @@ -427,7 +442,7 @@ def SlaScala(g) { reductionTable_tmp; /* extern WeightOfSweyl; */ ww = WeightOfSweyl; - Print("WeghtOfSweyl="); Println(WeightOfSweyl); + Print("WeightOfSweyl="); Println(WeightOfSweyl); rf = SresolutionFrameWithTower(g); redundant_seq = 1; redundant_seq_ordinary = 1; tower = rf[1]; @@ -1024,22 +1039,71 @@ def Sannfs(f,v) { def Sannfs2(f) { local p,pp; p = Sannfs(f,"x,y"); + sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set "); /* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1], ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */ - Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); + /* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1]]); */ + + Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); + pp = Map(p,"Spoly"); + return(Sminimal_v(pp)); + /* return(Sminimal(pp)); */ +} + +HelpAdd(["Sannfs2", +["Sannfs2(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.", + "See also Sminimal_v, Sannfs3.", + "Example: a=Sannfs2(\"x^3-y^2\");", + " b=a[0]; sm1_pmat(b);", + " b[1]*b[0]:", + "Example: a=Sannfs2(\"x*y*(x-y)*(x+y)\");", + " b=a[0]; sm1_pmat(b);", + " b[1]*b[0]:" +]]); + +/* 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[0],"Spoly"); - return(Sminimal(pp)); + pp = Map(p,"Spoly"); + return(Sminimal_v(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.", + "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. @@ -1048,15 +1112,24 @@ def Sannfs3(f) { */ +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 is under construction. */ +/* 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; + reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase; /* extern WeightOfSweyl; */ ww = WeightOfSweyl; Print("WeghtOfSweyl="); Println(WeightOfSweyl); @@ -1121,19 +1194,49 @@ def Sschreyer(g) { /* i must be equal to f[2], I think. Double check. */ /* Correction Of Constant */ - c2 = f[6]; + /* Correction of syzygy */ + c2 = f[6]; /* or -f[6]? Double check. */ + Print("c2="); Println(c2); nn = Length(bases); for (ii=0; ii<nn;ii++) { - if (ii != place) { - bases[ii] = bases[ii]*c2; + if ((ii != i) && (! IsNull(bases[ii]))) { + m = Length(bases[ii]); + for (jj=0; jj<m; jj++) { + if (jj != place) { + bases[ii,jj] = bases[ii,jj]*c2; + } + } } } + Print("Old freeRes[level] = "); sm1_pmat(freeRes[level]); freeRes[level] = bases; - /* bases = freeRes[level-1]; - bases[place] = f[0]; - freeRes[level-1] = bases; It is already set. */ - reducer[level-1,place] = f[1]; + Print("New freeRes[level] = "); sm1_pmat(freeRes[level]); + + /* Update the freeRes[level-1] */ + Print("Old freeRes[level-1] = "); sm1_pmat(freeRes[level-1]); + bases = freeRes[level-1]; + bases[place] = f[0]; + freeRes[level-1] = bases; + Print("New freeRes[level-1] = "); sm1_pmat(freeRes[level-1]); + + reducer[level-1,place] = f[1]-SunitOfFormat(place,f[1]); + /* This reducer is different from that of SlaScala(). */ + + reducerBasis = reducer[level-1]; + nn = Length(reducerBasis); + for (ii=0; ii<nn;ii++) { + if ((ii != place) && (! IsNull(reducerBasis[ii]))) { + m = Length(reducerBasis[ii]); + for (jj=0; jj<m; jj++) { + if (jj != place) { + reducerBasis[ii,jj] = reducerBasis[ii,jj]*c2; + } + } + } + } + reducer[level-1] = reducerBasis; + }else{ /* redundantTable[level,i] = 0; */ bases = freeRes[level]; @@ -1143,6 +1246,33 @@ def Sschreyer(g) { } /* end of level >= 1 */ } } /* i loop */ + + /* Triangulate reducer */ + if (level >= 1) { + Println(" "); + Print("Triangulating reducer at level "); Println(level-1); + Println("freeRes[level]="); sm1_pmat(freeRes[level]); + reducerBase = reducer[level-1]; + Print("reducerBase="); Println(reducerBase); + Println("Compare freeRes[level] and reducerBase (put -1)"); + m = Length(reducerBase); + for (ii=m-1; 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); @@ -1151,6 +1281,18 @@ def Sschreyer(g) { bases = Sbases_to_vec(bases,bettiTable[i]); freeResV[i] = bases; } + + /* Mark the non-redundant elements. */ + for (i=0; i<n; i++) { + m = Length(redundantTable[i]); + for (jj=0; jj<m; jj++) { + if (IsNull(redundantTable[i,jj])) { + redundantTable[i,jj] = 0; + } + } + } + + return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]); } @@ -1158,7 +1300,7 @@ def SpairAndReduction2(skel,level,ii,freeRes,tower,ww, local i, j, myindex, p, bases, tower2, gi, gj, si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2, vdeg,vdeg_reduced,n,c2; - Println("SpairAndReduction2:"); + Println("SpairAndReduction2 : -------------------------"); if (level < 1) Error("level should be >= 1 in SpairAndReduction."); p = skel[level,ii]; @@ -1193,6 +1335,11 @@ def SpairAndReduction2(skel,level,ii,freeRes,tower,ww, 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]; @@ -1203,15 +1350,18 @@ def SpairAndReduction2(skel,level,ii,freeRes,tower,ww, /* tmp[0] must be zero */ n = Length(t_syz); for (i=0; i<n; i++) { - if (IsConstant(t_syz[i])) { + if (IsConstant(t_syz[i])){ + if (!IsZero(t_syz[i])) { if (IsNull(redundantTable[level-1,i])) { /* i must equal to pos2 below. */ c2 = -t_syz[i]; - tmp[0] = freeRes[level-1,i]; + tmp[0] = c2*Stoes_vec(freeRes[level-1,i]); t_syz[i] = 0; + /* tmp[0] = t_syz . g */ /* break; does not work. Use */ i = n; } + } } } @@ -1222,11 +1372,317 @@ def SpairAndReduction2(skel,level,ii,freeRes,tower,ww, Print("vdegree of the original = "); Println(vdeg); Print("vdegree of the remainder = "); Println(vdeg_reduced); + if (!IsNull(vdeg_reduced)) { + if (vdeg_reduced < vdeg) { + Println("--- Special in V-minimal!"); + Println(tmp[0]); + Println("syzygy="); sm1_pmat(t_syz); + Print("[vdeg, vdeg_reduced] = "); Println([vdeg,vdeg_reduced]); + } + } + + pos = SwhereInTower(syzHead,tower[level]); pos2 = SwhereInTower(tmp[0],tower[level-1]); ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2]; /* pos is the place to put syzygy at level. */ /* pos2 is the place to put a new GB at level-1. */ Println(ans); + Println(" "); return(ans); +} + +HelpAdd(["Sminimal_v", +["It constructs the V-minimal free resolution from the Schreyer resolution", + "step by step.", + "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],", + " [0, x^2+y^2],", + " [0, x*y]];", + " a=Sminimal_v(v);", + " sm1_pmat(a[0]); b=a[0]; b[1]*b[0]:", + "Note: a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]); + + +def Sminimal_v(g) { + local r, freeRes, redundantTable, reducer, maxLevel, + minRes, seq, maxSeq, level, betti, q, bases, dr, + betti_levelplus, newbases, i, j,qq,tminRes; + r = Sschreyer(g); + sm1_pmat(r); + Debug_Sminimal_v = r; + Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v"); + /* Should I turn off the tower?? */ + freeRes = r[0]; + redundantTable = r[1]; + reducer = r[2]; + minRes = SnewArrayOfFormat(freeRes); + seq = 0; + maxSeq = SgetMaxSeq(redundantTable); + maxLevel = Length(freeRes); + for (level = 0; level < maxLevel; level++) { + minRes[level] = freeRes[level]; + } + for (level = 0; level < maxLevel; level++) { + betti = Length(freeRes[level]); + for (q = betti-1; q>=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<betti_levelplus; i++) { + newbases[i] = dr[q]*bases[i] - bases[i,q]*dr; + } + Println(["level, q =", level,q]); + Println("bases="); sm1_pmat(bases); + Println("dr="); sm1_pmat(dr); + Println("newbases="); sm1_pmat(newbases); + minRes[level+1] = newbases; + freeRes = minRes; +#ifdef DEBUG + for (qq=q; qq<betti; qq++) { + for (i=0; i<betti_levelplus; i++) { + if ((!IsZero(newbases[i,qq])) && (redundantTable[level,qq] >0)) { + Println(["[i,qq]=",[i,qq]," is not zero in newbases."]); + Print("redundantTable ="); sm1_pmat(redundantTable[level]); + Error("Stop in Sminimal for debugging."); + } + } + } +#endif + } + } + } + } + tminRes = Stetris(minRes,redundantTable); + return([SpruneZeroRow(tminRes), tminRes, + [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]); + /* r[4] is the redundantTable_ordinary */ + /* r[0] is the freeResolution */ +} + +/* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */ +/* x y (x+y-1)(x-2), x^3-y^2, x^3 - y^2 z^2, + x y z (x+y+z-1) seems to be interesting, because the first syzygy + contains 1. +*/ + +def CopyArray(m) { + local ans,i,n; + if (IsArray(m)) { + n = Length(m); + ans = NewArray(n); + for (i=0; i<n; i++) { + ans[i] = CopyArray(m[i]); + } + return(ans); + }else{ + return(m); + } +} +HelpAdd(["CopyArray", +["It duplicates the argument array recursively.", + "Example: m=[1,[2,3]];", + " a=CopyArray(m); a[1] = \"Hello\";", + " Println(m); Println(a);"]]); + +def IsZeroVector(m) { + local n,i; + n = Length(m); + for (i=0; i<n; i++) { + if (!IsZero(m[i])) { + return(false); + } + } + return(true); +} + +def SpruneZeroRow(res) { + local minRes, n,i,j,m, base,base2,newbase,newbase2, newMinRes; + + minRes = CopyArray(res); + n = Length(minRes); + for (i=0; i<n; i++) { + base = minRes[i]; + m = Length(base); + if (i != n-1) { + base2 = minRes[i+1]; + base2 = Transpose(base2); + } + newbase = [ ]; + newbase2 = [ ]; + for (j=0; j<m; j++) { + if (!IsZeroVector(base[j])) { + newbase = Append(newbase,base[j]); + if (i != n-1) { + newbase2 = Append(newbase2,base2[j]); + } + } + } + minRes[i] = newbase; + if (i != n-1) { + if (newbase2 == [ ]) { + minRes[i+1] = [ ]; + }else{ + minRes[i+1] = Transpose(newbase2); + } + } + } + + newMinRes = [ ]; + n = Length(minRes); + i = 0; + while (i < n ) { + base = minRes[i]; + if (base == [ ]) { + i = n; /* break; */ + }else{ + newMinRes = Append(newMinRes,base); + } + i++; + } + return(newMinRes); +} + +def testAnnfs2(f) { + local a,i,n; + a = Sannfs2(f); + b=a[0]; + n = Length(b); + Println("------ V-minimal free resolution -----"); + sm1_pmat(b); + Println("----- Is it complex? ---------------"); + for (i=0; i<n-1; i++) { + Println(b[i+1]*b[i]); + } + return(a); +} +def testAnnfs3(f) { + local a,i,n; + a = Sannfs3(f); + b=a[0]; + n = Length(b); + Println("------ V-minimal free resolution -----"); + sm1_pmat(b); + Println("----- Is it complex? ---------------"); + for (i=0; i<n-1; i++) { + Println(b[i+1]*b[i]); + } + return(a); +} + +def ToString_array(p) { + local ans; + if (IsArray(p)) { + ans = Map(p,"ToString_array"); + }else{ + ans = ToString(p); + } + return(ans); +} + +/* sm1_res_div([[x],[y]],[[x^2],[x*y],[y^2]],[x,y]): */ + +def sm1_res_div(I,J,V) { + I = ToString_array(I); + J = ToString_array(J); + V = ToString_array(V); + sm1(" [[ I J] V ] res*div /FunctionValue set "); +} + +/* It has not yet been working */ +def sm1_res_kernel_image(m,n,v) { + m = ToString_array(m); + n = ToString_array(n); + v = ToString_array(v); + sm1(" [m n v] res-kernel-image /FunctionValue set "); +} +def Skernel(m,v) { + m = ToString_array(m); + v = ToString_array(v); + sm1(" [ m v ] syz /FunctionValue set "); +} + +def test3() { + local a1,a2,b1,b2; + a1 = Sannfs3("x^3-y^2*z^2"); + a1 = a1[0]; + a2 = Sannfs3_laScala2("x^3-y^2*z^2"); + a2 = a2[0]; + b1 = a1[1]; + b2 = a2[1]; + sm1_pmat(b2); + Println(" OVER "); + sm1_pmat(b1); + return([sm1_res_div(b2,b1,["x","y","z"]),b2,b1,a2,a1]); +} + +def test4() { + local a,b; + a = Sannfs3_laScala2("x^3-y^2*z^2"); + b = a[0]; + sm1_pmat( sm1_res_kernel_image(b[0],b[1],[x,y,z])); + sm1_pmat( sm1_res_kernel_image(b[1],b[2],[x,y,z])); + return(a); +} + +def sm1_gb(f,v) { + f =ToString_array(f); + v = ToString_array(v); + sm1(" [f v] gb /FunctionValue set "); +} + +def test5() { + local a,b,c,cc,v; + a = Sannfs3_laScala2("x^3-y^2*z^2"); + b = a[0]; + v = [x,y,z]; + c = Skernel(b[0],v); + c = c[0]; + sm1_pmat([c,b[1],v]); + Println("-----------------------------------"); + cc = sm1_res_div(c,b[1],v); + sm1_pmat(sm1_gb(cc,v)); + c = Skernel(b[1],v); + c = c[0]; + cc = sm1_res_div(c,b[2],v); + sm1_pmat(sm1_gb(cc,v)); + return(a); +} +def test6() { + local a,b,c,cc,v; + a = Sannfs3("x^3-y^2*z^2"); + b = a[0]; + v = [x,y,z]; + c = Skernel(b[0],v); + c = c[0]; + sm1_pmat([c,b[1],v]); + Println("-------ker = im for minimal ?---------------------"); + cc = sm1_res_div(c,b[1],v); + sm1_pmat(sm1_gb(cc,v)); + c = Skernel(b[1],v); + c = c[0]; + cc = sm1_res_div(c,b[2],v); + sm1_pmat(sm1_gb(cc,v)); + Println("------ ker=im for Schreyer ?------------------"); + b = a[3]; + c = Skernel(b[0],v); + c = c[0]; + sm1_pmat([c,b[1],v]); + cc = sm1_res_div(c,b[1],v); + sm1_pmat(sm1_gb(cc,v)); + c = Skernel(b[1],v); + c = c[0]; + cc = sm1_res_div(c,b[2],v); + sm1_pmat(sm1_gb(cc,v)); + return(a); }