OpenXM/src/k097/lib/minimal/minimal.k - diff

Return to minimal.k CVS log

Up to [local] / OpenXM / src / k097 / lib / minimal

Diff for /OpenXM/src/k097/lib/minimal/minimal.k between version 1.2 and 1.26

-version 1.2, 2000/05/03 07:50:38
+version 1.26, 2000/08/10 02:59:08
 Line 1
 Line 1
 Line 1
- /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.1 2000/05/03 06:42:07 takayama Exp $ */
+ /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.25 2000/08/02 05:14:31 takayama Exp $ */
  #define DEBUG 1
- /* #define ORDINARY 1 */
+ Sordinary = false;
+ /* 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 OFFSET 20*/
  /* Test sequences.
     Use load["minimal.k"];;
-Line 28  def load_tower() {
+Line 34  def load_tower() {
 Line 28  def load_tower() {
 Line 34  def load_tower() {
    if (Boundp("k0-tower.sm1.loaded")) {
    }else{
      sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
+     sm1(" [(parse) (new.sm1) pushfile ] extension ");
      sm1(" /k0-tower.sm1.loaded 1 def ");
    }
+   sm1(" oxNoX ");
  }
  load_tower();
  SonAutoReduce = true;
-Line 42  def Reverse(f) {
+Line 50  def Reverse(f) {
 Line 42  def Reverse(f) {
 Line 50  def Reverse(f) {
  def Sgroebner(f) {
     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."]]);
+ /*  End of standard functions that should be moved to standard libraries. */
  def test0() {
    local f;
    Sweyl("x,y,z");
-Line 124  sm1(" [(AvoidTheSameRing)] pushEnv
+Line 212  sm1(" [(AvoidTheSameRing)] pushEnv
 Line 124  sm1(" [(AvoidTheSameRing)] pushEnv
 Line 212  sm1(" [(AvoidTheSameRing)] pushEnv
        [ [(AvoidTheSameRing) 0] system_variable
          [(gbListTower) tower (list) dc] system_variable
        ] pop popEnv ");
+       /* sm1("(hoge) message show_ring "); */
  }
  def SresolutionFrameWithTower(g,opt) {
    local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
-         gbasis;
+         gbasis, nohomog,i,n;
+   /* extern Sordinary */
+   nohomog = false;
+   count = -1;  Sordinary = false; /* default value for options. */
    if (Length(Arglist) >= 2) {
-     if (IsInteger(opt)) count = opt;
+     if (IsArray(opt)) {
-   }else{
+       n = Length(opt);
-     count = -1;
+       for (i=0; i<n; i++) {
+         if (IsInteger(opt[i])) {
+           count = opt[i];
+         }
+         if (IsString(opt[i])) {
+           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("--------------------------------------------");
+     }
    }
    sm1(" setupEnvForResolution ");
-Line 144  def SresolutionFrameWithTower(g,opt) {
+Line 256  def SresolutionFrameWithTower(g,opt) {
 Line 144  def SresolutionFrameWithTower(g,opt) {
 Line 256  def SresolutionFrameWithTower(g,opt) {
    */
    sm1(" (mmLarger) (matrix) switch_function ");
-   g = Map(g,"Shomogenize");
+   if (! nohomog) {
+     Println("Automatic homogenization.");
+     g = Map(g,"Shomogenize");
+   }else{
+     Println("No automatic homogenization.");
+   }
    if (SonAutoReduce) {
      sm1("[ (AutoReduce) ] system_variable /autof set ");
      sm1("[ (AutoReduce) 1 ] system_variable ");
-Line 184  def SresolutionFrameWithTower(g,opt) {
+Line 301  def SresolutionFrameWithTower(g,opt) {
 Line 184  def SresolutionFrameWithTower(g,opt) {
 Line 301  def SresolutionFrameWithTower(g,opt) {
  }
  HelpAdd(["SresolutionFrameWithTower",
  ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
+  "option: \"homogenized\" (no automatic homogenization) ",
   "Example: Sweyl(\"x,y\");",
   "         a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
  def SresolutionFrame(f,opt) {
    local ans;
-   ans = SresolutionFrameWithTower(f);
+   ans = SresolutionFrameWithTower(f,opt);
    return(ans[0]);
  }
  /* ---------------------------- */
-Line 283  def Sres0FrameWithSkelton(g) {
+Line 401  def Sres0FrameWithSkelton(g) {
 Line 283  def Sres0FrameWithSkelton(g) {
 Line 401  def Sres0FrameWithSkelton(g) {
  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);
  }
+ 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]); */
  def Sord_w(f,w) {
    local neww,i,n;
-Line 332  def test_SinitOfArray() {
+Line 456  def test_SinitOfArray() {
 Line 332  def test_SinitOfArray() {
 Line 456  def test_SinitOfArray() {
  /* f is assumed to be a monomial with toes. */
  def Sdegree(f,tower,level) {
-   local i;
+   local i,ww, wd;
+   /* extern WeightOfSweyl; */
+   ww = WeightOfSweyl;
+   f = Init(f);
    if (level <= 1) return(StotalDegree(f));
    i = Degree(f,es);
    return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
  }
  def SgenerateTable(tower) {
    local height, n,i,j, ans, ans_at_each_floor;
+   /*
+   Print("SgenerateTable: tower=");Println(tower);
+   sm1(" print_switch_status "); */
    height = Length(tower);
    ans = NewArray(height);
    for (i=0; i<height; i++) {
      n = Length(tower[i]);
      ans_at_each_floor=NewArray(n);
      for (j=0; j<n; j++) {
-       ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1);
+       ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1)
+                             + OFFSET;
        /* Println([i,j,ans_at_each_floor[j]]); */
      }
      ans[i] = ans_at_each_floor;
-Line 367  def SnewArrayOfFormat(p) {
+Line 500  def SnewArrayOfFormat(p) {
 Line 367  def SnewArrayOfFormat(p) {
 Line 500  def SnewArrayOfFormat(p) {
       return(null);
    }
  }
+ def ScopyArray(a) {
+   local n, i,ans;
+   n = Length(a);
+   ans = NewArray(n);
+   for (i=0; i<n; i++) {
+     ans[i] = a[i];
+   }
+   return(ans);
+ }
  def SminOfStrategy(a) {
    local n,i,ans,tt;
    ans = 100000; /* very big number */
-Line 405  def SmaxOfStrategy(a) {
+Line 547  def SmaxOfStrategy(a) {
 Line 405  def SmaxOfStrategy(a) {
 Line 547  def SmaxOfStrategy(a) {
  }
- def SlaScala(g) {
+ def SlaScala(g,opt) {
    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;
+         redundantTable_ordinary, redundant_seq_ordinary,
+         reductionTable_tmp;
    /* extern WeightOfSweyl; */
    ww = WeightOfSweyl;
-   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
+   Print("WeightOfSweyl="); Println(WeightOfSweyl);
-   rf = SresolutionFrameWithTower(g);
+   rf = SresolutionFrameWithTower(g,opt);
+   Print("rf="); sm1_pmat(rf);
    redundant_seq = 1;   redundant_seq_ordinary = 1;
    tower = rf[1];
+   Println("Generating reduction table which gives an order of reduction.");
+   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
+   Print("tower"); Println(tower);
    reductionTable = SgenerateTable(tower);
+   Print("reductionTable="); sm1_pmat(reductionTable);
    skel = rf[2];
    redundantTable = SnewArrayOfFormat(rf[1]);
    redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
-Line 430  def SlaScala(g) {
+Line 580  def SlaScala(g) {
 Line 430  def SlaScala(g) {
 Line 580  def SlaScala(g) {
    while (strategy <= maxOfStrategy) {
      for (level = 0; level < height; 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) {
-            Print("Processing "); Print([level,i]);
+            Print("Processing [level,i]= "); Print([level,i]);
             Print("   Strategy = "); Println(strategy);
             if (level == 0) {
               if (IsNull(redundantTable[level,i])) {
-Line 453  def SlaScala(g) {
+Line 608  def SlaScala(g) {
 Line 453  def SlaScala(g) {
 Line 608  def SlaScala(g) {
                    place = f[3];
                    /* (level-1, place) is the place for f[0],
                       which is a newly obtained  GB. */
- #ifdef ORDINARY
+ if (Sordinary) {
                    redundantTable[level-1,place] = redundant_seq;
                    redundant_seq++;
- #else
+ }else{
                    if (f[4] > f[5]) {
                      /* Zero in the gr-module */
                      Print("v-degree of [org,remainder] = ");
-Line 467  def SlaScala(g) {
+Line 622  def SlaScala(g) {
 Line 467  def SlaScala(g) {
 Line 622  def SlaScala(g) {
                      redundantTable[level-1,place] = redundant_seq;
                      redundant_seq++;
                    }
- #endif
+ }
                    redundantTable_ordinary[level-1,place]
                       =redundant_seq_ordinary;
                    redundant_seq_ordinary++;
-Line 500  def SlaScala(g) {
+Line 655  def SlaScala(g) {
 Line 500  def SlaScala(g) {
 Line 655  def SlaScala(g) {
      bases = Sbases_to_vec(bases,bettiTable[i]);
      freeResV[i] = bases;
    }
-   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
+   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary,rf]);
  }
+ def SthereIs(reductionTable_tmp,strategy) {
+   local n,i;
+   n = Length(reductionTable_tmp);
+   for (i=0; i<n; i++) {
+     if (reductionTable_tmp[i] == strategy) {
+       return(true);
+     }
+   }
+   return(false);
+ }
+ def SnextI(reductionTable_tmp,strategy,redundantTable,
+                                   skel,level,freeRes)
+ {
+    local ii,n,p,myindex,i,j,bases;
+    n = Length(reductionTable_tmp);
+    if (level == 0) {
+      for (ii=0; ii<n; ii++) {
+        if (reductionTable_tmp[ii] == strategy) {
+           return(ii);
+         }
+       }
+    }else{
+      for (ii=0; ii<n; ii++) {
+        if (reductionTable_tmp[ii] == strategy) {
+          p = skel[level,ii];
+          myindex = p[0];
+          i = myindex[0]; j = myindex[1];
+          bases = freeRes[level-1];
+          if (IsNull(bases[i]) || IsNull(bases[j])) {
+          }else{
+            return(ii);
+          }
+        }
+      }
+    }
+    Print("reductionTable_tmp=");
+    Println(reductionTable_tmp);
+    Println("See also reductionTable, strategy, level,i");
+    Error("SnextI: bases[i] or bases[j] is null for all combinations.");
+ }
  def SsetBettiTable(freeRes,g) {
    local level,i, n,bases,ans;
    ans = NewArray(Length(freeRes)+1);
-Line 548  def SunitOfFormat(pos,forms) {
+Line 748  def SunitOfFormat(pos,forms) {
 Line 548  def SunitOfFormat(pos,forms) {
 Line 748  def SunitOfFormat(pos,forms) {
    return(ans);
  }
- def Error(s) {
-   sm1(" s error ");
- }
- def IsNull(s) {
-   if (Stag(s) == 0) return(true);
-   else return(false);
- }
  def StowerOf(tower,level) {
    local ans,i;
    ans = [ ];
-Line 577  def Sspolynomial(f,g) {
+Line 769  def Sspolynomial(f,g) {
 Line 577  def Sspolynomial(f,g) {
 Line 769  def Sspolynomial(f,g) {
    sm1("f g spol /FunctionValue set");
  }
- def MonomialPart(f) {
-   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) {
    local i,n,p,q;
    if (f == Poly("0")) return(-1);
-Line 617  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
+Line 812  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 617  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 812  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
    tower2 = StowerOf(tower,level-1);
    SsetTower(tower2);
+   Println(["level=",level]);
+   Println(["tower2=",tower2]);
    /** sm1(" show_ring ");   */
    gi = Stoes_vec(bases[i]);
-Line 638  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
+Line 835  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 638  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 835  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
    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);
    Print("vdegree of the original = "); Println(vdeg);
    Print("vdegree of the remainder = "); Println(vdeg_reduced);
-Line 648  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
+Line 847  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 648  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
 Line 847  def SpairAndReduction(skel,level,ii,freeRes,tower,ww)
    sj = sj*tmp[1]+t_syz[j];
    t_syz[i] = si;
    t_syz[j] = sj;
+   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];
    /* pos is the place to put syzygy at level. */
-Line 686  def Sreduction(f,myset) {
+Line 889  def Sreduction(f,myset) {
 Line 686  def Sreduction(f,myset) {
 Line 889  def Sreduction(f,myset) {
    return([tmp[0],tmp[1],t_syz]);
  }
- def Warning(s) {
-   Print("Warning: ");
-   Println(s);
- }
- def RingOf(f) {
-   local r;
-   if (IsPolynomial(f)) {
-     if (f != Poly("0")) {
-       sm1(f," (ring) dc /r set ");
-     }else{
-       sm1(" [(CurrentRingp)] system_variable /r set ");
-     }
-   }else{
-     Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
-     sm1(" [(CurrentRingp)] system_variable /r set ");
-   }
-   return(r);
- }
  def Sfrom_es(f,size) {
    local c,ans, i, d, myes, myee, j,n,r,ans2;
-Line 761  def Sbases_to_vec(bases,size) {
+Line 946  def Sbases_to_vec(bases,size) {
 Line 761  def Sbases_to_vec(bases,size) {
 Line 946  def Sbases_to_vec(bases,size) {
    return(newbases);
  }
- def Sminimal(g) {
+ HelpAdd(["Sminimal",
+ ["It constructs the V-minimal free resolution by LaScala's algorithm",
+  "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]]);",
+  "          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,opt) {
    local r, freeRes, redundantTable, reducer, maxLevel,
          minRes, seq, maxSeq, level, betti, q, bases, dr,
-         betti_levelplus, newbases, i, j,qq;
+         betti_levelplus, newbases, i, j,qq, tminRes,bettiTable, ansSminimal;
-   r = SlaScala(g);
+   if (Length(Arglist) < 2) {
+      opt = null;
+   }
+   /* Sordinary is set in SlaScala(g,opt) --> SresolutionFrameWithTower */
+   ScheckIfSchreyer("Sminimal:0");
+   r = SlaScala(g,opt);
    /* Should I turn off the tower?? */
+   ScheckIfSchreyer("Sminimal:1");
    freeRes = r[0];
    redundantTable = r[1];
    reducer = r[2];
+   bettiTable = SbettiTable(redundantTable);
+   Println("Betti numbers are ------");
+   sm1_pmat(bettiTable);
    minRes = SnewArrayOfFormat(freeRes);
    seq = 0;
    maxSeq = SgetMaxSeq(redundantTable);
-Line 822  def Sminimal(g) {
+Line 1033  def Sminimal(g) {
 Line 822  def Sminimal(g) {
 Line 1033  def Sminimal(g) {
        }
      }
     }
-    return([Stetris(minRes,redundantTable),
+    tminRes = Stetris(minRes,redundantTable);
-           [ minRes, redundantTable, reducer,r[3],r[4]]]);
+    ansSminimal = [SpruneZeroRow(tminRes), tminRes,
+                   [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]];
+    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(Map(ansSminimal[0],"Length"));
+    Print("Betti numbers of the Schreyer frame: ");
+    Println(Map(ansSminimal[3],"Length"));
+    Println("-----------------------------------------------");
+    sm1(" restoreEnvAfterResolution ");
+    Sordinary = false;
+    return(ansSminimal);
    /* r[4] is the redundantTable_ordinary */
+   /* r[0] is the freeResolution */
+   /* r[5] is the skelton */
  }
-Line 933  In(20)=SvDegree(x,tt,2,ww):
+Line 1162  In(20)=SvDegree(x,tt,2,ww):
 Line 933  In(20)=SvDegree(x,tt,2,ww):
 Line 1162  In(20)=SvDegree(x,tt,2,ww):
  def SvDegree(f,tower,level,w) {
    local i,ans;
    if (IsZero(f)) return(null);
+   f = Init(f);
    if (level <= 0) {
      return(Sord_w(f,w));
    }
-Line 955  def Sannfs(f,v) {
+Line 1185  def Sannfs(f,v) {
 Line 955  def Sannfs(f,v) {
 Line 1185  def Sannfs(f,v) {
  def Sannfs2(f) {
    local p,pp;
    p = Sannfs(f,"x,y");
-   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
+   sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
-   pp = Map(p[0],"Spoly");
+   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
-   return(Sminimal(pp));
+   pp = Map(p,"Spoly");
+   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, 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]:"
+ ]]);
+ /* 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");    : it causes an error. It should be fixed.
+   a=Sannfs2("x^3-y^2-x");
+   a=Sannfs2("x*y*(x-y)");
  */
+ 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(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, Sannfs2.",
+  "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
+  "         b=a[0]; sm1_pmat(b);",
+  "         b[1]*b[0]: b[2]*b[1]:"]]);
+ /* 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 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 IsSameIdeal_h(ii,jj,v) {
+   local a;
+   v = ToString_array(v);
+   a = [ii,jj,v];
+   sm1(a," isSameIdeal_h /FunctionValue set ");
+ }
+ HelpAdd(["IsSameIdeal_h",
+ ["IsSameIdeal_h(ii,jj,var): bool",
+  "It checks the given ideals are the same 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"
+ ]]);
+ def ScheckIfSchreyer(s) {
+   local ss;
+   sm1(" (report) (grade) switch_function /ss set ");
+   if (ss != "module1v") {
+      Print("ScheckIfSchreyer: from "); Println(s);
+      Error("grade is not module1v");
+   }
+   /*
+   sm1(" (report) (mmLarger) switch_function /ss set ");
+   if (ss != "tower") {
+      Print("ScheckIfSchreyer: from "); Println(s);
+      Error("mmLarger is not tower");
+   }
+   */
+   sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set ");
+   if (ss != 1) {
+      Print("ScheckIfSchreyer: from "); Println(s);
+      Error("Schreyer order is not set.");
+   }
+   /* More check will be necessary. */
+   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);
+   m0 = NewArray(Length(resmat[0,0]));
+   ans[0] = m0;
+   for (i=0; i<n-1; 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);
+ }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>