[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.5 and 1.14

version 1.5, 2000/05/05 08:13:49 version 1.14, 2000/06/09 08:04:54
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.4 2000/05/04 11:05:20 takayama Exp $ */  /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.13 2000/06/08 08:37:53 takayama Exp $ */
 #define DEBUG 1  #define DEBUG 1
 /* #define ORDINARY 1 */  /* #define ORDINARY 1 */
 /* 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 TOTAL_STRATEGY 1
   /* #define OFFSET 20*/
 /* Test sequences.  /* Test sequences.
    Use load["minimal.k"];;     Use load["minimal.k"];;
   
Line 34  def load_tower() {
Line 37  def load_tower() {
     sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");      sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
     sm1(" /k0-tower.sm1.loaded 1 def ");      sm1(" /k0-tower.sm1.loaded 1 def ");
   }    }
     sm1(" oxNoX ");
 }  }
 load_tower();  load_tower();
 SonAutoReduce = true;  SonAutoReduce = true;
Line 128  sm1(" [(AvoidTheSameRing)] pushEnv 
Line 132  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 "); */
 }  }
   
 def SresolutionFrameWithTower(g,opt) {  def SresolutionFrameWithTower(g,opt) {
Line 287  def Sres0FrameWithSkelton(g) {
Line 292  def Sres0FrameWithSkelton(g) {
   
   
 def StotalDegree(f) {  def StotalDegree(f) {
   sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");    local d0;
     sm1(" [(grade) f] gbext (universalNumber) dc /d0 set ");
     /* Print("degree of "); Print(f); Print(" is "); Println(d0); */
     return(d0);
 }  }
   
 /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */  /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
Line 336  def test_SinitOfArray() {
Line 344  def test_SinitOfArray() {
   
 /* f is assumed to be a monomial with toes. */  /* f is assumed to be a monomial with toes. */
 def Sdegree(f,tower,level) {  def Sdegree(f,tower,level) {
   local i;    local i,ww, wd;
     /* extern WeightOfSweyl; */
     ww = WeightOfSweyl;
   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);
   return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));  #ifdef TOTAL_STRATEGY
     return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
   #endif
     /* Strategy must be compatible with ordering.  */
     /* Weight vector must be non-negative, too.  */
     /* See Sdegree, SgenerateTable, reductionTable. */
     wd = Sord_w(f,ww);
     return(wd+Sdegree(tower[level-2,i],tower,level-1));
   
 }  }
   
 def SgenerateTable(tower) {  def SgenerateTable(tower) {
Line 351  def SgenerateTable(tower) {
Line 369  def SgenerateTable(tower) {
     n = Length(tower[i]);      n = Length(tower[i]);
     ans_at_each_floor=NewArray(n);      ans_at_each_floor=NewArray(n);
     for (j=0; j<n; j++) {      for (j=0; j<n; j++) {
       ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1);        ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1)
                               + OFFSET;
       /* Println([i,j,ans_at_each_floor[j]]); */        /* Println([i,j,ans_at_each_floor[j]]); */
     }      }
     ans[i] = ans_at_each_floor;      ans[i] = ans_at_each_floor;
Line 427  def SlaScala(g) {
Line 446  def SlaScala(g) {
         reductionTable_tmp;          reductionTable_tmp;
   /* extern WeightOfSweyl; */    /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;    ww = WeightOfSweyl;
   Print("WeghtOfSweyl="); Println(WeightOfSweyl);    Print("WeightOfSweyl="); Println(WeightOfSweyl);
   rf = SresolutionFrameWithTower(g);    rf = SresolutionFrameWithTower(g);
     Print("rf="); sm1_pmat(rf);
   redundant_seq = 1;   redundant_seq_ordinary = 1;    redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];    tower = rf[1];
   reductionTable = SgenerateTable(tower);    reductionTable = SgenerateTable(tower);
Line 646  def MonomialPart(f) {
Line 666  def MonomialPart(f) {
   sm1(" [(lmonom) f] gbext /FunctionValue set ");    sm1(" [(lmonom) f] gbext /FunctionValue set ");
 }  }
   
   /* WARNING:
     When you use SwhereInTower, you have to change gbList
     as below. Ofcourse, you should restrore the gbList
     SsetTower(StowerOf(tower,level));
     pos = SwhereInTower(syzHead,tower[level]);
   */
 def SwhereInTower(f,tower) {  def SwhereInTower(f,tower) {
   local i,n,p,q;    local i,n,p,q;
   if (f == Poly("0")) return(-1);    if (f == Poly("0")) return(-1);
Line 682  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 708  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   
   tower2 = StowerOf(tower,level-1);    tower2 = StowerOf(tower,level-1);
   SsetTower(tower2);    SsetTower(tower2);
     Println(["level=",level]);
     Println(["tower2=",tower2]);
   /** sm1(" show_ring ");   */    /** sm1(" show_ring ");   */
   
   gi = Stoes_vec(bases[i]);    gi = Stoes_vec(bases[i]);
Line 715  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 743  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   sj = sj*tmp[1]+t_syz[j];    sj = sj*tmp[1]+t_syz[j];
   t_syz[i] = si;    t_syz[i] = si;
   t_syz[j] = sj;    t_syz[j] = sj;
   
     SsetTower(StowerOf(tower,level));
   pos = SwhereInTower(syzHead,tower[level]);    pos = SwhereInTower(syzHead,tower[level]);
   
     SsetTower(StowerOf(tower,level-1));
   pos2 = SwhereInTower(tmp[0],tower[level-1]);    pos2 = SwhereInTower(tmp[0],tower[level-1]);
   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. */
Line 828  def Sbases_to_vec(bases,size) {
Line 860  def Sbases_to_vec(bases,size) {
   return(newbases);    return(newbases);
 }  }
   
   HelpAdd(["Sminimal",
   ["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm",
    "Example:  Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
    "          v=[[2*x*Dx + 3*y*Dy+6, 0],",
    "             [3*x^2*Dy + 2*y*Dx, 0],",
    "             [0,  x^2+y^2],",
    "             [0,  x*y]];",
    "         a=Sminimal(v);",
    "         Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
    "         b = ReParse(a[0]); sm1_pmat(b); ",
    "         IsExact_h(b,[x,y]):",
    "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
   
 def Sminimal(g) {  def Sminimal(g) {
   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;          betti_levelplus, newbases, i, j,qq, tminRes;
   r = SlaScala(g);    r = SlaScala(g);
   /* Should I turn off the tower?? */    /* Should I turn off the tower?? */
   freeRes = r[0];    freeRes = r[0];
Line 889  def Sminimal(g) {
Line 934  def Sminimal(g) {
       }        }
     }      }
    }     }
    return([Stetris(minRes,redundantTable),     tminRes = Stetris(minRes,redundantTable);
      return([SpruneZeroRow(tminRes), tminRes,
           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);            [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
   /* r[4] is the redundantTable_ordinary */    /* r[4] is the redundantTable_ordinary */
   /* r[0] is the freeResolution */    /* r[0] is the freeResolution */
Line 1024  def Sannfs(f,v) {
Line 1070  def Sannfs(f,v) {
 def Sannfs2(f) {  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 ");
 /*  /*
   Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],    Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
                ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */                 ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */
   Sweyl("x,y",[["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");    pp = Map(p[0],"Spoly");
   return(Sminimal(pp));    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 ");
   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[0],"Spoly");    pp = Map(p,"Spoly");
   return(Sminimal(pp));    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).    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*y*(x+y-1)"); ==> The betti numbers are 3, 2.
Line 1048  def Sannfs3(f) {
Line 1143  def Sannfs3(f) {
   
 */  */
   
   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 is under construction. */  /*  The below does not use LaScala-Stillman's algorithm. */
 def Sschreyer(g) {  def Sschreyer(g) {
   local rf, tower, reductionTable, skel, redundantTable, bases,    local rf, tower, reductionTable, skel, redundantTable, bases,
         strategy, maxOfStrategy, height, level, n, i,          strategy, maxOfStrategy, height, level, n, i,
         freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,          freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
         redundantTable_ordinary, redundant_seq_ordinary,          redundantTable_ordinary, redundant_seq_ordinary,
         reductionTable_tmp,c2,ii,nn;          reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
   /* extern WeightOfSweyl; */    /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;    ww = WeightOfSweyl;
   Print("WeghtOfSweyl="); Println(WeightOfSweyl);    Print("WeghtOfSweyl="); Println(WeightOfSweyl);
Line 1121  def Sschreyer(g) {
Line 1225  def Sschreyer(g) {
                   /* i must be equal to f[2], I think. Double check. */                    /* i must be equal to f[2], I think. Double check. */
   
                   /* Correction Of Constant */                    /* Correction Of Constant */
                   c2 = f[6];                    /* Correction of syzygy */
                     c2 = f[6];  /* or -f[6]?  Double check. */
                     Print("c2="); Println(c2);
                   nn = Length(bases);                    nn = Length(bases);
                   for (ii=0; ii<nn;ii++) {                    for (ii=0; ii<nn;ii++) {
                      if (ii != place) {                       if ((ii != i) && (! IsNull(bases[ii]))) {
                        bases[ii] = bases[ii]*c2;                         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;                    freeRes[level] = bases;
                   /* bases = freeRes[level-1];                    Print("New freeRes[level] = "); sm1_pmat(freeRes[level]);
                      bases[place] = f[0];  
                      freeRes[level-1] = bases;  It is already set. */                   /* Update the freeRes[level-1] */
                   reducer[level-1,place] = f[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{                 }else{
                   /* redundantTable[level,i] = 0; */                    /* redundantTable[level,i] = 0; */
                   bases = freeRes[level];                    bases = freeRes[level];
Line 1143  def Sschreyer(g) {
Line 1277  def Sschreyer(g) {
              }  /* end of level >= 1 */               }  /* end of level >= 1 */
           }            }
     } /* i loop */      } /* 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 */    } /* level loop */
   n = Length(freeRes);    n = Length(freeRes);
   freeResV = SnewArrayOfFormat(freeRes);    freeResV = SnewArrayOfFormat(freeRes);
Line 1151  def Sschreyer(g) {
Line 1312  def Sschreyer(g) {
     bases = Sbases_to_vec(bases,bettiTable[i]);      bases = Sbases_to_vec(bases,bettiTable[i]);
     freeResV[i] = bases;      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]);    return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
 }  }
   
Line 1158  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
Line 1331  def SpairAndReduction2(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,n,c2;         vdeg,vdeg_reduced,n,c2;
   Println("SpairAndReduction2:");    Println("SpairAndReduction2 : -------------------------");
   
   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];
Line 1173  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
Line 1346  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
   
   tower2 = StowerOf(tower,level-1);    tower2 = StowerOf(tower,level-1);
   SsetTower(tower2);    SsetTower(tower2);
     Println(["level=",level]);
     Println(["tower2=",tower2]);
   /** sm1(" show_ring ");   */    /** sm1(" show_ring ");   */
   
   gi = Stoes_vec(bases[i]);    gi = Stoes_vec(bases[i]);
Line 1193  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
Line 1368  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
   tmp = Sreduction(si*gi+sj*gj, bases);    tmp = Sreduction(si*gi+sj*gj, bases);
   
   Print("result is "); Println(tmp);    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];    t_syz = tmp[2];
   si = si*tmp[1]+t_syz[i];    si = si*tmp[1]+t_syz[i];
   sj = sj*tmp[1]+t_syz[j];    sj = sj*tmp[1]+t_syz[j];
Line 1203  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
Line 1383  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
   /* tmp[0] must be zero */    /* tmp[0] must be zero */
   n = Length(t_syz);    n = Length(t_syz);
   for (i=0; i<n; i++) {    for (i=0; i<n; i++) {
      if (IsConstant(t_syz[i])) {       if (IsConstant(t_syz[i])){
         if (!IsZero(t_syz[i])) {
        if (IsNull(redundantTable[level-1,i])) {         if (IsNull(redundantTable[level-1,i])) {
          /* i must equal to pos2 below. */           /* i must equal to pos2 below. */
          c2 = -t_syz[i];           c2 = -t_syz[i];
          tmp[0] = freeRes[level-1,i];           tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
          t_syz[i] = 0;           t_syz[i] = 0;
            /* tmp[0] = t_syz . g */
          /* break; does not work. Use */           /* break; does not work. Use */
          i = n;           i = n;
        }         }
         }
      }       }
   }    }
   
Line 1222  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
Line 1405  def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,
   Print("vdegree of the original = "); Println(vdeg);    Print("vdegree of the original = "); Println(vdeg);
   Print("vdegree of the remainder = "); Println(vdeg_reduced);    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]);    pos = SwhereInTower(syzHead,tower[level]);
   
     SsetTower(StowerOf(tower,level-1));
   pos2 = SwhereInTower(tmp[0],tower[level-1]);    pos2 = SwhereInTower(tmp[0],tower[level-1]);
   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];    ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
   /* 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);    Println(ans);
     Println("--- end of SpairAndReduction2  ");
   return(ans);    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 */
   /* 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 SisComplex(a) {
     local n,i,j,k,b,p,q;
     n = Length(a);
     for (i=0; i<n-1; i++) {
       if (Length(a[i+1]) != 0) {
         b = a[i+1]*a[i];
         p = Length(b); q = Length(b[0]);
         for (j=0; j<p; j++) {
           for (k=0; k<q; k++) {
             if (!IsZero(b[j,k])) {
                Print("Is is not complex at ");
                Println([i,j,k]);
                return(false);
             }
           }
         }
       }
     }
     return(true);
   }
   
   def IsExact_h(c,v) {
     local a;
     v = ToString_array(v);
     a = [c,v];
     sm1(a," isExact_h /FunctionValue set ");
   }
   HelpAdd(["IsExact_h",
   ["IsExact_h(complex,var): bool",
    "It checks the given complex is exact or not in D<h> (homogenized Weyl algebra)",
    "cf. ReParse"
   ]]);
   
   def ReParse(a) {
     local c;
     if (IsArray(a)) {
       c = Map(a,"ReParse");
     }else{
       sm1(a," toString . /c set");
     }
     return(c);
   }
   HelpAdd(["ReParse",
   ["Reparse(obj): obj",
    "It parses the given object in the current ring.",
    "Outputs from SlaScala, Sschreyer may cause a trouble in other functions,",
    "because it uses the Schreyer order.",
    "In this case, ReParse the outputs from these functions.",
    "cf. IsExaxt_h"
   ]]);
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.14
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>