/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.11 2000/05/19 11:16:51 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"];;

   a=Sminimal(v);
   b=a[0];
   b[1]*b[0]:
   b[2]*b[1]:

   a = test0();
   b = a[0];
   b[1]*b[0]:
   b[2]*b[1]:
   a = Sminimal(b[0]);

   a = test1();
   b=a[0];
   b[1]*b[0]:
   b[2]*b[1]:
     
*/


load("cohom.k");
def load_tower() {
  if (Boundp("k0-tower.sm1.loaded")) {
  }else{
    sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
    sm1(" /k0-tower.sm1.loaded 1 def ");
  }
  sm1(" oxNoX ");
}
load_tower();
SonAutoReduce = true;
def Factor(f) {
   sm1(f, " fctr /FunctionValue set");
}
def Reverse(f) {
   sm1(f," reverse /FunctionValue set");
}
def Sgroebner(f) {
   sm1(" [f] groebner /FunctionValue set");
}
def test0() {
  local f;
  Sweyl("x,y,z");
  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];
  frame=SresolutionFrame(f);
  Println(frame);
  /* return(frame); */
  return(SlaScala(f));
}
def test1() {
  local f;
  Sweyl("x,y,z");
  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];
  return(Sminimal(f));
}

  

def Sweyl(v,w) {
  /* extern WeightOfSweyl ; */
  local ww,i,n;
  if(Length(Arglist) == 1) {
    sm1(" [v s_ring_of_differential_operators 0 [(schreyer) 1]] define_ring ");
    sm1(" define_ring_variables ");

    sm1(" [ v to_records pop ] /ww set ");
    n = Length(ww);
    WeightOfSweyl = NewArray(n*4);
    for (i=0; i< n; i++) {
      WeightOfSweyl[2*i] = ww[i];
      WeightOfSweyl[2*i+1] = 1;
    }
    for (i=0; i< n; i++) {
      WeightOfSweyl[2*n+2*i] = AddString(["D",ww[i]]);
      WeightOfSweyl[2*n+2*i+1] = 1;
    }
      
  }else{
    sm1(" [v s_ring_of_differential_operators w s_weight_vector 0 [(schreyer) 1]] define_ring ");
    sm1(" define_ring_variables ");
    WeightOfSweyl = w[0];
  }
}


def Spoly(f) {
  sm1(f, " toString tparse /FunctionValue set ");
}

def SreplaceZeroByZeroPoly(f) {
  if (IsArray(f)) {
     return(Map(f,"SreplaceZeroByZeroPoly"));
  }else{
     if (IsInteger(f)) {
       return(Poly(ToString(f)));
     }else{
       return(f);
     }
  }
}
def Shomogenize(f) {
  f = SreplaceZeroByZeroPoly(f);
  if (IsArray(f)) {
    sm1(f," sHomogenize2  /FunctionValue set "); 
    /* sm1(f," {sHomogenize2} map  /FunctionValue set ");  */
    /* Is it correct? Double check.*/
  }else{
    sm1(f, " sHomogenize /FunctionValue set ");
  }
}

def StoTower() {
  sm1("  [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv ");
}

