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