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