def SsetTower(tower) {
sm1(" [(AvoidTheSameRing)] pushEnv 
      [ [(AvoidTheSameRing) 0] system_variable 
        [(gbListTower) tower (list) dc] system_variable
      ] pop popEnv ");
}

def SresolutionFrameWithTower(g,opt) {
  local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
        gbasis;
  if (Length(Arglist) >= 2) {
    if (IsInteger(opt)) count = opt;
  }else{
    count = -1;
  }

  sm1(" setupEnvForResolution ");
  /* If I do not put this macro, homogenization
     make a strange behavior. For example,
     [(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
     [(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
     4/19, 2000.
  */

  sm1(" (mmLarger) (matrix) switch_function ");
  g = Map(g,"Shomogenize");
  if (SonAutoReduce) {
    sm1("[ (AutoReduce) ] system_variable /autof set ");
    sm1("[ (AutoReduce) 1 ] system_variable ");
  }
  gbasis = Sgroebner(g);
  g = gbasis[0];
  if (SonAutoReduce) {
    sm1("[ (AutoReduce) autof] system_variable  "); 
  }

  g = Init(g);

/*  sm1(" setupEnvForResolution-sugar "); */
  /* -sugar is fine? */
  sm1(" setupEnvForResolution "); 

  Println(g);
  startingGB = g;
  /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */
  ans = [ ];
  gbTower = [ ];
  skelton = [ ];
  while (true) {
    /* sm1(g," res0Frame /ff set "); */
    withSkel = Sres0FrameWithSkelton(g);
    ff = withSkel[0];
    ans = Append(ans, ff[0]);
    gbTower = Join([ ff[1] ], gbTower);
    skelton = Join([ withSkel[1] ], skelton);
    g = ff[0];
    if (Length(g) == 0) break;
    SsetTower( gbTower );
    if (count == 0) break;
    count = count - 1;
  }
  return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
}
HelpAdd(["SresolutionFrameWithTower",
["It returs [resolution of the initial, gbTower, skelton, gbasis]",
 "Example: Sweyl(\"x,y\");",
 "         a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);

def SresolutionFrame(f,opt) {
  local ans;
  ans = SresolutionFrameWithTower(f);
  return(ans[0]);
}
/* ---------------------------- */
def ToGradedPolySet(g) {
  sm1(g," (gradedPolySet) dc /FunctionValue set ");
}

def NewPolynomialVector(size) {
  sm1(size," (integer) dc newPolyVector /FunctionValue set ");
}

def  SturnOffHomogenization() {
  sm1("
    [(Homogenize)] system_variable 1 eq
    { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
      [(Homogenize) 0] system_variable
      [(ReduceLowerTerms) 0] system_variable
    } {  } ifelse
  ");
}
def  SturnOnHomogenization() {
  sm1("
    [(Homogenize)] system_variable 0 eq
    { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
      [(Homogenize) 1] system_variable
      [(ReduceLowerTerms) 1] system_variable
    } {  } ifelse
  ");
}

def SschreyerSkelton(g) {
  sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
}
def Stoes(g) {
  if (IsArray(g)) {
    sm1(g," {toes} map /FunctionValue set ");
  }else{
    sm1(g," toes /FunctionValue set ");
  }
}
def Stoes_vec(g) {
    sm1(g," toes /FunctionValue set ");
}

def Sres0Frame(g) {
  local ans;
  ans = Sres0FrameWithSkelton(g);
  return(ans[0]);
}
def Sres0FrameWithSkelton(g) {
  local t_syz, nexttower, m, t_gb, skel, betti,
        gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;

  SturnOffHomogenization();

  g = Stoes(g);
  skel = SschreyerSkelton(g);
  /* Print("Skelton is ");
  sm1_pmat(skel); */
  betti = Length(skel);
  
  gLength = Length(g);
  grG = ToGradedPolySet(g);
  syzAll = NewPolynomialVector(betti);
  for (k=0; k<betti; k++) {
    pair = skel[k];
    i = pair[0,0];
    j = pair[0,1];
    si = pair[1,0];
    sj = pair[1,1];
    /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
    Print(".");

    t_syz = NewPolynomialVector(gLength);
    t_syz[i] = si;
    t_syz[j] = sj;
    syzAll[k] = t_syz;
  }
  t_syz = syzAll;
  Print("Done. betti="); Println(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 , ...]
    [y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
    [3*es^3*y*Dy-es^5*x ] 
  */
  nexttower = Init(g);
  SturnOnHomogenization();
  return([[t_syz, nexttower],skel]);
}


def StotalDegree(f) {
  sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");
}

/* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
def Sord_w(f,w) {
  local neww,i,n;
  n = Length(w);
  neww = NewArray(n);
  for (i=0; i<n; i=i+2) {
    neww[i] = ToString(w[i]);
  }
  for (i=1; i<n; i=i+2) {
    neww[i] = IntegerToSm1Integer(w[i]);
  }
  sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
}


/* This is not satisfactory. */
def SinitOfArray(f) {
  local p,pos,top;
  if (IsArray(f)) {
     sm1(f," toes init /p set ");
     sm1(p," (es). degree (universalNumber) dc /pos set ");
     return([Init(f[pos]),pos]);
  } else {
     return(Init(f));
  }
}

def test_SinitOfArray() {
  local f, frame,p,tower,i,j,k;
  Sweyl("x,y,z");
  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]));
  return(p);
  frame = p[0];
  sm1_pmat(p[1]);
  sm1_pmat(frame);
  sm1_pmat(Map(frame[0],"SinitOfArray"));
  sm1_pmat(Map(frame[1],"SinitOfArray"));
  return(p);
}

/* f is assumed to be a monomial with toes. */
def Sdegree(f,tower,level) {
  local i,ww, wd;
  /* extern WeightOfSweyl; */
  ww = WeightOfSweyl;
  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;
  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)
                            + OFFSET;
      /* Println([i,j,ans_at_each_floor[j]]); */
    }
    ans[i] = ans_at_each_floor;
  }
  return(ans);
}
Sweyl("x,y,z");
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]];
/*  SresolutionFrameWithTower(v); */

