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