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