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