Annotation of OpenXM/src/kan96xx/Kan/gbGM.c, Revision 1.4
1.4 ! takayama 1: /* $OpenXM: OpenXM/src/kan96xx/Kan/gbGM.c,v 1.3 2001/05/04 01:06:23 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.4 ! takayama 343: struct gradedPolySet *groebner_gm(f,needBack,needSyz,grP,countDown,forceReduction,reduceOnly)
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.1 maekawa 351: {
352: int r;
353: struct pair_gm top;
354: int tindex;
355: POLY h;
356: int i,j;
357: int ss;
358: struct polySet_gm g;
359: struct pairSet d;
360: int t;
361: struct spValue sv;
362: struct gradedPolySet *ans;
363: int grade,indx;
364: extern int ReduceLowerTerms;
365:
366: if (needBack || needSyz) {
367: fprintf(stderr,"Warning: groebner_gm() does not compute the backward transformation and syzygies.\n");
1.4 ! takayama 368: }
! 369: if (reduceOnly) {
! 370: fprintf(stderr,"Warning: groebner_gm() does not implement reduceOnly.\n");
1.1 maekawa 371: }
372:
373: #ifdef STATISTICS
374: CountE = 0;
375: for (i=0; i<2*N; i++) CountI[i]=0;
376: #endif
377: /* taking the statistics. */
378: Criterion1 = Criterion2B = Criterion2M = Criterion2F = 0;
379: Spairs = 0;
380:
381: r = f->n;
382: g = newPolySet_gm(r + GLIMIT);
383: d = newPairSet(1);
384:
385: g.size = 0;
386: g.g[g.size] = getArrayOfPOLY(f,g.size);
387: g.del[g.size] = 0; g.size ++;
388:
389: for (t=1; t<r; t++) {
390: d = updatePair_gm(d,t,g,getArrayOfPOLY(f,t));
391: g.g[t] = getArrayOfPOLY(f,t); g.del[t] = 0; g.size ++ ;
392: g = markRedundant_gm(g,t);
393: }
394:
395: #ifdef DEBUG
396: outputPolySet_gm(g);
397: #endif
398:
399: while (d.size > 0) {
400: tindex = minPair_gm(d);
401: top = d.p[tindex];
402: if (Verbose) {
403: printf("\nRemaining pairs = %d : ",d.size);
404: }
405: i = top.i;
406: j = top.j;
407: if (Verbose >=2) printf("sp(%d,%d) ",i,j);
408: Spairs++;
409: sv = (*sp)(g.g[i],g.g[j]);
410: h = ppAddv(ppMult(sv.a,g.g[i]),ppMult(sv.b,g.g[j]));
411: h = mReductionBySet(h,g.g,g.size);
412: if (Verbose >=2) printf("->%s ",POLYToString(h,' ',0));
413:
414: if (!Verbose) {
415: printf(" <%d>",d.size);
416: fflush(stdout);
417: }
418:
419: if (!(h ISZERO)) {
420: if (ReduceLowerTerms) {
1.3 takayama 421: h = mReductionCdrBySet(h,g.g,g.size);
1.1 maekawa 422: }
423: d = updatePair_gm(d,r,g,h);
424: g.g[r] = h; g.del[r] = 0;
425: r++; g.size++;
426: g = markRedundant_gm(g,r-1);
427: }
428: d = deletePair_gm(d,tindex);
429: #ifdef DEBUG
430: outputPolySet_gm(g);
431: #endif
432: }
433:
434:
435: if (Verbose || Statistics) {
436: printf("\n\nThe number of s-pairs is %d.\n",Spairs);
437: printf("Criterion1 is applied %d times.\n",Criterion1);
438: printf("Criterions M,F and B are applied M=%d, F=%d, B=%d times.\n",Criterion2M,Criterion2F,Criterion2B);
439: Criterion1 = Criterion2M = Criterion2F = Criterion2B = 0;
440: Spairs = 0;
441: }
442: #ifdef STATISTICS
443: printf("\n\nEqual : %d\n",CountE);
444: for (i=0; i<2*N; i++) {
445: printf("%d times of loop : %d\n",i,CountI[i]);
446: }
447: #endif
448:
449: ans = newGradedPolySet(1);
450: for (i=0; i<ans->lim; i++) {
451: ans->polys[i] = newPolySet(1);
452: }
453:
454: for (i=0; i<g.size; i++) {
455: if (g.del[i] == 0) {
456: grade=-1;whereInG(ans,g.g[i],&grade,&indx,0); /* we do not use sugar. */
457: ans = putPolyInG(ans,g.g[i],grade,indx,(struct syz0 *)NULL,0,-1);
458: }
459: }
460:
461: printf("\n");
462: return(ans);
463: }
464:
465:
466: void outputPairSet(d)
1.3 takayama 467: struct pairSet d;
1.1 maekawa 468: { int i;
1.3 takayama 469: printf("\nOutput struct pairSet. ");
470: printf(".size=%d, .lim=%d :\n",d.size,d.lim);
471: for (i=0; i<d.size; i++) {
472: printf("%d[%d] i=%d,j=%d, lcm=%s ",i,d.del[i],d.p[i].i,d.p[i].j,
473: POLYToString(d.p[i].lcm,' ',0));
474: }
475: printf("\n"); fflush(stdout);
1.1 maekawa 476: }
477:
478: void outputPolySet_gm(g)
1.3 takayama 479: struct polySet_gm g;
1.1 maekawa 480: { int i;
1.3 takayama 481: printf("\nOutput struct polySet_gm. ");
482: printf(".size=%d, .lim=%d :\n",g.size,g.lim);
483: for (i=0; i<g.size; i++) {
484: printf("%d[%d] %s\n",i,g.del[i],POLYToString(g.g[i],' ',0));
485: }
486: printf("\n"); fflush(stdout);
1.1 maekawa 487: }
488:
489: int criterion1(f0,g0,lc0)
1.3 takayama 490: POLY f0;
491: POLY g0;
492: POLY lc0;
1.1 maekawa 493: {
494: /* This is used only for commutative case. */
495: register int i;
496: MONOMIAL f,g,lc;
497: int n;
498: f = f0->m; g = g0->m; lc = lc0->m;
499: if (f->ringp != g->ringp) return(0);
500: if (f->ringp != lc->ringp) return(0);
501: n = f->ringp->n;
502: for (i = 0; i<n ; i++) {
503: if (( f->e[i].x + g->e[i].x) != (lc->e[i].x)) return(0);
504: }
505: for (i = 0; i<n ; i++) {
506: if (( f->e[i].D + g->e[i].D) != (lc->e[i].D)) return(0);
507: }
508: return(1);
509: }
510:
511:
512: POLY mReductionBySet(f,s,size)
1.3 takayama 513: POLY f;
514: POLY *s;
515: int size;
1.1 maekawa 516: {
517: int reduced1;
518: int i;
519: POLY cc,cg;
520: do {
521: reduced1 = 0;
522: for (i=0; i<size; i++) {
523: if (f ISZERO) goto ss;
524: if ((*isReducible)(f,s[i])) {
1.3 takayama 525: f = (*reduction1)(f,s[i],0,&cc,&cg);
526: reduced1 = 1;
1.1 maekawa 527: }
528: }
529: } while (reduced1 != 0);
530: ss: ;
531: return(f);
532: }
533:
534: POLY mReductionCdrBySet(f,s,size)
1.3 takayama 535: POLY f;
536: POLY *s;
537: int size;
1.1 maekawa 538: {
539: int reduced1;
540: int i;
541: POLY cc,cg;
542: POLY fs;
543: do {
544: reduced1 = 0;
545: for (i=0; i<size; i++) {
546: if (f ISZERO) goto ss;
547: if ((fs=(*isCdrReducible)(f,s[i])) != ZERO) {
1.3 takayama 548: f = (*reduction1Cdr)(f,fs,s[i],0,&cc,&cg);
549: reduced1 = 1;
1.1 maekawa 550: }
551: }
552: } while (reduced1 != 0);
553: ss: ;
554: return(f);
555: }
556:
557:
558: void errorGBGM(s)
1.3 takayama 559: char *s;
1.1 maekawa 560: {
561: fprintf(stderr,"Fatal error in yaGbasis.c: %s \n",s);
562: exit(10);
563: }
564:
565:
566:
567:
568:
569:
570:
571:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>