[BACK]Return to gbGM.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / kan96xx / Kan

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>