Annotation of OpenXM/src/kan96xx/Kan/gbGM.c, Revision 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