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