File: [local] / OpenXM / src / k097 / lib / minimal / minimal.k (download)
Revision 1.36, Tue Jul 3 22:28:11 2007 UTC (17 years ago) by takayama
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_3_1_13b, RELEASE_1_2_3_12, HEAD, DEB_REL_1_2_3-9
Changes since 1.35: +5 -2
lines
Added a help message for Sminimal.
All free modules in the complex are sets of the row vectors.
|
/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.36 2007/07/03 22:28:11 takayama Exp $ */
#define DEBUG 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*/
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"];;
a=Sminimal(v);
b=a[0];
b[1]*b[0]:
b[2]*b[1]:
a = test0();
b = a[0];
b[1]*b[0]:
b[2]*b[1]:
a = Sminimal(b[0]);
a = test1();
b=a[0];
b[1]*b[0]:
b[2]*b[1]:
*/
/* 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 ");
SonAutoReduce = true;
def Factor(f) {
sm1(f, " fctr /FunctionValue set");
}
def Reverse(f) {
sm1(f," reverse /FunctionValue set");
}
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");
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];
frame=SresolutionFrame(f);
Println(frame);
/* return(frame); */
return(SlaScala(f));
}
def test1() {
local f;
Sweyl("x,y,z");
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];
return(Sminimal(f));
}
def Sweyl(v,w) {
/* extern WeightOfSweyl ; */
local ww,i,n;
if(Length(Arglist) == 1) {
sm1(" [v s_ring_of_differential_operators 0 [(schreyer) 1]] define_ring ");
sm1(" define_ring_variables ");
sm1(" [ v to_records pop ] /ww set ");
n = Length(ww);
WeightOfSweyl = NewArray(n*4);
for (i=0; i< n; i++) {
WeightOfSweyl[2*i] = ww[i];
WeightOfSweyl[2*i+1] = 1;
}
for (i=0; i< n; i++) {
WeightOfSweyl[2*n+2*i] = AddString(["D",ww[i]]);
WeightOfSweyl[2*n+2*i+1] = 1;
}
}else{
sm1(" [v s_ring_of_differential_operators w s_weight_vector 0 [(schreyer) 1]] define_ring ");
sm1(" define_ring_variables ");
WeightOfSweyl = w[0];
}
}
def Spoly(f) {
sm1(f, " toString tparse /FunctionValue set ");
}
def SreplaceZeroByZeroPoly(f) {
if (IsArray(f)) {
return(Map(f,"SreplaceZeroByZeroPoly"));
}else{
if (IsInteger(f)) {
return(Poly(ToString(f)));
}else{
return(f);
}
}
}
def Shomogenize(f) {
f = SreplaceZeroByZeroPoly(f);
if (IsArray(f)) {
sm1(f," sHomogenize2 /FunctionValue set ");
/* sm1(f," {sHomogenize2} map /FunctionValue set "); */
/* Is it correct? Double check.*/
}else{
sm1(f, " sHomogenize /FunctionValue set ");
}
}
def StoTower() {
sm1(" [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv ");
}
def SsetTower(tower) {
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, nohomog,i,n;
/* extern Sordinary */
nohomog = false;
count = -1; Sordinary = false; /* default value for options. */
if (Length(Arglist) >= 2) {
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 ");
/* If I do not put this macro, homogenization
make a strange behavior. For example,
[(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
[(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
4/19, 2000.
*/
sm1(" (mmLarger) (matrix) switch_function ");
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 ");
}
gbasis = Sgroebner(g);
g = gbasis[0];
if (SonAutoReduce) {
sm1("[ (AutoReduce) autof] system_variable ");
}
g = Init(g);
/* sm1(" setupEnvForResolution-sugar "); */
/* -sugar is fine? */
sm1(" setupEnvForResolution ");
Sprintln(g);
startingGB = g;
/* ans = [ SzeroMap(g) ]; It has not been implemented. see resol1.withZeroMap */
ans = [ ];
gbTower = [ ];
skelton = [ ];
while (true) {
/* sm1(g," res0Frame /ff set "); */
withSkel = Sres0FrameWithSkelton(g);
ff = withSkel[0];
ans = Append(ans, ff[0]);
gbTower = Join([ ff[1] ], gbTower);
skelton = Join([ withSkel[1] ], skelton);
g = ff[0];
if (Length(g) == 0) break;
SsetTower( gbTower );
if (count == 0) break;
count = count - 1;
}
return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
}
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,opt);
return(ans[0]);
}
/* ---------------------------- */
def ToGradedPolySet(g) {
sm1(g," (gradedPolySet) dc /FunctionValue set ");
}
def NewPolynomialVector(size) {
sm1(size," (integer) dc newPolyVector /FunctionValue set ");
}
def SturnOffHomogenization() {
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 \
{ Sverbose { \
(Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message } { } ifelse \
[(Homogenize) 1] system_variable \
[(ReduceLowerTerms) 1] system_variable \
} { } ifelse \
");
}
def SschreyerSkelton(g) {
sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
}
def Stoes(g) {
if (IsArray(g)) {
sm1(g," {toes} map /FunctionValue set ");
}else{
sm1(g," toes /FunctionValue set ");
}
}
def Stoes_vec(g) {
sm1(g," toes /FunctionValue set ");
}
def Sres0Frame(g) {
local ans;
ans = Sres0FrameWithSkelton(g);
return(ans[0]);
}
def Sres0FrameWithSkelton(g) {
local t_syz, nexttower, m, t_gb, skel, betti,
gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;
SturnOffHomogenization();
g = Stoes(g);
skel = SschreyerSkelton(g);
/* Print("Skelton is ");
sm1_pmat(skel); */
betti = Length(skel);
gLength = Length(g);
grG = ToGradedPolySet(g);
syzAll = NewPolynomialVector(betti);
for (k=0; k<betti; k++) {
pair = skel[k];
i = pair[0,0];
j = pair[0,1];
si = pair[1,0];
sj = pair[1,1];
/* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
Sprint(".");
t_syz = NewPolynomialVector(gLength);
t_syz[i] = si;
t_syz[j] = sj;
syzAll[k] = t_syz;
}
t_syz = syzAll;
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 , ...]
[y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
[3*es^3*y*Dy-es^5*x ]
*/
nexttower = Init(g);
SturnOnHomogenization();
return([[t_syz, nexttower],skel]);
}
def StotalDegree(f) {
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;
n = Length(w);
neww = NewArray(n);
for (i=0; i<n; i=i+2) {
neww[i] = ToString(w[i]);
}
for (i=1; i<n; i=i+2) {
neww[i] = IntegerToSm1Integer(w[i]);
}
sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
}
/* This is not satisfactory. */
def SinitOfArray(f) {
local p,pos,top;
if (IsArray(f)) {
sm1(f," toes init /p set ");
sm1(p," (es). degree (universalNumber) dc /pos set ");
return([Init(f[pos]),pos]);
} else {
return(Init(f));
}
}
def test_SinitOfArray() {
local f, frame,p,tower,i,j,k;
Sweyl("x,y,z");
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);
if (Sverbose) {
sm1_pmat(p);
sm1_pmat(SgenerateTable(p[1]));
}
return(p);
frame = p[0];
sm1_pmat(p[1]);
sm1_pmat(frame);
sm1_pmat(Map(frame[0],"SinitOfArray"));
sm1_pmat(Map(frame[1],"SinitOfArray"));
return(p);
}
/* f is assumed to be a monomial with toes. */
def Sdegree(f,tower,level) {
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;
/*
Sprint("SgenerateTable: tower=");Sprintln(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)
+ OFFSET;
/* Println([i,j,ans_at_each_floor[j]]); */
}
ans[i] = ans_at_each_floor;
}
return(ans);
}
Sweyl("x,y,z");
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]];
/* SresolutionFrameWithTower(v); */
def SnewArrayOfFormat(p) {
if (IsArray(p)) {
return(Map(p,"SnewArrayOfFormat"));
}else{
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 */
if (IsArray(a)) {
n = Length(a);
for (i=0; i<n; i++) {
if (IsArray(a[i])) {
tt = SminOfStrategy(a[i]);
if (tt < ans) ans = tt;
}else{
if (a[i] < ans) ans = a[i];
}
}
}else{
if (a < ans) ans = a;
}
return(ans);
}
def SmaxOfStrategy(a) {
local n,i,ans,tt;
ans = -100000; /* very small number */
if (IsArray(a)) {
n = Length(a);
for (i=0; i<n; i++) {
if (IsArray(a[i])) {
tt = SmaxOfStrategy(a[i]);
if (tt > ans) ans = tt;
}else{
if (a[i] > ans) ans = a[i];
}
}
}else{
if (a > ans) ans = a;
}
return(ans);
}
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,
reductionTable_tmp;
/* extern WeightOfSweyl; */
ww = WeightOfSweyl;
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]);
reducer = SnewArrayOfFormat(rf[1]);
freeRes = SnewArrayOfFormat(rf[1]);
bettiTable = SsetBettiTable(rf[1],g);
strategy = SminOfStrategy( reductionTable );
maxOfStrategy = SmaxOfStrategy( reductionTable );
height = Length(reductionTable);
while (strategy <= maxOfStrategy) {
for (level = 0; level < height; level++) {
n = Length(reductionTable[level]);
reductionTable_tmp = ScopyArray(reductionTable[level]);
while (SthereIs(reductionTable_tmp,strategy)) {
i = SnextI(reductionTable_tmp,strategy,redundantTable,
skel,level,freeRes);
Sprintln([level,i]);
reductionTable_tmp[i] = -200000;
if (reductionTable[level,i] == 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];
/* 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 = SpairAndReduction(skel,level,i,freeRes,tower,ww);
if (f[0] != Poly("0")) {
place = f[3];
/* (level-1, place) is the place for f[0],
which is a newly obtained GB. */
if (Sordinary) {
redundantTable[level-1,place] = redundant_seq;
redundant_seq++;
}else{
if (f[4] > f[5]) {
/* Zero in the gr-module */
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++;
}
}
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. */
freeRes[level] = bases;
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 */
}
}
}
strategy++;
}
Sprintln2(" ");
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;
}
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);
}
}
}
}
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.");
}
def SsetBettiTable(freeRes,g) {
local level,i, n,bases,ans;
ans = NewArray(Length(freeRes)+1);
n = Length(freeRes);
if (IsArray(g[0])) {
ans[0] = Length(g[0]);
}else{
ans[0] = 1;
}
for (level=0; level<n; level++) {
bases = freeRes[level];
if (IsArray(bases)) {
ans[level+1] = Length(bases);
}else{
ans[level+1] = 1;
}
}
return(ans);
}
def SwhereInGB(f,tower) {
local i,n,p,q;
n = Length(tower);
for (i=0; i<n; i++) {
p = MonomialPart(tower[i]);
q = MonomialPart(f);
if (p == q) return(i);
}
Sprintln([f,tower]);
Error("whereInGB : [f,myset]: f could not be found in the myset.");
}
def SunitOfFormat(pos,forms) {
local ans,i,n;
n = Length(forms);
ans = NewArray(n);
for (i=0; i<n; i++) {
if (i != pos) {
ans[i] = Poly("0");
}else{
ans[i] = Poly("1");
}
}
return(ans);
}
def StowerOf(tower,level) {
local ans,i;
ans = [ ];
if (level == 0) return([[]]);
for (i=0; i<level; i++) {
ans = Append(ans,tower[i]);
}
return(Reverse(ans));
}
def Sspolynomial(f,g) {
if (IsArray(f)) {
f = Stoes_vec(f);
}
if (IsArray(g)) {
g = Stoes_vec(g);
}
sm1("f g spol /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);
n = Length(tower);
for (i=0; i<n; i++) {
p = MonomialPart(tower[i]);
q = MonomialPart(f);
if (p == q) return(i);
}
Sprintln([f,tower]);
Error("[f,tower]: f could not be found in the tower.");
}
def Stag(f) {
sm1(f," tag (universalNumber) dc /FunctionValue set");
}
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;
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];
Sprintln(["p and bases ",p,bases]);
if (IsNull(bases[i]) || IsNull(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]);
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. */
Sprintln([si*es^i,sj*es^j]);
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);
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);
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. */
Sprintln(ans);
return(ans);
}
def Sreduction(f,myset) {
local n, indexTable, set2, i, j, tmp, t_syz;
n = Length(myset);
indexTable = NewArray(n);
set2 = [ ];
j = 0;
for (i=0; i<n; i++) {
if (IsNull(myset[i])) {
indexTable[i] = -1;
/* }else if (myset[i] == Poly("0")) {
indexTable[i] = -1; */
}else{
set2 = Append(set2,Stoes_vec(myset[i]));
indexTable[i] = j;
j++;
}
}
sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
t_syz = NewArray(n);
for (i=0; i<n; i++) {
if (indexTable[i] != -1) {
t_syz[i] = tmp[2, indexTable[i]];
}else{
t_syz[i] = Poly("0");
}
}
return([tmp[0],tmp[1],t_syz]);
}
def Sfrom_es(f,size) {
local c,ans, i, d, myes, myee, j,n,r,ans2;
if (Length(Arglist) < 2) size = -1;
if (IsArray(f)) return(f);
r = RingOf(f);
myes = PolyR("es",r);
myee = PolyR("e_",r);
if (Degree(f,myee) > 0 && size == -1) {
if (size == -1) {
sm1(f," (array) dc /ans set");
return(ans);
}
}
/*
Coefficients(x^2-1,x):
[ [ 2 , 0 ] , [ 1 , -1 ] ]
*/
if (Degree(f,myee) > 0) {
c = Coefficients(f,myee);
}else{
c = Coefficients(f,myes);
}
if (size < 0) {
size = c[0,0]+1;
}
ans = NewArray(size);
for (i=0; i<size; i++) {ans[i] = 0;}
n = Length(c[0]);
for (j=0; j<n; j++) {
d = c[0,j];
ans[d] = c[1,j];
}
return(ans);
}
def Sbases_to_vec(bases,size) {
local n, giveSize, newbases,i;
/* bases = [1+es*x, [1,2,3*x]] */
if (Length(Arglist) > 1) {
giveSize = true;
}else{
giveSize = false;
}
n = Length(bases);
newbases = NewArray(n);
for (i=0; i<n; i++) {
if (giveSize) {
newbases[i] = Sfrom_es(bases[i], size);
}else{
newbases[i] = Sfrom_es(bases[i]);
}
}
return(newbases);
}
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.",
" ---> D^{m_3} --b[2]--> D^{m_2} --b[1]--> D^{m_1} --b[0]--> D^{m_0} ",
" Here D^{m_i} are the set of row vectors. "
]]);
def Sminimal(g,opt) {
local r, freeRes, redundantTable, reducer, maxLevel,
minRes, seq, maxSeq, level, betti, q, bases, dr,
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);
maxLevel = Length(freeRes);
for (level = 0; level < maxLevel; level++) {
minRes[level] = freeRes[level];
}
seq=maxSeq+1;
while (seq > 1) {
seq--; Sprint2(seq);
for (level = 0; level < maxLevel; level++) {
betti = Length(freeRes[level]);
for (q = 0; q<betti; q++) {
if (redundantTable[level,q] == seq) {
Sprint("[seq,level,q]="); Sprintln([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;
}
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
for (qq=0; qq<betti; qq++) {
if ((redundantTable[level,qq] >= seq) &&
(redundantTable[level,qq] <= maxSeq)) {
for (i=0; i<betti_levelplus; i++) {
if (!IsZero(newbases[i,qq])) {
Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
Sprint("redundantTable ="); sm1_pmat(redundantTable[level]);
Error("Stop in Sminimal for debugging.");
}
}
}
}
#endif
}
}
}
}
}
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 */
}
def IsZero(f) {
if (IsPolynomial(f)) {
return( f == Poly("0"));
}else if (IsInteger(f)) {
return( f == 0);
}else if (IsSm1Integer(f)) {
return( f == true );
}else if (IsDouble(f)) {
return( f == 0.0 );
}else if (IsRational(f)) {
return(IsZero(Denominator(f)));
}else{
Error("IsZero: cannot deal with this data type.");
}
}
def SgetMaxSeq(redundantTable) {
local level,i,n,ans, levelMax,bases;
levelMax = Length( redundantTable );
ans = 0;
for (level = 0; level < levelMax; level++) {
bases = redundantTable[level];
n = Length(bases);
for (i=0; i<n; i++) {
if (IsInteger( bases[i] )) {
if (bases[i] > ans) {
ans = bases[i];
}
}
}
}
return(ans);
}
def Stetris(freeRes,redundantTable) {
local level, i, j, resLength, minRes,
bases, newbases, newbases2;
minRes = SnewArrayOfFormat(freeRes);
resLength = Length( freeRes );
for (level=0; level<resLength; level++) {
bases = freeRes[level];
newbases = SnewArrayOfFormat(bases);
betti = Length(bases); j = 0;
/* Delete rows */
for (i=0; i<betti; i++) {
if (redundantTable[level,i] < 1) {
newbases[j] = bases[i];
j++;
}
}
bases = SfirstN(newbases,j);
if (level > 0) {
/* Delete columns */
newbases = Transpose(bases);
betti = Length(newbases); j = 0;
newbases2 = SnewArrayOfFormat(newbases);
for (i=0; i<betti; i++) {
if (redundantTable[level-1,i] < 1) {
newbases2[j] = newbases[i];
j++;
}
}
newbases = Transpose(SfirstN(newbases2,j));
}else{
newbases = bases;
}
Sprintln(["level=", level]);
if (Sverbose){
sm1_pmat(bases);
sm1_pmat(newbases);
}
minRes[level] = newbases;
}
return(minRes);
}
def SfirstN(bases,k) {
local ans,i;
ans = NewArray(k);
for (i=0; i<k; i++) {
ans[i] = bases[i];
}
return(ans);
}
/* usage: tt is tower. ww is weight.
a = SresolutionFrameWithTower(v);
tt = a[1];
ww = [x,1,y,1,Dx,1,Dy,1];
SvDegree(x*es,tt,1,ww):
In(17)=tt:
[[2*x*Dx , e_*x^2 , e_*x*y , 3*x^2*Dy , e_*y^3 , 9*x*y*Dy^2 , 27*y^2*Dy^3 ] ,
[es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ] ,
[3*es^3*y*Dy ] ]
In(18)=SvDegree(x*es,tt,1,ww):
3
In(19)=SvDegree(x*es^3,tt,1,ww):
4
In(20)=SvDegree(x,tt,2,ww):
4
*/
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));
}
i = Degree(f,es);
ans = Sord_w(f,w) +
SvDegree(tower[level-1,i],tower,level-1,w);
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]]);
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");
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"
]]);
/*
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]);"
]]);
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to minimal.k CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / k097 / lib / minimal