[BACK]Return to minimal.k CVS log [TXT][DIR] Up to [local] / OpenXM / src / k097 / lib / minimal

Diff for /OpenXM/src/k097/lib/minimal/minimal.k between version 1.16 and 1.35

version 1.16, 2000/06/15 07:38:36 version 1.35, 2007/07/03 22:05:46
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.15 2000/06/14 07:44:05 takayama Exp $ */  /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.34 2001/01/05 11:14:28 takayama Exp $ */
 #define DEBUG 1  #define DEBUG 1
 /* #define ORDINARY 1 */  Sordinary = false;
 /* If you run this program on openxm version 1.1.2 (FreeBSD),  /* If you run this program on openxm version 1.1.2 (FreeBSD),
    make a symbolic link by the command     make a symbolic link by the command
    ln -s /usr/bin/cpp /lib/cpp     ln -s /usr/bin/cpp /lib/cpp
 */  */
 #define OFFSET 0  #define OFFSET 0
 #define TOTAL_STRATEGY 1  
 /* #define OFFSET 20*/  /* #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.  /* Test sequences.
    Use load["minimal.k"];;     Use load["minimal.k"];;
   
Line 29 
Line 44 
   
 */  */
   
   /* We cannot use load command in the if statement. */
   load("lib/minimal/cohom.k");
   Load_sm1(["k0-tower.sm1","lib/minimal/k0-tower.sm1"],"k0-tower.sm1.loaded");
   Load_sm1(["new.sm1","lib/minimal/new.sm1"],"new.sm1.loaded");
   sm1(" oxNoX ");
   
 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;  SonAutoReduce = true;
 def Factor(f) {  def Factor(f) {
    sm1(f, " fctr /FunctionValue set");     sm1(f, " fctr /FunctionValue set");
Line 50  def Reverse(f) {
Line 60  def Reverse(f) {
 def Sgroebner(f) {  def Sgroebner(f) {
    sm1(" [f] groebner /FunctionValue set");     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."]]);
   
   def Kernel(f,v) {
     local ans;
     /* v :  string or ring */
     if (Length(Arglist) < 2) {
       sm1(" [f] syz /ans set ");
     }else{
       sm1(" [f v] syz /ans set ");
     }
     return(ans);
   }
   def Syz(f) {
     sm1(" [f] syz /FunctionValue set ");
   }
   HelpAdd(["Kernel",
   ["Kernel(f) returns the syzygy of f.",
    "Return value [b, c]: b is a set of generators of the syzygies of f",
    "                   : c=[gb, backward transformation, syzygy without",
    "                                                   dehomogenization",
    "Example:  Weyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ",
    "          s=Kernel([x*Dx+1,Dx^2+x^5]); s[0]:"]]);
   /* cf. sm1_syz in cohom.k */
   def Gb(f) {
     sm1(" [f] gb /FunctionValue set ");
   }
   HelpAdd(["Gb",
   ["Gb(f) returns the Groebner basis of f.",
    "cf. Kernel, Weyl."]]);
   
   
   /*  End of standard functions that should be moved to standard libraries. */
 def test0() {  def test0() {
   local f;    local f;
   Sweyl("x,y,z");    Sweyl("x,y,z");
Line 69  def test1() {
Line 188  def test1() {
 }  }
   
   
   
 def Sweyl(v,w) {  def Sweyl(v,w) {
   /* extern WeightOfSweyl ; */    /* extern WeightOfSweyl ; */
   local ww,i,n;    local ww,i,n;
Line 128  def StoTower() {
Line 246  def StoTower() {
 }  }
   
 def SsetTower(tower) {  def SsetTower(tower) {
 sm1(" [(AvoidTheSameRing)] pushEnv  sm1(" [(AvoidTheSameRing)] pushEnv \
       [ [(AvoidTheSameRing) 0] system_variable        [ [(AvoidTheSameRing) 0] system_variable \
         [(gbListTower) tower (list) dc] system_variable          [(gbListTower) tower (list) dc] system_variable \
       ] pop popEnv ");        ] pop popEnv ");
       /* sm1("(hoge) message show_ring "); */        /* sm1("(hoge) message show_ring "); */
 }  }
   
 def SresolutionFrameWithTower(g,opt) {  def SresolutionFrameWithTower(g,opt) {
   local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,    local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
         gbasis, nohomog;          gbasis, nohomog,i,n;
     /* extern Sordinary */
   nohomog = false;    nohomog = false;
   count = -1;    count = -1;  Sordinary = false; /* default value for options. */
   if (Length(Arglist) >= 2) {    if (Length(Arglist) >= 2) {
     if (IsInteger(opt)) {      if (IsArray(opt)) {
       count = opt;        n = Length(opt);
     }else if (IsString(opt)) {        for (i=0; i<n; i++) {
       if (opt == "homogenized") {          if (IsInteger(opt[i])) {
          nohomog = true;            count = opt[i];
       }else{          }
          Println("Warning: unknown option");          if (IsString(opt[i])) {
          Println(opt);            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("--------------------------------------------");
     }      }
   }else{  
     count = -1;  
   }    }
   
   sm1(" setupEnvForResolution ");    sm1(" setupEnvForResolution ");
Line 186  def SresolutionFrameWithTower(g,opt) {
Line 316  def SresolutionFrameWithTower(g,opt) {
   /* -sugar is fine? */    /* -sugar is fine? */
   sm1(" setupEnvForResolution ");    sm1(" setupEnvForResolution ");
   
   Println(g);    Sprintln(g);
   startingGB = g;    startingGB = g;
   /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */    /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */
   ans = [ ];    ans = [ ];
Line 228  def NewPolynomialVector(size) {
Line 358  def NewPolynomialVector(size) {
 }  }
   
 def  SturnOffHomogenization() {  def  SturnOffHomogenization() {
   sm1("    sm1(" \
     [(Homogenize)] system_variable 1 eq      [(Homogenize)] system_variable 1 eq \
     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message      { Sverbose { \
       [(Homogenize) 0] system_variable        (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message } { } ifelse \
       [(ReduceLowerTerms) 0] system_variable        [(Homogenize) 0] system_variable \
     } {  } ifelse        [(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() {  def  SturnOnHomogenization() {
   sm1("    sm1(" \
     [(Homogenize)] system_variable 0 eq      [(Homogenize)] system_variable 0 eq \
     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message      { Sverbose { \
       [(Homogenize) 1] system_variable          (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message } {  } ifelse \
       [(ReduceLowerTerms) 1] system_variable        [(Homogenize) 1] system_variable \
     } {  } ifelse        [(ReduceLowerTerms) 1] system_variable \
       } {  } ifelse \
   ");    ");
 }  }
   
Line 287  def Sres0FrameWithSkelton(g) {
Line 423  def Sres0FrameWithSkelton(g) {
     si = pair[1,0];      si = pair[1,0];
     sj = pair[1,1];      sj = pair[1,1];
     /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */      /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
     Print(".");      Sprint(".");
   
     t_syz = NewPolynomialVector(gLength);      t_syz = NewPolynomialVector(gLength);
     t_syz[i] = si;      t_syz[i] = si;
Line 295  def Sres0FrameWithSkelton(g) {
Line 431  def Sres0FrameWithSkelton(g) {
     syzAll[k] = t_syz;      syzAll[k] = t_syz;
   }    }
   t_syz = syzAll;    t_syz = syzAll;
   Print("Done. betti="); Println(betti);    Sprint("Done. betti="); Sprintln(betti);
   /* Println(g);  g is in a format such as    /* 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 , ...]
     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]      [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
Line 315  def StotalDegree(f) {
Line 451  def StotalDegree(f) {
   return(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]); */  /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
 def Sord_w(f,w) {  def Sord_w(f,w) {
   local neww,i,n;    local neww,i,n;
Line 348  def test_SinitOfArray() {
Line 487  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,    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];         -y^2*z^2 + x*z^3 + y*z^3, -z^4];
   p=SresolutionFrameWithTower(f);    p=SresolutionFrameWithTower(f);
   sm1_pmat(p);    if (Sverbose) {
   sm1_pmat(SgenerateTable(p[1]));      sm1_pmat(p);
       sm1_pmat(SgenerateTable(p[1]));
     }
   return(p);    return(p);
   frame = p[0];    frame = p[0];
   sm1_pmat(p[1]);    sm1_pmat(p[1]);
Line 367  def Sdegree(f,tower,level) {
Line 508  def Sdegree(f,tower,level) {
   f = Init(f);    f = Init(f);
   if (level <= 1) return(StotalDegree(f));    if (level <= 1) return(StotalDegree(f));
   i = Degree(f,es);    i = Degree(f,es);
 #ifdef TOTAL_STRATEGY  
   return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));    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));  
   
 }  }
   
Line 382  def SgenerateTable(tower) {
Line 516  def SgenerateTable(tower) {
   local height, n,i,j, ans, ans_at_each_floor;    local height, n,i,j, ans, ans_at_each_floor;
   
   /*    /*
   Print("SgenerateTable: tower=");Println(tower);    Sprint("SgenerateTable: tower=");Sprintln(tower);
   sm1(" print_switch_status "); */    sm1(" print_switch_status "); */
   height = Length(tower);    height = Length(tower);
   ans = NewArray(height);    ans = NewArray(height);
Line 467  def SlaScala(g,opt) {
Line 601  def SlaScala(g,opt) {
         reductionTable_tmp;          reductionTable_tmp;
   /* extern WeightOfSweyl; */    /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;    ww = WeightOfSweyl;
   Print("WeightOfSweyl="); Println(WeightOfSweyl);    Sprint("WeightOfSweyl="); Sprintln(WeightOfSweyl);
   rf = SresolutionFrameWithTower(g,opt);    rf = SresolutionFrameWithTower(g,opt);
   Print("rf="); sm1_pmat(rf);    Sprint("rf="); if (Sverbose) {sm1_pmat(rf);}
   redundant_seq = 1;   redundant_seq_ordinary = 1;    redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];    tower = rf[1];
   
   Println("Generating reduction table which gives an order of reduction.");    Sprintln("Generating reduction table which gives an order of reduction.");
   Print("WeghtOfSweyl="); Println(WeightOfSweyl);    Sprint("WeghtOfSweyl="); Sprintln(WeightOfSweyl);
   Print("tower"); Println(tower);    Sprint2("tower="); Sprintln2(tower);
   reductionTable = SgenerateTable(tower);    reductionTable = SgenerateTable(tower);
   Print("reductionTable="); sm1_pmat(reductionTable);    Sprint2("reductionTable=");
     if (Sverbose || Sverbose2) {sm1_pmat(reductionTable);}
   
   skel = rf[2];    skel = rf[2];
   redundantTable = SnewArrayOfFormat(rf[1]);    redundantTable = SnewArrayOfFormat(rf[1]);
Line 496  def SlaScala(g,opt) {
Line 631  def SlaScala(g,opt) {
       while (SthereIs(reductionTable_tmp,strategy)) {        while (SthereIs(reductionTable_tmp,strategy)) {
         i = SnextI(reductionTable_tmp,strategy,redundantTable,          i = SnextI(reductionTable_tmp,strategy,redundantTable,
                    skel,level,freeRes);                     skel,level,freeRes);
         Println([level,i]);          Sprintln([level,i]);
         reductionTable_tmp[i] = -200000;          reductionTable_tmp[i] = -200000;
         if (reductionTable[level,i] == strategy) {          if (reductionTable[level,i] == strategy) {
            Print("Processing [level,i]= "); Print([level,i]);             Sprint("Processing [level,i]= "); Sprint([level,i]);
            Print("   Strategy = "); Println(strategy);             Sprint("   Strategy = "); Sprintln(strategy);
              Sprint2(strategy);
            if (level == 0) {             if (level == 0) {
              if (IsNull(redundantTable[level,i])) {               if (IsNull(redundantTable[level,i])) {
                bases = freeRes[level];                 bases = freeRes[level];
Line 520  def SlaScala(g,opt) {
Line 656  def SlaScala(g,opt) {
                   place = f[3];                    place = f[3];
                   /* (level-1, place) is the place for f[0],                    /* (level-1, place) is the place for f[0],
                      which is a newly obtained  GB. */                       which is a newly obtained  GB. */
 #ifdef ORDINARY  if (Sordinary) {
                   redundantTable[level-1,place] = redundant_seq;                    redundantTable[level-1,place] = redundant_seq;
                   redundant_seq++;                    redundant_seq++;
 #else  }else{
                   if (f[4] > f[5]) {                    if (f[4] > f[5]) {
                     /* Zero in the gr-module */                      /* Zero in the gr-module */
                     Print("v-degree of [org,remainder] = ");                      Sprint("v-degree of [org,remainder] = ");
                     Println([f[4],f[5]]);                      Sprintln([f[4],f[5]]);
                     Print("[level,i] = "); Println([level,i]);                      Sprint("[level,i] = "); Sprintln([level,i]);
                     redundantTable[level-1,place] = 0;                      redundantTable[level-1,place] = 0;
                   }else{                    }else{
                     redundantTable[level-1,place] = redundant_seq;                      redundantTable[level-1,place] = redundant_seq;
                     redundant_seq++;                      redundant_seq++;
                   }                    }
 #endif  }
                   redundantTable_ordinary[level-1,place]                    redundantTable_ordinary[level-1,place]
                      =redundant_seq_ordinary;                       =redundant_seq_ordinary;
                   redundant_seq_ordinary++;                    redundant_seq_ordinary++;
Line 560  def SlaScala(g,opt) {
Line 696  def SlaScala(g,opt) {
     }      }
     strategy++;      strategy++;
   }    }
     Sprintln2(" ");
   n = Length(freeRes);    n = Length(freeRes);
   freeResV = SnewArrayOfFormat(freeRes);    freeResV = SnewArrayOfFormat(freeRes);
   for (i=0; i<n; i++) {    for (i=0; i<n; i++) {
Line 567  def SlaScala(g,opt) {
Line 704  def SlaScala(g,opt) {
     bases = Sbases_to_vec(bases,bettiTable[i]);      bases = Sbases_to_vec(bases,bettiTable[i]);
     freeResV[i] = bases;      freeResV[i] = bases;
   }    }
   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);    return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary,rf]);
 }  }
   
 def SthereIs(reductionTable_tmp,strategy) {  def SthereIs(reductionTable_tmp,strategy) {
Line 607  def SnextI(reductionTable_tmp,strategy,redundantTable,
Line 744  def SnextI(reductionTable_tmp,strategy,redundantTable,
        }         }
      }       }
    }     }
    Print("reductionTable_tmp=");     Sprint("reductionTable_tmp=");
    Println(reductionTable_tmp);     Sprintln(reductionTable_tmp);
    Println("See also reductionTable, strategy, level,i");     Sprintln("See also reductionTable, strategy, level,i");
    Error("SnextI: bases[i] or bases[j] is null for all combinations.");     Error("SnextI: bases[i] or bases[j] is null for all combinations.");
 }  }
   
Line 643  def SwhereInGB(f,tower) {
Line 780  def SwhereInGB(f,tower) {
     q = MonomialPart(f);      q = MonomialPart(f);
     if (p == q) return(i);      if (p == q) return(i);
   }    }
   Println([f,tower]);    Sprintln([f,tower]);
   Error("whereInGB : [f,myset]: f could not be found in the myset.");    Error("whereInGB : [f,myset]: f could not be found in the myset.");
 }  }
 def SunitOfFormat(pos,forms) {  def SunitOfFormat(pos,forms) {
Line 660  def SunitOfFormat(pos,forms) {
Line 797  def SunitOfFormat(pos,forms) {
   return(ans);    return(ans);
 }  }
   
 def Error(s) {  
   sm1(" s error ");  
 }  
   
 def IsNull(s) {  
   if (Stag(s) == 0) return(true);  
   else return(false);  
 }  
   
 def StowerOf(tower,level) {  def StowerOf(tower,level) {
   local ans,i;    local ans,i;
   ans = [ ];    ans = [ ];
Line 689  def Sspolynomial(f,g) {
Line 818  def Sspolynomial(f,g) {
   sm1("f g spol /FunctionValue set");    sm1("f g spol /FunctionValue set");
 }  }
   
 def MonomialPart(f) {  
   sm1(" [(lmonom) f] gbext /FunctionValue set ");  
 }  
   
 /* WARNING:  /* WARNING:
   When you use SwhereInTower, you have to change gbList    When you use SwhereInTower, you have to change gbList
Line 708  def SwhereInTower(f,tower) {
Line 834  def SwhereInTower(f,tower) {
     q = MonomialPart(f);      q = MonomialPart(f);
     if (p == q) return(i);      if (p == q) return(i);
   }    }
   Println([f,tower]);    Sprintln([f,tower]);
   Error("[f,tower]: f could not be found in the tower.");    Error("[f,tower]: f could not be found in the tower.");
 }  }
   
Line 720  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 846  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   local i, j, myindex, p, bases, tower2, gi, gj,    local i, j, myindex, p, bases, tower2, gi, gj,
        si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,         si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
        vdeg,vdeg_reduced;         vdeg,vdeg_reduced;
   Println("SpairAndReduction:");    Sprintln("SpairAndReduction:");
   
   if (level < 1) Error("level should be >= 1 in SpairAndReduction.");    if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
   p = skel[level,ii];    p = skel[level,ii];
   myindex = p[0];    myindex = p[0];
   i = myindex[0]; j = myindex[1];    i = myindex[0]; j = myindex[1];
   bases = freeRes[level-1];    bases = freeRes[level-1];
   Println(["p and bases ",p,bases]);    Sprintln(["p and bases ",p,bases]);
   if (IsNull(bases[i]) || IsNull(bases[j])) {    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.");      Error("level, i, j : bases[i], bases[j]  must not be NULL.");
   }    }
   
   tower2 = StowerOf(tower,level-1);    tower2 = StowerOf(tower,level-1);
   SsetTower(tower2);    SsetTower(tower2);
   Println(["level=",level]);    Sprintln(["level=",level]);
   Println(["tower2=",tower2]);    Sprintln(["tower2=",tower2]);
   /** sm1(" show_ring ");   */    /** sm1(" show_ring ");   */
   
   gi = Stoes_vec(bases[i]);    gi = Stoes_vec(bases[i]);
Line 747  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 873  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   sj = ssp[0,1];    sj = ssp[0,1];
   syzHead = si*es^i;    syzHead = si*es^i;
   /* This will be the head term, I think. But, double check. */    /* 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]);    Sprint("[gi, gj] = "); Sprintln([gi,gj]);
   sm1(" [(Homogenize)] system_variable message ");    sm1(" [(Homogenize)] system_variable  ");
   Print("Reduce the element "); Println(si*gi+sj*gj);    Sprint("Reduce the element "); Sprintln(si*gi+sj*gj);
   Print("by  "); Println(bases);    Sprint("by  "); Sprintln(bases);
   
   tmp = Sreduction(si*gi+sj*gj, 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. */    /* This is essential part for V-minimal resolution. */
   /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */    /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
   vdeg = SvDegree(si*gi,tower,level-1,ww);    vdeg = SvDegree(si*gi,tower,level-1,ww);
   vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);    vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
   Print("vdegree of the original = "); Println(vdeg);    Sprint("vdegree of the original = "); Sprintln(vdeg);
   Print("vdegree of the remainder = "); Println(vdeg_reduced);    Sprint("vdegree of the remainder = "); Sprintln(vdeg_reduced);
   
   t_syz = tmp[2];    t_syz = tmp[2];
   si = si*tmp[1]+t_syz[i];    si = si*tmp[1]+t_syz[i];
Line 779  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 905  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];    ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
   /* pos is the place to put syzygy at level. */    /* pos is the place to put syzygy at level. */
   /* pos2 is the place to put a new GB at level-1. */    /* pos2 is the place to put a new GB at level-1. */
   Println(ans);    Sprintln(ans);
   return(ans);    return(ans);
 }  }
   
Line 812  def Sreduction(f,myset) {
Line 938  def Sreduction(f,myset) {
   return([tmp[0],tmp[1],t_syz]);    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) {  def Sfrom_es(f,size) {
   local c,ans, i, d, myes, myee, j,n,r,ans2;    local c,ans, i, d, myes, myee, j,n,r,ans2;
Line 888  def Sbases_to_vec(bases,size) {
Line 996  def Sbases_to_vec(bases,size) {
 }  }
   
 HelpAdd(["Sminimal",  HelpAdd(["Sminimal",
 ["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm",  ["It constructs the V-minimal free resolution by LaScala's algorithm",
  "option: \"homogenized\" (no automatic homogenization ",   "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]]);",   "Example:  Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
  "          v=[[2*x*Dx + 3*y*Dy+6, 0],",   "          v=[[2*x*Dx + 3*y*Dy+6, 0],",
  "             [3*x^2*Dy + 2*y*Dx, 0],",   "             [3*x^2*Dy + 2*y*Dx, 0],",
Line 904  HelpAdd(["Sminimal",
Line 1014  HelpAdd(["Sminimal",
 def Sminimal(g,opt) {  def Sminimal(g,opt) {
   local r, freeRes, redundantTable, reducer, maxLevel,    local r, freeRes, redundantTable, reducer, maxLevel,
         minRes, seq, maxSeq, level, betti, q, bases, dr,          minRes, seq, maxSeq, level, betti, q, bases, dr,
         betti_levelplus, newbases, i, j,qq, tminRes;          betti_levelplus, newbases, i, j,qq, tminRes,bettiTable, ansSminimal;
   if (Length(Arglist) < 2) {    if (Length(Arglist) < 2) {
      opt = null;       opt = null;
   }    }
     /* Sordinary is set in SlaScala(g,opt) --> SresolutionFrameWithTower */
   
   ScheckIfSchreyer("Sminimal:0");    ScheckIfSchreyer("Sminimal:0");
   r = SlaScala(g,opt);    r = SlaScala(g,opt);
   /* Should I turn off the tower?? */    /* Should I turn off the tower?? */
Line 915  def Sminimal(g,opt) {
Line 1027  def Sminimal(g,opt) {
   freeRes = r[0];    freeRes = r[0];
   redundantTable = r[1];    redundantTable = r[1];
   reducer = r[2];    reducer = r[2];
     bettiTable = SbettiTable(redundantTable);
     Sprintln2("BettiTable ------");
     if (Sverbose || Sverbose2) {sm1_pmat(bettiTable);}
   minRes = SnewArrayOfFormat(freeRes);    minRes = SnewArrayOfFormat(freeRes);
   seq = 0;    seq = 0;
   maxSeq = SgetMaxSeq(redundantTable);    maxSeq = SgetMaxSeq(redundantTable);
Line 924  def Sminimal(g,opt) {
Line 1039  def Sminimal(g,opt) {
   }    }
   seq=maxSeq+1;    seq=maxSeq+1;
   while (seq > 1) {    while (seq > 1) {
     seq--;      seq--;  Sprint2(seq);
     for (level = 0; level < maxLevel; level++) {      for (level = 0; level < maxLevel; level++) {
       betti = Length(freeRes[level]);        betti = Length(freeRes[level]);
       for (q = 0; q<betti; q++) {        for (q = 0; q<betti; q++) {
         if (redundantTable[level,q] == seq) {          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) {            if (level < maxLevel-1) {
             bases = freeRes[level+1];              bases = freeRes[level+1];
             dr = reducer[level,q];              dr = reducer[level,q];
Line 942  def Sminimal(g,opt) {
Line 1057  def Sminimal(g,opt) {
             for (i=0; i<betti_levelplus; i++) {              for (i=0; i<betti_levelplus; i++) {
               newbases[i] = bases[i] + bases[i,q]*dr;                newbases[i] = bases[i] + bases[i,q]*dr;
             }              }
             Println(["level, q =", level,q]);              Sprintln(["level, q =", level,q]);
             Println("bases="); sm1_pmat(bases);              Sprintln("bases="); if (Sverbose) {sm1_pmat(bases); }
             Println("dr="); sm1_pmat(dr);              Sprintln("dr="); if (Sverbose) {sm1_pmat(dr);}
             Println("newbases="); sm1_pmat(newbases);              Sprintln("newbases="); if (Sverbose) {sm1_pmat(newbases);}
             minRes[level+1] = newbases;              minRes[level+1] = newbases;
             freeRes = minRes;              freeRes = minRes;
 #ifdef DEBUG  #ifdef DEBUG
Line 955  def Sminimal(g,opt) {
Line 1070  def Sminimal(g,opt) {
                 for (i=0; i<betti_levelplus; i++) {                  for (i=0; i<betti_levelplus; i++) {
                   if (!IsZero(newbases[i,qq])) {                    if (!IsZero(newbases[i,qq])) {
                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);                      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.");                      Error("Stop in Sminimal for debugging.");
                   }                    }
                 }                  }
Line 968  def Sminimal(g,opt) {
Line 1083  def Sminimal(g,opt) {
     }      }
    }     }
    tminRes = Stetris(minRes,redundantTable);     tminRes = Stetris(minRes,redundantTable);
    return([SpruneZeroRow(tminRes), tminRes,     ansSminimal = [SpruneZeroRow(tminRes), tminRes,
           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);                    [ 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(Join([Length(ansSminimal[0,0,0])],Map(ansSminimal[0],"Length")));
      Print("Betti numbers of the Schreyer frame: ");
      Println(Join([Length(ansSminimal[3,0,0])],Map(ansSminimal[3],"Length")));
      Println("-----------------------------------------------");
   
      sm1(" restoreEnvAfterResolution ");
      Sordinary = false;
   
      return(ansSminimal);
   /* r[4] is the redundantTable_ordinary */    /* r[4] is the redundantTable_ordinary */
   /* r[0] is the freeResolution */    /* r[0] is the freeResolution */
     /* r[5] is the skelton */
 }  }
   
   
Line 1040  def Stetris(freeRes,redundantTable) {
Line 1172  def Stetris(freeRes,redundantTable) {
     }else{      }else{
       newbases = bases;        newbases = bases;
     }      }
     Println(["level=", level]);      Sprintln(["level=", level]);
     sm1_pmat(bases);      if (Sverbose){
     sm1_pmat(newbases);        sm1_pmat(bases);
         sm1_pmat(newbases);
       }
   
     minRes[level] = newbases;      minRes[level] = newbases;
   }    }
Line 1104  def Sannfs2(f) {
Line 1238  def Sannfs2(f) {
   local p,pp;    local p,pp;
   p = Sannfs(f,"x,y");    p = Sannfs(f,"x,y");
   sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");    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]]);    Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
   pp = Map(p,"Spoly");    pp = Map(p,"Spoly");
   return(Sminimal_v(pp));    return(Sminimal(pp));
   /* return(Sminimal(pp)); */  
 }  }
   
 HelpAdd(["Sannfs2",  HelpAdd(["Sannfs2",
 ["Sannfs2(f) constructs the V-minimal free resolution for the weight (-1,1)",  ["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.",   "of the Laplace transform of the annihilating ideal of the polynomial f in x,y.",
  "See also Sminimal_v, Sannfs3.",   "See also Sminimal, Sannfs3.",
  "Example: a=Sannfs2(\"x^3-y^2\");",   "Example: a=Sannfs2(\"x^3-y^2\");",
  "         b=a[0]; sm1_pmat(b);",   "         b=a[0]; sm1_pmat(b);",
  "         b[1]*b[0]:",   "         b[1]*b[0]:",
Line 1126  HelpAdd(["Sannfs2",
Line 1254  HelpAdd(["Sannfs2",
  "         b=a[0]; sm1_pmat(b);",   "         b=a[0]; sm1_pmat(b);",
  "         b[1]*b[0]:"   "         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)");
   */
   
 /* 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) {  def Sannfs3(f) {
   local p,pp;    local p,pp;
   p = Sannfs(f,"x,y,z");    p = Sannfs(f,"x,y,z");
   sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");    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]]);    Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
   pp = Map(p,"Spoly");    pp = Map(p,"Spoly");
   return(Sminimal_v(pp));    return(Sminimal(pp));
 }  }
   
 HelpAdd(["Sannfs3",  HelpAdd(["Sannfs3",
 ["Sannfs3(f) constructs the V-minimal free resolution for the weight (-1,1)",  ["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.",   "of the Laplace transform of the annihilating ideal of the polynomial f in x,y,z.",
  "See also Sminimal_v, Sannfs2.",   "See also Sminimal, Sannfs2.",
  "Example: a=Sannfs3(\"x^3-y^2*z^2\");",   "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
  "         b=a[0]; sm1_pmat(b);",   "         b=a[0]; sm1_pmat(b);",
  "         b[1]*b[0]: b[2]*b[1]:"]]);   "         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];  
   Println("Generating reduction table which gives an order of reduction.");  
   Println("But, you are in Sschreyer...., you may not use LaScala-Stillman");  
   Print("WeghtOfSweyl="); Println(WeightOfSweyl);  
   Print("tower"); Println(tower);  
   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);  
   Println(["level=",level]);  
   Println(["tower2=",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]);  
     }  
   }  
   
   SsetTower(StowerOf(tower,level));  
   pos = SwhereInTower(syzHead,tower[level]);  
   
   SsetTower(StowerOf(tower,level-1));  
   pos2 = SwhereInTower(tmp[0],tower[level-1]);  
   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];  
   /* pos is the place to put syzygy at level. */  
   /* pos2 is the place to put a new GB at level-1. */  
   Println(ans);  
   Println("--- end of SpairAndReduction2  ");  
   return(ans);  
 }  
   
 HelpAdd(["Sminimal_v",  
 ["It constructs the V-minimal free resolution from the Schreyer resolution",  
  "step by step.",  
  "This code still contains bugs. It sometimes outputs wrong answer.",  
  "Example:   Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",  
  "          v=[[2*x*Dx + 3*y*Dy+6, 0],",  
  "             [3*x^2*Dy + 2*y*Dx, 0],",  
  "             [0,  x^2+y^2],",  
  "             [0,  x*y]];",  
  "         a=Sminimal_v(v);",  
  "         sm1_pmat(a[0]); b=a[0]; b[1]*b[0]:",  
  "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);  
   
 /* This code still contains bugs. It sometimes outputs wrong answer. */  
 /* See test12() in minimal-test.k.  */  
 /* There may be remaining 1, too */  
 def Sminimal_v(g) {  
   local r, freeRes, redundantTable, reducer, maxLevel,  
         minRes, seq, maxSeq, level, betti, q, bases, dr,  
         betti_levelplus, newbases, i, j,qq,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 */  /* 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 (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     x y z (x+y+z-1) seems to be interesting, because the first syzygy
Line 1687  def Skernel(m,v) {
Line 1423  def Skernel(m,v) {
   sm1(" [ m v ] syz /FunctionValue set ");    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) {  def sm1_gb(f,v) {
   f =ToString_array(f);    f =ToString_array(f);
   v = ToString_array(v);    v = ToString_array(v);
Line 1750  HelpAdd(["IsExact_h",
Line 1464  HelpAdd(["IsExact_h",
  "cf. ReParse"   "cf. ReParse"
 ]]);  ]]);
   
 def ReParse(a) {  def IsSameIdeal_h(ii,jj,v) {
   local c;    local a;
   if (IsArray(a)) {    v = ToString_array(v);
     c = Map(a,"ReParse");    a = [ii,jj,v];
   }else{    sm1(a," isSameIdeal_h /FunctionValue set ");
     sm1(a," toString . /c set");  
   }  
   return(c);  
 }  }
 HelpAdd(["ReParse",  HelpAdd(["IsSameIdeal_h",
 ["Reparse(obj): obj",  ["IsSameIdeal_h(ii,jj,var): bool",
  "It parses the given object in the current ring.",   "It checks the given ideals are the same or not in D<h> (homogenized Weyl algebra)",
  "Outputs from SlaScala, Sschreyer may cause a trouble in other functions,",   "cf. ReParse"
  "because it uses the Schreyer order.",  
  "In this case, ReParse the outputs from these functions.",  
  "cf. IsExaxt_h"  
 ]]);  ]]);
   
   /*
     Output of S* functions may cause a trouble because it uses Schreyer orders.
     In this case, use ReParse().
   */
   
 def ScheckIfSchreyer(s) {  def ScheckIfSchreyer(s) {
   local ss;    local ss;
   sm1(" (report) (grade) switch_function /ss set ");    sm1(" (report) (grade) switch_function /ss set ");
Line 1784  def ScheckIfSchreyer(s) {
Line 1497  def ScheckIfSchreyer(s) {
   */    */
   sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set ");    sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set ");
   if (ss != 1) {    if (ss != 1) {
      Print("ScheckIfSchreyer: from "); Println(s);       Print("ScheckIfSchreyer: from "); Printl(s);
      Error("Schreyer order is not set.");       Error("Schreyer order is not set.");
   }    }
   /* More check will be necessary. */    /* More check will be necessary. */
   return(true);    return(true);
 }  
   
   }
   
   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]):"]]);
   
   def SgetShifts(resmat,w) {
     local i,n,ans,m0;
     n = Length(resmat);
     ans = NewArray(n+1);
     m0 = NewArray(Length(resmat[0,0]));
     ans[0] = m0;
     for (i=0; i<n; 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 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;
     }
     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:"
   ]]);
   
   /* 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);
   }
   
   def BfRoots1(G,V) {
      local bb,ans;
      sm1(" /BFparlist [ ] def ");
      if (IsString(V)) {
         sm1(" [ V to_records pop ] /V set ");
      }else {
        sm1(" V { toString } map /V set ");
      }
      sm1(" /BFvarlist V def ");
   
      sm1(" G flatten { toString } map  /G set ");
      sm1(" G V bfm /bb set ");
      if (IsSm1Integer(bb)) {
        return([ ]);
      }
      sm1(" bb 0 get findIntegralRoots { (universalNumber) dc } map /ans set ");
      return([ans, bb]);
   }
   
   HelpAdd(["BfRoots1",
   ["BfRoots1(g,v) returns the integral roots of g with respect to the weight",
    "vector (1,1,...,1) and the b-function itself",
    "Example:  BfRoots1([x*Dx-2, y*Dy-3],[x,y]);"
   ]]);
   
   
   
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.35
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>