===================================================================
RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v
retrieving revision 1.9
retrieving revision 1.27
diff -u -p -r1.9 -r1.27
--- OpenXM/src/k097/lib/minimal/minimal.k	2000/05/06 13:41:12	1.9
+++ OpenXM/src/k097/lib/minimal/minimal.k	2000/08/16 22:38:52	1.27
@@ -1,13 +1,28 @@
-/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.8 2000/05/06 10:45:43 takayama Exp $ */
+/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.26 2000/08/10 02:59:08 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
 /* #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,86 @@ 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<m>} 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<m> { (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<m> /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<n; i++) {
+    if (v[i] > t) { t = v[i];}
+  }
+  return(t);
+}
+HelpAdd(["Max",
+["Max(v) returns the maximal element in v."]]);
+
+/*  End of standard functions that should be moved to standard libraries. */
 def test0() {
   local f;
   Sweyl("x,y,z");
@@ -132,15 +228,39 @@ 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,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<n; i++) {
+        if (IsInteger(opt[i])) {
+          count = opt[i];
+        }
+        if (IsString(opt[i])) {
+          if (opt[i] == "homogenized") {
+            nohomog = true;
+          }else if (opt[i] == "Sordinary") {
+            Sordinary = true;
+          }else{
+            Println("Warning: unknown option");
+            Println(opt);
+          }
+        }
+      }
+    } else if (IsNull(opt)){
+    } else {
+      Println("Warning: option should be given by an array.");
+      Println(opt);
+      Println("--------------------------------------------");
+    }
   }
 
   sm1(" setupEnvForResolution ");
@@ -152,7 +272,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 ");
@@ -169,7 +294,7 @@ def SresolutionFrameWithTower(g,opt) {
   /* -sugar is fine? */
   sm1(" setupEnvForResolution "); 
 
-  Println(g);
+  Sprintln(g);
   startingGB = g;
   /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */
   ans = [ ];
@@ -192,12 +317,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]);
 }
 /* ---------------------------- */