def SnewArrayOfFormat(p) {
  if (IsArray(p)) {
     return(Map(p,"SnewArrayOfFormat"));
  }else{
     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 */
  if (IsArray(a)) {
    n = Length(a);
    for (i=0; i<n; i++) {
      if (IsArray(a[i])) {
        tt = SminOfStrategy(a[i]);
        if (tt < ans) ans = tt;
      }else{
        if (a[i] < ans) ans = a[i];
      }
    }
  }else{
     if (a < ans) ans = a;
  }
  return(ans);
}  
def SmaxOfStrategy(a) {
  local n,i,ans,tt;
  ans = -100000; /* very small number */
  if (IsArray(a)) {
    n = Length(a);
    for (i=0; i<n; i++) {
      if (IsArray(a[i])) {
        tt = SmaxOfStrategy(a[i]);
        if (tt > ans) ans = tt;
      }else{
        if (a[i] > ans) ans = a[i];
      }
    }
  }else{
     if (a > ans) ans = a;
  }
  return(ans);
}  


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,
        reductionTable_tmp;
  /* extern WeightOfSweyl; */
  ww = WeightOfSweyl;
  Print("WeightOfSweyl="); 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);

  strategy = SminOfStrategy( reductionTable );
  maxOfStrategy = SmaxOfStrategy( reductionTable );
  height = Length(reductionTable);
  while (strategy <= maxOfStrategy) {
    for (level = 0; level < height; level++) {
      n = Length(reductionTable[level]);
      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);
           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 = SpairAndReduction(skel,level,i,freeRes,tower,ww);
               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. */
                  freeRes[level] = bases;
                  bases = freeRes[level-1];
                  bases[place] = f[0];
                  freeRes[level-1] = bases;
                  reducer[level-1,place] = f[1];
               }else{
                  redundantTable[level,i] = 0;
                  bases = freeRes[level];
                  bases[i] = f[1];  /* Put the syzygy. */
                  freeRes[level] = bases;
               }
             }
           } /* end of level >= 1 */
        }
      }
    }
    strategy++;
  }
  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;
  }
  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);
  n = Length(freeRes);
  if (IsArray(g[0])) {
    ans[0] = Length(g[0]);
  }else{
    ans[0] = 1;
  }
  for (level=0; level<n; level++) {
    bases = freeRes[level];
    if (IsArray(bases)) {
      ans[level+1] = Length(bases);
    }else{
      ans[level+1] = 1;
    }
  }
  return(ans);
}

