=================================================================== RCS file: /home/cvs/OpenXM/src/k097/lib/minimal/minimal.k,v retrieving revision 1.13 retrieving revision 1.35 diff -u -p -r1.13 -r1.35 --- OpenXM/src/k097/lib/minimal/minimal.k 2000/06/08 08:37:53 1.13 +++ OpenXM/src/k097/lib/minimal/minimal.k 2007/07/03 22:05:46 1.35 @@ -1,13 +1,28 @@ -/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.12 2000/05/24 15:24:54 takayama Exp $ */ +/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.34 2001/01/05 11:14:28 takayama Exp $ */ #define DEBUG 1 -/* #define 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 TOTAL_STRATEGY */ /* #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. Use load["minimal.k"];; @@ -29,17 +44,12 @@ */ +/* We cannot use load command in the if statement. */ +load("lib/minimal/cohom.k"); +Load_sm1(["k0-tower.sm1","lib/minimal/k0-tower.sm1"],"k0-tower.sm1.loaded"); +Load_sm1(["new.sm1","lib/minimal/new.sm1"],"new.sm1.loaded"); +sm1(" oxNoX "); -load("cohom.k"); -def load_tower() { - if (Boundp("k0-tower.sm1.loaded")) { - }else{ - sm1(" [(parse) (k0-tower.sm1) pushfile ] extension "); - sm1(" /k0-tower.sm1.loaded 1 def "); - } - sm1(" oxNoX "); -} -load_tower(); SonAutoReduce = true; def Factor(f) { sm1(f, " fctr /FunctionValue set"); @@ -50,6 +60,115 @@ 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."]]); + +def Kernel(f,v) { + local ans; + /* v : string or ring */ + if (Length(Arglist) < 2) { + sm1(" [f] syz /ans set "); + }else{ + sm1(" [f v] syz /ans set "); + } + return(ans); +} +def Syz(f) { + sm1(" [f] syz /FunctionValue set "); +} +HelpAdd(["Kernel", +["Kernel(f) returns the syzygy of f.", + "Return value [b, c]: b is a set of generators of the syzygies of f", + " : c=[gb, backward transformation, syzygy without", + " dehomogenization", + "Example: Weyl(\"x,y\",[[\"x\",-1,\"Dx\",1]]); ", + " s=Kernel([x*Dx+1,Dx^2+x^5]); s[0]:"]]); +/* cf. sm1_syz in cohom.k */ +def Gb(f) { + sm1(" [f] gb /FunctionValue set "); +} +HelpAdd(["Gb", +["Gb(f) returns the Groebner basis of f.", + "cf. Kernel, Weyl."]]); + + +/* End of standard functions that should be moved to standard libraries. */ def test0() { local f; Sweyl("x,y,z"); @@ -69,7 +188,6 @@ def test1() { } - def Sweyl(v,w) { /* extern WeightOfSweyl ; */ local ww,i,n; @@ -128,19 +246,43 @@ def StoTower() { } def SsetTower(tower) { -sm1(" [(AvoidTheSameRing)] pushEnv - [ [(AvoidTheSameRing) 0] system_variable - [(gbListTower) tower (list) dc] system_variable +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; - }else{ - count = -1; + if (IsArray(opt)) { + n = Length(opt); + 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 "); @@ -152,7 +294,12 @@ 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 "); @@ -169,7 +316,7 @@ def SresolutionFrameWithTower(g,opt) { /* -sugar is fine? */ sm1(" setupEnvForResolution "); - Println(g); + Sprintln(g); startingGB = g; /* ans = [ SzeroMap(g) ]; It has not been implemented. see resol1.withZeroMap */ ans = [ ]; @@ -192,12 +339,13 @@ 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]); } /* ---------------------------- */ @@ -210,21 +358,27 @@ def NewPolynomialVector(size) { } def SturnOffHomogenization() { - sm1(" - [(Homogenize)] system_variable 1 eq - { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message - [(Homogenize) 0] system_variable - [(ReduceLowerTerms) 0] system_variable - } { } ifelse + sm1(" \ + [(Homogenize)] system_variable 1 eq \ + { Sverbose { \ + (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message } { } ifelse \ + [(Homogenize) 0] system_variable \ + [(ReduceLowerTerms) 0] system_variable \ + } { } ifelse \ "); } +/* NOTE!!! Be careful these changes of global environmental variables. + We should make a standard set of values and restore these values + after computation and interruption. August 15, 2000. +*/ def SturnOnHomogenization() { - sm1(" - [(Homogenize)] system_variable 0 eq - { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message - [(Homogenize) 1] system_variable - [(ReduceLowerTerms) 1] system_variable - } { } ifelse + sm1(" \ + [(Homogenize)] system_variable 0 eq \ + { Sverbose { \ + (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message } { } ifelse \ + [(Homogenize) 1] system_variable \ + [(ReduceLowerTerms) 1] system_variable \ + } { } ifelse \ "); } @@ -269,7 +423,7 @@ def Sres0FrameWithSkelton(g) { si = pair[1,0]; sj = pair[1,1]; /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */ - Print("."); + Sprint("."); t_syz = NewPolynomialVector(gLength); t_syz[i] = si; @@ -277,7 +431,7 @@ def Sres0FrameWithSkelton(g) { syzAll[k] = t_syz; } t_syz = syzAll; - Print("Done. betti="); Println(betti); + Sprint("Done. betti="); Sprintln(betti); /* 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 , ...] @@ -291,9 +445,15 @@ 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; @@ -327,8 +487,10 @@ 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, -y^2*z^2 + x*z^3 + y*z^3, -z^4]; p=SresolutionFrameWithTower(f); - sm1_pmat(p); - sm1_pmat(SgenerateTable(p[1])); + if (Sverbose) { + sm1_pmat(p); + sm1_pmat(SgenerateTable(p[1])); + } return(p); frame = p[0]; sm1_pmat(p[1]); @@ -346,19 +508,16 @@ def Sdegree(f,tower,level) { f = Init(f); if (level <= 1) return(StotalDegree(f)); i = Degree(f,es); -#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) { local height, n,i,j, ans, ans_at_each_floor; + + /* + Sprint("SgenerateTable: tower=");Sprintln(tower); + sm1(" print_switch_status "); */ height = Length(tower); ans = NewArray(height); for (i=0; i<height; i++) { @@ -434,7 +593,7 @@ 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, @@ -442,11 +601,19 @@ def SlaScala(g) { reductionTable_tmp; /* extern WeightOfSweyl; */ ww = WeightOfSweyl; - Print("WeightOfSweyl="); Println(WeightOfSweyl); - rf = SresolutionFrameWithTower(g); + Sprint("WeightOfSweyl="); Sprintln(WeightOfSweyl); + rf = SresolutionFrameWithTower(g,opt); + Sprint("rf="); if (Sverbose) {sm1_pmat(rf);} redundant_seq = 1; redundant_seq_ordinary = 1; tower = rf[1]; + + Sprintln("Generating reduction table which gives an order of reduction."); + Sprint("WeghtOfSweyl="); Sprintln(WeightOfSweyl); + Sprint2("tower="); Sprintln2(tower); reductionTable = SgenerateTable(tower); + Sprint2("reductionTable="); + if (Sverbose || Sverbose2) {sm1_pmat(reductionTable);} + skel = rf[2]; redundantTable = SnewArrayOfFormat(rf[1]); redundantTable_ordinary = SnewArrayOfFormat(rf[1]); @@ -464,11 +631,12 @@ def SlaScala(g) { while (SthereIs(reductionTable_tmp,strategy)) { i = SnextI(reductionTable_tmp,strategy,redundantTable, skel,level,freeRes); - Println([level,i]); + Sprintln([level,i]); reductionTable_tmp[i] = -200000; if (reductionTable[level,i] == strategy) { - Print("Processing "); Print([level,i]); - Print(" Strategy = "); Println(strategy); + Sprint("Processing [level,i]= "); Sprint([level,i]); + Sprint(" Strategy = "); Sprintln(strategy); + Sprint2(strategy); if (level == 0) { if (IsNull(redundantTable[level,i])) { bases = freeRes[level]; @@ -488,21 +656,21 @@ 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] = "); - Println([f[4],f[5]]); - Print("[level,i] = "); Println([level,i]); + Sprint("v-degree of [org,remainder] = "); + Sprintln([f[4],f[5]]); + Sprint("[level,i] = "); Sprintln([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++; @@ -528,6 +696,7 @@ def SlaScala(g) { } strategy++; } + Sprintln2(" "); n = Length(freeRes); freeResV = SnewArrayOfFormat(freeRes); for (i=0; i<n; i++) { @@ -535,7 +704,7 @@ 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) { @@ -575,9 +744,9 @@ def SnextI(reductionTable_tmp,strategy,redundantTable, } } } - Print("reductionTable_tmp="); - Println(reductionTable_tmp); - Println("See also reductionTable, strategy, level,i"); + Sprint("reductionTable_tmp="); + Sprintln(reductionTable_tmp); + Sprintln("See also reductionTable, strategy, level,i"); Error("SnextI: bases[i] or bases[j] is null for all combinations."); } @@ -611,7 +780,7 @@ def SwhereInGB(f,tower) { q = MonomialPart(f); if (p == q) return(i); } - Println([f,tower]); + Sprintln([f,tower]); Error("whereInGB : [f,myset]: f could not be found in the myset."); } def SunitOfFormat(pos,forms) { @@ -628,15 +797,7 @@ 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 = [ ]; @@ -657,10 +818,13 @@ 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); @@ -670,7 +834,7 @@ def SwhereInTower(f,tower) { q = MonomialPart(f); if (p == q) return(i); } - Println([f,tower]); + Sprintln([f,tower]); Error("[f,tower]: f could not be found in the tower."); } @@ -682,21 +846,23 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) local i, j, myindex, p, bases, tower2, gi, gj, si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2, vdeg,vdeg_reduced; - Println("SpairAndReduction:"); + Sprintln("SpairAndReduction:"); 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]); + Sprintln(["p and bases ",p,bases]); 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."); } tower2 = StowerOf(tower,level-1); SsetTower(tower2); + Sprintln(["level=",level]); + Sprintln(["tower2=",tower2]); /** sm1(" show_ring "); */ gi = Stoes_vec(bases[i]); @@ -707,35 +873,39 @@ def SpairAndReduction(skel,level,ii,freeRes,tower,ww) 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]); + Sprintln([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); + Sprint("[gi, gj] = "); Sprintln([gi,gj]); + sm1(" [(Homogenize)] system_variable "); + Sprint("Reduce the element "); Sprintln(si*gi+sj*gj); + Sprint("by "); Sprintln(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. */ /* 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); + Sprint("vdegree of the original = "); Sprintln(vdeg); + Sprint("vdegree of the remainder = "); Sprintln(vdeg_reduced); 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; + + 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. */ /* pos2 is the place to put a new GB at level-1. */ - Println(ans); + Sprintln(ans); return(ans); } @@ -768,24 +938,6 @@ 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; @@ -843,15 +995,41 @@ 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; - r = SlaScala(g); + betti_levelplus, newbases, i, j,qq, tminRes,bettiTable, ansSminimal; + 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); + Sprintln2("BettiTable ------"); + if (Sverbose || Sverbose2) {sm1_pmat(bettiTable);} minRes = SnewArrayOfFormat(freeRes); seq = 0; maxSeq = SgetMaxSeq(redundantTable); @@ -861,12 +1039,12 @@ def Sminimal(g) { } seq=maxSeq+1; while (seq > 1) { - seq--; + seq--; Sprint2(seq); for (level = 0; level < maxLevel; level++) { betti = Length(freeRes[level]); for (q = 0; q<betti; q++) { 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) { bases = freeRes[level+1]; dr = reducer[level,q]; @@ -879,10 +1057,10 @@ def Sminimal(g) { 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); + Sprintln(["level, q =", level,q]); + Sprintln("bases="); if (Sverbose) {sm1_pmat(bases); } + Sprintln("dr="); if (Sverbose) {sm1_pmat(dr);} + Sprintln("newbases="); if (Sverbose) {sm1_pmat(newbases);} minRes[level+1] = newbases; freeRes = minRes; #ifdef DEBUG @@ -892,7 +1070,7 @@ def Sminimal(g) { 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]); + Sprint("redundantTable ="); sm1_pmat(redundantTable[level]); Error("Stop in Sminimal for debugging."); } } @@ -904,10 +1082,28 @@ def Sminimal(g) { } } } - return([Stetris(minRes,redundantTable), - [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]); + tminRes = Stetris(minRes,redundantTable); + ansSminimal = [SpruneZeroRow(tminRes), tminRes, + [ minRes, redundantTable, reducer,r[3],r[4]],r[0],r[5]]; + Sprintln2(" "); + Println("------------ Note -----------------------------"); + Println("To get shift vectors, use Reparse and SgetShifts(resmat,w)"); + Println("To get initial of the complex, use Reparse and Sinit_w(resmat,w)"); + Println("0: minimal resolution, 3: Schreyer resolution "); + Println("------------ Resolution Summary --------------"); + Print("Betti numbers : "); + Println(Join([Length(ansSminimal[0,0,0])],Map(ansSminimal[0],"Length"))); + Print("Betti numbers of the Schreyer frame: "); + Println(Join([Length(ansSminimal[3,0,0])],Map(ansSminimal[3],"Length"))); + Println("-----------------------------------------------"); + + sm1(" restoreEnvAfterResolution "); + Sordinary = false; + + return(ansSminimal); /* r[4] is the redundantTable_ordinary */ /* r[0] is the freeResolution */ + /* r[5] is the skelton */ } @@ -976,9 +1172,11 @@ def Stetris(freeRes,redundantTable) { }else{ newbases = bases; } - Println(["level=", level]); - sm1_pmat(bases); - sm1_pmat(newbases); + Sprintln(["level=", level]); + if (Sverbose){ + sm1_pmat(bases); + sm1_pmat(newbases); + } minRes[level] = newbases; } @@ -1040,21 +1238,15 @@ 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)); */ + 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.", + "See also Sminimal, Sannfs3.", "Example: a=Sannfs2(\"x^3-y^2\");", " b=a[0]; sm1_pmat(b);", " b[1]*b[0]:", @@ -1062,414 +1254,33 @@ HelpAdd(["Sannfs2", " 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"); + a=Sannfs2("x*y*(x-y)"); +*/ -/* Do not forget to turn on TOTAL_STRATEGY */ -def Sannfs2_laScala(f) { - local p,pp; - p = Sannfs(f,"x,y"); - /* Do not make laplace transform. - sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set "); - p = [p]; - */ - Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]); - pp = Map(p[0],"Spoly"); - return(Sminimal(pp)); -} -def Sannfs2_laScala2(f) { - local p,pp; - p = Sannfs(f,"x,y"); - sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set "); - p = [p]; - Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1], - ["x",-1,"y",-1,"Dx",1,"Dy",1]]); - pp = Map(p[0],"Spoly"); - return(Sminimal(pp)); -} - def Sannfs3(f) { 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)); + 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_v, Sannfs2.", + "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]:"]]); -/* - The betti numbers of most examples are 2,1. (0-th and 1-th). - a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2. - a=Sannfs2("x^3-y^2-x"); : it causes an error. It should be fixed. - a=Sannfs2("x*y*(x-y)"); : it causes an error. It should be fixed. - -*/ -def Sannfs3_laScala2(f) { - local p,pp; - p = Sannfs(f,"x,y,z"); - sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set "); - Sweyl("x,y,z",[["x",1,"y",1,"z",1,"Dx",1,"Dy",1,"Dz",1,"h",1], - ["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]); - pp = Map(p,"Spoly"); - return(Sminimal(pp)); -} - -/* The below does not use LaScala-Stillman's algorithm. */ -def Sschreyer(g) { - local rf, tower, reductionTable, skel, redundantTable, bases, - strategy, maxOfStrategy, height, level, n, i, - freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww, - redundantTable_ordinary, redundant_seq_ordinary, - reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase; - /* extern WeightOfSweyl; */ - ww = WeightOfSweyl; - Print("WeghtOfSweyl="); Println(WeightOfSweyl); - rf = SresolutionFrameWithTower(g); - redundant_seq = 1; redundant_seq_ordinary = 1; - tower = rf[1]; - reductionTable = SgenerateTable(tower); - skel = rf[2]; - redundantTable = SnewArrayOfFormat(rf[1]); - redundantTable_ordinary = SnewArrayOfFormat(rf[1]); - reducer = SnewArrayOfFormat(rf[1]); - freeRes = SnewArrayOfFormat(rf[1]); - bettiTable = SsetBettiTable(rf[1],g); - - height = Length(reductionTable); - for (level = 0; level < height; level++) { - n = Length(reductionTable[level]); - for (i=0; i<n; i++) { - Println([level,i]); - Print("Processing "); Print([level,i]); - if (level == 0) { - if (IsNull(redundantTable[level,i])) { - bases = freeRes[level]; - /* Println(["At floor : GB=",i,bases,tower[0,i]]); */ - pos = SwhereInGB(tower[0,i],rf[3,0]); - bases[i] = rf[3,0,pos]; - /* redundantTable[level,i] = 0; - redundantTable_ordinary[level,i] = 0; */ - freeRes[level] = bases; - /* Println(["GB=",i,bases,tower[0,i]]); */ - } - }else{ /* level >= 1 */ - if (IsNull(redundantTable[level,i])) { - bases = freeRes[level]; - f = SpairAndReduction2(skel,level,i,freeRes,tower, - ww,redundantTable); - if (f[0] != Poly("0")) { - place = f[3]; - /* (level-1, place) is the place for f[0], - which is a newly obtained GB. */ -#ifdef ORDINARY - redundantTable[level-1,place] = redundant_seq; - redundant_seq++; -#else - if (f[4] > f[5]) { - /* Zero in the gr-module */ - Print("v-degree of [org,remainder] = "); - Println([f[4],f[5]]); - Print("[level,i] = "); Println([level,i]); - redundantTable[level-1,place] = 0; - }else{ - redundantTable[level-1,place] = redundant_seq; - redundant_seq++; - } -#endif - redundantTable_ordinary[level-1,place] - =redundant_seq_ordinary; - redundant_seq_ordinary++; - bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */ - /* redundantTable[level,i] = 0; - redundantTable_ordinary[level,i] = 0; */ - /* i must be equal to f[2], I think. Double check. */ - - /* Correction Of Constant */ - /* Correction of syzygy */ - c2 = f[6]; /* or -f[6]? Double check. */ - Print("c2="); Println(c2); - nn = Length(bases); - for (ii=0; ii<nn;ii++) { - if ((ii != i) && (! IsNull(bases[ii]))) { - m = Length(bases[ii]); - for (jj=0; jj<m; jj++) { - if (jj != place) { - bases[ii,jj] = bases[ii,jj]*c2; - } - } - } - } - - Print("Old freeRes[level] = "); sm1_pmat(freeRes[level]); - freeRes[level] = bases; - Print("New freeRes[level] = "); sm1_pmat(freeRes[level]); - - /* Update the freeRes[level-1] */ - Print("Old freeRes[level-1] = "); sm1_pmat(freeRes[level-1]); - bases = freeRes[level-1]; - bases[place] = f[0]; - freeRes[level-1] = bases; - Print("New freeRes[level-1] = "); sm1_pmat(freeRes[level-1]); - - reducer[level-1,place] = f[1]-SunitOfFormat(place,f[1]); - /* This reducer is different from that of SlaScala(). */ - - reducerBasis = reducer[level-1]; - nn = Length(reducerBasis); - for (ii=0; ii<nn;ii++) { - if ((ii != place) && (! IsNull(reducerBasis[ii]))) { - m = Length(reducerBasis[ii]); - for (jj=0; jj<m; jj++) { - if (jj != place) { - reducerBasis[ii,jj] = reducerBasis[ii,jj]*c2; - } - } - } - } - reducer[level-1] = reducerBasis; - - }else{ - /* redundantTable[level,i] = 0; */ - bases = freeRes[level]; - bases[i] = f[1]; /* Put the syzygy. */ - freeRes[level] = bases; - } - } /* end of level >= 1 */ - } - } /* i loop */ - - /* Triangulate reducer */ - if (level >= 1) { - Println(" "); - Print("Triangulating reducer at level "); Println(level-1); - Println("freeRes[level]="); sm1_pmat(freeRes[level]); - reducerBase = reducer[level-1]; - Print("reducerBase="); Println(reducerBase); - Println("Compare freeRes[level] and reducerBase (put -1)"); - m = Length(reducerBase); - for (ii=m-1; ii>=0; ii--) { - if (!IsNull(reducerBase[ii])) { - for (jj=ii-1; jj>=0; jj--) { - if (!IsNull(reducerBase[jj])) { - if (!IsZero(reducerBase[jj,ii])) { - /* reducerBase[ii,ii] should be always constant. */ - reducerBase[jj] = reducerBase[ii,ii]*reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii]; - } - } - } - } - } - Println("New reducer"); - sm1_pmat(reducerBase); - reducer[level-1] = reducerBase; - } - - } /* level loop */ - n = Length(freeRes); - freeResV = SnewArrayOfFormat(freeRes); - for (i=0; i<n; i++) { - bases = freeRes[i]; - bases = Sbases_to_vec(bases,bettiTable[i]); - freeResV[i] = bases; - } - - /* Mark the non-redundant elements. */ - for (i=0; i<n; i++) { - m = Length(redundantTable[i]); - for (jj=0; jj<m; jj++) { - if (IsNull(redundantTable[i,jj])) { - redundantTable[i,jj] = 0; - } - } - } - - - return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]); -} - -def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) { - local i, j, myindex, p, bases, tower2, gi, gj, - si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2, - vdeg,vdeg_reduced,n,c2; - Println("SpairAndReduction2 : -------------------------"); - - if (level < 1) Error("level should be >= 1 in SpairAndReduction."); - p = skel[level,ii]; - myindex = p[0]; - i = myindex[0]; j = myindex[1]; - bases = freeRes[level-1]; - Println(["p and bases ",p,bases]); - if (IsNull(bases[i]) || IsNull(bases[j])) { - Println([level,i,j,bases[i],bases[j]]); - Error("level, i, j : bases[i], bases[j] must not be NULL."); - } - - tower2 = StowerOf(tower,level-1); - SsetTower(tower2); - /** sm1(" show_ring "); */ - - gi = Stoes_vec(bases[i]); - gj = Stoes_vec(bases[j]); - - ssp = Sspolynomial(gi,gj); - si = ssp[0,0]; - sj = ssp[0,1]; - syzHead = si*es^i; - /* This will be the head term, I think. But, double check. */ - Println([si*es^i,sj*es^j]); - - Print("[gi, gj] = "); Println([gi,gj]); - sm1(" [(Homogenize)] system_variable message "); - Print("Reduce the element "); Println(si*gi+sj*gj); - Print("by "); Println(bases); - - tmp = Sreduction(si*gi+sj*gj, bases); - - Print("result is "); Println(tmp); - if (!IsZero(tmp[0])) { - Print("Error: base = "); - Println(Map(bases,"Stoes_vec")); - Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring."); - } - t_syz = tmp[2]; - si = si*tmp[1]+t_syz[i]; - sj = sj*tmp[1]+t_syz[j]; - t_syz[i] = si; - t_syz[j] = sj; - - c2 = null; - /* tmp[0] must be zero */ - n = Length(t_syz); - for (i=0; i<n; i++) { - if (IsConstant(t_syz[i])){ - if (!IsZero(t_syz[i])) { - if (IsNull(redundantTable[level-1,i])) { - /* i must equal to pos2 below. */ - c2 = -t_syz[i]; - tmp[0] = c2*Stoes_vec(freeRes[level-1,i]); - t_syz[i] = 0; - /* tmp[0] = t_syz . g */ - /* break; does not work. Use */ - i = n; - } - } - } - } - - /* This is essential part for V-minimal resolution. */ - /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */ - vdeg = SvDegree(si*gi,tower,level-1,ww); - vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww); - Print("vdegree of the original = "); Println(vdeg); - Print("vdegree of the remainder = "); Println(vdeg_reduced); - - if (!IsNull(vdeg_reduced)) { - if (vdeg_reduced < vdeg) { - Println("--- Special in V-minimal!"); - Println(tmp[0]); - Println("syzygy="); sm1_pmat(t_syz); - Print("[vdeg, vdeg_reduced] = "); Println([vdeg,vdeg_reduced]); - } - } - - - 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); -} - -HelpAdd(["Sminimal_v", -["It constructs the V-minimal free resolution from the Schreyer resolution", - "step by step.", - "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."]]); - - -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 @@ -1612,41 +1423,13 @@ def Skernel(m,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 SisExact_h(c,v) { - local a; - v = ToString_array(v); - a = [c,v]; - sm1(a," isExact /FunctionValue set "); -} def SisComplex(a) { local n,i,j,k,b,p,q; @@ -1667,4 +1450,155 @@ def SisComplex(a) { } } return(true); -} \ No newline at end of file +} + +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" +]]); + +/* + Output of S* functions may cause a trouble because it uses Schreyer orders. + In this case, use ReParse(). +*/ + +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 "); Printl(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+1); + m0 = NewArray(Length(resmat[0,0])); + ans[0] = m0; + for (i=0; i<n; i++) { + ans[i+1] = SgetShift(resmat[i],w,m0); + m0 = ans[i+1]; + } + return(ans); +} +HelpAdd(["SgetShifts", +["SgetShifts(resmat,w) returns the shift vectors of the resolution resmat", + " with respect to w with the shift m.", + "Note that the order of the ring and the weight w must be the same.", + "Zero row is not allowed.", + "Example: a=Sannfs2(\"x^3-y^2\");", + " b=a[0]; w = [\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1];", + " Sweyl(\"x,y\",[w]); b = Reparse(b);", + " SgetShifts(b,w):"]]); + +def Sinit_w(resmat,w) { + local shifts,ans,n,i,m,mat,j; + shifts = SgetShifts(resmat,w); + n = Length(resmat); + ans = NewArray(n); + for (i=0; i<n; i++) { + m = shifts[i]; + mat = ScopyArray(resmat[i]); + for (j=0; j<Length(mat); j++) { + mat[j] = Init_w_m(mat[j],w,m); + } + ans[i] = mat; + } + return(ans); +} +HelpAdd(["Sinit_w", +["Sinit_w(resmat,w) returns the initial of the complex resmat with respect to the weight w.", + "Example: a=Sannfs2(\"x^3-y^2\");", + " b=a[0]; w = [\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1];", + " Sweyl(\"x,y\",[w]); b = Reparse(b);", + " c=Sinit_w(b,w); c:" +]]); + +/* This method does not work, because we have zero rows. + Think about it later. */ +def SbettiTable(rtable) { + local ans,i,j,pp; + ans = SnewArrayOfFormat(rtable); + for (i=0; i<Length(rtable); i++) { + pp = 0; + for (j=0; j<Length(rtable[i]); j++) { + if (rtable[i,j] != 0) {pp = pp+1;} + } + ans[i] = pp; + } + return(ans); +} + +def BfRoots1(G,V) { + local bb,ans; + sm1(" /BFparlist [ ] def "); + if (IsString(V)) { + sm1(" [ V to_records pop ] /V set "); + }else { + sm1(" V { toString } map /V set "); + } + sm1(" /BFvarlist V def "); + + sm1(" G flatten { toString } map /G set "); + sm1(" G V bfm /bb set "); + if (IsSm1Integer(bb)) { + return([ ]); + } + sm1(" bb 0 get findIntegralRoots { (universalNumber) dc } map /ans set "); + return([ans, bb]); +} + +HelpAdd(["BfRoots1", +["BfRoots1(g,v) returns the integral roots of g with respect to the weight", + "vector (1,1,...,1) and the b-function itself", + "Example: BfRoots1([x*Dx-2, y*Dy-3],[x,y]);" +]]); + + +