===================================================================
RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v
retrieving revision 1.6
retrieving revision 1.14
diff -u -p -r1.6 -r1.14
--- OpenXM/src/k097/lib/minimal/minimal.k	2000/05/06 07:58:37	1.6
+++ OpenXM/src/k097/lib/minimal/minimal.k	2000/06/09 08:04:54	1.14
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.5 2000/05/05 08:13:49 takayama Exp $ */
+/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.13 2000/06/08 08:37:53 takayama Exp $ */
 #define DEBUG 1 
 /* #define ORDINARY 1 */
 /* If you run this program on openxm version 1.1.2 (FreeBSD),
@@ -6,7 +6,7 @@
    ln -s /usr/bin/cpp /lib/cpp
 */
 #define OFFSET 0 
-#define TOTAL_STRATEGY
+#define TOTAL_STRATEGY 1
 /* #define OFFSET 20*/
 /* Test sequences. 
    Use load["minimal.k"];;
@@ -37,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;
@@ -131,6 +132,7 @@ sm1(" [(AvoidTheSameRing)] pushEnv 
       [ [(AvoidTheSameRing) 0] system_variable 
         [(gbListTower) tower (list) dc] system_variable
       ] pop popEnv ");
+      /* sm1("(hoge) message show_ring "); */
 }
 
 def SresolutionFrameWithTower(g,opt) {
@@ -290,7 +292,10 @@ def Sres0FrameWithSkelton(g) {
 
 
 def StotalDegree(f) {
-  sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");
+  local d0;
+  sm1(" [(grade) f] gbext (universalNumber) dc /d0 set ");
+  /* Print("degree of "); Print(f); Print(" is "); Println(d0); */
+  return(d0);
 }
 
 /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
@@ -443,6 +448,7 @@ def SlaScala(g) {
   ww = WeightOfSweyl;
   Print("WeightOfSweyl="); Println(WeightOfSweyl);
   rf = SresolutionFrameWithTower(g);
+  Print("rf="); sm1_pmat(rf);
   redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];
   reductionTable = SgenerateTable(tower);
@@ -660,6 +666,12 @@ def MonomialPart(f) {
   sm1(" [(lmonom) f] gbext /FunctionValue set ");
 }
 
+/* WARNING:
+  When you use SwhereInTower, you have to change gbList
+  as below. Ofcourse, you should restrore the gbList
+  SsetTower(StowerOf(tower,level));
+  pos = SwhereInTower(syzHead,tower[level]);
+*/
 def SwhereInTower(f,tower) {
   local i,n,p,q;
   if (f == Poly("0")) return(-1);
@@ -696,6 +708,8 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
 
   tower2 = StowerOf(tower,level-1);
   SsetTower(tower2);  
+  Println(["level=",level]);
+  Println(["tower2=",tower2]);
   /** sm1(" show_ring ");   */
 
   gi = Stoes_vec(bases[i]);
@@ -729,7 +743,11 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   sj = sj*tmp[1]+t_syz[j];
   t_syz[i] = si; 
   t_syz[j] = sj;
+
+  SsetTower(StowerOf(tower,level));
   pos = SwhereInTower(syzHead,tower[level]);
+
+  SsetTower(StowerOf(tower,level-1));
   pos2 = SwhereInTower(tmp[0],tower[level-1]);
   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
   /* pos is the place to put syzygy at level. */
@@ -842,10 +860,23 @@ def Sbases_to_vec(bases,size) {
   return(newbases);
 }
 
+HelpAdd(["Sminimal",
+["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm",
+ "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);",
+ "         Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);", 
+ "         b = ReParse(a[0]); sm1_pmat(b); ",
+ "         IsExact_h(b,[x,y]):",
+ "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
+
 def Sminimal(g) {
   local r, freeRes, redundantTable, reducer, maxLevel,
         minRes, seq, maxSeq, level, betti, q, bases, dr,
-        betti_levelplus, newbases, i, j,qq;
+        betti_levelplus, newbases, i, j,qq, tminRes;
   r = SlaScala(g);
   /* Should I turn off the tower?? */
   freeRes = r[0];
@@ -903,7 +934,8 @@ def Sminimal(g) {
       }
     }
    }
-   return([Stetris(minRes,redundantTable),
+   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 */
@@ -1043,12 +1075,25 @@ def Sannfs2(f) {
   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,"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;
@@ -1062,6 +1107,17 @@ def Sannfs2_laScala(f) {
   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");
@@ -1071,6 +1127,14 @@ def Sannfs3(f) {
   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. 
@@ -1079,6 +1143,15 @@ 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 does not use LaScala-Stillman's algorithm. */
@@ -1152,22 +1225,49 @@ def Sschreyer(g) {
                   /* i must be equal to f[2], I think. Double check. */
 
                   /* Correction Of Constant */
-                  c2 = -f[6];  /* or f[6]?  Double check. */
+                  /* 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) && (! IsNull(bases[ii]))) {
-                       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;
+                  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];
+                  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];
@@ -1182,15 +1282,18 @@ def Sschreyer(g) {
     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[jj] = reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
+                /* reducerBase[ii,ii] should be always constant. */  
+                reducerBase[jj] = reducerBase[ii,ii]*reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
               }
              }
            }
@@ -1243,6 +1346,8 @@ def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
 
   tower2 = StowerOf(tower,level-1);
   SsetTower(tower2);  
+  Println(["level=",level]);
+  Println(["tower2=",tower2]);
   /** sm1(" show_ring ");   */
 
   gi = Stoes_vec(bases[i]);
@@ -1300,20 +1405,48 @@ 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]);
+    }
+  }
+
+  SsetTower(StowerOf(tower,level));
   pos = SwhereInTower(syzHead,tower[level]);
+
+  SsetTower(StowerOf(tower,level-1));
   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("  ");
+  Println("--- end of SpairAndReduction2  ");
   return(ans);
 }
 
+HelpAdd(["Sminimal_v",
+["It constructs the V-minimal free resolution from the Schreyer resolution",
+ "step by step.",
+ "This code still contains bugs. It sometimes outputs wrong answer.",
+ "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."]]);
+
+/* This code still contains bugs. It sometimes outputs wrong answer. */
+/* See test12() in minimal-test.k.  */
+/* There may be remaining 1, too */
 def Sminimal_v(g) {
   local r, freeRes, redundantTable, reducer, maxLevel,
         minRes, seq, maxSeq, level, betti, q, bases, dr,
-        betti_levelplus, newbases, i, j,qq;
+        betti_levelplus, newbases, i, j,qq,tminRes;
   r = Sschreyer(g);
   sm1_pmat(r);
   Debug_Sminimal_v = r;
@@ -1337,14 +1470,15 @@ def Sminimal_v(g) {
           if (level < maxLevel-1) {
             bases = freeRes[level+1];
             dr = reducer[level,q];
-            dr[q] = -1;
+            /* 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] = bases[i] + bases[i,q]*dr;
+              newbases[i] = dr[q]*bases[i] - bases[i,q]*dr;
             }
             Println(["level, q =", level,q]);
             Println("bases="); sm1_pmat(bases); 
@@ -1353,26 +1487,246 @@ def Sminimal_v(g) {
             minRes[level+1] = newbases;
             freeRes = minRes;
 #ifdef DEBUG
-/*  Do it later. 
-            for (qq=0; qq<betti; qq++) {
+            for (qq=q; qq<betti; qq++) {
                 for (i=0; i<betti_levelplus; i++) {
-                  if (!IsZero(newbases[i,qq])) {
+                  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
           }
         }
       }
    }
-   return([Stetris(minRes,redundantTable),
+   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 SisComplex(a) {
+  local n,i,j,k,b,p,q;
+  n = Length(a);
+  for (i=0; i<n-1; i++) {
+    if (Length(a[i+1]) != 0) {
+      b = a[i+1]*a[i];
+      p = Length(b); q = Length(b[0]);
+      for (j=0; j<p; j++) {
+        for (k=0; k<q; k++) {
+          if (!IsZero(b[j,k])) {
+             Print("Is is not complex at ");
+             Println([i,j,k]);
+             return(false);
+          }
+        }
+      }
+    }
+  }
+  return(true);
+}
+
+def IsExact_h(c,v) {
+  local a;
+  v = ToString_array(v);
+  a = [c,v];
+  sm1(a," isExact_h /FunctionValue set ");
+}
+HelpAdd(["IsExact_h",
+["IsExact_h(complex,var): bool",
+ "It checks the given complex is exact or not in D<h> (homogenized Weyl algebra)",
+ "cf. ReParse"
+]]);
+
+def ReParse(a) {
+  local c;
+  if (IsArray(a)) {
+    c = Map(a,"ReParse");
+  }else{
+    sm1(a," toString . /c set");
+  }
+  return(c);
+}
+HelpAdd(["ReParse",
+["Reparse(obj): obj",
+ "It parses the given object in the current ring.",
+ "Outputs from SlaScala, Sschreyer may cause a trouble in other functions,",
+ "because it uses the Schreyer order.",
+ "In this case, ReParse the outputs from these functions.",
+ "cf. IsExaxt_h"
+]]);