File: [local] / OpenXM / src / k097 / lib / minimal / minimal.k (download)
Revision 1.7, Sat May 6 10:35:33 2000 UTC (24 years, 2 months ago) by takayama
Branch: MAIN
Changes since 1.6: +14 -3
lines
Temporary change of repository to work on a different computer.
|
/* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.7 2000/05/06 10:35:33 takayama Exp $ */
#define DEBUG 1
/* #define ORDINARY 1 */
/* If you run this program on openxm version 1.1.2 (FreeBSD),
make a symbolic link by the command
ln -s /usr/bin/cpp /lib/cpp
*/
#define OFFSET 0
#define TOTAL_STRATEGY
/* #define OFFSET 20*/
/* Test sequences.
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]:
*/
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");
}
def Reverse(f) {
sm1(f," reverse /FunctionValue set");
}
def Sgroebner(f) {
sm1(" [f] groebner /FunctionValue set");
}
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 ");
}
def SresolutionFrameWithTower(g,opt) {
local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
gbasis;
if (Length(Arglist) >= 2) {
if (IsInteger(opt)) count = opt;
}else{
count = -1;
}
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 ");
g = Map(g,"Shomogenize");
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 ");
Println(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]",
"Example: Sweyl(\"x,y\");",
" a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
def SresolutionFrame(f,opt) {
local ans;
ans = SresolutionFrameWithTower(f);
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
{ (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
[(Homogenize) 0] system_variable
[(ReduceLowerTerms) 0] system_variable
} { } ifelse
");
}
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
");
}
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 */
Print(".");
t_syz = NewPolynomialVector(gLength);
t_syz[i] = si;
t_syz[j] = sj;
syzAll[k] = t_syz;
}
t_syz = syzAll;
Print("Done. betti="); Println(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) {
sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");
}
/* 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);
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);
#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;
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) {
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;
Print("WeightOfSweyl="); 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);
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);
Println([level,i]);
reductionTable_tmp[i] = -200000;
if (reductionTable[level,i] == strategy) {
Print("Processing "); Print([level,i]);
Print(" Strategy = "); Println(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. */
#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. */
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++;
}
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]);
}
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);
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);
}
Println([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 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 = [ ];
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");
}
def MonomialPart(f) {
sm1(" [(lmonom) f] gbext /FunctionValue set ");
}
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);
}
Println([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;
Println("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]);
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);
/* 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);
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;
pos = SwhereInTower(syzHead,tower[level]);
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);
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 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;
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);
}
def Sminimal(g) {
local r, freeRes, redundantTable, reducer, maxLevel,
minRes, seq, maxSeq, level, betti, q, bases, dr,
betti_levelplus, newbases, i, j,qq;
r = SlaScala(g);
/* 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];
}
seq=maxSeq+1;
while (seq > 1) {
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]);
if (level < maxLevel-1) {
bases = freeRes[level+1];
dr = reducer[level,q];
dr[q] = -1;
newbases = SnewArrayOfFormat(bases);
betti_levelplus = Length(bases);
/*
bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
*/
for (i=0; i<betti_levelplus; i++) {
newbases[i] = bases[i] + bases[i,q]*dr;
}
Println(["level, q =", level,q]);
Println("bases="); sm1_pmat(bases);
Println("dr="); sm1_pmat(dr);
Println("newbases="); sm1_pmat(newbases);
minRes[level+1] = newbases;
freeRes = minRes;
#ifdef DEBUG
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."]);
Print("redundantTable ="); sm1_pmat(redundantTable[level]);
Error("Stop in Sminimal for debugging.");
}
}
}
}
#endif
}
}
}
}
}
return([Stetris(minRes,redundantTable),
[ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
/* r[4] is the redundantTable_ordinary */
/* r[0] is the freeResolution */
}
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;
}
Println(["level=", level]);
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,"h",1],
["x",-1,"y",-1,"Dx",1,"Dy",1]]); */
/* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1]]); */
Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
pp = Map(p,"Spoly");
return(Sminimal_v(pp));
/* return(Sminimal(pp)); */
}
/* Do not forget to turn on TOTAL_STRATEGY */
def Sannfs2_laScala(f) {
local p,pp;
p = Sannfs(f,"x,y");
/* Do not make laplace transform.
sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
p = [p];
*/
Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
pp = Map(p[0],"Spoly");
return(Sminimal(pp));
}
def Sannfs3(f) {
local p,pp;
p = Sannfs(f,"x,y,z");
sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
pp = Map(p,"Spoly");
return(Sminimal_v(pp));
}
/*
The betti numbers of most examples are 2,1. (0-th and 1-th).
a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2.
a=Sannfs2("x^3-y^2-x"); : it causes an error. It should be fixed.
a=Sannfs2("x*y*(x-y)"); : it causes an error. It should be fixed.
*/
/* The below does not use LaScala-Stillman's algorithm. */
def Sschreyer(g) {
local rf, tower, reductionTable, skel, redundantTable, bases,
strategy, maxOfStrategy, height, level, n, i,
freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
redundantTable_ordinary, redundant_seq_ordinary,
reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
/* extern WeightOfSweyl; */
ww = WeightOfSweyl;
Print("WeghtOfSweyl="); Println(WeightOfSweyl);
rf = SresolutionFrameWithTower(g);
redundant_seq = 1; redundant_seq_ordinary = 1;
tower = rf[1];
reductionTable = SgenerateTable(tower);
skel = rf[2];
redundantTable = SnewArrayOfFormat(rf[1]);
redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
reducer = SnewArrayOfFormat(rf[1]);
freeRes = SnewArrayOfFormat(rf[1]);
bettiTable = SsetBettiTable(rf[1],g);
height = Length(reductionTable);
for (level = 0; level < height; level++) {
n = Length(reductionTable[level]);
for (i=0; i<n; i++) {
Println([level,i]);
Print("Processing "); Print([level,i]);
if (level == 0) {
if (IsNull(redundantTable[level,i])) {
bases = freeRes[level];
/* Println(["At floor : GB=",i,bases,tower[0,i]]); */
pos = SwhereInGB(tower[0,i],rf[3,0]);
bases[i] = rf[3,0,pos];
/* redundantTable[level,i] = 0;
redundantTable_ordinary[level,i] = 0; */
freeRes[level] = bases;
/* Println(["GB=",i,bases,tower[0,i]]); */
}
}else{ /* level >= 1 */
if (IsNull(redundantTable[level,i])) {
bases = freeRes[level];
f = SpairAndReduction2(skel,level,i,freeRes,tower,
ww,redundantTable);
if (f[0] != Poly("0")) {
place = f[3];
/* (level-1, place) is the place for f[0],
which is a newly obtained GB. */
#ifdef ORDINARY
redundantTable[level-1,place] = redundant_seq;
redundant_seq++;
#else
if (f[4] > f[5]) {
/* Zero in the gr-module */
Print("v-degree of [org,remainder] = ");
Println([f[4],f[5]]);
Print("[level,i] = "); Println([level,i]);
redundantTable[level-1,place] = 0;
}else{
redundantTable[level-1,place] = redundant_seq;
redundant_seq++;
}
#endif
redundantTable_ordinary[level-1,place]
=redundant_seq_ordinary;
redundant_seq_ordinary++;
bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */
/* redundantTable[level,i] = 0;
redundantTable_ordinary[level,i] = 0; */
/* i must be equal to f[2], I think. Double check. */
/* Correction Of Constant */
c2 = f[6]; /* or -f[6]? Double check. */
Print("c2="); Println(c2);
nn = Length(bases);
for (ii=0; ii<nn;ii++) {
if ((ii != place) && (! 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];
}else{
/* redundantTable[level,i] = 0; */
bases = freeRes[level];
bases[i] = f[1]; /* Put the syzygy. */
freeRes[level] = bases;
}
} /* end of level >= 1 */
}
} /* i loop */
/* Triangulate reducer */
if (level >= 1) {
Println(" ");
Print("Triangulating reducer at level "); Println(level-1);
reducerBase = reducer[level-1];
Print("reducerBase="); Println(reducerBase);
m = Length(reducerBase);
for (ii=m-1; ii>=0; ii--) {
if (!IsNull(reducerBase[ii])) {
for (jj=ii-1; jj>=0; jj--) {
if (!IsNull(reducerBase[jj])) {
if (!IsZero(reducerBase[jj,ii])) {
reducerBase[jj] = reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
}
}
}
}
}
Println("New reducer");
sm1_pmat(reducerBase);
reducer[level-1] = reducerBase;
}
} /* level loop */
n = Length(freeRes);
freeResV = SnewArrayOfFormat(freeRes);
for (i=0; i<n; i++) {
bases = freeRes[i];
bases = Sbases_to_vec(bases,bettiTable[i]);
freeResV[i] = bases;
}
/* Mark the non-redundant elements. */
for (i=0; i<n; i++) {
m = Length(redundantTable[i]);
for (jj=0; jj<m; jj++) {
if (IsNull(redundantTable[i,jj])) {
redundantTable[i,jj] = 0;
}
}
}
return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
}
def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
local i, j, myindex, p, bases, tower2, gi, gj,
si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
vdeg,vdeg_reduced,n,c2;
Println("SpairAndReduction2 : -------------------------");
if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
p = skel[level,ii];
myindex = p[0];
i = myindex[0]; j = myindex[1];
bases = freeRes[level-1];
Println(["p and bases ",p,bases]);
if (IsNull(bases[i]) || IsNull(bases[j])) {
Println([level,i,j,bases[i],bases[j]]);
Error("level, i, j : bases[i], bases[j] must not be NULL.");
}
tower2 = StowerOf(tower,level-1);
SsetTower(tower2);
/** sm1(" show_ring "); */
gi = Stoes_vec(bases[i]);
gj = Stoes_vec(bases[j]);
ssp = Sspolynomial(gi,gj);
si = ssp[0,0];
sj = ssp[0,1];
syzHead = si*es^i;
/* This will be the head term, I think. But, double check. */
Println([si*es^i,sj*es^j]);
Print("[gi, gj] = "); Println([gi,gj]);
sm1(" [(Homogenize)] system_variable message ");
Print("Reduce the element "); Println(si*gi+sj*gj);
Print("by "); Println(bases);
tmp = Sreduction(si*gi+sj*gj, bases);
Print("result is "); Println(tmp);
if (!IsZero(tmp[0])) {
Print("Error: base = ");
Println(Map(bases,"Stoes_vec"));
Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
}
t_syz = tmp[2];
si = si*tmp[1]+t_syz[i];
sj = sj*tmp[1]+t_syz[j];
t_syz[i] = si;
t_syz[j] = sj;
c2 = null;
/* tmp[0] must be zero */
n = Length(t_syz);
for (i=0; i<n; i++) {
if (IsConstant(t_syz[i])){
if (!IsZero(t_syz[i])) {
if (IsNull(redundantTable[level-1,i])) {
/* i must equal to pos2 below. */
c2 = -t_syz[i];
tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
t_syz[i] = 0;
/* tmp[0] = t_syz . g */
/* break; does not work. Use */
i = n;
}
}
}
}
/* This is essential part for V-minimal resolution. */
/* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
vdeg = SvDegree(si*gi,tower,level-1,ww);
vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
Print("vdegree of the original = "); Println(vdeg);
Print("vdegree of the remainder = "); Println(vdeg_reduced);
pos = SwhereInTower(syzHead,tower[level]);
pos2 = SwhereInTower(tmp[0],tower[level-1]);
ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
/* pos is the place to put syzygy at level. */
/* pos2 is the place to put a new GB at level-1. */
Println(ans);
Println(" ");
return(ans);
}
def Sminimal_v(g) {
local r, freeRes, redundantTable, reducer, maxLevel,
minRes, seq, maxSeq, level, betti, q, bases, dr,
betti_levelplus, newbases, i, j,qq;
r = Sschreyer(g);
sm1_pmat(r);
Debug_Sminimal_v = r;
Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
/* Should I turn off the tower?? */
freeRes = r[0];
redundantTable = r[1];
reducer = r[2];
minRes = SnewArrayOfFormat(freeRes);
seq = 0;
maxSeq = SgetMaxSeq(redundantTable);
maxLevel = Length(freeRes);
for (level = 0; level < maxLevel; level++) {
minRes[level] = freeRes[level];
}
for (level = 0; level < maxLevel; level++) {
betti = Length(freeRes[level]);
for (q = betti-1; q>=0; q--) {
if (redundantTable[level,q] > 0) {
Print("[seq,level,q]="); Println([seq,level,q]);
if (level < maxLevel-1) {
bases = freeRes[level+1];
dr = reducer[level,q];
dr[q] = -1;
newbases = SnewArrayOfFormat(bases);
betti_levelplus = Length(bases);
/*
bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
*/
for (i=0; i<betti_levelplus; i++) {
newbases[i] = bases[i] + bases[i,q]*dr;
}
Println(["level, q =", level,q]);
Println("bases="); sm1_pmat(bases);
Println("dr="); sm1_pmat(dr);
Println("newbases="); sm1_pmat(newbases);
minRes[level+1] = newbases;
freeRes = minRes;
#ifdef DEBUG
/* Do it later.
for (qq=0; qq<betti; qq++) {
for (i=0; i<betti_levelplus; i++) {
if (!IsZero(newbases[i,qq])) {
Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
Print("redundantTable ="); sm1_pmat(redundantTable[level]);
Error("Stop in Sminimal for debugging.");
}
}
}
*/
#endif
}
}
}
}
return([Stetris(minRes,redundantTable),
[ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
/* r[4] is the redundantTable_ordinary */
/* r[0] is the freeResolution */
}
/* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */