[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.26 and 1.27

version 1.26, 2000/08/10 02:59:08 version 1.27, 2000/08/16 22:38:52
Line 1 
Line 1 
 /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.25 2000/08/02 05:14:31 takayama Exp $ */  /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.26 2000/08/10 02:59:08 takayama Exp $ */
 #define DEBUG 1  #define DEBUG 1
 Sordinary = false;  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),
Line 7  Sordinary = false;
Line 7  Sordinary = false;
 */  */
 #define OFFSET 0  #define OFFSET 0
 /* #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 278  def SresolutionFrameWithTower(g,opt) {
Line 294  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 322  def NewPolynomialVector(size) {
Line 338  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 {
         (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message } { } ifelse
       [(Homogenize) 0] system_variable        [(Homogenize) 0] system_variable
       [(ReduceLowerTerms) 0] system_variable        [(ReduceLowerTerms) 0] system_variable
     } {  } ifelse      } {  } 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 {
           (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message } {  } ifelse
       [(Homogenize) 1] system_variable        [(Homogenize) 1] system_variable
       [(ReduceLowerTerms) 1] system_variable        [(ReduceLowerTerms) 1] system_variable
     } {  } ifelse      } {  } ifelse
Line 379  def Sres0FrameWithSkelton(g) {
Line 401  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 387  def Sres0FrameWithSkelton(g) {
Line 409  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 443  def test_SinitOfArray() {
Line 465  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 470  def SgenerateTable(tower) {
Line 494  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 555  def SlaScala(g,opt) {
Line 579  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 584  def SlaScala(g,opt) {
Line 609  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 614  if (Sordinary) {
Line 640  if (Sordinary) {
 }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;
Line 648  if (Sordinary) {
Line 674  if (Sordinary) {
     }      }
     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 695  def SnextI(reductionTable_tmp,strategy,redundantTable,
Line 722  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 731  def SwhereInGB(f,tower) {
Line 758  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 785  def SwhereInTower(f,tower) {
Line 812  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 797  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 824  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 824  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 851  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 856  def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 
Line 883  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 948  def Sbases_to_vec(bases,size) {
Line 975  def Sbases_to_vec(bases,size) {
   
 HelpAdd(["Sminimal",  HelpAdd(["Sminimal",
 ["It constructs the V-minimal free resolution by LaScala'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)",   "      : \"Sordinary\"   (no (u,v)-minimal resolution)",
  "Options should be given as an array.",   "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]]);",
Line 979  def Sminimal(g,opt) {
Line 1006  def Sminimal(g,opt) {
   redundantTable = r[1];    redundantTable = r[1];
   reducer = r[2];    reducer = r[2];
   bettiTable = SbettiTable(redundantTable);    bettiTable = SbettiTable(redundantTable);
   Println("Betti numbers are ------");    Sprintln2("Betti numbers are ------");
   sm1_pmat(bettiTable);    if (Sverbose || Sverbose2) {sm1_pmat(bettiTable);}
   minRes = SnewArrayOfFormat(freeRes);    minRes = SnewArrayOfFormat(freeRes);
   seq = 0;    seq = 0;
   maxSeq = SgetMaxSeq(redundantTable);    maxSeq = SgetMaxSeq(redundantTable);
Line 990  def Sminimal(g,opt) {
Line 1017  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 1008  def Sminimal(g,opt) {
Line 1035  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 1021  def Sminimal(g,opt) {
Line 1048  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 1036  def Sminimal(g,opt) {
Line 1063  def Sminimal(g,opt) {
    tminRes = Stetris(minRes,redundantTable);     tminRes = Stetris(minRes,redundantTable);
    ansSminimal = [SpruneZeroRow(tminRes), tminRes,     ansSminimal = [SpruneZeroRow(tminRes), tminRes,
                   [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]];                    [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]];
      Sprintln2(" ");
    Println("------------ Note -----------------------------");     Println("------------ Note -----------------------------");
    Println("To get shift vectors, use Reparse and SgetShifts(resmat,w)");     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("To get initial of the complex, use Reparse and Sinit_w(resmat,w)");
Line 1122  def Stetris(freeRes,redundantTable) {
Line 1150  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 1458  def ScheckIfSchreyer(s) {
Line 1488  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. */
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.27
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>