[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.1 and 1.6

version 1.1, 2000/05/03 06:42:07 version 1.6, 2000/05/06 07:58:37
Line 1 
Line 1 
 /* $OpenXM$ */  /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.5 2000/05/05 08:13:49 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),
      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.  /* Test sequences.
    Use load["minimal.k"];;     Use load["minimal.k"];;
   
Line 332  def test_SinitOfArray() {
Line 339  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);
   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 346  def SgenerateTable(tower) {
Line 364  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 367  def SnewArrayOfFormat(p) {
Line 386  def SnewArrayOfFormat(p) {
      return(null);       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) {  def SminOfStrategy(a) {
   local n,i,ans,tt;    local n,i,ans,tt;
   ans = 100000; /* very big number */    ans = 100000; /* very big number */
Line 409  def SlaScala(g) {
Line 437  def SlaScala(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;
   /* extern WeightOfSweyl; */    /* extern WeightOfSweyl; */
   ww = WeightOfSweyl;    ww = WeightOfSweyl;
   Print("WeghtOfSweyl="); Println(WeightOfSweyl);    Print("WeightOfSweyl="); Println(WeightOfSweyl);
   rf = SresolutionFrameWithTower(g);    rf = SresolutionFrameWithTower(g);
   redundant_seq = 1;   redundant_seq_ordinary = 1;    redundant_seq = 1;   redundant_seq_ordinary = 1;
   tower = rf[1];    tower = rf[1];
Line 430  def SlaScala(g) {
Line 459  def SlaScala(g) {
   while (strategy <= maxOfStrategy) {    while (strategy <= maxOfStrategy) {
     for (level = 0; level < height; level++) {      for (level = 0; level < height; level++) {
       n = Length(reductionTable[level]);        n = Length(reductionTable[level]);
       for (i=0; i<n; i++) {        reductionTable_tmp = ScopyArray(reductionTable[level]);
         while (SthereIs(reductionTable_tmp,strategy)) {
           i = SnextI(reductionTable_tmp,strategy,redundantTable,
                      skel,level,freeRes);
           Println([level,i]);
           reductionTable_tmp[i] = -200000;
         if (reductionTable[level,i] == strategy) {          if (reductionTable[level,i] == strategy) {
            Print("Processing "); Print([level,i]);             Print("Processing "); Print([level,i]);
            Print("   Strategy = "); Println(strategy);             Print("   Strategy = "); Println(strategy);
Line 503  def SlaScala(g) {
Line 537  def SlaScala(g) {
   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);    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) {  def SsetBettiTable(freeRes,g) {
   local level,i, n,bases,ans;    local level,i, n,bases,ans;
   ans = NewArray(Length(freeRes)+1);    ans = NewArray(Length(freeRes)+1);
Line 638  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 717  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
   
   Print("result is "); Println(tmp);    Print("result is "); Println(tmp);
   
   vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww);    /* 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);    vdeg_reduced = SvDegree(tmp[0],tower,level-1,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);
Line 823  def Sminimal(g) {
Line 904  def Sminimal(g) {
     }      }
    }     }
    return([Stetris(minRes,redundantTable),     return([Stetris(minRes,redundantTable),
           [ minRes, redundantTable, reducer,r[3],r[4]]]);            [ 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 */
 }  }
   
   
Line 933  In(20)=SvDegree(x,tt,2,ww):
Line 1015  In(20)=SvDegree(x,tt,2,ww):
 def SvDegree(f,tower,level,w) {  def SvDegree(f,tower,level,w) {
   local i,ans;    local i,ans;
   if (IsZero(f)) return(null);    if (IsZero(f)) return(null);
     f = Init(f);
   if (level <= 0) {    if (level <= 0) {
     return(Sord_w(f,w));      return(Sord_w(f,w));
   }    }
Line 942  def SvDegree(f,tower,level,w) {
Line 1025  def SvDegree(f,tower,level,w) {
   return(ans);    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)); */
   }
   
   /* 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 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));
   }
   
   /*
     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.
   
   */
   
   
   
   /*  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 */
                     c2 = -f[6];  /* or f[6]?  Double check. */
                     nn = Length(bases);
                     for (ii=0; ii<nn;ii++) {
                        if ((ii != place) && (! IsNull(bases[ii]))) {
                          bases[ii] = bases[ii]*c2;
                        }
                     }
   
                     freeRes[level] = bases;
   
                    /* Update the freeRes[level-1] */
                     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 */
             }
       } /* i loop */
   
       /* Triangulate reducer */
       if (level >= 1) {
         Println(" ");
         Print("Triangulating reducer at level "); Println(level-1);
         reducerBase = reducer[level-1];
         Print("reducerBase=");  Println(reducerBase);
         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[jj] = reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
                 }
                }
              }
            }
          }
          Println("New reducer");
          sm1_pmat(reducerBase);
          reducer[level-1] = reducerBase;
       }
   
     } /* level loop */
     n = Length(freeRes);
     freeResV = SnewArrayOfFormat(freeRes);
     for (i=0; i<n; i++) {
       bases = freeRes[i];
       bases = Sbases_to_vec(bases,bettiTable[i]);
       freeResV[i] = bases;
     }
   
     /* Mark the non-redundant elements. */
     for (i=0; i<n; i++) {
       m = Length(redundantTable[i]);
       for (jj=0; jj<m; jj++) {
         if (IsNull(redundantTable[i,jj])) {
           redundantTable[i,jj] = 0;
         }
       }
     }
   
   
     return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
   }
   
   def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
     local i, j, myindex, p, bases, tower2, gi, gj,
          si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
          vdeg,vdeg_reduced,n,c2;
     Println("SpairAndReduction2 : -------------------------");
   
     if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
     p = skel[level,ii];
     myindex = p[0];
     i = myindex[0]; j = myindex[1];
     bases = freeRes[level-1];
     Println(["p and bases ",p,bases]);
     if (IsNull(bases[i]) || IsNull(bases[j])) {
       Println([level,i,j,bases[i],bases[j]]);
       Error("level, i, j : bases[i], bases[j]  must not be NULL.");
     }
   
     tower2 = StowerOf(tower,level-1);
     SsetTower(tower2);
     /** sm1(" show_ring ");   */
   
     gi = Stoes_vec(bases[i]);
     gj = Stoes_vec(bases[j]);
   
     ssp = Sspolynomial(gi,gj);
     si = ssp[0,0];
     sj = ssp[0,1];
     syzHead = si*es^i;
     /* This will be the head term, I think. But, double check. */
     Println([si*es^i,sj*es^j]);
   
     Print("[gi, gj] = "); Println([gi,gj]);
     sm1(" [(Homogenize)] system_variable message ");
     Print("Reduce the element "); Println(si*gi+sj*gj);
     Print("by  "); Println(bases);
   
     tmp = Sreduction(si*gi+sj*gj, bases);
   
     Print("result is "); Println(tmp);
     if (!IsZero(tmp[0])) {
       Print("Error: base = ");
       Println(Map(bases,"Stoes_vec"));
       Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
     }
     t_syz = tmp[2];
     si = si*tmp[1]+t_syz[i];
     sj = sj*tmp[1]+t_syz[j];
     t_syz[i] = si;
     t_syz[j] = sj;
   
     c2 = null;
     /* tmp[0] must be zero */
     n = Length(t_syz);
     for (i=0; i<n; i++) {
        if (IsConstant(t_syz[i])){
         if (!IsZero(t_syz[i])) {
          if (IsNull(redundantTable[level-1,i])) {
            /* i must equal to pos2 below. */
            c2 = -t_syz[i];
            tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
            t_syz[i] = 0;
            /* tmp[0] = t_syz . g */
            /* break; does not work. Use */
            i = n;
          }
         }
        }
     }
   
     /* This is essential part for V-minimal resolution. */
     /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
     vdeg = SvDegree(si*gi,tower,level-1,ww);
     vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
     Print("vdegree of the original = "); Println(vdeg);
     Print("vdegree of the remainder = "); Println(vdeg_reduced);
   
     pos = SwhereInTower(syzHead,tower[level]);
     pos2 = SwhereInTower(tmp[0],tower[level-1]);
     ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
     /* pos is the place to put syzygy at level. */
     /* pos2 is the place to put a new GB at level-1. */
     Println(ans);
     Println("  ");
     return(ans);
   }
   
   def Sminimal_v(g) {
     local r, freeRes, redundantTable, reducer, maxLevel,
           minRes, seq, maxSeq, level, betti, q, bases, dr,
           betti_levelplus, newbases, i, j,qq;
     r = Sschreyer(g);
     sm1_pmat(r);
     Debug_Sminimal_v = r;
     Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
     /* Should I turn off the tower?? */
     freeRes = r[0];
     redundantTable = r[1];
     reducer = r[2];
     minRes = SnewArrayOfFormat(freeRes);
     seq = 0;
     maxSeq = SgetMaxSeq(redundantTable);
     maxLevel = Length(freeRes);
     for (level = 0; level < maxLevel; level++) {
       minRes[level] = freeRes[level];
     }
     for (level = 0; level < maxLevel; level++) {
         betti = Length(freeRes[level]);
         for (q = betti-1; q>=0; q--) {
           if (redundantTable[level,q] > 0) {
             Print("[seq,level,q]="); Println([seq,level,q]);
             if (level < maxLevel-1) {
               bases = freeRes[level+1];
               dr = reducer[level,q];
               dr[q] = -1;
               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
   /*  Do it later.
               for (qq=0; qq<betti; qq++) {
                   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 */
   }
   
   /* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>