===================================================================
RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v
retrieving revision 1.3
retrieving revision 1.9
diff -u -p -r1.3 -r1.9
--- OpenXM/src/k097/lib/minimal/minimal.k	2000/05/04 06:55:28	1.3
+++ OpenXM/src/k097/lib/minimal/minimal.k	2000/05/06 13:41:12	1.9
@@ -1,6 +1,13 @@
-/* $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.8 2000/05/06 10:45:43 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"];;
 
@@ -30,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;
@@ -332,10 +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) {
@@ -346,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;
@@ -367,6 +387,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 */
@@ -409,10 +438,11 @@ def SlaScala(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;
+        redundantTable_ordinary, redundant_seq_ordinary,
+        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];
@@ -430,7 +460,12 @@ 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("   Strategy = "); Println(strategy);
@@ -503,6 +538,51 @@ def SlaScala(g) {
   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
 }
 
+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);
@@ -959,7 +1039,26 @@ 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]]);
+  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)); */
+}
+
+/* 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));
 }
@@ -967,9 +1066,10 @@ def Sannfs2(f) {
 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));
 }
 
 /*
@@ -980,4 +1080,329 @@ def Sannfs3(f) {
 
 */
      
- 
+
+
+/*  The below does not use LaScala-Stillman's algorithm. */
+def Sschreyer(g) {
+  local rf, tower, reductionTable, skel, redundantTable, bases,
+        strategy, maxOfStrategy, height, level, n, i,
+        freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
+        redundantTable_ordinary, redundant_seq_ordinary,
+        reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
+  /* extern WeightOfSweyl; */
+  ww = WeightOfSweyl;
+  Print("WeghtOfSweyl="); Println(WeightOfSweyl);
+  rf = SresolutionFrameWithTower(g);
+  redundant_seq = 1;   redundant_seq_ordinary = 1;
+  tower = rf[1];
+  reductionTable = SgenerateTable(tower);
+  skel = rf[2];
+  redundantTable = SnewArrayOfFormat(rf[1]);
+  redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
+  reducer = SnewArrayOfFormat(rf[1]);
+  freeRes = SnewArrayOfFormat(rf[1]);
+  bettiTable = SsetBettiTable(rf[1],g);
+
+  height = Length(reductionTable);
+  for (level = 0; level < height; level++) {
+      n = Length(reductionTable[level]);
+      for (i=0; i<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);
+  for (i=0; i<n; i++) {
+    bases = freeRes[i];
+    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]);
+}
+
+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 : -------------------------");
+
+  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.");
+  }
+
+  tower2 = StowerOf(tower,level-1);
+  SsetTower(tower2);  
+  /** sm1(" show_ring ");   */
+
+  gi = Stoes_vec(bases[i]);
+  gj = Stoes_vec(bases[j]);
+
+  ssp = Sspolynomial(gi,gj); 
+  si = ssp[0,0];  
+  sj = ssp[0,1];
+  syzHead = si*es^i;
+  /* This will be the head term, I think. But, double check. */
+  Println([si*es^i,sj*es^j]);
+
+  Print("[gi, gj] = "); Println([gi,gj]);
+  sm1(" [(Homogenize)] system_variable message ");
+  Print("Reduce the element "); Println(si*gi+sj*gj);
+  Print("by  "); Println(bases);
+
+  tmp = Sreduction(si*gi+sj*gj, bases);
+
+  Print("result is "); Println(tmp);
+  if (!IsZero(tmp[0])) {
+    Print("Error: base = ");
+    Println(Map(bases,"Stoes_vec"));
+    Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
+  }
+  t_syz = tmp[2];
+  si = si*tmp[1]+t_syz[i];
+  sj = sj*tmp[1]+t_syz[j];
+  t_syz[i] = si; 
+  t_syz[j] = sj;
+
+  c2 = null;
+  /* tmp[0] must be zero */
+  n = Length(t_syz);
+  for (i=0; i<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;
+       }
+      }
+     }
+  }
+
+  /* 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);
+
+  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);
+}
+
+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];
+  }
+  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 */
+}  
+
+/* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */