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