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