def SwhereInGB(f,tower) {
  local i,n,p,q;
  n = Length(tower);
  for (i=0; i<n; i++) {
    p = MonomialPart(tower[i]);
    q = MonomialPart(f);
    if (p == q) return(i);
  }
  Println([f,tower]);
  Error("whereInGB : [f,myset]: f could not be found in the myset.");      
}
def SunitOfFormat(pos,forms) {
  local ans,i,n;
  n = Length(forms);
  ans = NewArray(n);
  for (i=0; i<n; i++) {
    if (i != pos) {
      ans[i] = Poly("0");
    }else{
      ans[i] = Poly("1");
    }
  }
  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 = [ ];
  if (level == 0) return([[]]);
  for (i=0; i<level; i++) {
    ans = Append(ans,tower[i]);
  } 
  return(Reverse(ans));
}

def Sspolynomial(f,g) {
  if (IsArray(f)) {
    f = Stoes_vec(f);
  }
  if (IsArray(g)) {
    g = Stoes_vec(g);
  }
  sm1("f g spol /FunctionValue set");
}

def MonomialPart(f) {
  sm1(" [(lmonom) f] gbext /FunctionValue set ");
}

def SwhereInTower(f,tower) {
  local i,n,p,q;
  if (f == Poly("0")) return(-1);
  n = Length(tower);
  for (i=0; i<n; i++) {
    p = MonomialPart(tower[i]);
    q = MonomialPart(f);
    if (p == q) return(i);
  }
  Println([f,tower]);
  Error("[f,tower]: f could not be found in the tower.");      
}

def Stag(f) {
  sm1(f," tag (universalNumber) dc /FunctionValue set");
}

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:");

  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);

  /* 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);

  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;
  pos = SwhereInTower(syzHead,tower[level]);
  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);
  return(ans);
}

def Sreduction(f,myset) {
  local n, indexTable, set2, i, j, tmp, t_syz;
  n = Length(myset);
  indexTable = NewArray(n);
  set2 = [ ];
  j = 0;
  for (i=0; i<n; i++) {
    if (IsNull(myset[i])) {
      indexTable[i] = -1;
/*    }else if (myset[i] == Poly("0")) {
      indexTable[i] = -1;  */
    }else{
      set2 = Append(set2,Stoes_vec(myset[i]));
      indexTable[i] = j;
      j++;
    }
  }
  sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
  t_syz = NewArray(n);
  for (i=0; i<n; i++) {
    if (indexTable[i] != -1) {
      t_syz[i] = tmp[2, indexTable[i]];
    }else{
      t_syz[i] = Poly("0");
    }
  }
  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;
  if (Length(Arglist) < 2) size = -1;
  if (IsArray(f)) return(f);
  r = RingOf(f);
  myes = PolyR("es",r);
  myee = PolyR("e_",r);
  if (Degree(f,myee) > 0 && size == -1) {
    if (size == -1) {
       sm1(f," (array) dc /ans set");
       return(ans);
    }
  }

/*
    Coefficients(x^2-1,x):
    [    [    2 , 0 ]  , [    1 , -1 ]  ] 
*/
  if (Degree(f,myee) > 0) {
    c = Coefficients(f,myee);
  }else{
    c = Coefficients(f,myes);
  }
  if (size < 0) {
    size = c[0,0]+1;
  }
  ans = NewArray(size);
  for (i=0; i<size; i++) {ans[i] = 0;}
  n = Length(c[0]);
  for (j=0; j<n; j++) {
    d = c[0,j];
    ans[d] = c[1,j];
  }
  return(ans);
}

def Sbases_to_vec(bases,size) {
  local n, giveSize, newbases,i;
  /*  bases = [1+es*x, [1,2,3*x]] */
  if (Length(Arglist) > 1) {
    giveSize = true;
  }else{
    giveSize = false;
  }
  n = Length(bases);
  newbases = NewArray(n);
  for (i=0; i<n; i++) {
     if (giveSize) {
       newbases[i] = Sfrom_es(bases[i], size);
     }else{
       newbases[i] = Sfrom_es(bases[i]);
     }
  }
  return(newbases);
}

