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