===================================================================
RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v
retrieving revision 1.3
retrieving revision 1.18
diff -u -p -r1.3 -r1.18
--- OpenXM/src/k097/lib/minimal/minimal.k	2000/05/04 06:55:28	1.3
+++ OpenXM/src/k097/lib/minimal/minimal.k	2000/07/30 02:26:25	1.18
@@ -1,6 +1,12 @@
-/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.2 2000/05/03 07:50:38 takayama Exp $ */
+/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.17 2000/07/26 12:56:36 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 OFFSET 20*/
 /* Test sequences. 
    Use load["minimal.k"];;
 
@@ -30,6 +36,7 @@ def load_tower() {
     sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
     sm1(" /k0-tower.sm1.loaded 1 def ");
   }
+  sm1(" oxNoX ");
 }
 load_tower();
 SonAutoReduce = true;
@@ -124,13 +131,25 @@ sm1(" [(AvoidTheSameRing)] pushEnv 
       [ [(AvoidTheSameRing) 0] system_variable 
         [(gbListTower) tower (list) dc] system_variable
       ] pop popEnv ");
+      /* sm1("(hoge) message show_ring "); */
 }
 
 def SresolutionFrameWithTower(g,opt) {
   local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
-        gbasis;
+        gbasis, nohomog;
+  nohomog = false;
+  count = -1;
   if (Length(Arglist) >= 2) {
-    if (IsInteger(opt)) count = opt;
+    if (IsInteger(opt)) {
+      count = opt;
+    }else if (IsString(opt)) {
+      if (opt == "homogenized") {
+         nohomog = true;
+      }else{
+         Println("Warning: unknown option");
+         Println(opt);
+      }
+    }
   }else{
     count = -1;
   }
@@ -144,7 +163,12 @@ def SresolutionFrameWithTower(g,opt) {
   */
 
   sm1(" (mmLarger) (matrix) switch_function ");
-  g = Map(g,"Shomogenize");
+  if (! nohomog) {
+    Println("Automatic homogenization.");
+    g = Map(g,"Shomogenize");
+  }else{
+    Println("No automatic homogenization.");
+  }
   if (SonAutoReduce) {
     sm1("[ (AutoReduce) ] system_variable /autof set ");
     sm1("[ (AutoReduce) 1 ] system_variable ");
@@ -184,12 +208,13 @@ def SresolutionFrameWithTower(g,opt) {
 }
 HelpAdd(["SresolutionFrameWithTower",
 ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
+ "option: \"homogenized\" (no automatic homogenization) ",
  "Example: Sweyl(\"x,y\");",
  "         a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
 
 def SresolutionFrame(f,opt) {
   local ans;
-  ans = SresolutionFrameWithTower(f);
+  ans = SresolutionFrameWithTower(f,opt);
   return(ans[0]);
 }
 /* ---------------------------- */
@@ -283,7 +308,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]); */
@@ -332,21 +360,30 @@ 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));
+  return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1)); 
+    
 }
 
 def SgenerateTable(tower) {
   local height, n,i,j, ans, ans_at_each_floor;
+
+  /*
+  Print("SgenerateTable: tower=");Println(tower);  
+  sm1(" print_switch_status "); */
   height = Length(tower);
   ans = NewArray(height);
   for (i=0; i<height; i++) {
     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;
@@ -367,6 +404,15 @@ def SnewArrayOfFormat(p) {
      return(null);
   }
 }
+def ScopyArray(a) {
+  local n, i,ans;
+  n = Length(a);
+  ans = NewArray(n);
+  for (i=0; i<n; i++) {
+    ans[i] = a[i];
+  }
+  return(ans);
+}
 def SminOfStrategy(a) {
   local n,i,ans,tt;
   ans = 100000; /* very big number */
@@ -405,18 +451,26 @@ def SmaxOfStrategy(a) {
 }  
 
 
-def SlaScala(g) {
+def SlaScala(g,opt) {
   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;
+        redundantTable_ordinary, redundant_seq_ordinary,
+        reductionTable_tmp;
   /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;
-  Print("WeghtOfSweyl="); Println(WeightOfSweyl);
-  rf = SresolutionFrameWithTower(g);
+  Print("WeightOfSweyl="); Println(WeightOfSweyl);
+  rf = SresolutionFrameWithTower(g,opt);
+  Print("rf="); sm1_pmat(rf);
   redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];
+
+  Println("Generating reduction table which gives an order of reduction.");
+  Print("WeghtOfSweyl="); Println(WeightOfSweyl);
+  Print("tower"); Println(tower);
   reductionTable = SgenerateTable(tower);
+  Print("reductionTable="); sm1_pmat(reductionTable);
+
   skel = rf[2];
   redundantTable = SnewArrayOfFormat(rf[1]);
   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
@@ -430,9 +484,14 @@ def SlaScala(g) {
   while (strategy <= maxOfStrategy) {
     for (level = 0; level < height; level++) {
       n = Length(reductionTable[level]);
-      for (i=0; i<n; i++) {
+      reductionTable_tmp = ScopyArray(reductionTable[level]); 
+      while (SthereIs(reductionTable_tmp,strategy)) {
+        i = SnextI(reductionTable_tmp,strategy,redundantTable,
+                   skel,level,freeRes);
+        Println([level,i]);
+        reductionTable_tmp[i] = -200000;
         if (reductionTable[level,i] == strategy) {
-           Print("Processing "); Print([level,i]);
+           Print("Processing [level,i]= "); Print([level,i]);
            Print("   Strategy = "); Println(strategy);
            if (level == 0) {
              if (IsNull(redundantTable[level,i])) {
@@ -500,9 +559,54 @@ def SlaScala(g) {
     bases = Sbases_to_vec(bases,bettiTable[i]);
     freeResV[i] = bases;
   }
-  return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
+  return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary,rf]);
 }
 
+def SthereIs(reductionTable_tmp,strategy) {
+  local n,i;
+  n = Length(reductionTable_tmp);
+  for (i=0; i<n; i++) {
+    if (reductionTable_tmp[i] == strategy) {
+      return(true);
+    }
+  }
+  return(false);
+}
+
+def SnextI(reductionTable_tmp,strategy,redundantTable,
+                                  skel,level,freeRes)
+{
+   local ii,n,p,myindex,i,j,bases;
+   n = Length(reductionTable_tmp);
+   if (level == 0) {
+     for (ii=0; ii<n; ii++) {
+       if (reductionTable_tmp[ii] == strategy) {
+          return(ii);
+        }
+      }
+   }else{
+     for (ii=0; ii<n; ii++) {
+       if (reductionTable_tmp[ii] == strategy) {
+         p = skel[level,ii];
+         myindex = p[0];
+         i = myindex[0]; j = myindex[1];
+         bases = freeRes[level-1];
+         if (IsNull(bases[i]) || IsNull(bases[j])) {
+            
+         }else{
+           return(ii);
+         }
+       }
+     }
+   }
+   Print("reductionTable_tmp=");
+   Println(reductionTable_tmp);
+   Println("See also reductionTable, strategy, level,i");
+   Error("SnextI: bases[i] or bases[j] is null for all combinations.");
+}
+
+
+
 def SsetBettiTable(freeRes,g) {
   local level,i, n,bases,ans;
   ans = NewArray(Length(freeRes)+1);
@@ -581,6 +685,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);
@@ -617,6 +727,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]);
@@ -650,7 +762,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. */
@@ -763,12 +879,31 @@ def Sbases_to_vec(bases,size) {
   return(newbases);
 }
 
-def Sminimal(g) {
+HelpAdd(["Sminimal",
+["It constructs the V-minimal free resolution by LaScala's algorithm",
+ "option: \"homogenized\" (no automatic homogenization ",
+ "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,opt) {
   local r, freeRes, redundantTable, reducer, maxLevel,
         minRes, seq, maxSeq, level, betti, q, bases, dr,
-        betti_levelplus, newbases, i, j,qq;
-  r = SlaScala(g);
+        betti_levelplus, newbases, i, j,qq, tminRes;
+  if (Length(Arglist) < 2) {
+     opt = null;
+  }
+  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];
@@ -824,10 +959,12 @@ def Sminimal(g) {
       }
     }
    }
-   return([Stetris(minRes,redundantTable),
-          [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
+   tminRes = Stetris(minRes,redundantTable);
+   return([SpruneZeroRow(tminRes), tminRes,
+          [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]]);
   /* r[4] is the redundantTable_ordinary */
   /* r[0] is the freeResolution */
+  /* r[5] is the skelton */
 }  
 
 
@@ -959,25 +1096,270 @@ def Sannfs(f,v) {
 def Sannfs2(f) {
   local p,pp;
   p = Sannfs(f,"x,y");
-  Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
-  pp = Map(p[0],"Spoly");
-  return(Sminimal(pp));
+  sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
+  Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); 
+  pp = Map(p,"Spoly");
+  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, 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]:"
+]]);
+/* Some samples.
+  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");    
+  a=Sannfs2("x*y*(x-y)");    
+*/
+
+
 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");
+  pp = Map(p,"Spoly");
   return(Sminimal(pp));
 }
 
-/*
-  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.
+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, Sannfs2.",
+ "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
+ "         b=a[0]; sm1_pmat(b);",
+ "         b[1]*b[0]: b[2]*b[1]:"]]);
 
-*/
      
- 
+
+/* 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 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"
+]]);
+
+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 "); Println(s);
+     Error("Schreyer order is not set.");
+  }
+  /* More check will be necessary. */
+  return(true);
+}
\ No newline at end of file