def Sminimal(g) {
  local r, freeRes, redundantTable, reducer, maxLevel,
        minRes, seq, maxSeq, level, betti, q, bases, dr,
        betti_levelplus, newbases, i, j,qq;
  r = SlaScala(g);
  /* 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];
  }
  seq=maxSeq+1;
  while (seq > 1) {
    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]);
          if (level < maxLevel-1) {
            bases = freeRes[level+1];
            dr = reducer[level,q];
            dr[q] = -1;
            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;
            }
            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=0; qq<betti; qq++) {
              if ((redundantTable[level,qq] >= seq) && 
                  (redundantTable[level,qq] <= maxSeq)) {
                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]);
                    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 */
}  


def IsZero(f) {
  if (IsPolynomial(f)) {
    return( f == Poly("0"));
  }else if (IsInteger(f)) {
    return( f == 0);
  }else if (IsSm1Integer(f)) {
    return( f == true );
  }else if (IsDouble(f)) {
    return( f == 0.0 );
  }else if (IsRational(f)) {
    return(IsZero(Denominator(f)));
  }else{
    Error("IsZero: cannot deal with this data type.");
  }
}
def SgetMaxSeq(redundantTable) {
   local level,i,n,ans, levelMax,bases;
   levelMax = Length( redundantTable );
   ans = 0;
   for (level = 0; level < levelMax; level++) {
     bases = redundantTable[level];
     n = Length(bases);
     for (i=0; i<n; i++) {
       if (IsInteger( bases[i] )) {
          if (bases[i] > ans) {
             ans = bases[i];
          }
       }
     }
   }
   return(ans);
}

def Stetris(freeRes,redundantTable) {
  local level, i, j, resLength, minRes,
        bases, newbases, newbases2;
  minRes = SnewArrayOfFormat(freeRes);
  resLength = Length( freeRes );
  for (level=0; level<resLength; level++) {
    bases = freeRes[level];
    newbases = SnewArrayOfFormat(bases);
    betti = Length(bases); j = 0;
    /* Delete rows */
    for (i=0; i<betti; i++) {
      if (redundantTable[level,i] < 1) {
         newbases[j] = bases[i];
         j++;
      }
    }
    bases = SfirstN(newbases,j);
    if (level > 0) {
      /* Delete columns */
      newbases = Transpose(bases);
      betti = Length(newbases); j = 0;
      newbases2 = SnewArrayOfFormat(newbases);
      for (i=0; i<betti; i++) {
        if (redundantTable[level-1,i] < 1) {
           newbases2[j] = newbases[i];
           j++;
        }
      }
      newbases = Transpose(SfirstN(newbases2,j));
    }else{
      newbases = bases;
    }
    Println(["level=", level]);
    sm1_pmat(bases);
    sm1_pmat(newbases);

    minRes[level] = newbases;
  }
  return(minRes);    
}

def SfirstN(bases,k) {
   local ans,i;
   ans = NewArray(k);
   for (i=0; i<k; i++) {
     ans[i] = bases[i];
   }
   return(ans);
}


/* usage:  tt is tower. ww is weight.
    a = SresolutionFrameWithTower(v);
    tt = a[1];    
    ww = [x,1,y,1,Dx,1,Dy,1];
    SvDegree(x*es,tt,1,ww):

In(17)=tt:
[[2*x*Dx , e_*x^2 , e_*x*y , 3*x^2*Dy , e_*y^3 , 9*x*y*Dy^2 , 27*y^2*Dy^3 ]  ,
 [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ]  ,
 [3*es^3*y*Dy ]  ] 
In(18)=SvDegree(x*es,tt,1,ww):
3
In(19)=SvDegree(x*es^3,tt,1,ww):
4
In(20)=SvDegree(x,tt,2,ww):
4

*/
def SvDegree(f,tower,level,w) {
  local i,ans;
  if (IsZero(f)) return(null);
  f = Init(f);
  if (level <= 0) {
    return(Sord_w(f,w));
  }
  i = Degree(f,es);
  ans = Sord_w(f,w) +
        SvDegree(tower[level-1,i],tower,level-1,w);
  return(ans);
}

def Sannfs(f,v) {
  local f2;
  f2 = ToString(f);
  if (IsArray(v)) {
     v = Map(v,"ToString");
  }
  sm1(" [f2 v] annfs /FunctionValue set ");
}

