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