[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.5

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>