/* Sannfs2("x^3-y^2"); */
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)); */
}

HelpAdd(["Sannfs2",
["Sannfs2(f) constructs the V-minimal free resolution for the weight (-1,1)",
 "of the Laplace transform of the annihilating ideal of the polynomial f in x,y.",
 "See also Sminimal_v, Sannfs3.",
 "Example: a=Sannfs2(\"x^3-y^2\");",
 "         b=a[0]; sm1_pmat(b);",
 "         b[1]*b[0]:",
 "Example: a=Sannfs2(\"x*y*(x-y)*(x+y)\");",
 "         b=a[0]; sm1_pmat(b);",
 "         b[1]*b[0]:"
]]);

/* Do not forget to turn on TOTAL_STRATEGY */ 
def Sannfs2_laScala(f) {
  local p,pp;
  p = Sannfs(f,"x,y");
  /*   Do not make laplace transform. 
    sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
    p = [p];
  */
  Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); 
  pp = Map(p[0],"Spoly");
  return(Sminimal(pp));
}

def Sannfs2_laScala2(f) {
  local p,pp;
  p = Sannfs(f,"x,y");
  sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
  p = [p];
  Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
               ["x",-1,"y",-1,"Dx",1,"Dy",1]]); 
  pp = Map(p[0],"Spoly");
  return(Sminimal(pp));
}

def Sannfs3(f) {
  local p,pp;
  p = Sannfs(f,"x,y,z");
  sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
  Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
  pp = Map(p,"Spoly");
  return(Sminimal_v(pp));
}

HelpAdd(["Sannfs3",
["Sannfs3(f) constructs the V-minimal free resolution for the weight (-1,1)",
 "of the Laplace transform of the annihilating ideal of the polynomial f in x,y,z.",
 "See also Sminimal_v, Sannfs2.",
 "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
 "         b=a[0]; sm1_pmat(b);",
 "         b[1]*b[0]: b[2]*b[1]:"]]);

/*
  The betti numbers of most examples are 2,1. (0-th and 1-th).
  a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2. 
  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.

*/
     
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. */
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);

  if (!IsNull(vdeg_reduced)) {
    if (vdeg_reduced < vdeg) {
      Println("--- Special in V-minimal!");
      Println(tmp[0]);
      Println("syzygy="); sm1_pmat(t_syz);
      Print("[vdeg, vdeg_reduced] = "); Println([vdeg,vdeg_reduced]);
    }
  }


  pos = SwhereInTower(syzHead,tower[level]);
  pos2 = SwhereInTower(tmp[0],tower[level-1]);
  ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
  /* pos is the place to put syzygy at level. */
  /* pos2 is the place to put a new GB at level-1. */   
  Println(ans);
  Println("  ");
  return(ans);
}

