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