Annotation of OpenXM/src/k097/lib/minimal/minimal.k, Revision 1.1
1.1 ! takayama 1: /* $OpenXM$ */
! 2: #define DEBUG 1
! 3: /* #define ORDINARY 1 */
! 4: /* Test sequences.
! 5: Use load["minimal.k"];;
! 6:
! 7: a=Sminimal(v);
! 8: b=a[0];
! 9: b[1]*b[0]:
! 10: b[2]*b[1]:
! 11:
! 12: a = test0();
! 13: b = a[0];
! 14: b[1]*b[0]:
! 15: b[2]*b[1]:
! 16: a = Sminimal(b[0]);
! 17:
! 18: a = test1();
! 19: b=a[0];
! 20: b[1]*b[0]:
! 21: b[2]*b[1]:
! 22:
! 23: */
! 24:
! 25:
! 26: load("cohom.k");
! 27: def load_tower() {
! 28: if (Boundp("k0-tower.sm1.loaded")) {
! 29: }else{
! 30: sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
! 31: sm1(" /k0-tower.sm1.loaded 1 def ");
! 32: }
! 33: }
! 34: load_tower();
! 35: SonAutoReduce = true;
! 36: def Factor(f) {
! 37: sm1(f, " fctr /FunctionValue set");
! 38: }
! 39: def Reverse(f) {
! 40: sm1(f," reverse /FunctionValue set");
! 41: }
! 42: def Sgroebner(f) {
! 43: sm1(" [f] groebner /FunctionValue set");
! 44: }
! 45: def test0() {
! 46: local f;
! 47: Sweyl("x,y,z");
! 48: 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,
! 49: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
! 50: frame=SresolutionFrame(f);
! 51: Println(frame);
! 52: /* return(frame); */
! 53: return(SlaScala(f));
! 54: }
! 55: def test1() {
! 56: local f;
! 57: Sweyl("x,y,z");
! 58: 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,
! 59: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
! 60: return(Sminimal(f));
! 61: }
! 62:
! 63:
! 64:
! 65: def Sweyl(v,w) {
! 66: /* extern WeightOfSweyl ; */
! 67: local ww,i,n;
! 68: if(Length(Arglist) == 1) {
! 69: sm1(" [v s_ring_of_differential_operators 0 [(schreyer) 1]] define_ring ");
! 70: sm1(" define_ring_variables ");
! 71:
! 72: sm1(" [ v to_records pop ] /ww set ");
! 73: n = Length(ww);
! 74: WeightOfSweyl = NewArray(n*4);
! 75: for (i=0; i< n; i++) {
! 76: WeightOfSweyl[2*i] = ww[i];
! 77: WeightOfSweyl[2*i+1] = 1;
! 78: }
! 79: for (i=0; i< n; i++) {
! 80: WeightOfSweyl[2*n+2*i] = AddString(["D",ww[i]]);
! 81: WeightOfSweyl[2*n+2*i+1] = 1;
! 82: }
! 83:
! 84: }else{
! 85: sm1(" [v s_ring_of_differential_operators w s_weight_vector 0 [(schreyer) 1]] define_ring ");
! 86: sm1(" define_ring_variables ");
! 87: WeightOfSweyl = w[0];
! 88: }
! 89: }
! 90:
! 91:
! 92: def Spoly(f) {
! 93: sm1(f, " toString tparse /FunctionValue set ");
! 94: }
! 95:
! 96: def SreplaceZeroByZeroPoly(f) {
! 97: if (IsArray(f)) {
! 98: return(Map(f,"SreplaceZeroByZeroPoly"));
! 99: }else{
! 100: if (IsInteger(f)) {
! 101: return(Poly(ToString(f)));
! 102: }else{
! 103: return(f);
! 104: }
! 105: }
! 106: }
! 107: def Shomogenize(f) {
! 108: f = SreplaceZeroByZeroPoly(f);
! 109: if (IsArray(f)) {
! 110: sm1(f," sHomogenize2 /FunctionValue set ");
! 111: /* sm1(f," {sHomogenize2} map /FunctionValue set "); */
! 112: /* Is it correct? Double check.*/
! 113: }else{
! 114: sm1(f, " sHomogenize /FunctionValue set ");
! 115: }
! 116: }
! 117:
! 118: def StoTower() {
! 119: sm1(" [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv ");
! 120: }
! 121:
! 122: def SsetTower(tower) {
! 123: sm1(" [(AvoidTheSameRing)] pushEnv
! 124: [ [(AvoidTheSameRing) 0] system_variable
! 125: [(gbListTower) tower (list) dc] system_variable
! 126: ] pop popEnv ");
! 127: }
! 128:
! 129: def SresolutionFrameWithTower(g,opt) {
! 130: local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
! 131: gbasis;
! 132: if (Length(Arglist) >= 2) {
! 133: if (IsInteger(opt)) count = opt;
! 134: }else{
! 135: count = -1;
! 136: }
! 137:
! 138: sm1(" setupEnvForResolution ");
! 139: /* If I do not put this macro, homogenization
! 140: make a strange behavior. For example,
! 141: [(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
! 142: [(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
! 143: 4/19, 2000.
! 144: */
! 145:
! 146: sm1(" (mmLarger) (matrix) switch_function ");
! 147: g = Map(g,"Shomogenize");
! 148: if (SonAutoReduce) {
! 149: sm1("[ (AutoReduce) ] system_variable /autof set ");
! 150: sm1("[ (AutoReduce) 1 ] system_variable ");
! 151: }
! 152: gbasis = Sgroebner(g);
! 153: g = gbasis[0];
! 154: if (SonAutoReduce) {
! 155: sm1("[ (AutoReduce) autof] system_variable ");
! 156: }
! 157:
! 158: g = Init(g);
! 159:
! 160: /* sm1(" setupEnvForResolution-sugar "); */
! 161: /* -sugar is fine? */
! 162: sm1(" setupEnvForResolution ");
! 163:
! 164: Println(g);
! 165: startingGB = g;
! 166: /* ans = [ SzeroMap(g) ]; It has not been implemented. see resol1.withZeroMap */
! 167: ans = [ ];
! 168: gbTower = [ ];
! 169: skelton = [ ];
! 170: while (true) {
! 171: /* sm1(g," res0Frame /ff set "); */
! 172: withSkel = Sres0FrameWithSkelton(g);
! 173: ff = withSkel[0];
! 174: ans = Append(ans, ff[0]);
! 175: gbTower = Join([ ff[1] ], gbTower);
! 176: skelton = Join([ withSkel[1] ], skelton);
! 177: g = ff[0];
! 178: if (Length(g) == 0) break;
! 179: SsetTower( gbTower );
! 180: if (count == 0) break;
! 181: count = count - 1;
! 182: }
! 183: return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
! 184: }
! 185: HelpAdd(["SresolutionFrameWithTower",
! 186: ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
! 187: "Example: Sweyl(\"x,y\");",
! 188: " a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
! 189:
! 190: def SresolutionFrame(f,opt) {
! 191: local ans;
! 192: ans = SresolutionFrameWithTower(f);
! 193: return(ans[0]);
! 194: }
! 195: /* ---------------------------- */
! 196: def ToGradedPolySet(g) {
! 197: sm1(g," (gradedPolySet) dc /FunctionValue set ");
! 198: }
! 199:
! 200: def NewPolynomialVector(size) {
! 201: sm1(size," (integer) dc newPolyVector /FunctionValue set ");
! 202: }
! 203:
! 204: def SturnOffHomogenization() {
! 205: sm1("
! 206: [(Homogenize)] system_variable 1 eq
! 207: { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
! 208: [(Homogenize) 0] system_variable
! 209: [(ReduceLowerTerms) 0] system_variable
! 210: } { } ifelse
! 211: ");
! 212: }
! 213: def SturnOnHomogenization() {
! 214: sm1("
! 215: [(Homogenize)] system_variable 0 eq
! 216: { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
! 217: [(Homogenize) 1] system_variable
! 218: [(ReduceLowerTerms) 1] system_variable
! 219: } { } ifelse
! 220: ");
! 221: }
! 222:
! 223: def SschreyerSkelton(g) {
! 224: sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
! 225: }
! 226: def Stoes(g) {
! 227: if (IsArray(g)) {
! 228: sm1(g," {toes} map /FunctionValue set ");
! 229: }else{
! 230: sm1(g," toes /FunctionValue set ");
! 231: }
! 232: }
! 233: def Stoes_vec(g) {
! 234: sm1(g," toes /FunctionValue set ");
! 235: }
! 236:
! 237: def Sres0Frame(g) {
! 238: local ans;
! 239: ans = Sres0FrameWithSkelton(g);
! 240: return(ans[0]);
! 241: }
! 242: def Sres0FrameWithSkelton(g) {
! 243: local t_syz, nexttower, m, t_gb, skel, betti,
! 244: gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;
! 245:
! 246: SturnOffHomogenization();
! 247:
! 248: g = Stoes(g);
! 249: skel = SschreyerSkelton(g);
! 250: /* Print("Skelton is ");
! 251: sm1_pmat(skel); */
! 252: betti = Length(skel);
! 253:
! 254: gLength = Length(g);
! 255: grG = ToGradedPolySet(g);
! 256: syzAll = NewPolynomialVector(betti);
! 257: for (k=0; k<betti; k++) {
! 258: pair = skel[k];
! 259: i = pair[0,0];
! 260: j = pair[0,1];
! 261: si = pair[1,0];
! 262: sj = pair[1,1];
! 263: /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
! 264: Print(".");
! 265:
! 266: t_syz = NewPolynomialVector(gLength);
! 267: t_syz[i] = si;
! 268: t_syz[j] = sj;
! 269: syzAll[k] = t_syz;
! 270: }
! 271: t_syz = syzAll;
! 272: Print("Done. betti="); Println(betti);
! 273: /* Println(g); g is in a format such as
! 274: [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
! 275: [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
! 276: [y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
! 277: [3*es^3*y*Dy-es^5*x ]
! 278: */
! 279: nexttower = Init(g);
! 280: SturnOnHomogenization();
! 281: return([[t_syz, nexttower],skel]);
! 282: }
! 283:
! 284:
! 285: def StotalDegree(f) {
! 286: sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");
! 287: }
! 288:
! 289: /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
! 290: def Sord_w(f,w) {
! 291: local neww,i,n;
! 292: n = Length(w);
! 293: neww = NewArray(n);
! 294: for (i=0; i<n; i=i+2) {
! 295: neww[i] = ToString(w[i]);
! 296: }
! 297: for (i=1; i<n; i=i+2) {
! 298: neww[i] = IntegerToSm1Integer(w[i]);
! 299: }
! 300: sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
! 301: }
! 302:
! 303:
! 304: /* This is not satisfactory. */
! 305: def SinitOfArray(f) {
! 306: local p,pos,top;
! 307: if (IsArray(f)) {
! 308: sm1(f," toes init /p set ");
! 309: sm1(p," (es). degree (universalNumber) dc /pos set ");
! 310: return([Init(f[pos]),pos]);
! 311: } else {
! 312: return(Init(f));
! 313: }
! 314: }
! 315:
! 316: def test_SinitOfArray() {
! 317: local f, frame,p,tower,i,j,k;
! 318: Sweyl("x,y,z");
! 319: 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,
! 320: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
! 321: p=SresolutionFrameWithTower(f);
! 322: sm1_pmat(p);
! 323: sm1_pmat(SgenerateTable(p[1]));
! 324: return(p);
! 325: frame = p[0];
! 326: sm1_pmat(p[1]);
! 327: sm1_pmat(frame);
! 328: sm1_pmat(Map(frame[0],"SinitOfArray"));
! 329: sm1_pmat(Map(frame[1],"SinitOfArray"));
! 330: return(p);
! 331: }
! 332:
! 333: /* f is assumed to be a monomial with toes. */
! 334: def Sdegree(f,tower,level) {
! 335: local i;
! 336: if (level <= 1) return(StotalDegree(f));
! 337: i = Degree(f,es);
! 338: return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
! 339: }
! 340:
! 341: def SgenerateTable(tower) {
! 342: local height, n,i,j, ans, ans_at_each_floor;
! 343: height = Length(tower);
! 344: ans = NewArray(height);
! 345: for (i=0; i<height; i++) {
! 346: n = Length(tower[i]);
! 347: ans_at_each_floor=NewArray(n);
! 348: for (j=0; j<n; j++) {
! 349: ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1);
! 350: /* Println([i,j,ans_at_each_floor[j]]); */
! 351: }
! 352: ans[i] = ans_at_each_floor;
! 353: }
! 354: return(ans);
! 355: }
! 356: Sweyl("x,y,z");
! 357: v=[[2*x*Dx + 3*y*Dy+6, 0],
! 358: [3*x^2*Dy + 2*y*Dx, 0],
! 359: [0, x^2+y^2],
! 360: [0, x*y]];
! 361: /* SresolutionFrameWithTower(v); */
! 362:
! 363: def SnewArrayOfFormat(p) {
! 364: if (IsArray(p)) {
! 365: return(Map(p,"SnewArrayOfFormat"));
! 366: }else{
! 367: return(null);
! 368: }
! 369: }
! 370: def SminOfStrategy(a) {
! 371: local n,i,ans,tt;
! 372: ans = 100000; /* very big number */
! 373: if (IsArray(a)) {
! 374: n = Length(a);
! 375: for (i=0; i<n; i++) {
! 376: if (IsArray(a[i])) {
! 377: tt = SminOfStrategy(a[i]);
! 378: if (tt < ans) ans = tt;
! 379: }else{
! 380: if (a[i] < ans) ans = a[i];
! 381: }
! 382: }
! 383: }else{
! 384: if (a < ans) ans = a;
! 385: }
! 386: return(ans);
! 387: }
! 388: def SmaxOfStrategy(a) {
! 389: local n,i,ans,tt;
! 390: ans = -100000; /* very small number */
! 391: if (IsArray(a)) {
! 392: n = Length(a);
! 393: for (i=0; i<n; i++) {
! 394: if (IsArray(a[i])) {
! 395: tt = SmaxOfStrategy(a[i]);
! 396: if (tt > ans) ans = tt;
! 397: }else{
! 398: if (a[i] > ans) ans = a[i];
! 399: }
! 400: }
! 401: }else{
! 402: if (a > ans) ans = a;
! 403: }
! 404: return(ans);
! 405: }
! 406:
! 407:
! 408: def SlaScala(g) {
! 409: local rf, tower, reductionTable, skel, redundantTable, bases,
! 410: strategy, maxOfStrategy, height, level, n, i,
! 411: freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
! 412: redundantTable_ordinary, redundant_seq_ordinary;
! 413: /* extern WeightOfSweyl; */
! 414: ww = WeightOfSweyl;
! 415: Print("WeghtOfSweyl="); Println(WeightOfSweyl);
! 416: rf = SresolutionFrameWithTower(g);
! 417: redundant_seq = 1; redundant_seq_ordinary = 1;
! 418: tower = rf[1];
! 419: reductionTable = SgenerateTable(tower);
! 420: skel = rf[2];
! 421: redundantTable = SnewArrayOfFormat(rf[1]);
! 422: redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
! 423: reducer = SnewArrayOfFormat(rf[1]);
! 424: freeRes = SnewArrayOfFormat(rf[1]);
! 425: bettiTable = SsetBettiTable(rf[1],g);
! 426:
! 427: strategy = SminOfStrategy( reductionTable );
! 428: maxOfStrategy = SmaxOfStrategy( reductionTable );
! 429: height = Length(reductionTable);
! 430: while (strategy <= maxOfStrategy) {
! 431: for (level = 0; level < height; level++) {
! 432: n = Length(reductionTable[level]);
! 433: for (i=0; i<n; i++) {
! 434: if (reductionTable[level,i] == strategy) {
! 435: Print("Processing "); Print([level,i]);
! 436: Print(" Strategy = "); Println(strategy);
! 437: if (level == 0) {
! 438: if (IsNull(redundantTable[level,i])) {
! 439: bases = freeRes[level];
! 440: /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
! 441: pos = SwhereInGB(tower[0,i],rf[3,0]);
! 442: bases[i] = rf[3,0,pos];
! 443: redundantTable[level,i] = 0;
! 444: redundantTable_ordinary[level,i] = 0;
! 445: freeRes[level] = bases;
! 446: /* Println(["GB=",i,bases,tower[0,i]]); */
! 447: }
! 448: }else{ /* level >= 1 */
! 449: if (IsNull(redundantTable[level,i])) {
! 450: bases = freeRes[level];
! 451: f = SpairAndReduction(skel,level,i,freeRes,tower,ww);
! 452: if (f[0] != Poly("0")) {
! 453: place = f[3];
! 454: /* (level-1, place) is the place for f[0],
! 455: which is a newly obtained GB. */
! 456: #ifdef ORDINARY
! 457: redundantTable[level-1,place] = redundant_seq;
! 458: redundant_seq++;
! 459: #else
! 460: if (f[4] > f[5]) {
! 461: /* Zero in the gr-module */
! 462: Print("v-degree of [org,remainder] = ");
! 463: Println([f[4],f[5]]);
! 464: Print("[level,i] = "); Println([level,i]);
! 465: redundantTable[level-1,place] = 0;
! 466: }else{
! 467: redundantTable[level-1,place] = redundant_seq;
! 468: redundant_seq++;
! 469: }
! 470: #endif
! 471: redundantTable_ordinary[level-1,place]
! 472: =redundant_seq_ordinary;
! 473: redundant_seq_ordinary++;
! 474: bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */
! 475: redundantTable[level,i] = 0;
! 476: redundantTable_ordinary[level,i] = 0;
! 477: /* i must be equal to f[2], I think. Double check. */
! 478: freeRes[level] = bases;
! 479: bases = freeRes[level-1];
! 480: bases[place] = f[0];
! 481: freeRes[level-1] = bases;
! 482: reducer[level-1,place] = f[1];
! 483: }else{
! 484: redundantTable[level,i] = 0;
! 485: bases = freeRes[level];
! 486: bases[i] = f[1]; /* Put the syzygy. */
! 487: freeRes[level] = bases;
! 488: }
! 489: }
! 490: } /* end of level >= 1 */
! 491: }
! 492: }
! 493: }
! 494: strategy++;
! 495: }
! 496: n = Length(freeRes);
! 497: freeResV = SnewArrayOfFormat(freeRes);
! 498: for (i=0; i<n; i++) {
! 499: bases = freeRes[i];
! 500: bases = Sbases_to_vec(bases,bettiTable[i]);
! 501: freeResV[i] = bases;
! 502: }
! 503: return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
! 504: }
! 505:
! 506: def SsetBettiTable(freeRes,g) {
! 507: local level,i, n,bases,ans;
! 508: ans = NewArray(Length(freeRes)+1);
! 509: n = Length(freeRes);
! 510: if (IsArray(g[0])) {
! 511: ans[0] = Length(g[0]);
! 512: }else{
! 513: ans[0] = 1;
! 514: }
! 515: for (level=0; level<n; level++) {
! 516: bases = freeRes[level];
! 517: if (IsArray(bases)) {
! 518: ans[level+1] = Length(bases);
! 519: }else{
! 520: ans[level+1] = 1;
! 521: }
! 522: }
! 523: return(ans);
! 524: }
! 525:
! 526: def SwhereInGB(f,tower) {
! 527: local i,n,p,q;
! 528: n = Length(tower);
! 529: for (i=0; i<n; i++) {
! 530: p = MonomialPart(tower[i]);
! 531: q = MonomialPart(f);
! 532: if (p == q) return(i);
! 533: }
! 534: Println([f,tower]);
! 535: Error("whereInGB : [f,myset]: f could not be found in the myset.");
! 536: }
! 537: def SunitOfFormat(pos,forms) {
! 538: local ans,i,n;
! 539: n = Length(forms);
! 540: ans = NewArray(n);
! 541: for (i=0; i<n; i++) {
! 542: if (i != pos) {
! 543: ans[i] = Poly("0");
! 544: }else{
! 545: ans[i] = Poly("1");
! 546: }
! 547: }
! 548: return(ans);
! 549: }
! 550:
! 551: def Error(s) {
! 552: sm1(" s error ");
! 553: }
! 554:
! 555: def IsNull(s) {
! 556: if (Stag(s) == 0) return(true);
! 557: else return(false);
! 558: }
! 559:
! 560: def StowerOf(tower,level) {
! 561: local ans,i;
! 562: ans = [ ];
! 563: if (level == 0) return([[]]);
! 564: for (i=0; i<level; i++) {
! 565: ans = Append(ans,tower[i]);
! 566: }
! 567: return(Reverse(ans));
! 568: }
! 569:
! 570: def Sspolynomial(f,g) {
! 571: if (IsArray(f)) {
! 572: f = Stoes_vec(f);
! 573: }
! 574: if (IsArray(g)) {
! 575: g = Stoes_vec(g);
! 576: }
! 577: sm1("f g spol /FunctionValue set");
! 578: }
! 579:
! 580: def MonomialPart(f) {
! 581: sm1(" [(lmonom) f] gbext /FunctionValue set ");
! 582: }
! 583:
! 584: def SwhereInTower(f,tower) {
! 585: local i,n,p,q;
! 586: if (f == Poly("0")) return(-1);
! 587: n = Length(tower);
! 588: for (i=0; i<n; i++) {
! 589: p = MonomialPart(tower[i]);
! 590: q = MonomialPart(f);
! 591: if (p == q) return(i);
! 592: }
! 593: Println([f,tower]);
! 594: Error("[f,tower]: f could not be found in the tower.");
! 595: }
! 596:
! 597: def Stag(f) {
! 598: sm1(f," tag (universalNumber) dc /FunctionValue set");
! 599: }
! 600:
! 601: def SpairAndReduction(skel,level,ii,freeRes,tower,ww) {
! 602: local i, j, myindex, p, bases, tower2, gi, gj,
! 603: si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
! 604: vdeg,vdeg_reduced;
! 605: Println("SpairAndReduction:");
! 606:
! 607: if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
! 608: p = skel[level,ii];
! 609: myindex = p[0];
! 610: i = myindex[0]; j = myindex[1];
! 611: bases = freeRes[level-1];
! 612: Println(["p and bases ",p,bases]);
! 613: if (IsNull(bases[i]) || IsNull(bases[j])) {
! 614: Println([level,i,j,bases[i],bases[j]]);
! 615: Error("level, i, j : bases[i], bases[j] must not be NULL.");
! 616: }
! 617:
! 618: tower2 = StowerOf(tower,level-1);
! 619: SsetTower(tower2);
! 620: /** sm1(" show_ring "); */
! 621:
! 622: gi = Stoes_vec(bases[i]);
! 623: gj = Stoes_vec(bases[j]);
! 624:
! 625: ssp = Sspolynomial(gi,gj);
! 626: si = ssp[0,0];
! 627: sj = ssp[0,1];
! 628: syzHead = si*es^i;
! 629: /* This will be the head term, I think. But, double check. */
! 630: Println([si*es^i,sj*es^j]);
! 631:
! 632: Print("[gi, gj] = "); Println([gi,gj]);
! 633: sm1(" [(Homogenize)] system_variable message ");
! 634: Print("Reduce the element "); Println(si*gi+sj*gj);
! 635: Print("by "); Println(bases);
! 636:
! 637: tmp = Sreduction(si*gi+sj*gj, bases);
! 638:
! 639: Print("result is "); Println(tmp);
! 640:
! 641: vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww);
! 642: vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
! 643: Print("vdegree of the original = "); Println(vdeg);
! 644: Print("vdegree of the remainder = "); Println(vdeg_reduced);
! 645:
! 646: t_syz = tmp[2];
! 647: si = si*tmp[1]+t_syz[i];
! 648: sj = sj*tmp[1]+t_syz[j];
! 649: t_syz[i] = si;
! 650: t_syz[j] = sj;
! 651: pos = SwhereInTower(syzHead,tower[level]);
! 652: pos2 = SwhereInTower(tmp[0],tower[level-1]);
! 653: ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
! 654: /* pos is the place to put syzygy at level. */
! 655: /* pos2 is the place to put a new GB at level-1. */
! 656: Println(ans);
! 657: return(ans);
! 658: }
! 659:
! 660: def Sreduction(f,myset) {
! 661: local n, indexTable, set2, i, j, tmp, t_syz;
! 662: n = Length(myset);
! 663: indexTable = NewArray(n);
! 664: set2 = [ ];
! 665: j = 0;
! 666: for (i=0; i<n; i++) {
! 667: if (IsNull(myset[i])) {
! 668: indexTable[i] = -1;
! 669: /* }else if (myset[i] == Poly("0")) {
! 670: indexTable[i] = -1; */
! 671: }else{
! 672: set2 = Append(set2,Stoes_vec(myset[i]));
! 673: indexTable[i] = j;
! 674: j++;
! 675: }
! 676: }
! 677: sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
! 678: t_syz = NewArray(n);
! 679: for (i=0; i<n; i++) {
! 680: if (indexTable[i] != -1) {
! 681: t_syz[i] = tmp[2, indexTable[i]];
! 682: }else{
! 683: t_syz[i] = Poly("0");
! 684: }
! 685: }
! 686: return([tmp[0],tmp[1],t_syz]);
! 687: }
! 688:
! 689: def Warning(s) {
! 690: Print("Warning: ");
! 691: Println(s);
! 692: }
! 693: def RingOf(f) {
! 694: local r;
! 695: if (IsPolynomial(f)) {
! 696: if (f != Poly("0")) {
! 697: sm1(f," (ring) dc /r set ");
! 698: }else{
! 699: sm1(" [(CurrentRingp)] system_variable /r set ");
! 700: }
! 701: }else{
! 702: Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
! 703: sm1(" [(CurrentRingp)] system_variable /r set ");
! 704: }
! 705: return(r);
! 706: }
! 707:
! 708: def Sfrom_es(f,size) {
! 709: local c,ans, i, d, myes, myee, j,n,r,ans2;
! 710: if (Length(Arglist) < 2) size = -1;
! 711: if (IsArray(f)) return(f);
! 712: r = RingOf(f);
! 713: myes = PolyR("es",r);
! 714: myee = PolyR("e_",r);
! 715: if (Degree(f,myee) > 0 && size == -1) {
! 716: if (size == -1) {
! 717: sm1(f," (array) dc /ans set");
! 718: return(ans);
! 719: }
! 720: }
! 721:
! 722: /*
! 723: Coefficients(x^2-1,x):
! 724: [ [ 2 , 0 ] , [ 1 , -1 ] ]
! 725: */
! 726: if (Degree(f,myee) > 0) {
! 727: c = Coefficients(f,myee);
! 728: }else{
! 729: c = Coefficients(f,myes);
! 730: }
! 731: if (size < 0) {
! 732: size = c[0,0]+1;
! 733: }
! 734: ans = NewArray(size);
! 735: for (i=0; i<size; i++) {ans[i] = 0;}
! 736: n = Length(c[0]);
! 737: for (j=0; j<n; j++) {
! 738: d = c[0,j];
! 739: ans[d] = c[1,j];
! 740: }
! 741: return(ans);
! 742: }
! 743:
! 744: def Sbases_to_vec(bases,size) {
! 745: local n, giveSize, newbases,i;
! 746: /* bases = [1+es*x, [1,2,3*x]] */
! 747: if (Length(Arglist) > 1) {
! 748: giveSize = true;
! 749: }else{
! 750: giveSize = false;
! 751: }
! 752: n = Length(bases);
! 753: newbases = NewArray(n);
! 754: for (i=0; i<n; i++) {
! 755: if (giveSize) {
! 756: newbases[i] = Sfrom_es(bases[i], size);
! 757: }else{
! 758: newbases[i] = Sfrom_es(bases[i]);
! 759: }
! 760: }
! 761: return(newbases);
! 762: }
! 763:
! 764: def Sminimal(g) {
! 765: local r, freeRes, redundantTable, reducer, maxLevel,
! 766: minRes, seq, maxSeq, level, betti, q, bases, dr,
! 767: betti_levelplus, newbases, i, j,qq;
! 768: r = SlaScala(g);
! 769: /* Should I turn off the tower?? */
! 770: freeRes = r[0];
! 771: redundantTable = r[1];
! 772: reducer = r[2];
! 773: minRes = SnewArrayOfFormat(freeRes);
! 774: seq = 0;
! 775: maxSeq = SgetMaxSeq(redundantTable);
! 776: maxLevel = Length(freeRes);
! 777: for (level = 0; level < maxLevel; level++) {
! 778: minRes[level] = freeRes[level];
! 779: }
! 780: seq=maxSeq+1;
! 781: while (seq > 1) {
! 782: seq--;
! 783: for (level = 0; level < maxLevel; level++) {
! 784: betti = Length(freeRes[level]);
! 785: for (q = 0; q<betti; q++) {
! 786: if (redundantTable[level,q] == seq) {
! 787: Print("[seq,level,q]="); Println([seq,level,q]);
! 788: if (level < maxLevel-1) {
! 789: bases = freeRes[level+1];
! 790: dr = reducer[level,q];
! 791: dr[q] = -1;
! 792: newbases = SnewArrayOfFormat(bases);
! 793: betti_levelplus = Length(bases);
! 794: /*
! 795: bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
! 796: */
! 797: for (i=0; i<betti_levelplus; i++) {
! 798: newbases[i] = bases[i] + bases[i,q]*dr;
! 799: }
! 800: Println(["level, q =", level,q]);
! 801: Println("bases="); sm1_pmat(bases);
! 802: Println("dr="); sm1_pmat(dr);
! 803: Println("newbases="); sm1_pmat(newbases);
! 804: minRes[level+1] = newbases;
! 805: freeRes = minRes;
! 806: #ifdef DEBUG
! 807: for (qq=0; qq<betti; qq++) {
! 808: if ((redundantTable[level,qq] >= seq) &&
! 809: (redundantTable[level,qq] <= maxSeq)) {
! 810: for (i=0; i<betti_levelplus; i++) {
! 811: if (!IsZero(newbases[i,qq])) {
! 812: Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
! 813: Print("redundantTable ="); sm1_pmat(redundantTable[level]);
! 814: Error("Stop in Sminimal for debugging.");
! 815: }
! 816: }
! 817: }
! 818: }
! 819: #endif
! 820: }
! 821: }
! 822: }
! 823: }
! 824: }
! 825: return([Stetris(minRes,redundantTable),
! 826: [ minRes, redundantTable, reducer,r[3],r[4]]]);
! 827: /* r[4] is the redundantTable_ordinary */
! 828: }
! 829:
! 830:
! 831: def IsZero(f) {
! 832: if (IsPolynomial(f)) {
! 833: return( f == Poly("0"));
! 834: }else if (IsInteger(f)) {
! 835: return( f == 0);
! 836: }else if (IsSm1Integer(f)) {
! 837: return( f == true );
! 838: }else if (IsDouble(f)) {
! 839: return( f == 0.0 );
! 840: }else if (IsRational(f)) {
! 841: return(IsZero(Denominator(f)));
! 842: }else{
! 843: Error("IsZero: cannot deal with this data type.");
! 844: }
! 845: }
! 846: def SgetMaxSeq(redundantTable) {
! 847: local level,i,n,ans, levelMax,bases;
! 848: levelMax = Length( redundantTable );
! 849: ans = 0;
! 850: for (level = 0; level < levelMax; level++) {
! 851: bases = redundantTable[level];
! 852: n = Length(bases);
! 853: for (i=0; i<n; i++) {
! 854: if (IsInteger( bases[i] )) {
! 855: if (bases[i] > ans) {
! 856: ans = bases[i];
! 857: }
! 858: }
! 859: }
! 860: }
! 861: return(ans);
! 862: }
! 863:
! 864: def Stetris(freeRes,redundantTable) {
! 865: local level, i, j, resLength, minRes,
! 866: bases, newbases, newbases2;
! 867: minRes = SnewArrayOfFormat(freeRes);
! 868: resLength = Length( freeRes );
! 869: for (level=0; level<resLength; level++) {
! 870: bases = freeRes[level];
! 871: newbases = SnewArrayOfFormat(bases);
! 872: betti = Length(bases); j = 0;
! 873: /* Delete rows */
! 874: for (i=0; i<betti; i++) {
! 875: if (redundantTable[level,i] < 1) {
! 876: newbases[j] = bases[i];
! 877: j++;
! 878: }
! 879: }
! 880: bases = SfirstN(newbases,j);
! 881: if (level > 0) {
! 882: /* Delete columns */
! 883: newbases = Transpose(bases);
! 884: betti = Length(newbases); j = 0;
! 885: newbases2 = SnewArrayOfFormat(newbases);
! 886: for (i=0; i<betti; i++) {
! 887: if (redundantTable[level-1,i] < 1) {
! 888: newbases2[j] = newbases[i];
! 889: j++;
! 890: }
! 891: }
! 892: newbases = Transpose(SfirstN(newbases2,j));
! 893: }else{
! 894: newbases = bases;
! 895: }
! 896: Println(["level=", level]);
! 897: sm1_pmat(bases);
! 898: sm1_pmat(newbases);
! 899:
! 900: minRes[level] = newbases;
! 901: }
! 902: return(minRes);
! 903: }
! 904:
! 905: def SfirstN(bases,k) {
! 906: local ans,i;
! 907: ans = NewArray(k);
! 908: for (i=0; i<k; i++) {
! 909: ans[i] = bases[i];
! 910: }
! 911: return(ans);
! 912: }
! 913:
! 914:
! 915: /* usage: tt is tower. ww is weight.
! 916: a = SresolutionFrameWithTower(v);
! 917: tt = a[1];
! 918: ww = [x,1,y,1,Dx,1,Dy,1];
! 919: SvDegree(x*es,tt,1,ww):
! 920:
! 921: In(17)=tt:
! 922: [[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 ] ,
! 923: [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ] ,
! 924: [3*es^3*y*Dy ] ]
! 925: In(18)=SvDegree(x*es,tt,1,ww):
! 926: 3
! 927: In(19)=SvDegree(x*es^3,tt,1,ww):
! 928: 4
! 929: In(20)=SvDegree(x,tt,2,ww):
! 930: 4
! 931:
! 932: */
! 933: def SvDegree(f,tower,level,w) {
! 934: local i,ans;
! 935: if (IsZero(f)) return(null);
! 936: if (level <= 0) {
! 937: return(Sord_w(f,w));
! 938: }
! 939: i = Degree(f,es);
! 940: ans = Sord_w(f,w) +
! 941: SvDegree(tower[level-1,i],tower,level-1,w);
! 942: return(ans);
! 943: }
! 944:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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