Annotation of OpenXM/src/kan96xx/Kan/gbGM.c, Revision 1.5
1.5 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/gbGM.c,v 1.4 2005/06/09 04:09:22 takayama Exp $ */
1.1 maekawa 2: /* gbGM.c GM=Gebauer and Moller
1.3 takayama 3: */
1.1 maekawa 4:
5: #include <stdio.h>
6: #include "datatype.h"
7: #include "extern2.h"
8: #include "matrix.h"
9: #include "gradedset.h"
10:
11:
12: /*#define DEBUG*/
13: /* Take statistics of mmLarger_matrix in orderT.c */
14: /* #define STATISTICS */
15:
16: struct polySet_gm {
17: POLY *g;
18: int lim;
19: int size;
20: int *del;
21: };
22:
23: struct pair_gm {
24: POLY lcm;
25: int i;
26: int j;
27: };
28:
29: struct pairSet {
30: struct pair_gm *p;
31: int lim;
32: int size;
33: int *del;
34: };
35:
36: /* prototypes */
37: struct polySet_gm newPolySet_gm(int n);
38: struct pairSet newPairSet(int n);
39: int pairSetSize(struct pairSet d);
40: struct pairSet pairSetJoin(struct pairSet a,struct pairSet b);
41: struct polySet_gm enlargePolySet_gm(struct polySet_gm g);
42: struct pairSet deletePair_gm(struct pairSet d,int index);
43: int minPair_gm(struct pairSet d);
44: struct pairSet updatePair_gm(struct pairSet d,int t,
1.3 takayama 45: struct polySet_gm g,POLY gt);
1.1 maekawa 46: struct polySet_gm markRedundant_gm(struct polySet_gm g,int j);
47: /*
1.3 takayama 48: struct gradedPolySet *groebner_gm(struct arrayOfPOLY f,
49: int needBack,
50: int needSyz, struct pair **grP);
1.1 maekawa 51: */
52: void outputPairSet(struct pairSet d);
53: void outputPolySet_gm(struct polySet_gm g);
54: void errorGBGM(char *s);
55: POLY mReductionBySet(POLY f, POLY *s,int size);
56: POLY mReductionCdrBySet(POLY f, POLY *s,int size);
57: int criterion1(POLY f,POLY g,POLY lc);
58:
59: extern int UseCriterion1;
60: extern int Verbose;
61:
62: #define GLIMIT 100
63: int Criterion1,Criterion2; /* for taking a statistics. */
64: int Spairs;
65: extern int Criterion2B, Criterion2F, Criterion2M;
66: extern int Statistics;
67:
68: #ifdef STATISTICS
69: int CountE;
70: int CountI[2*N0];
71: #endif
72:
73: struct polySet_gm newPolySet_gm(n)
1.3 takayama 74: int n;
1.1 maekawa 75: {
76: struct polySet_gm g;
77: int i;
78: g.g = (POLY *) sGC_malloc(sizeof(POLY)*(n+1));
79: g.del = (int *) sGC_malloc(sizeof(int)*(n+1));
80: if ((g.g == (POLY *)NULL) ||
81: (g.del == (int *)NULL))
82: errorGBGM("No more memory.");
83: g.lim = n;
84: g.size = 0;
85: for (i=0; i<n; i++) g.del[i] = 0;
86: return(g);
87: }
88:
89: struct pairSet newPairSet(n)
1.3 takayama 90: int n;
1.1 maekawa 91: {
92: struct pairSet g;
93: int i;
94: g.p = (struct pair_gm *) sGC_malloc(sizeof(struct pair_gm)*(n+1));
95: g.del = (int *) sGC_malloc(sizeof(int)*(n+1));
96: if ((g.p == (struct pair_gm *)NULL) ||
97: (g.del == (int *)NULL))
98: errorGBGM("No more memory.");
99: g.lim = n;
100: g.size = 0;
101: for (i=0; i<n; i++) g.del[i] = 0;
102: return(g);
103: }
104:
105: int pairSetSize(d)
1.3 takayama 106: struct pairSet d;
1.1 maekawa 107: {
108: int s,i;
109: s = 0;
110: for (i=0; i< d.size; i++) {
111: if (d.del[i] == 0) s++;
112: }
113: return(s);
114: }
115:
116: struct pairSet pairSetJoin(a,b)
1.3 takayama 117: struct pairSet a,b;
1.1 maekawa 118: {
119: int m,n,k,i;
120: struct pairSet ans;
121: m = pairSetSize(a);
122: n = pairSetSize(b);
123: ans = newPairSet(m+n);
124: k = 0;
125: for (i=0 ; i<a.size; i++) {
126: if (a.del[i] == 0) {
127: ans.p[k] = a.p[i]; ans.del[k++] = 0;
128: }
129: }
130: for (i=0; i<b.size; i++) {
131: if (b.del[i] == 0) {
132: ans.p[k] = b.p[i]; ans.del[k++] = 0;
133: }
134: }
135: if (k != (m+n)) errorGBGM("Internal error in pairSetJoin().");
136: ans.size = m+n;
137: return(ans);
138: }
139:
140: struct polySet_gm enlargePolySet_gm(g)
1.3 takayama 141: struct polySet_gm g;
1.1 maekawa 142: {
143: int i;
144: struct polySet_gm ans;
145: if (g.lim < g.size+3) {
146: ans = newPolySet_gm((g.lim)*2);
147: for (i=0; i<g.size; i++) {
148: ans.g[i] = g.g[i];
149: ans.del[i] = g.del[i];
150: }
151: ans.size = g.size;
152: }else{ ans = g; }
153: return(ans);
154: }
155:
156: struct pairSet deletePair_gm(d,index)
1.3 takayama 157: struct pairSet d;
158: int index;
159: /* delete d[index] */
1.1 maekawa 160: {
161: int i;
162: d.del[index] = 0;
163: for (i= index; i<d.size-1; i++) {
164: d.p[i] = d.p[i+1];
165: }
166: d.size--;
167: return(d);
168: }
169:
170: int minPair_gm(d)
1.3 takayama 171: struct pairSet d;
1.1 maekawa 172: {
173: POLY min;
174: int index,i;
175: min = d.p[0].lcm;
176: index = 0;
177: for (i=0; i < d.size; i++) {
178: if ( (*mmLarger)(min, d.p[i].lcm)) {
179: index = i;
180: min = d.p[i].lcm;
181: }
182: }
183: return(index);
184: }
185:
186: struct pairSet updatePair_gm(d,t,g,gt)
1.3 takayama 187: struct pairSet d;
188: int t;
189: struct polySet_gm g;
190: POLY gt;
1.1 maekawa 191: {
192: int i,j,k;
193: struct pairSet new;
194: POLY it;
195: POLY jt;
196:
197: if (Verbose >= 2) {
198: printf("\nupdatepair_gm(): ");
199: printf("g.size=%d ",g.size);
200: printf("d.size=%d; ",d.size);
201: }
202:
203: /* make list (1,t), (2,t), .... ====> new ( D1 ) */
204: new = newPairSet(g.size); new.size = 0;
205: for (i=0; i<g.size; i++) {
206: if (g.del[i] == 0) {
207: new.p[i].i = i;
208: new.p[i].j = t;
209: new.p[i].lcm = (*lcm)(g.g[i],gt);
210: }else{
211: new.p[i].i = i;
212: new.p[i].j = t;
213: new.p[i].lcm = ZERO; /* Compute it when it is necessary. */
214: }
215: new.del[i] = g.del[i];
216: new.size++;
217: }
218: #ifdef DEBUG
219: /*printf("\nnew is ...");
1.3 takayama 220: outputPairSet(new);*/
1.1 maekawa 221: #endif
222:
223: /* Candel in D=d all (i,j) such that B_t(i,j) */
224:
225: for (k=0; k<d.size; k++) {
226: if ((d.del[k] == 0) && ((*isReducible)(d.p[k].lcm,gt)) ) {
227: /* check T(it) != T(i,j) != T(j,t) */
228: i = d.p[k].i; j = d.p[k].j;
229: if ((new.del[i] == 1) || (new.del[j] == 1)) {
1.3 takayama 230: /* fprintf(stderr,"Warning in updatepair_gm(): i=%d, j=%d; rewriting new.\n",i,j); */
231: if (new.del[i] == 1) {
232: new.p[i].lcm = (*lcm)(g.g[i],gt);
233: }
234: if (new.del[j] == 1) {
235: new.p[j].lcm = (*lcm)(g.g[j],gt);
236: }
1.1 maekawa 237: }
238: if (((*mmLarger)((new.p[i].lcm),(new.p[j].lcm)) != 2) &&
1.3 takayama 239: ((*mmLarger)((new.p[i].lcm),(d.p[k].lcm)) != 2) &&
240: ((*mmLarger)((new.p[j].lcm),(d.p[k].lcm)) != 2)) {
241: /* delete T(i,j) in d */
242: d.del[k] = 1;
243: if (Verbose >= 2) printf("B%d(%d,%d) ",t,i,j);
244: Criterion2B++;
1.1 maekawa 245: }
246: }
247: }
248:
249: /* Then, d is D' */
250: /* Next, check M(i,t) and F(i,t) */
251: for (i=0; i<g.size; i++) {
252: for (j=i+1; j<g.size; j++) {
253: if ((g.del[i] == 0) && (g.del[j] == 0)) {
1.3 takayama 254: it = new.p[i].lcm;
255: jt = new.p[j].lcm;
256: switch ( (*mmLarger)(it,jt)) {
257: case 2:
258: /* F(j,t), i<j */
259: new.del[j] = 1;
260: if (Verbose >= 2) printf("F%d(%d,%d) ",i,j,t);
261: Criterion2F++;
262: break;
263: case 1:
264: /* g[i] > g[j], M(i,t) */
265: if ((*isReducible)(it,jt)) {
266: new.del[i] = 1;
267: if (Verbose >=2) printf("M%d(%d,%d) ",j,i,t);
268: Criterion2M++;
269: }
270: break;
271: case 0: /* M(j,t) */
272: if ((*isReducible)(jt,it)) {
273: new.del[j] = 1;
274: if (Verbose >=2) printf("M%d(%d,%d) ",i,j,t);
275: Criterion2M++;
276: }
277: break;
1.1 maekawa 278:
1.3 takayama 279: }
1.1 maekawa 280: }
281: }
282: }
283:
284: /* Finally, check T(i) T(t) = T(i,t) */
285: if (UseCriterion1) {
286: for (i=0; i<g.size; i++) {
287: if (new.del[i] == 0) {
1.3 takayama 288: if (criterion1(g.g[i],gt,new.p[i].lcm)) {
289: new.del[i] = 1;
290: if (Verbose >=2) printf("1(%d,%d) ",i,t);
291: Criterion1++;
292: }
1.1 maekawa 293: }
294: }
295: }
296:
297: /* d <-- D' \cup new */
298: d = pairSetJoin(d,new);
299: if (Verbose >=2) printf("new d.size = %d.\n",d.size);
300: #ifdef DEBUG
301: outputPairSet(d);
302: #endif
303: return(d);
304: }
305:
306:
307: struct polySet_gm markRedundant_gm(g,j)
1.3 takayama 308: struct polySet_gm g;
309: int j;
310: /* compare only with g[j] */
1.1 maekawa 311: {
312: int i;
313: for (i=0; i<g.size; i++) {
314: if ((g.del[i] == 0) && (i != j)) {
315: if (g.del[j] == 0) {
1.3 takayama 316: switch((*mmLarger)((g.g[i]),(g.g[j]))) {
317: case 2:
318: g.del[j] = 1;
319: break;
320: case 1:
321: if ((*isReducible)((g.g[i]),(g.g[j]))) {
322: g.del[i] = 1;
323: }
324: break;
325: case 0:
326: if ((*isReducible)((g.g[j]),(g.g[i]))) {
327: g.del[j] = 1;
328: }
329: break;
330: }
1.1 maekawa 331: }
332: }
333: }
334: #ifdef DEBUG
335: /*printf("\nend of markRedundant_gm...");
1.3 takayama 336: outputPolySet_gm(g);*/
1.1 maekawa 337: #endif
338: return(enlargePolySet_gm(g));
339: }
340:
341:
342:
1.5 ! takayama 343: struct gradedPolySet *groebner_gm(f,needBack,needSyz,grP,countDown,forceReduction,reduceOnly,gbCheck)
1.3 takayama 344: struct arrayOfPOLY *f;
345: int needBack;
346: int needSyz;
347: struct pair **grP;
348: int countDown;
349: int forceReduction;
1.4 takayama 350: int reduceOnly;
1.5 ! takayama 351: int gbCheck;
1.1 maekawa 352: {
353: int r;
354: struct pair_gm top;
355: int tindex;
356: POLY h;
357: int i,j;
358: int ss;
359: struct polySet_gm g;
360: struct pairSet d;
361: int t;
362: struct spValue sv;
363: struct gradedPolySet *ans;
364: int grade,indx;
365: extern int ReduceLowerTerms;
366:
367: if (needBack || needSyz) {
368: fprintf(stderr,"Warning: groebner_gm() does not compute the backward transformation and syzygies.\n");
1.4 takayama 369: }
1.5 ! takayama 370: if (reduceOnly || gbCheck) {
! 371: fprintf(stderr,"Warning: groebner_gm() does not implement reduceOnly/gbCheck.\n");
1.1 maekawa 372: }
373:
374: #ifdef STATISTICS
375: CountE = 0;
376: for (i=0; i<2*N; i++) CountI[i]=0;
377: #endif
378: /* taking the statistics. */
379: Criterion1 = Criterion2B = Criterion2M = Criterion2F = 0;
380: Spairs = 0;
381:
382: r = f->n;
383: g = newPolySet_gm(r + GLIMIT);
384: d = newPairSet(1);
385:
386: g.size = 0;
387: g.g[g.size] = getArrayOfPOLY(f,g.size);
388: g.del[g.size] = 0; g.size ++;
389:
390: for (t=1; t<r; t++) {
391: d = updatePair_gm(d,t,g,getArrayOfPOLY(f,t));
392: g.g[t] = getArrayOfPOLY(f,t); g.del[t] = 0; g.size ++ ;
393: g = markRedundant_gm(g,t);
394: }
395:
396: #ifdef DEBUG
397: outputPolySet_gm(g);
398: #endif
399:
400: while (d.size > 0) {
401: tindex = minPair_gm(d);
402: top = d.p[tindex];
403: if (Verbose) {
404: printf("\nRemaining pairs = %d : ",d.size);
405: }
406: i = top.i;
407: j = top.j;
408: if (Verbose >=2) printf("sp(%d,%d) ",i,j);
409: Spairs++;
410: sv = (*sp)(g.g[i],g.g[j]);
411: h = ppAddv(ppMult(sv.a,g.g[i]),ppMult(sv.b,g.g[j]));
412: h = mReductionBySet(h,g.g,g.size);
413: if (Verbose >=2) printf("->%s ",POLYToString(h,' ',0));
414:
415: if (!Verbose) {
416: printf(" <%d>",d.size);
417: fflush(stdout);
418: }
419:
420: if (!(h ISZERO)) {
421: if (ReduceLowerTerms) {
1.3 takayama 422: h = mReductionCdrBySet(h,g.g,g.size);
1.1 maekawa 423: }
424: d = updatePair_gm(d,r,g,h);
425: g.g[r] = h; g.del[r] = 0;
426: r++; g.size++;
427: g = markRedundant_gm(g,r-1);
428: }
429: d = deletePair_gm(d,tindex);
430: #ifdef DEBUG
431: outputPolySet_gm(g);
432: #endif
433: }
434:
435:
436: if (Verbose || Statistics) {
437: printf("\n\nThe number of s-pairs is %d.\n",Spairs);
438: printf("Criterion1 is applied %d times.\n",Criterion1);
439: printf("Criterions M,F and B are applied M=%d, F=%d, B=%d times.\n",Criterion2M,Criterion2F,Criterion2B);
440: Criterion1 = Criterion2M = Criterion2F = Criterion2B = 0;
441: Spairs = 0;
442: }
443: #ifdef STATISTICS
444: printf("\n\nEqual : %d\n",CountE);
445: for (i=0; i<2*N; i++) {
446: printf("%d times of loop : %d\n",i,CountI[i]);
447: }
448: #endif
449:
450: ans = newGradedPolySet(1);
1.5 ! takayama 451: ans->gb = 1;
1.1 maekawa 452: for (i=0; i<ans->lim; i++) {
453: ans->polys[i] = newPolySet(1);
454: }
455:
456: for (i=0; i<g.size; i++) {
457: if (g.del[i] == 0) {
458: grade=-1;whereInG(ans,g.g[i],&grade,&indx,0); /* we do not use sugar. */
459: ans = putPolyInG(ans,g.g[i],grade,indx,(struct syz0 *)NULL,0,-1);
460: }
461: }
462:
463: printf("\n");
464: return(ans);
465: }
466:
467:
468: void outputPairSet(d)
1.3 takayama 469: struct pairSet d;
1.1 maekawa 470: { int i;
1.3 takayama 471: printf("\nOutput struct pairSet. ");
472: printf(".size=%d, .lim=%d :\n",d.size,d.lim);
473: for (i=0; i<d.size; i++) {
474: printf("%d[%d] i=%d,j=%d, lcm=%s ",i,d.del[i],d.p[i].i,d.p[i].j,
475: POLYToString(d.p[i].lcm,' ',0));
476: }
477: printf("\n"); fflush(stdout);
1.1 maekawa 478: }
479:
480: void outputPolySet_gm(g)
1.3 takayama 481: struct polySet_gm g;
1.1 maekawa 482: { int i;
1.3 takayama 483: printf("\nOutput struct polySet_gm. ");
484: printf(".size=%d, .lim=%d :\n",g.size,g.lim);
485: for (i=0; i<g.size; i++) {
486: printf("%d[%d] %s\n",i,g.del[i],POLYToString(g.g[i],' ',0));
487: }
488: printf("\n"); fflush(stdout);
1.1 maekawa 489: }
490:
491: int criterion1(f0,g0,lc0)
1.3 takayama 492: POLY f0;
493: POLY g0;
494: POLY lc0;
1.1 maekawa 495: {
496: /* This is used only for commutative case. */
497: register int i;
498: MONOMIAL f,g,lc;
499: int n;
500: f = f0->m; g = g0->m; lc = lc0->m;
501: if (f->ringp != g->ringp) return(0);
502: if (f->ringp != lc->ringp) return(0);
503: n = f->ringp->n;
504: for (i = 0; i<n ; i++) {
505: if (( f->e[i].x + g->e[i].x) != (lc->e[i].x)) return(0);
506: }
507: for (i = 0; i<n ; i++) {
508: if (( f->e[i].D + g->e[i].D) != (lc->e[i].D)) return(0);
509: }
510: return(1);
511: }
512:
513:
514: POLY mReductionBySet(f,s,size)
1.3 takayama 515: POLY f;
516: POLY *s;
517: int size;
1.1 maekawa 518: {
519: int reduced1;
520: int i;
521: POLY cc,cg;
522: do {
523: reduced1 = 0;
524: for (i=0; i<size; i++) {
525: if (f ISZERO) goto ss;
526: if ((*isReducible)(f,s[i])) {
1.3 takayama 527: f = (*reduction1)(f,s[i],0,&cc,&cg);
528: reduced1 = 1;
1.1 maekawa 529: }
530: }
531: } while (reduced1 != 0);
532: ss: ;
533: return(f);
534: }
535:
536: POLY mReductionCdrBySet(f,s,size)
1.3 takayama 537: POLY f;
538: POLY *s;
539: int size;
1.1 maekawa 540: {
541: int reduced1;
542: int i;
543: POLY cc,cg;
544: POLY fs;
545: do {
546: reduced1 = 0;
547: for (i=0; i<size; i++) {
548: if (f ISZERO) goto ss;
549: if ((fs=(*isCdrReducible)(f,s[i])) != ZERO) {
1.3 takayama 550: f = (*reduction1Cdr)(f,fs,s[i],0,&cc,&cg);
551: reduced1 = 1;
1.1 maekawa 552: }
553: }
554: } while (reduced1 != 0);
555: ss: ;
556: return(f);
557: }
558:
559:
560: void errorGBGM(s)
1.3 takayama 561: char *s;
1.1 maekawa 562: {
563: fprintf(stderr,"Fatal error in yaGbasis.c: %s \n",s);
564: exit(10);
565: }
566:
567:
568:
569:
570:
571:
572:
573:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>