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

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

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