Annotation of OpenXM/src/k097/lib/minimal/minimal.k, Revision 1.4
1.4 ! takayama 1: /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.3 2000/05/04 06:55:28 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;
340: if (level <= 1) return(StotalDegree(f));
341: i = Degree(f,es);
342: return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
343: }
344:
345: def SgenerateTable(tower) {
346: local height, n,i,j, ans, ans_at_each_floor;
347: height = Length(tower);
348: ans = NewArray(height);
349: for (i=0; i<height; i++) {
350: n = Length(tower[i]);
351: ans_at_each_floor=NewArray(n);
352: for (j=0; j<n; j++) {
353: ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1);
354: /* Println([i,j,ans_at_each_floor[j]]); */
355: }
356: ans[i] = ans_at_each_floor;
357: }
358: return(ans);
359: }
360: Sweyl("x,y,z");
361: v=[[2*x*Dx + 3*y*Dy+6, 0],
362: [3*x^2*Dy + 2*y*Dx, 0],
363: [0, x^2+y^2],
364: [0, x*y]];
365: /* SresolutionFrameWithTower(v); */
366:
367: def SnewArrayOfFormat(p) {
368: if (IsArray(p)) {
369: return(Map(p,"SnewArrayOfFormat"));
370: }else{
371: return(null);
372: }
373: }
1.4 ! takayama 374: def ScopyArray(a) {
! 375: local n, i,ans;
! 376: n = Length(a);
! 377: ans = NewArray(n);
! 378: for (i=0; i<n; i++) {
! 379: ans[i] = a[i];
! 380: }
! 381: return(ans);
! 382: }
1.1 takayama 383: def SminOfStrategy(a) {
384: local n,i,ans,tt;
385: ans = 100000; /* very big number */
386: if (IsArray(a)) {
387: n = Length(a);
388: for (i=0; i<n; i++) {
389: if (IsArray(a[i])) {
390: tt = SminOfStrategy(a[i]);
391: if (tt < ans) ans = tt;
392: }else{
393: if (a[i] < ans) ans = a[i];
394: }
395: }
396: }else{
397: if (a < ans) ans = a;
398: }
399: return(ans);
400: }
401: def SmaxOfStrategy(a) {
402: local n,i,ans,tt;
403: ans = -100000; /* very small number */
404: if (IsArray(a)) {
405: n = Length(a);
406: for (i=0; i<n; i++) {
407: if (IsArray(a[i])) {
408: tt = SmaxOfStrategy(a[i]);
409: if (tt > ans) ans = tt;
410: }else{
411: if (a[i] > ans) ans = a[i];
412: }
413: }
414: }else{
415: if (a > ans) ans = a;
416: }
417: return(ans);
418: }
419:
420:
421: def SlaScala(g) {
422: local rf, tower, reductionTable, skel, redundantTable, bases,
423: strategy, maxOfStrategy, height, level, n, i,
424: freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
1.4 ! takayama 425: redundantTable_ordinary, redundant_seq_ordinary,
! 426: reductionTable_tmp;
1.1 takayama 427: /* extern WeightOfSweyl; */
428: ww = WeightOfSweyl;
429: Print("WeghtOfSweyl="); Println(WeightOfSweyl);
430: rf = SresolutionFrameWithTower(g);
431: redundant_seq = 1; redundant_seq_ordinary = 1;
432: tower = rf[1];
433: reductionTable = SgenerateTable(tower);
434: skel = rf[2];
435: redundantTable = SnewArrayOfFormat(rf[1]);
436: redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
437: reducer = SnewArrayOfFormat(rf[1]);
438: freeRes = SnewArrayOfFormat(rf[1]);
439: bettiTable = SsetBettiTable(rf[1],g);
440:
441: strategy = SminOfStrategy( reductionTable );
442: maxOfStrategy = SmaxOfStrategy( reductionTable );
443: height = Length(reductionTable);
444: while (strategy <= maxOfStrategy) {
445: for (level = 0; level < height; level++) {
446: n = Length(reductionTable[level]);
1.4 ! takayama 447: reductionTable_tmp = ScopyArray(reductionTable[level]);
! 448: while (SthereIs(reductionTable_tmp,strategy)) {
! 449: i = SnextI(reductionTable_tmp,strategy,redundantTable,
! 450: skel,level,freeRes);
! 451: Println([level,i]);
! 452: reductionTable_tmp[i] = -200000;
1.1 takayama 453: if (reductionTable[level,i] == strategy) {
454: Print("Processing "); Print([level,i]);
455: Print(" Strategy = "); Println(strategy);
456: if (level == 0) {
457: if (IsNull(redundantTable[level,i])) {
458: bases = freeRes[level];
459: /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
460: pos = SwhereInGB(tower[0,i],rf[3,0]);
461: bases[i] = rf[3,0,pos];
462: redundantTable[level,i] = 0;
463: redundantTable_ordinary[level,i] = 0;
464: freeRes[level] = bases;
465: /* Println(["GB=",i,bases,tower[0,i]]); */
466: }
467: }else{ /* level >= 1 */
468: if (IsNull(redundantTable[level,i])) {
469: bases = freeRes[level];
470: f = SpairAndReduction(skel,level,i,freeRes,tower,ww);
471: if (f[0] != Poly("0")) {
472: place = f[3];
473: /* (level-1, place) is the place for f[0],
474: which is a newly obtained GB. */
475: #ifdef ORDINARY
476: redundantTable[level-1,place] = redundant_seq;
477: redundant_seq++;
478: #else
479: if (f[4] > f[5]) {
480: /* Zero in the gr-module */
481: Print("v-degree of [org,remainder] = ");
482: Println([f[4],f[5]]);
483: Print("[level,i] = "); Println([level,i]);
484: redundantTable[level-1,place] = 0;
485: }else{
486: redundantTable[level-1,place] = redundant_seq;
487: redundant_seq++;
488: }
489: #endif
490: redundantTable_ordinary[level-1,place]
491: =redundant_seq_ordinary;
492: redundant_seq_ordinary++;
493: bases[i] = SunitOfFormat(place,f[1])-f[1]; /* syzygy */
494: redundantTable[level,i] = 0;
495: redundantTable_ordinary[level,i] = 0;
496: /* i must be equal to f[2], I think. Double check. */
497: freeRes[level] = bases;
498: bases = freeRes[level-1];
499: bases[place] = f[0];
500: freeRes[level-1] = bases;
501: reducer[level-1,place] = f[1];
502: }else{
503: redundantTable[level,i] = 0;
504: bases = freeRes[level];
505: bases[i] = f[1]; /* Put the syzygy. */
506: freeRes[level] = bases;
507: }
508: }
509: } /* end of level >= 1 */
510: }
511: }
512: }
513: strategy++;
514: }
515: n = Length(freeRes);
516: freeResV = SnewArrayOfFormat(freeRes);
517: for (i=0; i<n; i++) {
518: bases = freeRes[i];
519: bases = Sbases_to_vec(bases,bettiTable[i]);
520: freeResV[i] = bases;
521: }
522: return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
523: }
1.4 ! takayama 524:
! 525: def SthereIs(reductionTable_tmp,strategy) {
! 526: local n,i;
! 527: n = Length(reductionTable_tmp);
! 528: for (i=0; i<n; i++) {
! 529: if (reductionTable_tmp[i] == strategy) {
! 530: return(true);
! 531: }
! 532: }
! 533: return(false);
! 534: }
! 535:
! 536: def SnextI(reductionTable_tmp,strategy,redundantTable,
! 537: skel,level,freeRes)
! 538: {
! 539: local ii,n,p,myindex,i,j,bases;
! 540: n = Length(reductionTable_tmp);
! 541: if (level == 0) {
! 542: for (ii=0; ii<n; ii++) {
! 543: if (reductionTable_tmp[ii] == strategy) {
! 544: return(ii);
! 545: }
! 546: }
! 547: }else{
! 548: for (ii=0; ii<n; ii++) {
! 549: if (reductionTable_tmp[ii] == strategy) {
! 550: p = skel[level,ii];
! 551: myindex = p[0];
! 552: i = myindex[0]; j = myindex[1];
! 553: bases = freeRes[level-1];
! 554: if (IsNull(bases[i]) || IsNull(bases[j])) {
! 555:
! 556: }else{
! 557: return(ii);
! 558: }
! 559: }
! 560: }
! 561: }
! 562: Println(reductionTable_tmp);
! 563: Error("SnextI: bases[i] or bases[j] is null for all combinations.");
! 564: }
! 565:
! 566:
1.1 takayama 567:
568: def SsetBettiTable(freeRes,g) {
569: local level,i, n,bases,ans;
570: ans = NewArray(Length(freeRes)+1);
571: n = Length(freeRes);
572: if (IsArray(g[0])) {
573: ans[0] = Length(g[0]);
574: }else{
575: ans[0] = 1;
576: }
577: for (level=0; level<n; level++) {
578: bases = freeRes[level];
579: if (IsArray(bases)) {
580: ans[level+1] = Length(bases);
581: }else{
582: ans[level+1] = 1;
583: }
584: }
585: return(ans);
586: }
587:
588: def SwhereInGB(f,tower) {
589: local i,n,p,q;
590: n = Length(tower);
591: for (i=0; i<n; i++) {
592: p = MonomialPart(tower[i]);
593: q = MonomialPart(f);
594: if (p == q) return(i);
595: }
596: Println([f,tower]);
597: Error("whereInGB : [f,myset]: f could not be found in the myset.");
598: }
599: def SunitOfFormat(pos,forms) {
600: local ans,i,n;
601: n = Length(forms);
602: ans = NewArray(n);
603: for (i=0; i<n; i++) {
604: if (i != pos) {
605: ans[i] = Poly("0");
606: }else{
607: ans[i] = Poly("1");
608: }
609: }
610: return(ans);
611: }
612:
613: def Error(s) {
614: sm1(" s error ");
615: }
616:
617: def IsNull(s) {
618: if (Stag(s) == 0) return(true);
619: else return(false);
620: }
621:
622: def StowerOf(tower,level) {
623: local ans,i;
624: ans = [ ];
625: if (level == 0) return([[]]);
626: for (i=0; i<level; i++) {
627: ans = Append(ans,tower[i]);
628: }
629: return(Reverse(ans));
630: }
631:
632: def Sspolynomial(f,g) {
633: if (IsArray(f)) {
634: f = Stoes_vec(f);
635: }
636: if (IsArray(g)) {
637: g = Stoes_vec(g);
638: }
639: sm1("f g spol /FunctionValue set");
640: }
641:
642: def MonomialPart(f) {
643: sm1(" [(lmonom) f] gbext /FunctionValue set ");
644: }
645:
646: def SwhereInTower(f,tower) {
647: local i,n,p,q;
648: if (f == Poly("0")) return(-1);
649: n = Length(tower);
650: for (i=0; i<n; i++) {
651: p = MonomialPart(tower[i]);
652: q = MonomialPart(f);
653: if (p == q) return(i);
654: }
655: Println([f,tower]);
656: Error("[f,tower]: f could not be found in the tower.");
657: }
658:
659: def Stag(f) {
660: sm1(f," tag (universalNumber) dc /FunctionValue set");
661: }
662:
663: def SpairAndReduction(skel,level,ii,freeRes,tower,ww) {
664: local i, j, myindex, p, bases, tower2, gi, gj,
665: si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
666: vdeg,vdeg_reduced;
667: Println("SpairAndReduction:");
668:
669: if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
670: p = skel[level,ii];
671: myindex = p[0];
672: i = myindex[0]; j = myindex[1];
673: bases = freeRes[level-1];
674: Println(["p and bases ",p,bases]);
675: if (IsNull(bases[i]) || IsNull(bases[j])) {
676: Println([level,i,j,bases[i],bases[j]]);
677: Error("level, i, j : bases[i], bases[j] must not be NULL.");
678: }
679:
680: tower2 = StowerOf(tower,level-1);
681: SsetTower(tower2);
682: /** sm1(" show_ring "); */
683:
684: gi = Stoes_vec(bases[i]);
685: gj = Stoes_vec(bases[j]);
686:
687: ssp = Sspolynomial(gi,gj);
688: si = ssp[0,0];
689: sj = ssp[0,1];
690: syzHead = si*es^i;
691: /* This will be the head term, I think. But, double check. */
692: Println([si*es^i,sj*es^j]);
693:
694: Print("[gi, gj] = "); Println([gi,gj]);
695: sm1(" [(Homogenize)] system_variable message ");
696: Print("Reduce the element "); Println(si*gi+sj*gj);
697: Print("by "); Println(bases);
698:
699: tmp = Sreduction(si*gi+sj*gj, bases);
700:
701: Print("result is "); Println(tmp);
702:
1.3 takayama 703: /* This is essential part for V-minimal resolution. */
704: /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
705: vdeg = SvDegree(si*gi,tower,level-1,ww);
1.1 takayama 706: vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
707: Print("vdegree of the original = "); Println(vdeg);
708: Print("vdegree of the remainder = "); Println(vdeg_reduced);
709:
710: t_syz = tmp[2];
711: si = si*tmp[1]+t_syz[i];
712: sj = sj*tmp[1]+t_syz[j];
713: t_syz[i] = si;
714: t_syz[j] = sj;
715: pos = SwhereInTower(syzHead,tower[level]);
716: pos2 = SwhereInTower(tmp[0],tower[level-1]);
717: ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
718: /* pos is the place to put syzygy at level. */
719: /* pos2 is the place to put a new GB at level-1. */
720: Println(ans);
721: return(ans);
722: }
723:
724: def Sreduction(f,myset) {
725: local n, indexTable, set2, i, j, tmp, t_syz;
726: n = Length(myset);
727: indexTable = NewArray(n);
728: set2 = [ ];
729: j = 0;
730: for (i=0; i<n; i++) {
731: if (IsNull(myset[i])) {
732: indexTable[i] = -1;
733: /* }else if (myset[i] == Poly("0")) {
734: indexTable[i] = -1; */
735: }else{
736: set2 = Append(set2,Stoes_vec(myset[i]));
737: indexTable[i] = j;
738: j++;
739: }
740: }
741: sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
742: t_syz = NewArray(n);
743: for (i=0; i<n; i++) {
744: if (indexTable[i] != -1) {
745: t_syz[i] = tmp[2, indexTable[i]];
746: }else{
747: t_syz[i] = Poly("0");
748: }
749: }
750: return([tmp[0],tmp[1],t_syz]);
751: }
752:
753: def Warning(s) {
754: Print("Warning: ");
755: Println(s);
756: }
757: def RingOf(f) {
758: local r;
759: if (IsPolynomial(f)) {
760: if (f != Poly("0")) {
761: sm1(f," (ring) dc /r set ");
762: }else{
763: sm1(" [(CurrentRingp)] system_variable /r set ");
764: }
765: }else{
766: Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
767: sm1(" [(CurrentRingp)] system_variable /r set ");
768: }
769: return(r);
770: }
771:
772: def Sfrom_es(f,size) {
773: local c,ans, i, d, myes, myee, j,n,r,ans2;
774: if (Length(Arglist) < 2) size = -1;
775: if (IsArray(f)) return(f);
776: r = RingOf(f);
777: myes = PolyR("es",r);
778: myee = PolyR("e_",r);
779: if (Degree(f,myee) > 0 && size == -1) {
780: if (size == -1) {
781: sm1(f," (array) dc /ans set");
782: return(ans);
783: }
784: }
785:
786: /*
787: Coefficients(x^2-1,x):
788: [ [ 2 , 0 ] , [ 1 , -1 ] ]
789: */
790: if (Degree(f,myee) > 0) {
791: c = Coefficients(f,myee);
792: }else{
793: c = Coefficients(f,myes);
794: }
795: if (size < 0) {
796: size = c[0,0]+1;
797: }
798: ans = NewArray(size);
799: for (i=0; i<size; i++) {ans[i] = 0;}
800: n = Length(c[0]);
801: for (j=0; j<n; j++) {
802: d = c[0,j];
803: ans[d] = c[1,j];
804: }
805: return(ans);
806: }
807:
808: def Sbases_to_vec(bases,size) {
809: local n, giveSize, newbases,i;
810: /* bases = [1+es*x, [1,2,3*x]] */
811: if (Length(Arglist) > 1) {
812: giveSize = true;
813: }else{
814: giveSize = false;
815: }
816: n = Length(bases);
817: newbases = NewArray(n);
818: for (i=0; i<n; i++) {
819: if (giveSize) {
820: newbases[i] = Sfrom_es(bases[i], size);
821: }else{
822: newbases[i] = Sfrom_es(bases[i]);
823: }
824: }
825: return(newbases);
826: }
827:
828: def Sminimal(g) {
829: local r, freeRes, redundantTable, reducer, maxLevel,
830: minRes, seq, maxSeq, level, betti, q, bases, dr,
831: betti_levelplus, newbases, i, j,qq;
832: r = SlaScala(g);
833: /* Should I turn off the tower?? */
834: freeRes = r[0];
835: redundantTable = r[1];
836: reducer = r[2];
837: minRes = SnewArrayOfFormat(freeRes);
838: seq = 0;
839: maxSeq = SgetMaxSeq(redundantTable);
840: maxLevel = Length(freeRes);
841: for (level = 0; level < maxLevel; level++) {
842: minRes[level] = freeRes[level];
843: }
844: seq=maxSeq+1;
845: while (seq > 1) {
846: seq--;
847: for (level = 0; level < maxLevel; level++) {
848: betti = Length(freeRes[level]);
849: for (q = 0; q<betti; q++) {
850: if (redundantTable[level,q] == seq) {
851: Print("[seq,level,q]="); Println([seq,level,q]);
852: if (level < maxLevel-1) {
853: bases = freeRes[level+1];
854: dr = reducer[level,q];
855: dr[q] = -1;
856: newbases = SnewArrayOfFormat(bases);
857: betti_levelplus = Length(bases);
858: /*
859: bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
860: */
861: for (i=0; i<betti_levelplus; i++) {
862: newbases[i] = bases[i] + bases[i,q]*dr;
863: }
864: Println(["level, q =", level,q]);
865: Println("bases="); sm1_pmat(bases);
866: Println("dr="); sm1_pmat(dr);
867: Println("newbases="); sm1_pmat(newbases);
868: minRes[level+1] = newbases;
869: freeRes = minRes;
870: #ifdef DEBUG
871: for (qq=0; qq<betti; qq++) {
872: if ((redundantTable[level,qq] >= seq) &&
873: (redundantTable[level,qq] <= maxSeq)) {
874: for (i=0; i<betti_levelplus; i++) {
875: if (!IsZero(newbases[i,qq])) {
876: Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
877: Print("redundantTable ="); sm1_pmat(redundantTable[level]);
878: Error("Stop in Sminimal for debugging.");
879: }
880: }
881: }
882: }
883: #endif
884: }
885: }
886: }
887: }
888: }
889: return([Stetris(minRes,redundantTable),
1.3 takayama 890: [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
1.1 takayama 891: /* r[4] is the redundantTable_ordinary */
1.3 takayama 892: /* r[0] is the freeResolution */
1.1 takayama 893: }
894:
895:
896: def IsZero(f) {
897: if (IsPolynomial(f)) {
898: return( f == Poly("0"));
899: }else if (IsInteger(f)) {
900: return( f == 0);
901: }else if (IsSm1Integer(f)) {
902: return( f == true );
903: }else if (IsDouble(f)) {
904: return( f == 0.0 );
905: }else if (IsRational(f)) {
906: return(IsZero(Denominator(f)));
907: }else{
908: Error("IsZero: cannot deal with this data type.");
909: }
910: }
911: def SgetMaxSeq(redundantTable) {
912: local level,i,n,ans, levelMax,bases;
913: levelMax = Length( redundantTable );
914: ans = 0;
915: for (level = 0; level < levelMax; level++) {
916: bases = redundantTable[level];
917: n = Length(bases);
918: for (i=0; i<n; i++) {
919: if (IsInteger( bases[i] )) {
920: if (bases[i] > ans) {
921: ans = bases[i];
922: }
923: }
924: }
925: }
926: return(ans);
927: }
928:
929: def Stetris(freeRes,redundantTable) {
930: local level, i, j, resLength, minRes,
931: bases, newbases, newbases2;
932: minRes = SnewArrayOfFormat(freeRes);
933: resLength = Length( freeRes );
934: for (level=0; level<resLength; level++) {
935: bases = freeRes[level];
936: newbases = SnewArrayOfFormat(bases);
937: betti = Length(bases); j = 0;
938: /* Delete rows */
939: for (i=0; i<betti; i++) {
940: if (redundantTable[level,i] < 1) {
941: newbases[j] = bases[i];
942: j++;
943: }
944: }
945: bases = SfirstN(newbases,j);
946: if (level > 0) {
947: /* Delete columns */
948: newbases = Transpose(bases);
949: betti = Length(newbases); j = 0;
950: newbases2 = SnewArrayOfFormat(newbases);
951: for (i=0; i<betti; i++) {
952: if (redundantTable[level-1,i] < 1) {
953: newbases2[j] = newbases[i];
954: j++;
955: }
956: }
957: newbases = Transpose(SfirstN(newbases2,j));
958: }else{
959: newbases = bases;
960: }
961: Println(["level=", level]);
962: sm1_pmat(bases);
963: sm1_pmat(newbases);
964:
965: minRes[level] = newbases;
966: }
967: return(minRes);
968: }
969:
970: def SfirstN(bases,k) {
971: local ans,i;
972: ans = NewArray(k);
973: for (i=0; i<k; i++) {
974: ans[i] = bases[i];
975: }
976: return(ans);
977: }
978:
979:
980: /* usage: tt is tower. ww is weight.
981: a = SresolutionFrameWithTower(v);
982: tt = a[1];
983: ww = [x,1,y,1,Dx,1,Dy,1];
984: SvDegree(x*es,tt,1,ww):
985:
986: In(17)=tt:
987: [[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 ] ,
988: [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ] ,
989: [3*es^3*y*Dy ] ]
990: In(18)=SvDegree(x*es,tt,1,ww):
991: 3
992: In(19)=SvDegree(x*es^3,tt,1,ww):
993: 4
994: In(20)=SvDegree(x,tt,2,ww):
995: 4
996:
997: */
998: def SvDegree(f,tower,level,w) {
999: local i,ans;
1000: if (IsZero(f)) return(null);
1.3 takayama 1001: f = Init(f);
1.1 takayama 1002: if (level <= 0) {
1003: return(Sord_w(f,w));
1004: }
1005: i = Degree(f,es);
1006: ans = Sord_w(f,w) +
1007: SvDegree(tower[level-1,i],tower,level-1,w);
1008: return(ans);
1009: }
1010:
1.2 takayama 1011: def Sannfs(f,v) {
1012: local f2;
1013: f2 = ToString(f);
1014: if (IsArray(v)) {
1015: v = Map(v,"ToString");
1016: }
1017: sm1(" [f2 v] annfs /FunctionValue set ");
1018: }
1019:
1020: /* Sannfs2("x^3-y^2"); */
1021: def Sannfs2(f) {
1022: local p,pp;
1023: p = Sannfs(f,"x,y");
1024: Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1025: pp = Map(p[0],"Spoly");
1026: return(Sminimal(pp));
1027: }
1028:
1.3 takayama 1029: def Sannfs3(f) {
1030: local p,pp;
1031: p = Sannfs(f,"x,y,z");
1032: Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
1033: pp = Map(p[0],"Spoly");
1034: return(Sminimal(pp));
1035: }
1036:
1.2 takayama 1037: /*
1038: The betti numbers of most examples are 2,1. (0-th and 1-th).
1039: a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2.
1040: a=Sannfs2("x^3-y^2-x"); : it causes an error. It should be fixed.
1.3 takayama 1041: a=Sannfs2("x*y*(x-y)"); : it causes an error. It should be fixed.
1.2 takayama 1042:
1043: */
1044:
1045:
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