Annotation of OpenXM/src/kan96xx/Kan/gb.c, Revision 1.14
1.14 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/gb.c,v 1.13 2005/07/07 02:59:47 takayama Exp $ */
1.5 takayama 2: /* untabify on May 4, 2001 */
1.1 maekawa 3: #include <stdio.h>
1.14 ! takayama 4: #include <string.h>
1.1 maekawa 5: #include "datatype.h"
6: #include "extern2.h"
7: #include "matrix.h"
8: #include "gradedset.h"
9:
10: #define INITGRADE 4
11: #define INITSIZE 2
1.6 takayama 12: #define DMAX 100
1.1 maekawa 13:
14: int KanGBmessage = 1;
15:
16: static int Debug = 0;
17: static int sugarGrade(struct pair *inode,struct gradedPolySet *grG,
1.5 takayama 18: POLY gt,int gtSugarGrade);
1.1 maekawa 19: extern int Sugar;
20: extern int Statistics;
21: extern int Criterion1;
22: extern int UseCriterion1;
23: extern int UseCriterion2B;
24: extern int Spairs;
25: extern int Criterion2B, Criterion2F, Criterion2M;
1.3 takayama 26: extern int AutoReduce;
1.8 takayama 27: extern int TraceLift;
28: extern struct ring *TraceLift_ringmod;
1.6 takayama 29: static int MaxLength[DMAX];
30: static int SpNumber[DMAX];
1.1 maekawa 31:
32: struct gradedPairs *updatePairs(grD,gt,gtGrade,t,grG)
1.5 takayama 33: struct gradedPairs *grD; /* set of pairs */
34: POLY gt; /* new polynomial */
35: int gtGrade;
36: int t;
37: struct gradedPolySet *grG; /* (f,gt), f \in grG, should be added to grD. */
1.1 maekawa 38: {
39: int gmax,newGrade;
40: struct pair *node,*new,*inode,*jnode;
41: int i,k;
42: struct polySet *g;
43: POLY lcmp;
44: POLY it,jt;
45: extern int Verbose;
46:
47:
48: if (Verbose) {
49: printf("\nupdatePairs(): ");
50: }
51:
52: gmax = grG->maxGrade;
53: node = new = newPair((struct pair *)NULL);
54: /* (f,gt) are stored in new as a list. */
55: for (k=0; k<gmax; k++) {
56: g = grG->polys[k];
57: for (i=0; i<g->size; i++) {
58: if (g->del[i] == 0 && (gtGrade != k || t != i)) {
1.5 takayama 59: lcmp = (*lcm)(g->g[i],gt);
60: if (lcmp == ZERO) {
1.1 maekawa 61:
1.5 takayama 62: }else {
63: new->del = 0;
64: new->next = newPair(new);
65: new = new->next;
66: new->lcm = lcmp;
67: new->ig = k; new->ii = i; /* g->g[i] */
68: new->jg = gtGrade; new->ji = t; /* gt */
69: new->grade = -1; /* one do not need to set grade here. */
70: }
1.1 maekawa 71: }
72: }
73: }
74:
75:
76: /* Check the criterion */
77: /* See Gebauer and Mora, 1988, Journal of Symbolic Computation
78: 279. We must not use the test T(i)T(j) = T(i,j) because it
79: does not hold in the ring of diff op.
1.5 takayama 80: */
1.1 maekawa 81:
82: inode = node->next;
83: while (inode != (struct pair *)NULL) {
84: jnode = inode->next;
85: while (jnode != (struct pair *)NULL) {
86: it = inode->lcm;
87: jt = jnode->lcm;
88: if ((*mmLarger)(it,jt) == 2) {
1.5 takayama 89: /* F(j,t), i<j */
90: jnode->del = 1;
91: if (Verbose) printf("F[%d,%d]([%d,%d],[%d,%d]) ",
92: inode->ig,inode->ii,
93: jnode->ig,jnode->ii,
94: gtGrade,t);
95: Criterion2F++;
1.1 maekawa 96: }else {
1.5 takayama 97: /* g[i] > g[j] ?, M(i,t) */
98: if ((*isReducible)(it,jt)) {
99: inode->del = 1;
100: if (Verbose) printf("M[%d,%d]([%d,%d],[%d,%d]) ",
101: jnode->ig,jnode->ii,
102: inode->ig,inode->ii,
103: gtGrade,t);
104: Criterion2M++;
105: }
1.1 maekawa 106: /* M(j,t) */
1.5 takayama 107: if ((*isReducible)(jt,it)) {
108: jnode->del = 1;
109: if (Verbose) printf("M[%d,%d]([%d,%d],[%d,%d]) ",
110: inode->ig,inode->ii,
111: jnode->ig,jnode->ii,
112: gtGrade,t);
113: Criterion2M++;
114: }
1.1 maekawa 115: }
116: jnode = jnode->next;
117: }
118: inode = inode->next;
119: }
1.5 takayama 120:
1.1 maekawa 121: if (UseCriterion1) {
122: inode = node->next;
123: while (inode != NULL) {
124: if (inode->del == 0) {
1.5 takayama 125: if (criterion1(gt,grG->polys[inode->ig]->g[inode->ii],inode->lcm)) {
126: inode->del = 1;
127: if (Verbose) printf("Criterion1([%d,%d],[%d,%d]) ",
128: inode->ig,inode->ii,
129: gtGrade,t);
130: Criterion1++;
131: }
1.1 maekawa 132: }
133: inode = inode->next;
134: }
135: }
136:
137: if (UseCriterion2B) {
138: Criterion2B += deletePairByCriterion2B(grD,gt,grG);
139: }
140:
141: /* Merge to grD */
142: inode = node->next;
143: if (Debug) outputNode(inode);
144: while (inode != (struct pair *)NULL) {
145: if (inode->del == 0) {
146: if (Sugar) {
1.5 takayama 147: inode->grade = sugarGrade(inode,grG,gt,gtGrade);
1.1 maekawa 148: }else{
1.5 takayama 149: inode->grade = (*grade)(inode->lcm);
1.1 maekawa 150: }
151: if (grD->lim <= inode->grade) {
1.5 takayama 152: grD = enlargeGradedPairs(2*(inode->grade)+1,grD);
1.1 maekawa 153: }
154: insertPair(inode,grD->pairs[inode->grade]);
155: grD->maxGrade = max(grD->maxGrade,inode->grade+1); /* don't forget */
156: }
157: inode = inode->next;
158: }
159: if (Debug) printf("OK.\n");
160: if (Verbose) {
161: printf(" Remaining pairs are %d. maxGrade=%d\n",countPairs(grD),grD->maxGrade);
162: }
163: return(grD);
164: }
165:
1.10 takayama 166: struct gradedPolySet *groebner_gen(f,needBack,needSyz,grP,countDown,forceReduction,reduceOnly,gbCheck)
1.5 takayama 167: struct arrayOfPOLY *f;
168: int needBack;
169: int needSyz;
170: struct pair **grP; /* if (needSyz), it is set. */
171: int countDown;
172: int forceReduction;
1.9 takayama 173: int reduceOnly;
1.10 takayama 174: int gbCheck;
1.1 maekawa 175: {
176: int r;
177: struct gradedPolySet *g;
178: struct gradedPairs *d;
179: int i;
180: int grade,indx;
181: POLY gt;
182: struct pair *top;
183: int ig,ii,jg,ji;
184: POLY gi,gj;
185: struct spValue h;
186: struct syz0 syz;
187: int pgrade = 0;
188: POLY rd;
189: POLY syzPoly;
190: POLY syzCf;
191: struct syz0 *syzp;
192: int serial;
193: struct pair *listP;
194: extern int Verbose;
195: extern int ReduceLowerTerms;
196: extern int StopDegree;
197: extern struct ring *CurrentRingp;
198: extern char *F_mpMult;
1.13 takayama 199: extern int Homogenize;
1.1 maekawa 200: struct ring *rp;
1.8 takayama 201: int first;
1.6 takayama 202: int statisticsPL, statisticsCount;
1.1 maekawa 203:
1.12 takayama 204: if (KanGBmessage) {
1.13 takayama 205: printf("Homogenize=%d, countDown=%d, forceReduction=%d, reduceOnly=%d, gbCheck=%d\n",
206: Homogenize,countDown, forceReduction, reduceOnly, gbCheck);
1.12 takayama 207: }
1.6 takayama 208: if (Statistics) {
209: for (i=0; i<DMAX; i++) MaxLength[i] = SpNumber[i] = 0;
210: }
1.1 maekawa 211: r = f->n;
212: if (r<=0) return((struct gradedPolySet *)NULL);
213: if (UseCriterion1) {
214: if (needSyz) {
215: warningGradedSet("You cannot use both the options needSyz and UseCriterion1.\nInput [(UseCriterion1) 0] system_variable to turn off the use of the Criterion 1.");
216: }
217: if (strcmp(F_mpMult,"poly")!=0) {
218: warningGradedSet("You can use the option UseCriterion1 only for the computation in the ring of polynomials.\nInput [(UseCriterion1) 0] system_variable to turn off the use of the Criterion 1.");
219: }
220: }
221:
222: Spairs = Criterion1 = Criterion2B = Criterion2F = Criterion2M = 0;
223:
1.10 takayama 224: g = newGradedPolySet(INITGRADE); g->gb = 1;
1.11 takayama 225: if (reduceOnly) g->gb = 2; /* unknown */
1.1 maekawa 226: d = newGradedPairs(INITGRADE*2);
227: for (i=0; i<g->lim; i++) {
228: g->polys[i] = newPolySet(INITSIZE);
229: }
230:
1.8 takayama 231: first = 1;
1.1 maekawa 232: for (i=0; i<r; i++) {
233: gt = getArrayOfPOLY(f,i);
234: if (gt ISZERO) { rp = CurrentRingp; } else { rp = gt->m->ringp; }
1.8 takayama 235: if (TraceLift && (!(gt ISZERO)) && first) {
236: TraceLift_ringmod = newRingOverFp(rp,getPrime(TraceLift)); first = 0;
237: if (KanGBmessage) printf("Prime number for the trace lift is %d.\n",
238: TraceLift_ringmod->p);
239: }
1.1 maekawa 240: grade = -1; whereInG(g,gt,&grade,&indx,Sugar);
241: if (KanGBmessage == 2) {
242: printf("init=%s, ",POLYToString(head(gt),'*',1));
243: printf("(gr,indx)=(%d,%d).\n",grade,indx);
244: }
245: d = updatePairs(d,gt,grade,indx,g);
246: serial = i;
247: if (!needBack) {
248: g = putPolyInG(g,gt,grade,indx,
1.5 takayama 249: (struct syz0 *)NULL,1,serial);
1.1 maekawa 250: }else {
251: syzp = newSyz0();
252: syzp->cf = cxx(1,0,0,rp);
253: syzp->syz = toSyzPoly(cxx(1,0,0,rp),grade,indx);
254: g = putPolyInG(g,gt,grade,indx,syzp,1,serial);
255: }
256:
1.7 takayama 257: /* markRedundant0(g,grade,indx); ?*/
258: markGeneratorInG(g,grade,indx); /*?*/
1.1 maekawa 259: if (Debug) {
260: outputGradedPairs(d); outputGradedPolySet(g,needSyz);
261: }
262: }
263: if (needSyz) {
264: *grP = newPair((struct pair *)NULL);
265: listP = *grP;
266: }
267:
268: while ((top = getPair(d)) != (struct pair *)NULL) {
269: ig = top->ig; ii = top->ii; /* [ig,ii] */
270: jg = top->jg; ji = top->ji; /* [jg,ji] */
1.7 takayama 271: /*
272: if (g->polys[ig]->del[ii] || g->polys[jg]->del[ji]) {
273: if (KanGBmessage) printf("p");
274: continue;
275: } Don't do this.
276: */
1.1 maekawa 277: gi = g->polys[ig]->g[ii];
278: gj = g->polys[jg]->g[ji];
279:
280: Spairs++;
1.9 takayama 281: if (reduceOnly && (!needSyz) && (!needBack)) h = spZero(); /* rd = 0 */
282: else h = (*sp)(gi,gj);
1.1 maekawa 283: rd = ppAddv(ppMult(h.a,gi),ppMult(h.b,gj));
1.6 takayama 284:
285: if (Statistics) {
286: if (top->grade >=0 && top->grade < DMAX) {
287: statisticsPL = pLength(rd);
288: SpNumber[top->grade]++;
289: if (MaxLength[top->grade] < statisticsPL) {
290: MaxLength[top->grade] = statisticsPL;
291: }
292: }
293: }
294:
1.1 maekawa 295: if (!Sugar || forceReduction) {
296: rd = (*reduction)(rd,g,needBack,&syz);
297: }else{
298: rd = reduction_sugar(rd,g,needBack,&syz,top->grade);
299: }
300: syzPoly = syz.syz;
301: syzCf = syz.cf;
302:
303: if (KanGBmessage) {
304: if (KanGBmessage == 2) {
1.5 takayama 305: printf("sp([%d,%d],[%d,%d]), ",ig,ii,jg,ji);
306: if (rd ISZERO) {
307: printf(" 0 \n");
308: } else{
309: printf(" terms=%d, grade=%d, ",pLength(rd),top->grade);
310: printf(" init=%s, ",POLYToString(head(rd),'*',1));
311: }
1.1 maekawa 312: }else{
1.5 takayama 313: if (pgrade != top->grade) {
314: pgrade = top->grade;
315: printf(" %d",pgrade);
316: fflush(stdout);
317: }
318: if (rd ISZERO) {
319: printf("o"); fflush(stdout);
320: }else{
321: printf("."); fflush(stdout);
322: }
1.1 maekawa 323: }
324: }
325:
326: if (!(rd ISZERO)) {
1.10 takayama 327: if (gbCheck) {
328: /* Abort the calculation. */
329: g->gb = 0;
330: if (KanGBmessage) {
331: printf("gbCheck failed. \n");
332: printf("Note that the result is NOT groebner basis.\n");
333: }
334: break;
335: }
1.1 maekawa 336: if (needBack || needSyz) {
1.5 takayama 337: syzp = newSyz0();
338: syzp->cf = syzCf; /* no meaning */
339: syzp->syz = ppAdd(toSyzPoly(h.a,ig,ii),toSyzPoly(h.b,jg,ji));
340: syzp->syz = cpMult(toSyzCoeff(syzCf),syzp->syz);
341: syzp->syz = ppAdd(syzp->syz,syzPoly);
1.1 maekawa 342: }
343:
344: if (ReduceLowerTerms && !(Sugar)) {
1.5 takayama 345: rd = (*reductionCdr)(rd,g,needBack,&syz);
346: if (needBack || needSyz) {
347: /* syzp->cf = ppMult(syz.cf,syzp->cf); no meaning */
348: syzp->syz = ppAdd(syz.syz,
349: cpMult(toSyzCoeff(syz.cf),syzp->syz));
350: }
1.1 maekawa 351: }
352:
353: if(Sugar && (!forceReduction)){grade=top->grade;}else{grade=-1;}whereInG(g,rd,&grade,&indx,Sugar);
354: if (KanGBmessage == 2) {
1.5 takayama 355: printf("(gr,indx)=(%d,%d).\n",grade,indx);
356: /*
357: printf("sp(%s,%s)-->%s\n",POLYToString(gi,' ',1),
358: POLYToString(gj,' ',1),
359: POLYToString(rd,' ',1));
360: */
1.1 maekawa 361: }
1.5 takayama 362:
1.1 maekawa 363: d = updatePairs(d,rd,grade,indx,g);
364: g = putPolyInG(g,rd,grade,indx,syzp,0,-1);
365: if (Sugar) { markRedundant0(g,grade,indx);}
366: else {markRedundant(g,rd,grade,indx,Sugar);}
367:
368: if (countDown) {
1.5 takayama 369: if (eliminated(rd) == 1) {
370: --countDown;
371: printf("x"); fflush(stdout);
372: if (countDown == 0) {
373: printf("\nThe computation of the Groebner basis is suspended because of countDown==0.\n");
374: printf("Note that the result is NOT groebner basis.\n");
1.10 takayama 375: g->gb = 0;
1.5 takayama 376: break;
377: }
378: }
1.1 maekawa 379: }
380: if (Debug) {
1.5 takayama 381: outputGradedPairs(d); outputGradedPolySet(g,needSyz);
1.1 maekawa 382: }
383: }else{
384: if (needSyz) {
1.5 takayama 385: top->syz = ppAdd(toSyzPoly(h.a,ig,ii),toSyzPoly(h.b,jg,ji));
386: top->syz = cpMult(toSyzCoeff(syzCf),top->syz);
387: top->syz = ppAdd(top->syz,syzPoly);
388: listP->next = top; top->prev = listP; listP = listP->next;
389: }
390: }
391: if (StopDegree < pgrade) {
392: fprintf(stderr,"Obtained a partial GB (StopDegree=%d)\n",StopDegree);
393: if (KanGBmessage) {
394: printf("Computation of the Groebner basis is suspended bacause of StopDegree < computing grade.\n");
395: printf("Note that the result is NOT groebner basis.\n");
1.1 maekawa 396: }
1.10 takayama 397: g->gb = 0;
1.5 takayama 398: break;
1.1 maekawa 399: }
400: }
401:
402: if (KanGBmessage == 2) {
403: outputGradedPolySet(g,needSyz);
404: }
405: if (KanGBmessage) {
406: if (Sugar) {
407: printf("\nCompleted (GB with sugar).\n");
408: }else{
409: printf("\nCompleted.\n");
410: }
411: }
412:
413: if (Statistics) {
414: printf("\nThe number of s-pairs is %d.\n",Spairs);
415: printf("Criterion1 is applied %d times.\n",Criterion1);
416: printf("Criterions M,F and B are applied M=%d, F=%d, B=%d times.\n",Criterion2M,Criterion2F,Criterion2B);
417: Spairs = Criterion1 = Criterion2M = Criterion2F = Criterion2B = 0;
1.6 takayama 418:
419: printf("degree(number of spolys): maximal polynomial size\n");
420: statisticsCount = 0;
421: for (i=0; i<DMAX; i++) {
422: if (SpNumber[i] > 0) {
423: printf("%3d(%3d): %5d, ",i,SpNumber[i],MaxLength[i]);
424: if (statisticsCount % 4 == 3) {
425: printf("\n");
426: statisticsCount = 0;
427: }else{ statisticsCount++; }
428: }
429: }
430: printf("\n");
1.1 maekawa 431: }
432:
1.9 takayama 433: if (AutoReduce || reduceOnly) {
1.3 takayama 434: toReducedBasis(g,needBack,needSyz);
435: }
1.1 maekawa 436:
437: return(g);
438: }
439:
440:
441: static int sugarGrade(struct pair *inode,struct gradedPolySet *grG,
1.5 takayama 442: POLY gt,int gtSugarGrade)
1.1 maekawa 443: {
444: int a,b,ans;
445: int debug = 0;
446: a = grade_gen(inode->lcm)-grade_gen(grG->polys[inode->ig]->g[inode->ii]);
447: b = grade_gen(inode->lcm)-grade_gen(gt);
448: /* inode = lcm(f_i, gt) = p f_i = q gt modulo lower order terms.
449: a = tdeg(p), b = tdeg(gt);
1.5 takayama 450: */
1.1 maekawa 451: if (debug) {
452: printf("Sugar grade of sp([%d,%d],[%d,%d]) ",inode->ig,inode->ii,
1.5 takayama 453: inode->jg,inode->ji);
1.1 maekawa 454: printf("is max(%d+%d,%d+%d)\n",a,inode->ig,b,gtSugarGrade);
455: }
456: a = a+(inode->ig); /* new sugar degree of p f_i. */
457: b = b + gtSugarGrade; /* new sugar degree of q gt */
458: return( a > b ? a : b);
459: }
460:
1.3 takayama 461: void toReducedBasis(struct gradedPolySet *grP,int needBack, int needSyz)
462: {
463: int changed, grd, i, reduced, grade,indx;
464: struct syz0 syz;
465: struct syz0 *syzp;
466: POLY f;
467: POLY rd;
468: struct polySet *set;
469:
470: /* KanGBmessage=1; */
471: do {
472: if (KanGBmessage) {
473: printf("s"); fflush(stdout);
474: }
475: changed = 0;
476: grd = 0;
477: while (grd < grP->maxGrade) {
478: set = grP->polys[grd];
479: for (i=0; i<set->size; i++) {
1.5 takayama 480: if (set->del[i] == 0) {
481: f = set->g[i];
482: if (KanGBmessage) {
483: /* printf("(%d,%d)",grd,i); */
484: fflush(stdout);
485: }
486: rd = reductionCdr_except_grd_i(f,grP,needBack,&syz,grd,i,&reduced);
487: if (KanGBmessage) {
488: if (reduced) {
489: printf(".");
490: }else{
491: printf("o");
492: }
493: fflush(stdout);
494: }
495: if (reduced) {
496: changed = 1;
497: set->del[i] = 1;
498: if (rd != ZERO) {
499: if (needSyz) {
500: syzp = newSyz0();
501: syzp->cf = syz.cf; /* no meaning */
502: syzp->syz = toSyzPoly(cxx(1,0,0,rd->m->ringp),grd,i);
503: syzp->syz = cpMult(toSyzCoeff(syz.cf),syzp->syz);
504: syzp->syz = ppAdd(syzp->syz,syz.syz);
505: /* rd = c*f + \sum c_{d,i} g_{d,i}
506: c : syz.cf, \sum c_{d,j} g_{d,j} : syz.syz.
507: c*grade^grd*index^i + \sum c_{d,j} grade^d*index^j
508: grP is a set of polynomials. Polynomials are indexed by
509: grade and index.
510: */
511: /* printf("%s, ",POLYToString(syzp->cf,' ',1));
512: printf("%s\n",POLYToString(syzp->syz,' ',1)); */
513: }else{
514: syzp = NULL;
515: }
516: grade = -1; whereInG(grP,rd,&grade,&indx,Sugar);
517: /* Do not forget to set grade to -1 */
518: /* printf("grade=%d, indx=%d, ",grade,indx); */
519: putPolyInG(grP,rd,grade,indx,syzp,0,-1);
520: }
521: }
522: }
1.3 takayama 523: }
524: grd++;
525: }
526: } while(changed);
527: if (KanGBmessage) {
528: printf("Done(reduced basis)\n");
529: }
1.10 takayama 530: grP->reduced = 1;
1.3 takayama 531: }
1.1 maekawa 532:
533:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gb.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / kan96xx / Kan