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