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