HelpAdd(["Sminimal_v",
["It constructs the V-minimal free resolution from the Schreyer resolution",
 "step by step.",
 "Example:   Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);", 
 "          v=[[2*x*Dx + 3*y*Dy+6, 0],",
 "             [3*x^2*Dy + 2*y*Dx, 0],",
 "             [0,  x^2+y^2],",
 "             [0,  x*y]];",
 "         a=Sminimal_v(v);",
 "         sm1_pmat(a[0]); b=a[0]; b[1]*b[0]:",
 "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);


def Sminimal_v(g) {
  local r, freeRes, redundantTable, reducer, maxLevel,
        minRes, seq, maxSeq, level, betti, q, bases, dr,
        betti_levelplus, newbases, i, j,qq,tminRes;
  r = Sschreyer(g);
  sm1_pmat(r);
  Debug_Sminimal_v = r;
  Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
  /* Should I turn off the tower?? */
  freeRes = r[0];
  redundantTable = r[1];
  reducer = r[2];
  minRes = SnewArrayOfFormat(freeRes);
  seq = 0;
  maxSeq = SgetMaxSeq(redundantTable);
  maxLevel = Length(freeRes);
  for (level = 0; level < maxLevel; level++) {
    minRes[level] = freeRes[level];
  }
  for (level = 0; level < maxLevel; level++) {
      betti = Length(freeRes[level]);
      for (q = betti-1; q>=0; q--) {
        if (redundantTable[level,q] > 0) {
          Print("[seq,level,q]="); Println([seq,level,q]);
          if (level < maxLevel-1) {
            bases = freeRes[level+1];
            dr = reducer[level,q];
            /* dr[q] = -1;  We do not need this in our reducer format. */
            /* dr[q] should be a non-zero constant. */
            newbases = SnewArrayOfFormat(bases);
            betti_levelplus = Length(bases);
            /*
               bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
            */
            for (i=0; i<betti_levelplus; i++) {
              newbases[i] = dr[q]*bases[i] - bases[i,q]*dr;
            }
            Println(["level, q =", level,q]);
            Println("bases="); sm1_pmat(bases); 
            Println("dr="); sm1_pmat(dr);
            Println("newbases="); sm1_pmat(newbases);
            minRes[level+1] = newbases;
            freeRes = minRes;
#ifdef DEBUG
            for (qq=q; qq<betti; qq++) {
                for (i=0; i<betti_levelplus; i++) {
                  if ((!IsZero(newbases[i,qq])) && (redundantTable[level,qq] >0)) {
                    Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
                    Print("redundantTable ="); sm1_pmat(redundantTable[level]);
                    Error("Stop in Sminimal for debugging.");
                  }
                }
            }
#endif
          }
        }
      }
   }
   tminRes = Stetris(minRes,redundantTable);
   return([SpruneZeroRow(tminRes), tminRes,
          [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
  /* r[4] is the redundantTable_ordinary */
  /* r[0] is the freeResolution */
}  

/* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */
/* x y (x+y-1)(x-2),  x^3-y^2, x^3 - y^2 z^2,
   x y z (x+y+z-1) seems to be interesting, because the first syzygy
  contains 1.
*/

def CopyArray(m) {
  local ans,i,n;
  if (IsArray(m)) {
     n = Length(m);
     ans = NewArray(n);
     for (i=0; i<n; i++) {
       ans[i] = CopyArray(m[i]);
     }
     return(ans);
  }else{
     return(m);
  }
}
HelpAdd(["CopyArray",
["It duplicates the argument array recursively.",
 "Example: m=[1,[2,3]];",
 "         a=CopyArray(m); a[1] = \"Hello\";",
 "         Println(m); Println(a);"]]);

def IsZeroVector(m) {
  local n,i;
  n = Length(m);
  for (i=0; i<n; i++) {
    if (!IsZero(m[i])) { 
      return(false);
    }
  }
  return(true);
}

def SpruneZeroRow(res) {
  local minRes, n,i,j,m, base,base2,newbase,newbase2, newMinRes;

  minRes = CopyArray(res);
  n = Length(minRes);
  for (i=0; i<n; i++) {
    base = minRes[i];
    m = Length(base);
    if (i != n-1) {
      base2 = minRes[i+1];
      base2 = Transpose(base2);
    }
    newbase = [ ];
    newbase2 = [ ];
    for (j=0; j<m; j++) {
      if (!IsZeroVector(base[j])) {
        newbase = Append(newbase,base[j]);
        if (i != n-1) {
          newbase2 = Append(newbase2,base2[j]);
        }
      }
    }
    minRes[i] = newbase;
    if (i != n-1) {
      if (newbase2 == [ ]) {
        minRes[i+1] = [ ];
      }else{
        minRes[i+1] = Transpose(newbase2);
      }
    }
  }
  
  newMinRes = [ ];
  n = Length(minRes);
  i = 0;
  while (i < n ) {
    base = minRes[i];
    if (base == [ ]) {
      i = n; /* break; */
    }else{
      newMinRes = Append(newMinRes,base);
    }
    i++;
  }
  return(newMinRes);
}

def testAnnfs2(f) {
  local a,i,n;
  a = Sannfs2(f);
  b=a[0];
  n = Length(b);
  Println("------ V-minimal free resolution -----");
  sm1_pmat(b);
  Println("----- Is it complex?  ---------------");
  for (i=0; i<n-1; i++) {
    Println(b[i+1]*b[i]);
  }
  return(a);
}
def testAnnfs3(f) {
  local a,i,n;
  a = Sannfs3(f);
  b=a[0];
  n = Length(b);
  Println("------ V-minimal free resolution -----");
  sm1_pmat(b);
  Println("----- Is it complex?  ---------------");
  for (i=0; i<n-1; i++) {
    Println(b[i+1]*b[i]);
  }
  return(a);
}

def ToString_array(p) {
  local ans;
  if (IsArray(p)) {
    ans = Map(p,"ToString_array");
  }else{
    ans = ToString(p);
  }
  return(ans);
}

/* sm1_res_div([[x],[y]],[[x^2],[x*y],[y^2]],[x,y]): */

def sm1_res_div(I,J,V) {
  I = ToString_array(I);
  J = ToString_array(J);
  V = ToString_array(V);
  sm1(" [[ I J]  V ] res*div /FunctionValue set ");
}

/* It has not yet been working */
def sm1_res_kernel_image(m,n,v) {
  m = ToString_array(m);
  n = ToString_array(n);
  v = ToString_array(v);
  sm1(" [m n v] res-kernel-image /FunctionValue set ");
}
def Skernel(m,v) {
  m = ToString_array(m);
  v = ToString_array(v);
  sm1(" [ m v ] syz /FunctionValue set ");
}

def test3() {
  local a1,a2,b1,b2;
  a1 = Sannfs3("x^3-y^2*z^2");
  a1 = a1[0];
  a2 = Sannfs3_laScala2("x^3-y^2*z^2");
  a2 = a2[0];
  b1 = a1[1];
  b2 = a2[1];
  sm1_pmat(b2);
  Println("  OVER ");
  sm1_pmat(b1);
  return([sm1_res_div(b2,b1,["x","y","z"]),b2,b1,a2,a1]);
}

def test4() {
  local a,b;
  a = Sannfs3_laScala2("x^3-y^2*z^2");
  b = a[0];
  sm1_pmat( sm1_res_kernel_image(b[0],b[1],[x,y,z]));
  sm1_pmat( sm1_res_kernel_image(b[1],b[2],[x,y,z]));
  return(a);
}

def sm1_gb(f,v) {
  f =ToString_array(f);
  v = ToString_array(v);
  sm1(" [f v] gb /FunctionValue set ");
}

def test5() {
  local a,b,c,cc,v;
  a = Sannfs3_laScala2("x^3-y^2*z^2");
  b = a[0];
  v = [x,y,z];
  c = Skernel(b[0],v);
  c = c[0];
  sm1_pmat([c,b[1],v]);
  Println("-----------------------------------");
  cc = sm1_res_div(c,b[1],v);
  sm1_pmat(sm1_gb(cc,v));
  c = Skernel(b[1],v);
  c = c[0];
  cc = sm1_res_div(c,b[2],v);
  sm1_pmat(sm1_gb(cc,v));
  return(a);
}
def test6() {
  local a,b,c,cc,v;
  a = Sannfs3("x^3-y^2*z^2");
  b = a[0];
  v = [x,y,z];
  c = Skernel(b[0],v);
  c = c[0];
  sm1_pmat([c,b[1],v]);
  Println("-------ker = im for minimal ?---------------------");
  cc = sm1_res_div(c,b[1],v);
  sm1_pmat(sm1_gb(cc,v));
  c = Skernel(b[1],v);
  c = c[0];
  cc = sm1_res_div(c,b[2],v);
  sm1_pmat(sm1_gb(cc,v));
  Println("------ ker=im for Schreyer ?------------------");
  b = a[3];
  c = Skernel(b[0],v);
  c = c[0];
  sm1_pmat([c,b[1],v]);
  cc = sm1_res_div(c,b[1],v);
  sm1_pmat(sm1_gb(cc,v));
  c = Skernel(b[1],v);
  c = c[0];
  cc = sm1_res_div(c,b[2],v);
  sm1_pmat(sm1_gb(cc,v));
  return(a);
}