Annotation of OpenXM/src/k097/lib/minimal/minimal.k, Revision 1.16
1.16 ! takayama 1: /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.15 2000/06/14 07:44:05 takayama Exp $ */
1.1 takayama 2: #define DEBUG 1
3: /* #define ORDINARY 1 */
1.4 takayama 4: /* If you run this program on openxm version 1.1.2 (FreeBSD),
5: make a symbolic link by the command
6: ln -s /usr/bin/cpp /lib/cpp
7: */
1.6 takayama 8: #define OFFSET 0
1.14 takayama 9: #define TOTAL_STRATEGY 1
1.6 takayama 10: /* #define OFFSET 20*/
1.1 takayama 11: /* Test sequences.
12: Use load["minimal.k"];;
13:
14: a=Sminimal(v);
15: b=a[0];
16: b[1]*b[0]:
17: b[2]*b[1]:
18:
19: a = test0();
20: b = a[0];
21: b[1]*b[0]:
22: b[2]*b[1]:
23: a = Sminimal(b[0]);
24:
25: a = test1();
26: b=a[0];
27: b[1]*b[0]:
28: b[2]*b[1]:
29:
30: */
31:
32:
33: load("cohom.k");
34: def load_tower() {
35: if (Boundp("k0-tower.sm1.loaded")) {
36: }else{
37: sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
38: sm1(" /k0-tower.sm1.loaded 1 def ");
39: }
1.7 takayama 40: sm1(" oxNoX ");
1.1 takayama 41: }
42: load_tower();
43: SonAutoReduce = true;
44: def Factor(f) {
45: sm1(f, " fctr /FunctionValue set");
46: }
47: def Reverse(f) {
48: sm1(f," reverse /FunctionValue set");
49: }
50: def Sgroebner(f) {
51: sm1(" [f] groebner /FunctionValue set");
52: }
53: def test0() {
54: local f;
55: Sweyl("x,y,z");
56: 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,
57: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
58: frame=SresolutionFrame(f);
59: Println(frame);
60: /* return(frame); */
61: return(SlaScala(f));
62: }
63: def test1() {
64: local f;
65: Sweyl("x,y,z");
66: 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,
67: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
68: return(Sminimal(f));
69: }
70:
71:
72:
73: def Sweyl(v,w) {
74: /* extern WeightOfSweyl ; */
75: local ww,i,n;
76: if(Length(Arglist) == 1) {
77: sm1(" [v s_ring_of_differential_operators 0 [(schreyer) 1]] define_ring ");
78: sm1(" define_ring_variables ");
79:
80: sm1(" [ v to_records pop ] /ww set ");
81: n = Length(ww);
82: WeightOfSweyl = NewArray(n*4);
83: for (i=0; i< n; i++) {
84: WeightOfSweyl[2*i] = ww[i];
85: WeightOfSweyl[2*i+1] = 1;
86: }
87: for (i=0; i< n; i++) {
88: WeightOfSweyl[2*n+2*i] = AddString(["D",ww[i]]);
89: WeightOfSweyl[2*n+2*i+1] = 1;
90: }
91:
92: }else{
93: sm1(" [v s_ring_of_differential_operators w s_weight_vector 0 [(schreyer) 1]] define_ring ");
94: sm1(" define_ring_variables ");
95: WeightOfSweyl = w[0];
96: }
97: }
98:
99:
100: def Spoly(f) {
101: sm1(f, " toString tparse /FunctionValue set ");
102: }
103:
104: def SreplaceZeroByZeroPoly(f) {
105: if (IsArray(f)) {
106: return(Map(f,"SreplaceZeroByZeroPoly"));
107: }else{
108: if (IsInteger(f)) {
109: return(Poly(ToString(f)));
110: }else{
111: return(f);
112: }
113: }
114: }
115: def Shomogenize(f) {
116: f = SreplaceZeroByZeroPoly(f);
117: if (IsArray(f)) {
118: sm1(f," sHomogenize2 /FunctionValue set ");
119: /* sm1(f," {sHomogenize2} map /FunctionValue set "); */
120: /* Is it correct? Double check.*/
121: }else{
122: sm1(f, " sHomogenize /FunctionValue set ");
123: }
124: }
125:
126: def StoTower() {
127: sm1(" [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv ");
128: }
129:
130: def SsetTower(tower) {
131: sm1(" [(AvoidTheSameRing)] pushEnv
132: [ [(AvoidTheSameRing) 0] system_variable
133: [(gbListTower) tower (list) dc] system_variable
134: ] pop popEnv ");
1.14 takayama 135: /* sm1("(hoge) message show_ring "); */
1.1 takayama 136: }
137:
138: def SresolutionFrameWithTower(g,opt) {
139: local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
1.15 takayama 140: gbasis, nohomog;
141: nohomog = false;
142: count = -1;
1.1 takayama 143: if (Length(Arglist) >= 2) {
1.15 takayama 144: if (IsInteger(opt)) {
145: count = opt;
146: }else if (IsString(opt)) {
147: if (opt == "homogenized") {
148: nohomog = true;
149: }else{
150: Println("Warning: unknown option");
151: Println(opt);
152: }
153: }
1.1 takayama 154: }else{
155: count = -1;
156: }
157:
158: sm1(" setupEnvForResolution ");
159: /* If I do not put this macro, homogenization
160: make a strange behavior. For example,
161: [(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
162: [(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
163: 4/19, 2000.
164: */
165:
166: sm1(" (mmLarger) (matrix) switch_function ");
1.15 takayama 167: if (! nohomog) {
168: Println("Automatic homogenization.");
169: g = Map(g,"Shomogenize");
170: }else{
171: Println("No automatic homogenization.");
172: }
1.1 takayama 173: if (SonAutoReduce) {
174: sm1("[ (AutoReduce) ] system_variable /autof set ");
175: sm1("[ (AutoReduce) 1 ] system_variable ");
176: }
177: gbasis = Sgroebner(g);
178: g = gbasis[0];
179: if (SonAutoReduce) {
180: sm1("[ (AutoReduce) autof] system_variable ");
181: }
182:
183: g = Init(g);
184:
185: /* sm1(" setupEnvForResolution-sugar "); */
186: /* -sugar is fine? */
187: sm1(" setupEnvForResolution ");
188:
189: Println(g);
190: startingGB = g;
191: /* ans = [ SzeroMap(g) ]; It has not been implemented. see resol1.withZeroMap */
192: ans = [ ];
193: gbTower = [ ];
194: skelton = [ ];
195: while (true) {
196: /* sm1(g," res0Frame /ff set "); */
197: withSkel = Sres0FrameWithSkelton(g);
198: ff = withSkel[0];
199: ans = Append(ans, ff[0]);
200: gbTower = Join([ ff[1] ], gbTower);
201: skelton = Join([ withSkel[1] ], skelton);
202: g = ff[0];
203: if (Length(g) == 0) break;
204: SsetTower( gbTower );
205: if (count == 0) break;
206: count = count - 1;
207: }
208: return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
209: }
210: HelpAdd(["SresolutionFrameWithTower",
211: ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
1.15 takayama 212: "option: \"homogenized\" (no automatic homogenization) ",
1.1 takayama 213: "Example: Sweyl(\"x,y\");",
214: " a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
215:
216: def SresolutionFrame(f,opt) {
217: local ans;
1.15 takayama 218: ans = SresolutionFrameWithTower(f,opt);
1.1 takayama 219: return(ans[0]);
220: }
221: /* ---------------------------- */
222: def ToGradedPolySet(g) {
223: sm1(g," (gradedPolySet) dc /FunctionValue set ");
224: }
225:
226: def NewPolynomialVector(size) {
227: sm1(size," (integer) dc newPolyVector /FunctionValue set ");
228: }
229:
230: def SturnOffHomogenization() {
231: sm1("
232: [(Homogenize)] system_variable 1 eq
233: { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
234: [(Homogenize) 0] system_variable
235: [(ReduceLowerTerms) 0] system_variable
236: } { } ifelse
237: ");
238: }
239: def SturnOnHomogenization() {
240: sm1("
241: [(Homogenize)] system_variable 0 eq
242: { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
243: [(Homogenize) 1] system_variable
244: [(ReduceLowerTerms) 1] system_variable
245: } { } ifelse
246: ");
247: }
248:
249: def SschreyerSkelton(g) {
250: sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
251: }
252: def Stoes(g) {
253: if (IsArray(g)) {
254: sm1(g," {toes} map /FunctionValue set ");
255: }else{
256: sm1(g," toes /FunctionValue set ");
257: }
258: }
259: def Stoes_vec(g) {
260: sm1(g," toes /FunctionValue set ");
261: }
262:
263: def Sres0Frame(g) {
264: local ans;
265: ans = Sres0FrameWithSkelton(g);
266: return(ans[0]);
267: }
268: def Sres0FrameWithSkelton(g) {
269: local t_syz, nexttower, m, t_gb, skel, betti,
270: gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;
271:
272: SturnOffHomogenization();
273:
274: g = Stoes(g);
275: skel = SschreyerSkelton(g);
276: /* Print("Skelton is ");
277: sm1_pmat(skel); */
278: betti = Length(skel);
279:
280: gLength = Length(g);
281: grG = ToGradedPolySet(g);
282: syzAll = NewPolynomialVector(betti);
283: for (k=0; k<betti; k++) {
284: pair = skel[k];
285: i = pair[0,0];
286: j = pair[0,1];
287: si = pair[1,0];
288: sj = pair[1,1];
289: /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
290: Print(".");
291:
292: t_syz = NewPolynomialVector(gLength);
293: t_syz[i] = si;
294: t_syz[j] = sj;
295: syzAll[k] = t_syz;
296: }
297: t_syz = syzAll;
298: Print("Done. betti="); Println(betti);
299: /* Println(g); g is in a format such as
300: [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
301: [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
302: [y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
303: [3*es^3*y*Dy-es^5*x ]
304: */
305: nexttower = Init(g);
306: SturnOnHomogenization();
307: return([[t_syz, nexttower],skel]);
308: }
309:
310:
311: def StotalDegree(f) {
1.14 takayama 312: local d0;
313: sm1(" [(grade) f] gbext (universalNumber) dc /d0 set ");
314: /* Print("degree of "); Print(f); Print(" is "); Println(d0); */
315: return(d0);
1.1 takayama 316: }
317:
318: /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
319: def Sord_w(f,w) {
320: local neww,i,n;
321: n = Length(w);
322: neww = NewArray(n);
323: for (i=0; i<n; i=i+2) {
324: neww[i] = ToString(w[i]);
325: }
326: for (i=1; i<n; i=i+2) {
327: neww[i] = IntegerToSm1Integer(w[i]);
328: }
329: sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
330: }
331:
332:
333: /* This is not satisfactory. */
334: def SinitOfArray(f) {
335: local p,pos,top;
336: if (IsArray(f)) {
337: sm1(f," toes init /p set ");
338: sm1(p," (es). degree (universalNumber) dc /pos set ");
339: return([Init(f[pos]),pos]);
340: } else {
341: return(Init(f));
342: }
343: }
344:
345: def test_SinitOfArray() {
346: local f, frame,p,tower,i,j,k;
347: Sweyl("x,y,z");
348: 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,
349: -y^2*z^2 + x*z^3 + y*z^3, -z^4];
350: p=SresolutionFrameWithTower(f);
351: sm1_pmat(p);
352: sm1_pmat(SgenerateTable(p[1]));
353: return(p);
354: frame = p[0];
355: sm1_pmat(p[1]);
356: sm1_pmat(frame);
357: sm1_pmat(Map(frame[0],"SinitOfArray"));
358: sm1_pmat(Map(frame[1],"SinitOfArray"));
359: return(p);
360: }
361:
362: /* f is assumed to be a monomial with toes. */
363: def Sdegree(f,tower,level) {
1.6 takayama 364: local i,ww, wd;
365: /* extern WeightOfSweyl; */
366: ww = WeightOfSweyl;
1.5 takayama 367: f = Init(f);
1.1 takayama 368: if (level <= 1) return(StotalDegree(f));
369: i = Degree(f,es);
1.6 takayama 370: #ifdef TOTAL_STRATEGY
371: return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
372: #endif
373: /* Strategy must be compatible with ordering. */
374: /* Weight vector must be non-negative, too. */
375: /* See Sdegree, SgenerateTable, reductionTable. */
376: wd = Sord_w(f,ww);
377: return(wd+Sdegree(tower[level-2,i],tower,level-1));
378:
1.1 takayama 379: }
380:
381: def SgenerateTable(tower) {
382: local height, n,i,j, ans, ans_at_each_floor;
1.16 ! takayama 383:
! 384: /*
! 385: Print("SgenerateTable: tower=");Println(tower);
! 386: sm1(" print_switch_status "); */
1.1 takayama 387: height = Length(tower);
388: ans = NewArray(height);
389: for (i=0; i<height; i++) {
390: n = Length(tower[i]);
391: ans_at_each_floor=NewArray(n);
392: for (j=0; j<n; j++) {
1.6 takayama 393: ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1)
394: + OFFSET;
1.1 takayama 395: /* Println([i,j,ans_at_each_floor[j]]); */
396: }
397: ans[i] = ans_at_each_floor;
398: }
399: return(ans);
400: }
401: Sweyl("x,y,z");
402: v=[[2*x*Dx + 3*y*Dy+6, 0],
403: [3*x^2*Dy + 2*y*Dx, 0],
404: [0, x^2+y^2],
405: [0, x*y]];
406: /* SresolutionFrameWithTower(v); */
407:
408: def SnewArrayOfFormat(p) {
409: if (IsArray(p)) {
410: return(Map(p,"SnewArrayOfFormat"));
411: }else{
412: return(null);
413: }
414: }
1.4 takayama 415: def ScopyArray(a) {
416: local n, i,ans;
417: n = Length(a);
418: ans = NewArray(n);
419: for (i=0; i<n; i++) {
420: ans[i] = a[i];
421: }
422: return(ans);
423: }
1.1 takayama 424: def SminOfStrategy(a) {
425: local n,i,ans,tt;
426: ans = 100000; /* very big number */
427: if (IsArray(a)) {
428: n = Length(a);
429: for (i=0; i<n; i++) {
430: if (IsArray(a[i])) {
431: tt = SminOfStrategy(a[i]);
432: if (tt < ans) ans = tt;
433: }else{
434: if (a[i] < ans) ans = a[i];
435: }
436: }
437: }else{
438: if (a < ans) ans = a;
439: }
440: return(ans);
441: }
442: def SmaxOfStrategy(a) {
443: local n,i,ans,tt;
444: ans = -100000; /* very small number */
445: if (IsArray(a)) {
446: n = Length(a);
447: for (i=0; i<n; i++) {
448: if (IsArray(a[i])) {
449: tt = SmaxOfStrategy(a[i]);
450: if (tt > ans) ans = tt;
451: }else{
452: if (a[i] > ans) ans = a[i];
453: }
454: }
455: }else{
456: if (a > ans) ans = a;
457: }
458: return(ans);
459: }
460:
461:
1.15 takayama 462: def SlaScala(g,opt) {
1.1 takayama 463: local rf, tower, reductionTable, skel, redundantTable, bases,
464: strategy, maxOfStrategy, height, level, n, i,
465: freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
1.4 takayama 466: redundantTable_ordinary, redundant_seq_ordinary,
467: reductionTable_tmp;
1.1 takayama 468: /* extern WeightOfSweyl; */
469: ww = WeightOfSweyl;
1.6 takayama 470: Print("WeightOfSweyl="); Println(WeightOfSweyl);
1.15 takayama 471: rf = SresolutionFrameWithTower(g,opt);
1.14 takayama 472: Print("rf="); sm1_pmat(rf);
1.1 takayama 473: redundant_seq = 1; redundant_seq_ordinary = 1;
474: tower = rf[1];
1.16 ! takayama 475:
! 476: Println("Generating reduction table which gives an order of reduction.");
! 477: Print("WeghtOfSweyl="); Println(WeightOfSweyl);
! 478: Print("tower"); Println(tower);
1.1 takayama 479: reductionTable = SgenerateTable(tower);
1.16 ! takayama 480: Print("reductionTable="); sm1_pmat(reductionTable);
! 481:
1.1 takayama 482: skel = rf[2];
483: redundantTable = SnewArrayOfFormat(rf[1]);
484: redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
485: reducer = SnewArrayOfFormat(rf[1]);
486: freeRes = SnewArrayOfFormat(rf[1]);
487: bettiTable = SsetBettiTable(rf[1],g);
488:
489: strategy = SminOfStrategy( reductionTable );
490: maxOfStrategy = SmaxOfStrategy( reductionTable );
491: height = Length(reductionTable);
492: while (strategy <= maxOfStrategy) {
493: for (level = 0; level < height; level++) {
494: n = Length(reductionTable[level]);
1.4 takayama 495: reductionTable_tmp = ScopyArray(reductionTable[level]);
496: while (SthereIs(reductionTable_tmp,strategy)) {
497: i = SnextI(reductionTable_tmp,strategy,redundantTable,
498: skel,level,freeRes);
499: Println([level,i]);
500: reductionTable_tmp[i] = -200000;
1.1 takayama 501: if (reductionTable[level,i] == strategy) {
1.16 ! takayama 502: Print("Processing [level,i]= "); Print([level,i]);
1.1 takayama 503: Print(" Strategy = "); Println(strategy);
504: if (level == 0) {
505: if (IsNull(redundantTable[level,i])) {
506: bases = freeRes[level];
507: /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
508: pos = SwhereInGB(tower[0,i],rf[3,0]);
509: bases[i] = rf[3,0,pos];
510: redundantTable[level,i] = 0;
511: redundantTable_ordinary[level,i] = 0;
512: freeRes[level] = bases;
513: /* Println(["GB=",i,bases,tower[0,i]]); */
514: }
515: }else{ /* level >= 1 */
516: if (IsNull(redundantTable[level,i])) {
517: bases = freeRes[level];
518: f = SpairAndReduction(skel,level,i,freeRes,tower,ww);
519: if (f[0] != Poly("0")) {
520: place = f[3];
521: /* (level-1, place) is the place for f[0],
522: which is a newly obtained GB. */
523: #ifdef ORDINARY
524: redundantTable[level-1,place] = redundant_seq;
525: redundant_seq++;
526: #else
527: if (f[4] > f[5]) {
528: /* Zero in the gr-module */
529: Print("v-degree of [org,remainder] = ");
530: Println([f[4],f[5]]);
531: Print("[level,i] = "); Println([level,i]);
532: redundantTable[level-1,place] = 0;
533: }else{
534: redundantTable[level-1,place] = redundant_seq;
535: redundant_seq++;
536: }
537: #endif
538: redundantTable_ordinary[level-1,place]
539: =redundant_seq_ordinary;
540: redundant_seq_ordinary++;
541: bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */
542: redundantTable[level,i] = 0;
543: redundantTable_ordinary[level,i] = 0;
544: /* i must be equal to f[2], I think. Double check. */
545: freeRes[level] = bases;
546: bases = freeRes[level-1];
547: bases[place] = f[0];
548: freeRes[level-1] = bases;
549: reducer[level-1,place] = f[1];
550: }else{
551: redundantTable[level,i] = 0;
552: bases = freeRes[level];
553: bases[i] = f[1]; /* Put the syzygy. */
554: freeRes[level] = bases;
555: }
556: }
557: } /* end of level >= 1 */
558: }
559: }
560: }
561: strategy++;
562: }
563: n = Length(freeRes);
564: freeResV = SnewArrayOfFormat(freeRes);
565: for (i=0; i<n; i++) {
566: bases = freeRes[i];
567: bases = Sbases_to_vec(bases,bettiTable[i]);
568: freeResV[i] = bases;
569: }
570: return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
571: }
1.4 takayama 572:
573: def SthereIs(reductionTable_tmp,strategy) {
574: local n,i;
575: n = Length(reductionTable_tmp);
576: for (i=0; i<n; i++) {
577: if (reductionTable_tmp[i] == strategy) {
578: return(true);
579: }
580: }
581: return(false);
582: }
583:
584: def SnextI(reductionTable_tmp,strategy,redundantTable,
585: skel,level,freeRes)
586: {
587: local ii,n,p,myindex,i,j,bases;
588: n = Length(reductionTable_tmp);
589: if (level == 0) {
590: for (ii=0; ii<n; ii++) {
591: if (reductionTable_tmp[ii] == strategy) {
592: return(ii);
593: }
594: }
595: }else{
596: for (ii=0; ii<n; ii++) {
597: if (reductionTable_tmp[ii] == strategy) {
598: p = skel[level,ii];
599: myindex = p[0];
600: i = myindex[0]; j = myindex[1];
601: bases = freeRes[level-1];
602: if (IsNull(bases[i]) || IsNull(bases[j])) {
603:
604: }else{
605: return(ii);
606: }
607: }
608: }
609: }
1.5 takayama 610: Print("reductionTable_tmp=");
1.4 takayama 611: Println(reductionTable_tmp);
1.5 takayama 612: Println("See also reductionTable, strategy, level,i");
1.4 takayama 613: Error("SnextI: bases[i] or bases[j] is null for all combinations.");
614: }
615:
616:
1.1 takayama 617:
618: def SsetBettiTable(freeRes,g) {
619: local level,i, n,bases,ans;
620: ans = NewArray(Length(freeRes)+1);
621: n = Length(freeRes);
622: if (IsArray(g[0])) {
623: ans[0] = Length(g[0]);
624: }else{
625: ans[0] = 1;
626: }
627: for (level=0; level<n; level++) {
628: bases = freeRes[level];
629: if (IsArray(bases)) {
630: ans[level+1] = Length(bases);
631: }else{
632: ans[level+1] = 1;
633: }
634: }
635: return(ans);
636: }
637:
638: def SwhereInGB(f,tower) {
639: local i,n,p,q;
640: n = Length(tower);
641: for (i=0; i<n; i++) {
642: p = MonomialPart(tower[i]);
643: q = MonomialPart(f);
644: if (p == q) return(i);
645: }
646: Println([f,tower]);
647: Error("whereInGB : [f,myset]: f could not be found in the myset.");
648: }
649: def SunitOfFormat(pos,forms) {
650: local ans,i,n;
651: n = Length(forms);
652: ans = NewArray(n);
653: for (i=0; i<n; i++) {
654: if (i != pos) {
655: ans[i] = Poly("0");
656: }else{
657: ans[i] = Poly("1");
658: }
659: }
660: return(ans);
661: }
662:
663: def Error(s) {
664: sm1(" s error ");
665: }
666:
667: def IsNull(s) {
668: if (Stag(s) == 0) return(true);
669: else return(false);
670: }
671:
672: def StowerOf(tower,level) {
673: local ans,i;
674: ans = [ ];
675: if (level == 0) return([[]]);
676: for (i=0; i<level; i++) {
677: ans = Append(ans,tower[i]);
678: }
679: return(Reverse(ans));
680: }
681:
682: def Sspolynomial(f,g) {
683: if (IsArray(f)) {
684: f = Stoes_vec(f);
685: }
686: if (IsArray(g)) {
687: g = Stoes_vec(g);
688: }
689: sm1("f g spol /FunctionValue set");
690: }
691:
692: def MonomialPart(f) {
693: sm1(" [(lmonom) f] gbext /FunctionValue set ");
694: }
695:
1.14 takayama 696: /* WARNING:
697: When you use SwhereInTower, you have to change gbList
698: as below. Ofcourse, you should restrore the gbList
699: SsetTower(StowerOf(tower,level));
700: pos = SwhereInTower(syzHead,tower[level]);
701: */
1.1 takayama 702: def SwhereInTower(f,tower) {
703: local i,n,p,q;
704: if (f == Poly("0")) return(-1);
705: n = Length(tower);
706: for (i=0; i<n; i++) {
707: p = MonomialPart(tower[i]);
708: q = MonomialPart(f);
709: if (p == q) return(i);
710: }
711: Println([f,tower]);
712: Error("[f,tower]: f could not be found in the tower.");
713: }
714:
715: def Stag(f) {
716: sm1(f," tag (universalNumber) dc /FunctionValue set");
717: }
718:
719: def SpairAndReduction(skel,level,ii,freeRes,tower,ww) {
720: local i, j, myindex, p, bases, tower2, gi, gj,
721: si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
722: vdeg,vdeg_reduced;
723: Println("SpairAndReduction:");
724:
725: if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
726: p = skel[level,ii];
727: myindex = p[0];
728: i = myindex[0]; j = myindex[1];
729: bases = freeRes[level-1];
730: Println(["p and bases ",p,bases]);
731: if (IsNull(bases[i]) || IsNull(bases[j])) {
732: Println([level,i,j,bases[i],bases[j]]);
733: Error("level, i, j : bases[i], bases[j] must not be NULL.");
734: }
735:
736: tower2 = StowerOf(tower,level-1);
737: SsetTower(tower2);
1.14 takayama 738: Println(["level=",level]);
739: Println(["tower2=",tower2]);
1.1 takayama 740: /** sm1(" show_ring "); */
741:
742: gi = Stoes_vec(bases[i]);
743: gj = Stoes_vec(bases[j]);
744:
745: ssp = Sspolynomial(gi,gj);
746: si = ssp[0,0];
747: sj = ssp[0,1];
748: syzHead = si*es^i;
749: /* This will be the head term, I think. But, double check. */
750: Println([si*es^i,sj*es^j]);
751:
752: Print("[gi, gj] = "); Println([gi,gj]);
753: sm1(" [(Homogenize)] system_variable message ");
754: Print("Reduce the element "); Println(si*gi+sj*gj);
755: Print("by "); Println(bases);
756:
757: tmp = Sreduction(si*gi+sj*gj, bases);
758:
759: Print("result is "); Println(tmp);
760:
1.3 takayama 761: /* This is essential part for V-minimal resolution. */
762: /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
763: vdeg = SvDegree(si*gi,tower,level-1,ww);
1.1 takayama 764: vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
765: Print("vdegree of the original = "); Println(vdeg);
766: Print("vdegree of the remainder = "); Println(vdeg_reduced);
767:
768: t_syz = tmp[2];
769: si = si*tmp[1]+t_syz[i];
770: sj = sj*tmp[1]+t_syz[j];
771: t_syz[i] = si;
772: t_syz[j] = sj;
1.14 takayama 773:
774: SsetTower(StowerOf(tower,level));
1.1 takayama 775: pos = SwhereInTower(syzHead,tower[level]);
1.14 takayama 776:
777: SsetTower(StowerOf(tower,level-1));
1.1 takayama 778: pos2 = SwhereInTower(tmp[0],tower[level-1]);
779: ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
780: /* pos is the place to put syzygy at level. */
781: /* pos2 is the place to put a new GB at level-1. */
782: Println(ans);
783: return(ans);
784: }
785:
786: def Sreduction(f,myset) {
787: local n, indexTable, set2, i, j, tmp, t_syz;
788: n = Length(myset);
789: indexTable = NewArray(n);
790: set2 = [ ];
791: j = 0;
792: for (i=0; i<n; i++) {
793: if (IsNull(myset[i])) {
794: indexTable[i] = -1;
795: /* }else if (myset[i] == Poly("0")) {
796: indexTable[i] = -1; */
797: }else{
798: set2 = Append(set2,Stoes_vec(myset[i]));
799: indexTable[i] = j;
800: j++;
801: }
802: }
803: sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
804: t_syz = NewArray(n);
805: for (i=0; i<n; i++) {
806: if (indexTable[i] != -1) {
807: t_syz[i] = tmp[2, indexTable[i]];
808: }else{
809: t_syz[i] = Poly("0");
810: }
811: }
812: return([tmp[0],tmp[1],t_syz]);
813: }
814:
815: def Warning(s) {
816: Print("Warning: ");
817: Println(s);
818: }
819: def RingOf(f) {
820: local r;
821: if (IsPolynomial(f)) {
822: if (f != Poly("0")) {
823: sm1(f," (ring) dc /r set ");
824: }else{
825: sm1(" [(CurrentRingp)] system_variable /r set ");
826: }
827: }else{
828: Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
829: sm1(" [(CurrentRingp)] system_variable /r set ");
830: }
831: return(r);
832: }
833:
834: def Sfrom_es(f,size) {
835: local c,ans, i, d, myes, myee, j,n,r,ans2;
836: if (Length(Arglist) < 2) size = -1;
837: if (IsArray(f)) return(f);
838: r = RingOf(f);
839: myes = PolyR("es",r);
840: myee = PolyR("e_",r);
841: if (Degree(f,myee) > 0 && size == -1) {
842: if (size == -1) {
843: sm1(f," (array) dc /ans set");
844: return(ans);
845: }
846: }
847:
848: /*
849: Coefficients(x^2-1,x):
850: [ [ 2 , 0 ] , [ 1 , -1 ] ]
851: */
852: if (Degree(f,myee) > 0) {
853: c = Coefficients(f,myee);
854: }else{
855: c = Coefficients(f,myes);
856: }
857: if (size < 0) {
858: size = c[0,0]+1;
859: }
860: ans = NewArray(size);
861: for (i=0; i<size; i++) {ans[i] = 0;}
862: n = Length(c[0]);
863: for (j=0; j<n; j++) {
864: d = c[0,j];
865: ans[d] = c[1,j];
866: }
867: return(ans);
868: }
869:
870: def Sbases_to_vec(bases,size) {
871: local n, giveSize, newbases,i;
872: /* bases = [1+es*x, [1,2,3*x]] */
873: if (Length(Arglist) > 1) {
874: giveSize = true;
875: }else{
876: giveSize = false;
877: }
878: n = Length(bases);
879: newbases = NewArray(n);
880: for (i=0; i<n; i++) {
881: if (giveSize) {
882: newbases[i] = Sfrom_es(bases[i], size);
883: }else{
884: newbases[i] = Sfrom_es(bases[i]);
885: }
886: }
887: return(newbases);
888: }
889:
1.14 takayama 890: HelpAdd(["Sminimal",
891: ["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm",
1.15 takayama 892: "option: \"homogenized\" (no automatic homogenization ",
1.14 takayama 893: "Example: Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
894: " v=[[2*x*Dx + 3*y*Dy+6, 0],",
895: " [3*x^2*Dy + 2*y*Dx, 0],",
896: " [0, x^2+y^2],",
897: " [0, x*y]];",
898: " a=Sminimal(v);",
899: " Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
900: " b = ReParse(a[0]); sm1_pmat(b); ",
901: " IsExact_h(b,[x,y]):",
902: "Note: a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
903:
1.15 takayama 904: def Sminimal(g,opt) {
1.1 takayama 905: local r, freeRes, redundantTable, reducer, maxLevel,
906: minRes, seq, maxSeq, level, betti, q, bases, dr,
1.14 takayama 907: betti_levelplus, newbases, i, j,qq, tminRes;
1.16 ! takayama 908: if (Length(Arglist) < 2) {
! 909: opt = null;
! 910: }
! 911: ScheckIfSchreyer("Sminimal:0");
1.15 takayama 912: r = SlaScala(g,opt);
1.1 takayama 913: /* Should I turn off the tower?? */
1.16 ! takayama 914: ScheckIfSchreyer("Sminimal:1");
1.1 takayama 915: freeRes = r[0];
916: redundantTable = r[1];
917: reducer = r[2];
918: minRes = SnewArrayOfFormat(freeRes);
919: seq = 0;
920: maxSeq = SgetMaxSeq(redundantTable);
921: maxLevel = Length(freeRes);
922: for (level = 0; level < maxLevel; level++) {
923: minRes[level] = freeRes[level];
924: }
925: seq=maxSeq+1;
926: while (seq > 1) {
927: seq--;
928: for (level = 0; level < maxLevel; level++) {
929: betti = Length(freeRes[level]);
930: for (q = 0; q<betti; q++) {
931: if (redundantTable[level,q] == seq) {
932: Print("[seq,level,q]="); Println([seq,level,q]);
933: if (level < maxLevel-1) {
934: bases = freeRes[level+1];
935: dr = reducer[level,q];
936: dr[q] = -1;
937: newbases = SnewArrayOfFormat(bases);
938: betti_levelplus = Length(bases);
939: /*
940: bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
941: */
942: for (i=0; i<betti_levelplus; i++) {
943: newbases[i] = bases[i] + bases[i,q]*dr;
944: }
945: Println(["level, q =", level,q]);
946: Println("bases="); sm1_pmat(bases);
947: Println("dr="); sm1_pmat(dr);
948: Println("newbases="); sm1_pmat(newbases);
949: minRes[level+1] = newbases;
950: freeRes = minRes;
951: #ifdef DEBUG
952: for (qq=0; qq<betti; qq++) {
953: if ((redundantTable[level,qq] >= seq) &&
954: (redundantTable[level,qq] <= maxSeq)) {
955: for (i=0; i<betti_levelplus; i++) {
956: if (!IsZero(newbases[i,qq])) {
957: Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
958: Print("redundantTable ="); sm1_pmat(redundantTable[level]);
959: Error("Stop in Sminimal for debugging.");
960: }
961: }
962: }
963: }
964: #endif
965: }
966: }
967: }
968: }
969: }
1.14 takayama 970: tminRes = Stetris(minRes,redundantTable);
971: return([SpruneZeroRow(tminRes), tminRes,
1.3 takayama 972: [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
1.1 takayama 973: /* r[4] is the redundantTable_ordinary */
1.3 takayama 974: /* r[0] is the freeResolution */
1.1 takayama 975: }
976:
977:
978: def IsZero(f) {
979: if (IsPolynomial(f)) {
980: return( f == Poly("0"));
981: }else if (IsInteger(f)) {
982: return( f == 0);
983: }else if (IsSm1Integer(f)) {
984: return( f == true );
985: }else if (IsDouble(f)) {
986: return( f == 0.0 );
987: }else if (IsRational(f)) {
988: return(IsZero(Denominator(f)));
989: }else{
990: Error("IsZero: cannot deal with this data type.");
991: }
992: }
993: def SgetMaxSeq(redundantTable) {
994: local level,i,n,ans, levelMax,bases;
995: levelMax = Length( redundantTable );
996: ans = 0;
997: for (level = 0; level < levelMax; level++) {
998: bases = redundantTable[level];
999: n = Length(bases);
1000: for (i=0; i<n; i++) {
1001: if (IsInteger( bases[i] )) {
1002: if (bases[i] > ans) {
1003: ans = bases[i];
1004: }
1005: }
1006: }
1007: }
1008: return(ans);
1009: }
1010:
1011: def Stetris(freeRes,redundantTable) {
1012: local level, i, j, resLength, minRes,
1013: bases, newbases, newbases2;
1014: minRes = SnewArrayOfFormat(freeRes);
1015: resLength = Length( freeRes );
1016: for (level=0; level<resLength; level++) {
1017: bases = freeRes[level];
1018: newbases = SnewArrayOfFormat(bases);
1019: betti = Length(bases); j = 0;
1020: /* Delete rows */
1021: for (i=0; i<betti; i++) {
1022: if (redundantTable[level,i] < 1) {
1023: newbases[j] = bases[i];
1024: j++;
1025: }
1026: }
1027: bases = SfirstN(newbases,j);
1028: if (level > 0) {
1029: /* Delete columns */
1030: newbases = Transpose(bases);
1031: betti = Length(newbases); j = 0;
1032: newbases2 = SnewArrayOfFormat(newbases);
1033: for (i=0; i<betti; i++) {
1034: if (redundantTable[level-1,i] < 1) {
1035: newbases2[j] = newbases[i];
1036: j++;
1037: }
1038: }
1039: newbases = Transpose(SfirstN(newbases2,j));
1040: }else{
1041: newbases = bases;
1042: }
1043: Println(["level=", level]);
1044: sm1_pmat(bases);
1045: sm1_pmat(newbases);
1046:
1047: minRes[level] = newbases;
1048: }
1049: return(minRes);
1050: }
1051:
1052: def SfirstN(bases,k) {
1053: local ans,i;
1054: ans = NewArray(k);
1055: for (i=0; i<k; i++) {
1056: ans[i] = bases[i];
1057: }
1058: return(ans);
1059: }
1060:
1061:
1062: /* usage: tt is tower. ww is weight.
1063: a = SresolutionFrameWithTower(v);
1064: tt = a[1];
1065: ww = [x,1,y,1,Dx,1,Dy,1];
1066: SvDegree(x*es,tt,1,ww):
1067:
1068: In(17)=tt:
1069: [[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 ] ,
1070: [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ] ,
1071: [3*es^3*y*Dy ] ]
1072: In(18)=SvDegree(x*es,tt,1,ww):
1073: 3
1074: In(19)=SvDegree(x*es^3,tt,1,ww):
1075: 4
1076: In(20)=SvDegree(x,tt,2,ww):
1077: 4
1078:
1079: */
1080: def SvDegree(f,tower,level,w) {
1081: local i,ans;
1082: if (IsZero(f)) return(null);
1.3 takayama 1083: f = Init(f);
1.1 takayama 1084: if (level <= 0) {
1085: return(Sord_w(f,w));
1086: }
1087: i = Degree(f,es);
1088: ans = Sord_w(f,w) +
1089: SvDegree(tower[level-1,i],tower,level-1,w);
1090: return(ans);
1091: }
1092:
1.2 takayama 1093: def Sannfs(f,v) {
1094: local f2;
1095: f2 = ToString(f);
1096: if (IsArray(v)) {
1097: v = Map(v,"ToString");
1098: }
1099: sm1(" [f2 v] annfs /FunctionValue set ");
1100: }
1101:
1102: /* Sannfs2("x^3-y^2"); */
1103: def Sannfs2(f) {
1104: local p,pp;
1105: p = Sannfs(f,"x,y");
1.6 takayama 1106: sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
1.5 takayama 1107: /*
1108: Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
1109: ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */
1.6 takayama 1110: /* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1]]); */
1.10 takayama 1111:
1.6 takayama 1112: Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1113: pp = Map(p,"Spoly");
1114: return(Sminimal_v(pp));
1115: /* return(Sminimal(pp)); */
1116: }
1117:
1.10 takayama 1118: HelpAdd(["Sannfs2",
1119: ["Sannfs2(f) constructs the V-minimal free resolution for the weight (-1,1)",
1120: "of the Laplace transform of the annihilating ideal of the polynomial f in x,y.",
1121: "See also Sminimal_v, Sannfs3.",
1122: "Example: a=Sannfs2(\"x^3-y^2\");",
1123: " b=a[0]; sm1_pmat(b);",
1124: " b[1]*b[0]:",
1125: "Example: a=Sannfs2(\"x*y*(x-y)*(x+y)\");",
1126: " b=a[0]; sm1_pmat(b);",
1127: " b[1]*b[0]:"
1128: ]]);
1129:
1.6 takayama 1130: /* Do not forget to turn on TOTAL_STRATEGY */
1131: def Sannfs2_laScala(f) {
1132: local p,pp;
1133: p = Sannfs(f,"x,y");
1134: /* Do not make laplace transform.
1135: sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
1136: p = [p];
1137: */
1138: Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1.2 takayama 1139: pp = Map(p[0],"Spoly");
1140: return(Sminimal(pp));
1141: }
1142:
1.11 takayama 1143: def Sannfs2_laScala2(f) {
1144: local p,pp;
1145: p = Sannfs(f,"x,y");
1146: sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
1147: p = [p];
1148: Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
1149: ["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1150: pp = Map(p[0],"Spoly");
1151: return(Sminimal(pp));
1152: }
1153:
1.3 takayama 1154: def Sannfs3(f) {
1155: local p,pp;
1156: p = Sannfs(f,"x,y,z");
1.6 takayama 1157: sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
1.3 takayama 1158: Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
1.6 takayama 1159: pp = Map(p,"Spoly");
1160: return(Sminimal_v(pp));
1.3 takayama 1161: }
1162:
1.10 takayama 1163: HelpAdd(["Sannfs3",
1164: ["Sannfs3(f) constructs the V-minimal free resolution for the weight (-1,1)",
1165: "of the Laplace transform of the annihilating ideal of the polynomial f in x,y,z.",
1166: "See also Sminimal_v, Sannfs2.",
1167: "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
1168: " b=a[0]; sm1_pmat(b);",
1169: " b[1]*b[0]: b[2]*b[1]:"]]);
1170:
1.2 takayama 1171: /*
1172: The betti numbers of most examples are 2,1. (0-th and 1-th).
1173: a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2.
1174: a=Sannfs2("x^3-y^2-x"); : it causes an error. It should be fixed.
1.3 takayama 1175: a=Sannfs2("x*y*(x-y)"); : it causes an error. It should be fixed.
1.2 takayama 1176:
1177: */
1178:
1.11 takayama 1179: def Sannfs3_laScala2(f) {
1180: local p,pp;
1181: p = Sannfs(f,"x,y,z");
1182: sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
1183: Sweyl("x,y,z",[["x",1,"y",1,"z",1,"Dx",1,"Dy",1,"Dz",1,"h",1],
1184: ["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
1185: pp = Map(p,"Spoly");
1186: return(Sminimal(pp));
1187: }
1.5 takayama 1188:
1189:
1.6 takayama 1190: /* The below does not use LaScala-Stillman's algorithm. */
1.5 takayama 1191: def Sschreyer(g) {
1192: local rf, tower, reductionTable, skel, redundantTable, bases,
1193: strategy, maxOfStrategy, height, level, n, i,
1194: freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
1195: redundantTable_ordinary, redundant_seq_ordinary,
1.6 takayama 1196: reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
1.5 takayama 1197: /* extern WeightOfSweyl; */
1198: ww = WeightOfSweyl;
1199: Print("WeghtOfSweyl="); Println(WeightOfSweyl);
1200: rf = SresolutionFrameWithTower(g);
1201: redundant_seq = 1; redundant_seq_ordinary = 1;
1202: tower = rf[1];
1.16 ! takayama 1203: Println("Generating reduction table which gives an order of reduction.");
! 1204: Println("But, you are in Sschreyer...., you may not use LaScala-Stillman");
! 1205: Print("WeghtOfSweyl="); Println(WeightOfSweyl);
! 1206: Print("tower"); Println(tower);
1.5 takayama 1207: reductionTable = SgenerateTable(tower);
1208: skel = rf[2];
1209: redundantTable = SnewArrayOfFormat(rf[1]);
1210: redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
1211: reducer = SnewArrayOfFormat(rf[1]);
1212: freeRes = SnewArrayOfFormat(rf[1]);
1213: bettiTable = SsetBettiTable(rf[1],g);
1214:
1215: height = Length(reductionTable);
1216: for (level = 0; level < height; level++) {
1217: n = Length(reductionTable[level]);
1218: for (i=0; i<n; i++) {
1219: Println([level,i]);
1220: Print("Processing "); Print([level,i]);
1221: if (level == 0) {
1222: if (IsNull(redundantTable[level,i])) {
1223: bases = freeRes[level];
1224: /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
1225: pos = SwhereInGB(tower[0,i],rf[3,0]);
1226: bases[i] = rf[3,0,pos];
1227: /* redundantTable[level,i] = 0;
1228: redundantTable_ordinary[level,i] = 0; */
1229: freeRes[level] = bases;
1230: /* Println(["GB=",i,bases,tower[0,i]]); */
1231: }
1232: }else{ /* level >= 1 */
1233: if (IsNull(redundantTable[level,i])) {
1234: bases = freeRes[level];
1235: f = SpairAndReduction2(skel,level,i,freeRes,tower,
1236: ww,redundantTable);
1237: if (f[0] != Poly("0")) {
1238: place = f[3];
1239: /* (level-1, place) is the place for f[0],
1240: which is a newly obtained GB. */
1241: #ifdef ORDINARY
1242: redundantTable[level-1,place] = redundant_seq;
1243: redundant_seq++;
1244: #else
1245: if (f[4] > f[5]) {
1246: /* Zero in the gr-module */
1247: Print("v-degree of [org,remainder] = ");
1248: Println([f[4],f[5]]);
1249: Print("[level,i] = "); Println([level,i]);
1250: redundantTable[level-1,place] = 0;
1251: }else{
1252: redundantTable[level-1,place] = redundant_seq;
1253: redundant_seq++;
1254: }
1255: #endif
1256: redundantTable_ordinary[level-1,place]
1257: =redundant_seq_ordinary;
1258: redundant_seq_ordinary++;
1259: bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */
1260: /* redundantTable[level,i] = 0;
1261: redundantTable_ordinary[level,i] = 0; */
1262: /* i must be equal to f[2], I think. Double check. */
1263:
1264: /* Correction Of Constant */
1.9 takayama 1265: /* Correction of syzygy */
1.7 takayama 1266: c2 = f[6]; /* or -f[6]? Double check. */
1267: Print("c2="); Println(c2);
1.5 takayama 1268: nn = Length(bases);
1269: for (ii=0; ii<nn;ii++) {
1.8 takayama 1270: if ((ii != i) && (! IsNull(bases[ii]))) {
1.7 takayama 1271: m = Length(bases[ii]);
1272: for (jj=0; jj<m; jj++) {
1273: if (jj != place) {
1274: bases[ii,jj] = bases[ii,jj]*c2;
1275: }
1276: }
1.5 takayama 1277: }
1278: }
1279:
1.7 takayama 1280: Print("Old freeRes[level] = "); sm1_pmat(freeRes[level]);
1.5 takayama 1281: freeRes[level] = bases;
1.7 takayama 1282: Print("New freeRes[level] = "); sm1_pmat(freeRes[level]);
1.6 takayama 1283:
1284: /* Update the freeRes[level-1] */
1.7 takayama 1285: Print("Old freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6 takayama 1286: bases = freeRes[level-1];
1287: bases[place] = f[0];
1288: freeRes[level-1] = bases;
1.7 takayama 1289: Print("New freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6 takayama 1290:
1.9 takayama 1291: reducer[level-1,place] = f[1]-SunitOfFormat(place,f[1]);
1292: /* This reducer is different from that of SlaScala(). */
1293:
1294: reducerBasis = reducer[level-1];
1295: nn = Length(reducerBasis);
1296: for (ii=0; ii<nn;ii++) {
1297: if ((ii != place) && (! IsNull(reducerBasis[ii]))) {
1298: m = Length(reducerBasis[ii]);
1299: for (jj=0; jj<m; jj++) {
1300: if (jj != place) {
1301: reducerBasis[ii,jj] = reducerBasis[ii,jj]*c2;
1302: }
1303: }
1304: }
1305: }
1306: reducer[level-1] = reducerBasis;
1307:
1.5 takayama 1308: }else{
1309: /* redundantTable[level,i] = 0; */
1310: bases = freeRes[level];
1311: bases[i] = f[1]; /* Put the syzygy. */
1312: freeRes[level] = bases;
1313: }
1314: } /* end of level >= 1 */
1315: }
1316: } /* i loop */
1.6 takayama 1317:
1318: /* Triangulate reducer */
1319: if (level >= 1) {
1320: Println(" ");
1321: Print("Triangulating reducer at level "); Println(level-1);
1.9 takayama 1322: Println("freeRes[level]="); sm1_pmat(freeRes[level]);
1.6 takayama 1323: reducerBase = reducer[level-1];
1324: Print("reducerBase="); Println(reducerBase);
1.9 takayama 1325: Println("Compare freeRes[level] and reducerBase (put -1)");
1.6 takayama 1326: m = Length(reducerBase);
1327: for (ii=m-1; ii>=0; ii--) {
1328: if (!IsNull(reducerBase[ii])) {
1329: for (jj=ii-1; jj>=0; jj--) {
1330: if (!IsNull(reducerBase[jj])) {
1331: if (!IsZero(reducerBase[jj,ii])) {
1.9 takayama 1332: /* reducerBase[ii,ii] should be always constant. */
1333: reducerBase[jj] = reducerBase[ii,ii]*reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
1.6 takayama 1334: }
1335: }
1336: }
1337: }
1338: }
1339: Println("New reducer");
1340: sm1_pmat(reducerBase);
1341: reducer[level-1] = reducerBase;
1342: }
1343:
1.5 takayama 1344: } /* level loop */
1345: n = Length(freeRes);
1346: freeResV = SnewArrayOfFormat(freeRes);
1347: for (i=0; i<n; i++) {
1348: bases = freeRes[i];
1349: bases = Sbases_to_vec(bases,bettiTable[i]);
1350: freeResV[i] = bases;
1351: }
1.6 takayama 1352:
1353: /* Mark the non-redundant elements. */
1354: for (i=0; i<n; i++) {
1355: m = Length(redundantTable[i]);
1356: for (jj=0; jj<m; jj++) {
1357: if (IsNull(redundantTable[i,jj])) {
1358: redundantTable[i,jj] = 0;
1359: }
1360: }
1361: }
1362:
1363:
1.5 takayama 1364: return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
1365: }
1366:
1367: def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
1368: local i, j, myindex, p, bases, tower2, gi, gj,
1369: si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
1370: vdeg,vdeg_reduced,n,c2;
1.6 takayama 1371: Println("SpairAndReduction2 : -------------------------");
1.5 takayama 1372:
1373: if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
1374: p = skel[level,ii];
1375: myindex = p[0];
1376: i = myindex[0]; j = myindex[1];
1377: bases = freeRes[level-1];
1378: Println(["p and bases ",p,bases]);
1379: if (IsNull(bases[i]) || IsNull(bases[j])) {
1380: Println([level,i,j,bases[i],bases[j]]);
1381: Error("level, i, j : bases[i], bases[j] must not be NULL.");
1382: }
1383:
1384: tower2 = StowerOf(tower,level-1);
1385: SsetTower(tower2);
1.14 takayama 1386: Println(["level=",level]);
1387: Println(["tower2=",tower2]);
1.5 takayama 1388: /** sm1(" show_ring "); */
1389:
1390: gi = Stoes_vec(bases[i]);
1391: gj = Stoes_vec(bases[j]);
1392:
1393: ssp = Sspolynomial(gi,gj);
1394: si = ssp[0,0];
1395: sj = ssp[0,1];
1396: syzHead = si*es^i;
1397: /* This will be the head term, I think. But, double check. */
1398: Println([si*es^i,sj*es^j]);
1399:
1400: Print("[gi, gj] = "); Println([gi,gj]);
1401: sm1(" [(Homogenize)] system_variable message ");
1402: Print("Reduce the element "); Println(si*gi+sj*gj);
1403: Print("by "); Println(bases);
1404:
1405: tmp = Sreduction(si*gi+sj*gj, bases);
1406:
1407: Print("result is "); Println(tmp);
1.6 takayama 1408: if (!IsZero(tmp[0])) {
1409: Print("Error: base = ");
1410: Println(Map(bases,"Stoes_vec"));
1411: Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
1412: }
1.5 takayama 1413: t_syz = tmp[2];
1414: si = si*tmp[1]+t_syz[i];
1415: sj = sj*tmp[1]+t_syz[j];
1416: t_syz[i] = si;
1417: t_syz[j] = sj;
1418:
1419: c2 = null;
1420: /* tmp[0] must be zero */
1421: n = Length(t_syz);
1422: for (i=0; i<n; i++) {
1.6 takayama 1423: if (IsConstant(t_syz[i])){
1424: if (!IsZero(t_syz[i])) {
1.5 takayama 1425: if (IsNull(redundantTable[level-1,i])) {
1426: /* i must equal to pos2 below. */
1427: c2 = -t_syz[i];
1.6 takayama 1428: tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
1.5 takayama 1429: t_syz[i] = 0;
1.6 takayama 1430: /* tmp[0] = t_syz . g */
1.5 takayama 1431: /* break; does not work. Use */
1432: i = n;
1433: }
1.6 takayama 1434: }
1.5 takayama 1435: }
1436: }
1437:
1438: /* This is essential part for V-minimal resolution. */
1439: /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
1440: vdeg = SvDegree(si*gi,tower,level-1,ww);
1441: vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
1442: Print("vdegree of the original = "); Println(vdeg);
1443: Print("vdegree of the remainder = "); Println(vdeg_reduced);
1444:
1.11 takayama 1445: if (!IsNull(vdeg_reduced)) {
1446: if (vdeg_reduced < vdeg) {
1447: Println("--- Special in V-minimal!");
1448: Println(tmp[0]);
1449: Println("syzygy="); sm1_pmat(t_syz);
1450: Print("[vdeg, vdeg_reduced] = "); Println([vdeg,vdeg_reduced]);
1451: }
1452: }
1453:
1.14 takayama 1454: SsetTower(StowerOf(tower,level));
1455: pos = SwhereInTower(syzHead,tower[level]);
1.11 takayama 1456:
1.14 takayama 1457: SsetTower(StowerOf(tower,level-1));
1.5 takayama 1458: pos2 = SwhereInTower(tmp[0],tower[level-1]);
1459: ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
1460: /* pos is the place to put syzygy at level. */
1461: /* pos2 is the place to put a new GB at level-1. */
1462: Println(ans);
1.14 takayama 1463: Println("--- end of SpairAndReduction2 ");
1.5 takayama 1464: return(ans);
1465: }
1.6 takayama 1466:
1.10 takayama 1467: HelpAdd(["Sminimal_v",
1468: ["It constructs the V-minimal free resolution from the Schreyer resolution",
1469: "step by step.",
1.14 takayama 1470: "This code still contains bugs. It sometimes outputs wrong answer.",
1.10 takayama 1471: "Example: Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
1472: " v=[[2*x*Dx + 3*y*Dy+6, 0],",
1473: " [3*x^2*Dy + 2*y*Dx, 0],",
1474: " [0, x^2+y^2],",
1475: " [0, x*y]];",
1476: " a=Sminimal_v(v);",
1477: " sm1_pmat(a[0]); b=a[0]; b[1]*b[0]:",
1478: "Note: a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
1479:
1.14 takayama 1480: /* This code still contains bugs. It sometimes outputs wrong answer. */
1481: /* See test12() in minimal-test.k. */
1482: /* There may be remaining 1, too */
1.6 takayama 1483: def Sminimal_v(g) {
1484: local r, freeRes, redundantTable, reducer, maxLevel,
1485: minRes, seq, maxSeq, level, betti, q, bases, dr,
1.10 takayama 1486: betti_levelplus, newbases, i, j,qq,tminRes;
1.6 takayama 1487: r = Sschreyer(g);
1488: sm1_pmat(r);
1489: Debug_Sminimal_v = r;
1490: Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
1491: /* Should I turn off the tower?? */
1492: freeRes = r[0];
1493: redundantTable = r[1];
1494: reducer = r[2];
1495: minRes = SnewArrayOfFormat(freeRes);
1496: seq = 0;
1497: maxSeq = SgetMaxSeq(redundantTable);
1498: maxLevel = Length(freeRes);
1499: for (level = 0; level < maxLevel; level++) {
1500: minRes[level] = freeRes[level];
1501: }
1502: for (level = 0; level < maxLevel; level++) {
1503: betti = Length(freeRes[level]);
1504: for (q = betti-1; q>=0; q--) {
1505: if (redundantTable[level,q] > 0) {
1506: Print("[seq,level,q]="); Println([seq,level,q]);
1507: if (level < maxLevel-1) {
1508: bases = freeRes[level+1];
1509: dr = reducer[level,q];
1.9 takayama 1510: /* dr[q] = -1; We do not need this in our reducer format. */
1511: /* dr[q] should be a non-zero constant. */
1.6 takayama 1512: newbases = SnewArrayOfFormat(bases);
1513: betti_levelplus = Length(bases);
1514: /*
1515: bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
1516: */
1517: for (i=0; i<betti_levelplus; i++) {
1.9 takayama 1518: newbases[i] = dr[q]*bases[i] - bases[i,q]*dr;
1.6 takayama 1519: }
1520: Println(["level, q =", level,q]);
1521: Println("bases="); sm1_pmat(bases);
1522: Println("dr="); sm1_pmat(dr);
1523: Println("newbases="); sm1_pmat(newbases);
1524: minRes[level+1] = newbases;
1525: freeRes = minRes;
1526: #ifdef DEBUG
1.9 takayama 1527: for (qq=q; qq<betti; qq++) {
1.6 takayama 1528: for (i=0; i<betti_levelplus; i++) {
1.9 takayama 1529: if ((!IsZero(newbases[i,qq])) && (redundantTable[level,qq] >0)) {
1.6 takayama 1530: Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
1531: Print("redundantTable ="); sm1_pmat(redundantTable[level]);
1532: Error("Stop in Sminimal for debugging.");
1533: }
1534: }
1535: }
1536: #endif
1537: }
1538: }
1539: }
1540: }
1.10 takayama 1541: tminRes = Stetris(minRes,redundantTable);
1542: return([SpruneZeroRow(tminRes), tminRes,
1.6 takayama 1543: [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
1544: /* r[4] is the redundantTable_ordinary */
1545: /* r[0] is the freeResolution */
1546: }
1547:
1548: /* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */
1.10 takayama 1549: /* x y (x+y-1)(x-2), x^3-y^2, x^3 - y^2 z^2,
1550: x y z (x+y+z-1) seems to be interesting, because the first syzygy
1551: contains 1.
1552: */
1553:
1554: def CopyArray(m) {
1555: local ans,i,n;
1556: if (IsArray(m)) {
1557: n = Length(m);
1558: ans = NewArray(n);
1559: for (i=0; i<n; i++) {
1560: ans[i] = CopyArray(m[i]);
1561: }
1562: return(ans);
1563: }else{
1564: return(m);
1565: }
1566: }
1567: HelpAdd(["CopyArray",
1568: ["It duplicates the argument array recursively.",
1569: "Example: m=[1,[2,3]];",
1570: " a=CopyArray(m); a[1] = \"Hello\";",
1571: " Println(m); Println(a);"]]);
1572:
1573: def IsZeroVector(m) {
1574: local n,i;
1575: n = Length(m);
1576: for (i=0; i<n; i++) {
1577: if (!IsZero(m[i])) {
1578: return(false);
1579: }
1580: }
1581: return(true);
1582: }
1583:
1584: def SpruneZeroRow(res) {
1585: local minRes, n,i,j,m, base,base2,newbase,newbase2, newMinRes;
1586:
1587: minRes = CopyArray(res);
1588: n = Length(minRes);
1589: for (i=0; i<n; i++) {
1590: base = minRes[i];
1591: m = Length(base);
1592: if (i != n-1) {
1593: base2 = minRes[i+1];
1594: base2 = Transpose(base2);
1595: }
1596: newbase = [ ];
1597: newbase2 = [ ];
1598: for (j=0; j<m; j++) {
1599: if (!IsZeroVector(base[j])) {
1600: newbase = Append(newbase,base[j]);
1601: if (i != n-1) {
1602: newbase2 = Append(newbase2,base2[j]);
1603: }
1604: }
1605: }
1606: minRes[i] = newbase;
1607: if (i != n-1) {
1608: if (newbase2 == [ ]) {
1609: minRes[i+1] = [ ];
1610: }else{
1611: minRes[i+1] = Transpose(newbase2);
1612: }
1613: }
1614: }
1615:
1616: newMinRes = [ ];
1617: n = Length(minRes);
1618: i = 0;
1619: while (i < n ) {
1620: base = minRes[i];
1621: if (base == [ ]) {
1622: i = n; /* break; */
1623: }else{
1624: newMinRes = Append(newMinRes,base);
1625: }
1626: i++;
1627: }
1628: return(newMinRes);
1629: }
1630:
1631: def testAnnfs2(f) {
1632: local a,i,n;
1633: a = Sannfs2(f);
1634: b=a[0];
1635: n = Length(b);
1636: Println("------ V-minimal free resolution -----");
1637: sm1_pmat(b);
1638: Println("----- Is it complex? ---------------");
1639: for (i=0; i<n-1; i++) {
1640: Println(b[i+1]*b[i]);
1641: }
1642: return(a);
1643: }
1644: def testAnnfs3(f) {
1645: local a,i,n;
1646: a = Sannfs3(f);
1647: b=a[0];
1648: n = Length(b);
1649: Println("------ V-minimal free resolution -----");
1650: sm1_pmat(b);
1651: Println("----- Is it complex? ---------------");
1652: for (i=0; i<n-1; i++) {
1653: Println(b[i+1]*b[i]);
1654: }
1.11 takayama 1655: return(a);
1656: }
1657:
1658: def ToString_array(p) {
1659: local ans;
1660: if (IsArray(p)) {
1661: ans = Map(p,"ToString_array");
1662: }else{
1663: ans = ToString(p);
1664: }
1665: return(ans);
1666: }
1667:
1668: /* sm1_res_div([[x],[y]],[[x^2],[x*y],[y^2]],[x,y]): */
1669:
1670: def sm1_res_div(I,J,V) {
1671: I = ToString_array(I);
1672: J = ToString_array(J);
1673: V = ToString_array(V);
1674: sm1(" [[ I J] V ] res*div /FunctionValue set ");
1675: }
1676:
1677: /* It has not yet been working */
1678: def sm1_res_kernel_image(m,n,v) {
1679: m = ToString_array(m);
1680: n = ToString_array(n);
1681: v = ToString_array(v);
1682: sm1(" [m n v] res-kernel-image /FunctionValue set ");
1683: }
1684: def Skernel(m,v) {
1685: m = ToString_array(m);
1686: v = ToString_array(v);
1687: sm1(" [ m v ] syz /FunctionValue set ");
1688: }
1689:
1690: def test3() {
1691: local a1,a2,b1,b2;
1692: a1 = Sannfs3("x^3-y^2*z^2");
1693: a1 = a1[0];
1694: a2 = Sannfs3_laScala2("x^3-y^2*z^2");
1695: a2 = a2[0];
1696: b1 = a1[1];
1697: b2 = a2[1];
1698: sm1_pmat(b2);
1699: Println(" OVER ");
1700: sm1_pmat(b1);
1701: return([sm1_res_div(b2,b1,["x","y","z"]),b2,b1,a2,a1]);
1702: }
1703:
1704: def test4() {
1705: local a,b;
1706: a = Sannfs3_laScala2("x^3-y^2*z^2");
1707: b = a[0];
1708: sm1_pmat( sm1_res_kernel_image(b[0],b[1],[x,y,z]));
1709: sm1_pmat( sm1_res_kernel_image(b[1],b[2],[x,y,z]));
1710: return(a);
1711: }
1712:
1713: def sm1_gb(f,v) {
1714: f =ToString_array(f);
1715: v = ToString_array(v);
1716: sm1(" [f v] gb /FunctionValue set ");
1.13 takayama 1717: }
1718:
1.11 takayama 1719:
1.12 takayama 1720: def SisComplex(a) {
1721: local n,i,j,k,b,p,q;
1722: n = Length(a);
1723: for (i=0; i<n-1; i++) {
1724: if (Length(a[i+1]) != 0) {
1725: b = a[i+1]*a[i];
1726: p = Length(b); q = Length(b[0]);
1727: for (j=0; j<p; j++) {
1728: for (k=0; k<q; k++) {
1729: if (!IsZero(b[j,k])) {
1730: Print("Is is not complex at ");
1731: Println([i,j,k]);
1732: return(false);
1733: }
1734: }
1735: }
1736: }
1737: }
1738: return(true);
1.14 takayama 1739: }
1740:
1741: def IsExact_h(c,v) {
1742: local a;
1743: v = ToString_array(v);
1744: a = [c,v];
1745: sm1(a," isExact_h /FunctionValue set ");
1746: }
1747: HelpAdd(["IsExact_h",
1748: ["IsExact_h(complex,var): bool",
1749: "It checks the given complex is exact or not in D<h> (homogenized Weyl algebra)",
1750: "cf. ReParse"
1751: ]]);
1752:
1753: def ReParse(a) {
1754: local c;
1755: if (IsArray(a)) {
1756: c = Map(a,"ReParse");
1757: }else{
1758: sm1(a," toString . /c set");
1759: }
1760: return(c);
1761: }
1762: HelpAdd(["ReParse",
1763: ["Reparse(obj): obj",
1764: "It parses the given object in the current ring.",
1765: "Outputs from SlaScala, Sschreyer may cause a trouble in other functions,",
1766: "because it uses the Schreyer order.",
1767: "In this case, ReParse the outputs from these functions.",
1768: "cf. IsExaxt_h"
1769: ]]);
1.16 ! takayama 1770:
! 1771: def ScheckIfSchreyer(s) {
! 1772: local ss;
! 1773: sm1(" (report) (grade) switch_function /ss set ");
! 1774: if (ss != "module1v") {
! 1775: Print("ScheckIfSchreyer: from "); Println(s);
! 1776: Error("grade is not module1v");
! 1777: }
! 1778: /*
! 1779: sm1(" (report) (mmLarger) switch_function /ss set ");
! 1780: if (ss != "tower") {
! 1781: Print("ScheckIfSchreyer: from "); Println(s);
! 1782: Error("mmLarger is not tower");
! 1783: }
! 1784: */
! 1785: sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set ");
! 1786: if (ss != 1) {
! 1787: Print("ScheckIfSchreyer: from "); Println(s);
! 1788: Error("Schreyer order is not set.");
! 1789: }
! 1790: /* More check will be necessary. */
! 1791: return(true);
! 1792: }
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