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