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