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