@@ -212,16 +338,22 @@ def NewPolynomialVector(size) {
 def  SturnOffHomogenization() {
   sm1("
     [(Homogenize)] system_variable 1 eq
-    { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
+    { Sverbose {
+      (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message } { } ifelse
       [(Homogenize) 0] system_variable
       [(ReduceLowerTerms) 0] system_variable
     } {  } ifelse
   ");
 }
+/* NOTE!!!  Be careful these changes of global environmental variables.
+   We should make a standard set of values and restore these values
+   after computation and interruption.  August 15, 2000.
+*/
 def  SturnOnHomogenization() {
   sm1("
     [(Homogenize)] system_variable 0 eq
-    { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
+    { Sverbose {
+        (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message } {  } ifelse
       [(Homogenize) 1] system_variable
       [(ReduceLowerTerms) 1] system_variable
     } {  } ifelse
@@ -269,7 +401,7 @@ def Sres0FrameWithSkelton(g) {
     si = pair[1,0];
     sj = pair[1,1];
     /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
-    Print(".");
+    Sprint(".");
 
     t_syz = NewPolynomialVector(gLength);
     t_syz[i] = si;
@@ -277,7 +409,7 @@ def Sres0FrameWithSkelton(g) {
     syzAll[k] = t_syz;
   }
   t_syz = syzAll;
-  Print("Done. betti="); Println(betti);
+  Sprint("Done. betti="); Sprintln(betti);
   /* Println(g);  g is in a format such as 
     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
@@ -291,9 +423,15 @@ 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);
 }
 
+HelpAdd(["Sord_w",
+["Sord_w(f,w) returns the w-order of f",
+ "Example: Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]):"]]);
 /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
 def Sord_w(f,w) {
   local neww,i,n;
@@ -327,8 +465,10 @@ def test_SinitOfArray() {
   f = [x^2+y^2+z^2, x*y+x*z+y*z, x*z^2+y*z^2, y^3-x^2*z - x*y*z+y*z^2,
        -y^2*z^2 + x*z^3 + y*z^3, -z^4];
   p=SresolutionFrameWithTower(f);
-  sm1_pmat(p); 
-  sm1_pmat(SgenerateTable(p[1]));
+  if (Sverbose) {
+    sm1_pmat(p); 
+    sm1_pmat(SgenerateTable(p[1]));
+  }
   return(p);
   frame = p[0];
   sm1_pmat(p[1]);
@@ -346,19 +486,16 @@ def Sdegree(f,tower,level) {
   f = Init(f);
   if (level <= 1) return(StotalDegree(f));
   i = Degree(f,es);
-#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) {
   local height, n,i,j, ans, ans_at_each_floor;
+
+  /*
+  Sprint("SgenerateTable: tower=");Sprintln(tower);  
+  sm1(" print_switch_status "); */
   height = Length(tower);
   ans = NewArray(height);
   for (i=0; i<height; i++) {
@@ -434,7 +571,7 @@ 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,
@@ -442,11 +579,19 @@ def SlaScala(g) {
         reductionTable_tmp;
   /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;
-  Print("WeightOfSweyl="); Println(WeightOfSweyl);
-  rf = SresolutionFrameWithTower(g);
+  Sprint("WeightOfSweyl="); Sprintln(WeightOfSweyl);
+  rf = SresolutionFrameWithTower(g,opt);  
+  Sprint("rf="); if (Sverbose) {sm1_pmat(rf);}
   redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];
+
+  Sprintln("Generating reduction table which gives an order of reduction.");
+  Sprint("WeghtOfSweyl="); Sprintln(WeightOfSweyl);
+  Sprint2("tower="); Sprintln2(tower);
   reductionTable = SgenerateTable(tower);
+  Sprint2("reductionTable="); 
+  if (Sverbose || Sverbose2) {sm1_pmat(reductionTable);}
+
   skel = rf[2];
   redundantTable = SnewArrayOfFormat(rf[1]);
   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
@@ -464,11 +609,12 @@ def SlaScala(g) {
       while (SthereIs(reductionTable_tmp,strategy)) {
         i = SnextI(reductionTable_tmp,strategy,redundantTable,
                    skel,level,freeRes);
-        Println([level,i]);
+        Sprintln([level,i]);
         reductionTable_tmp[i] = -200000;
         if (reductionTable[level,i] == strategy) {
-           Print("Processing "); Print([level,i]);
-           Print("   Strategy = "); Println(strategy);
+           Sprint("Processing [level,i]= "); Sprint([level,i]);
+           Sprint("   Strategy = "); Sprintln(strategy);
+           Sprint2(strategy);
            if (level == 0) {
              if (IsNull(redundantTable[level,i])) {
                bases = freeRes[level];
@@ -488,21 +634,21 @@ def SlaScala(g) {
                   place = f[3];
                   /* (level-1, place) is the place for f[0], 
                      which is a newly obtained  GB. */
-#ifdef ORDINARY
+if (Sordinary) {
                   redundantTable[level-1,place] = redundant_seq;
                   redundant_seq++;
-#else
+}else{
                   if (f[4] > f[5]) {
                     /* Zero in the gr-module */
-                    Print("v-degree of [org,remainder] = ");
-                    Println([f[4],f[5]]);
-                    Print("[level,i] = "); Println([level,i]);
+                    Sprint("v-degree of [org,remainder] = ");
+                    Sprintln([f[4],f[5]]);
+                    Sprint("[level,i] = "); Sprintln([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++;
@@ -528,6 +674,7 @@ def SlaScala(g) {
     }
     strategy++;
   }
+  Sprintln2(" ");
   n = Length(freeRes);
   freeResV = SnewArrayOfFormat(freeRes);
   for (i=0; i<n; i++) {
@@ -535,7 +682,7 @@ 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) {
@@ -575,9 +722,9 @@ def SnextI(reductionTable_tmp,strategy,redundantTable,
        }
      }
    }
-   Print("reductionTable_tmp=");
-   Println(reductionTable_tmp);
-   Println("See also reductionTable, strategy, level,i");
+   Sprint("reductionTable_tmp=");
+   Sprintln(reductionTable_tmp);
+   Sprintln("See also reductionTable, strategy, level,i");
    Error("SnextI: bases[i] or bases[j] is null for all combinations.");
 }
 
@@ -611,7 +758,7 @@ def SwhereInGB(f,tower) {
     q = MonomialPart(f);
     if (p == q) return(i);
   }
-  Println([f,tower]);
+  Sprintln([f,tower]);
   Error("whereInGB : [f,myset]: f could not be found in the myset.");      
 }
 def SunitOfFormat(pos,forms) {
@@ -628,15 +775,7 @@ def SunitOfFormat(pos,forms) {
   return(ans);
 }
 
-def Error(s) {
-  sm1(" s error ");
-}
 
-def IsNull(s) {
-  if (Stag(s) == 0) return(true);
-  else return(false);
-}
-
 def StowerOf(tower,level) {
   local ans,i;
   ans = [ ];
@@ -657,10 +796,13 @@ def Sspolynomial(f,g) {
   sm1("f g spol /FunctionValue set");
 }
 
-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);
@@ -670,7 +812,7 @@ def SwhereInTower(f,tower) {
     q = MonomialPart(f);
     if (p == q) return(i);
   }
-  Println([f,tower]);
+  Sprintln([f,tower]);
   Error("[f,tower]: f could not be found in the tower.");      
 }
 
@@ -682,21 +824,23 @@ def SpairAndReduction(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;
-  Println("SpairAndReduction:");
+  Sprintln("SpairAndReduction:");
 
   if (level < 1) Error("level should be >= 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);  
+  Sprintln(["level=",level]);
+  Sprintln(["tower2=",tower2]);
   /** sm1(" show_ring ");   */
 
   gi = Stoes_vec(bases[i]);
@@ -707,35 +851,39 @@ 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];
   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. */
   /* pos2 is the place to put a new GB at level-1. */   
-  Println(ans);
+  Sprintln(ans);
   return(ans);
 }
 
@@ -768,24 +916,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;
@@ -843,15 +973,41 @@ 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)",
+ "      : \"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],",
+ "             [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,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("Betti numbers are ------");
+  if (Sverbose || Sverbose2) {sm1_pmat(bettiTable);}
   minRes = SnewArrayOfFormat(freeRes);
   seq = 0;
   maxSeq = SgetMaxSeq(redundantTable);
@@ -861,12 +1017,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<betti; q++) {
         if (redundantTable[level,q] == seq) {
-          Print("[seq,level,q]="); Println([seq,level,q]);
+          Sprint("[seq,level,q]="); Sprintln([seq,level,q]);
           if (level < maxLevel-1) {
             bases = freeRes[level+1];
             dr = reducer[level,q];
@@ -879,10 +1035,10 @@ def Sminimal(g) {
             for (i=0; i<betti_levelplus; i++) {
               newbases[i] = 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);
+            Sprintln(["level, q =", level,q]);
+            Sprintln("bases="); if (Sverbose) {sm1_pmat(bases); }
+            Sprintln("dr="); if (Sverbose) {sm1_pmat(dr);}
+            Sprintln("newbases="); if (Sverbose) {sm1_pmat(newbases);}
             minRes[level+1] = newbases;
             freeRes = minRes;
 #ifdef DEBUG
@@ -892,7 +1048,7 @@ def Sminimal(g) {
                 for (i=0; i<betti_levelplus; i++) {
                   if (!IsZero(newbases[i,qq])) {
                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
-                    Print("redundantTable ="); sm1_pmat(redundantTable[level]);
+                    Sprint("redundantTable ="); sm1_pmat(redundantTable[level]);
                     Error("Stop in Sminimal for debugging.");
                   }
                 }
@@ -904,10 +1060,28 @@ def Sminimal(g) {
       }
     }
    }
-   return([Stetris(minRes,redundantTable),
-          [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
+   tminRes = Stetris(minRes,redundantTable);
+   ansSminimal = [SpruneZeroRow(tminRes), tminRes,
+                  [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]];
+   Sprintln2(" ");
+   Println("------------ Note -----------------------------");
+   Println("To get shift vectors, use Reparse and SgetShifts(resmat,w)");
+   Println("To get initial of the complex, use Reparse and Sinit_w(resmat,w)");
+   Println("0: minimal resolution, 3: Schreyer resolution ");
+   Println("------------ Resolution Summary  --------------");
+   Print("Betti numbers : "); 
+   Println(Map(ansSminimal[0],"Length"));
+   Print("Betti numbers of the Schreyer frame: ");
+   Println(Map(ansSminimal[3],"Length"));
+   Println("-----------------------------------------------");
+
+   sm1(" restoreEnvAfterResolution ");
+   Sordinary = false;
+
+   return(ansSminimal);
   /* r[4] is the redundantTable_ordinary */
   /* r[0] is the freeResolution */
+  /* r[5] is the skelton */
 }  
 
 
@@ -976,9 +1150,11 @@ def Stetris(freeRes,redundantTable) {
     }else{
       newbases = bases;
     }
-    Println(["level=", level]);
-    sm1_pmat(bases);
-    sm1_pmat(newbases);
+    Sprintln(["level=", level]);
+    if (Sverbose){
+      sm1_pmat(bases);
+      sm1_pmat(newbases);
+    }
 
     minRes[level] = newbases;
   }
@@ -1040,369 +1216,352 @@ 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,"h",1]]); */
   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); 
   pp = Map(p,"Spoly");
-  return(Sminimal_v(pp));
-  /* return(Sminimal(pp)); */
+  return(Sminimal(pp)); 
 }
 
-/* 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));
-}
+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,"Spoly");
-  return(Sminimal_v(pp));
+  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.
+*/
 
-/*  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);
+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);"]]);
 
-  height = Length(reductionTable);
-  for (level = 0; level < height; level++) {
-      n = Length(reductionTable[level]);
-      for (i=0; i<n; i++) {
-           Println([level,i]);
-           Print("Processing "); Print([level,i]);
-           if (level == 0) {
-             if (IsNull(redundantTable[level,i])) {
-               bases = freeRes[level];
-               /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
-               pos = SwhereInGB(tower[0,i],rf[3,0]);
-               bases[i] = rf[3,0,pos];
-               /* redundantTable[level,i] = 0;
-               redundantTable_ordinary[level,i] = 0; */
-               freeRes[level] = bases;
-               /* Println(["GB=",i,bases,tower[0,i]]); */
-             }
-           }else{ /* level >= 1 */
-             if (IsNull(redundantTable[level,i])) {
-               bases = freeRes[level];
-               f = SpairAndReduction2(skel,level,i,freeRes,tower,
-                                      ww,redundantTable);
-               if (f[0] != Poly("0")) {
-                  place = f[3];
-                  /* (level-1, place) is the place for f[0], 
-                     which is a newly obtained  GB. */
-#ifdef ORDINARY
-                  redundantTable[level-1,place] = redundant_seq;
-                  redundant_seq++;
-#else
-                  if (f[4] > 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<nn;ii++) {
-                     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]-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];
-                  bases[i] = f[1];  /* Put the syzygy. */
-                  freeRes[level] = bases;
-               }
-             }  /* 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);
+def IsZeroVector(m) {
+  local n,i;
+  n = Length(m);
   for (i=0; i<n; i++) {
-    bases = freeRes[i];
-    bases = Sbases_to_vec(bases,bettiTable[i]);
-    freeResV[i] = bases;
+    if (!IsZero(m[i])) { 
+      return(false);
+    }
   }
+  return(true);
+}
 
-  /* Mark the non-redundant elements. */
+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++) {
-    m = Length(redundantTable[i]);
-    for (jj=0; jj<m; jj++) {
-      if (IsNull(redundantTable[i,jj])) {
-        redundantTable[i,jj] = 0;
+    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);
+}
 
-  return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
+def ToString_array(p) {
+  local ans;
+  if (IsArray(p)) {
+    ans = Map(p,"ToString_array");
+  }else{
+    ans = ToString(p);
+  }
+  return(ans);
 }
 
-def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
-  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 : -------------------------");
+/* sm1_res_div([[x],[y]],[[x^2],[x*y],[y^2]],[x,y]): */
 
-  if (level < 1) Error("level should be >= 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.");
-  }
+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 ");
+}
 
-  tower2 = StowerOf(tower,level-1);
-  SsetTower(tower2);  
-  /** sm1(" show_ring ");   */
+/* 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 ");
+}
 
-  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]);
+def sm1_gb(f,v) {
+  f =ToString_array(f);
+  v = ToString_array(v);
+  sm1(" [f v] gb /FunctionValue set ");
+}
 
-  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);
+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);
+}
 
-  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.");
+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 IsSameIdeal_h(ii,jj,v) {
+  local a;
+  v = ToString_array(v);
+  a = [ii,jj,v];
+  sm1(a," isSameIdeal_h /FunctionValue set ");
+}
+HelpAdd(["IsSameIdeal_h",
+["IsSameIdeal_h(ii,jj,var): bool",
+ "It checks the given ideals are the same 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");
   }
-  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;
+  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"
+]]);
 
-  c2 = null;
-  /* tmp[0] must be zero */
-  n = Length(t_syz);
-  for (i=0; i<n; 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] = c2*Stoes_vec(freeRes[level-1,i]);
-         t_syz[i] = 0;
-         /* tmp[0] = t_syz . g */
-         /* break; does not work. Use */
-         i = n;
-       }
-      }
-     }
+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);
+}
 
-  /* 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);
+def SgetShift(mat,w,m) {
+  local omat;
+  sm1(" mat { w m ord_w<m> {(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]):"]]);
 
-  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("  ");
+def SgetShifts(resmat,w) {
+  local i,n,ans,m0;
+  n = Length(resmat);
+  ans = NewArray(n);
+  m0 = NewArray(Length(resmat[0,0]));
+  ans[0] = m0;
+  for (i=0; i<n-1; i++) {
+    ans[i+1] = SgetShift(resmat[i],w,m0);
+    m0 = ans[i+1];
+  }
   return(ans);
 }
+HelpAdd(["SgetShifts",
+["SgetShifts(resmat,w) returns the shift vectors of the resolution resmat",
+ " with respect to w with the shift m.",
+ "Note that the order of the ring and the weight w must be the same.",
+ "Zero row is not allowed.",
+ "Example:   a=Sannfs2(\"x^3-y^2\");",
+ "           b=a[0]; w = [\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1];",
+ "           Sweyl(\"x,y\",[w]); b = Reparse(b);",
+ "           SgetShifts(b,w):"]]);
 
-def Sminimal_v(g) {
-  local r, freeRes, redundantTable, reducer, maxLevel,
-        minRes, seq, maxSeq, level, betti, q, bases, dr,
-        betti_levelplus, newbases, i, j,qq;
-  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];
+def Sinit_w(resmat,w) {
+  local shifts,ans,n,i,m,mat,j;
+  shifts = SgetShifts(resmat,w);
+  n = Length(resmat);
+  ans = NewArray(n);
+  for (i=0; i<n; i++) {
+    m = shifts[i];
+    mat = ScopyArray(resmat[i]);
+    for (j=0; j<Length(mat); j++) {
+      mat[j] = Init_w_m(mat[j],w,m);
+    }
+    ans[i] = mat;
   }
-  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
-          }
-        }
-      }
-   }
-   return([Stetris(minRes,redundantTable),
-          [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
-  /* r[4] is the redundantTable_ordinary */
-  /* r[0] is the freeResolution */
-}  
+  return(ans);
+}
+HelpAdd(["Sinit_w",
+["Sinit_w(resmat,w) returns the initial of the complex resmat with respect to the weight w.",
+ "Example:   a=Sannfs2(\"x^3-y^2\");",
+ "           b=a[0]; w = [\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1];",
+ "           Sweyl(\"x,y\",[w]); b = Reparse(b);",
+ "           c=Sinit_w(b,w); c:"
+]]);
 
-/* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */
+/* This method does not work, because we have zero rows. 
+   Think about it later. */
+def SbettiTable(rtable) {
+  local ans,i,j,pp;
+  ans = SnewArrayOfFormat(rtable);
+  for (i=0; i<Length(rtable); i++) {
+    pp = 0;
+    for (j=0; j<Length(rtable[i]); j++) {
+       if (rtable[i,j] != 0) {pp = pp+1;}
+    }
+    ans[i] = pp;
+  }
+  return(ans);
+}
\ No newline at end of file