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