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