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

Annotation of OpenXM/src/kan96xx/Kan/ecart.c, Revision 1.6

1.6     ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/ecart.c,v 1.5 2003/07/19 06:03:57 takayama Exp $ */
1.1       takayama    2: #include <stdio.h>
                      3: #include "datatype.h"
                      4: #include "extern2.h"
                      5: #include "gradedset.h"
                      6:
                      7: #define outofmemory(a) if (a == NULL) {fprintf(stderr,"ecart.c: out of memory.\n"); exit(10);}
                      8: /* Local structures */
                      9: struct ecartReducer {
                     10:   int ell; /* s^\ell */
                     11:   int first; /* first =1 ==> gset else gg */
                     12:   int grade;  int gseti;  /* if first==1,  gradedPolySet */
                     13:   int ggi;                /* if first==0   ecartPoly [] */
                     14: };
                     15: struct ecartPolyArray {
                     16:   int size;
                     17:   int limit;
                     18:   POLY *pa;
1.5       takayama   19:   POLY *cf;
                     20:   POLY *syz;
1.1       takayama   21: };
                     22:
                     23: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,struct ecartPolyArray *epa);
                     24: static int ecartCheckPoly(POLY f);  /* check if it does not contain s-variables */
                     25: static int ecartCheckEnv();         /* check if the environment is OK for ecart div*/
1.5       takayama   26: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf, POLY syz);
1.1       takayama   27: static int ecartGetEll(POLY r,POLY g);
1.5       takayama   28: static POLY ecartDivideSv(POLY r,int *d);
                     29: /* No automatic homogenization and s is used as a standart var. */
                     30: static POLY reduction_ecart0(POLY r,struct gradedPolySet *gset,
                     31:                              int needSyz, struct syz0 *syzp);
                     32: /* Automatic homogenization and s->1 */
                     33: static POLY reduction_ecart1(POLY r,struct gradedPolySet *gset,
                     34:                              int needSyz, struct syz0 *syzp);
                     35: static POLY  ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
                     36:                             struct gradedPolySet *gg,POLY r);
1.1       takayama   37:
                     38: extern int DebugReductionRed;
1.6     ! takayama   39: int DebugReductionEcart = 0;
1.1       takayama   40:
                     41: /* This is used for goHomogenization */
                     42: extern int DegreeShifto_size;
                     43: extern int *DegreeShifto_vec;
                     44:
1.3       takayama   45: /* It is used reduction_ecart() and ecartFindReducer()
                     46:    to determine if we homogenize in this function */
                     47: extern int EcartAutomaticHomogenization;
                     48:
1.1       takayama   49: #define LARGE 0x7fffffff
                     50:
                     51:
1.5       takayama   52: static POLY ecartDivideSv(POLY r,int *d) {
1.1       takayama   53:   POLY f;
                     54:   int k;
1.5       takayama   55:   *d = 0;
1.1       takayama   56:   if (r == POLYNULL) return r;
                     57:   f = r; k = LARGE;
                     58:   while (r != POLYNULL) {
                     59:     if (r->m->e[0].x < k) {
                     60:       k = r->m->e[0].x;
                     61:     }
                     62:     r = r->next;
                     63:   }
                     64:   if (k > 0) {
1.5       takayama   65:     *d = k;
1.1       takayama   66:     f = r;
                     67:     while (r != POLYNULL) {
                     68:       r->m->e[0].x -= k;
                     69:       r = r->next;
                     70:     }
                     71:   }
                     72:   return f;
                     73: }
                     74:
                     75: static int ecartGetEll(POLY f,POLY g) {
                     76:   int n,i,p;
                     77:   MONOMIAL tf;
                     78:   MONOMIAL tg;
                     79:
                     80:   if (f ISZERO) return(-1);
                     81:   if (g ISZERO) return(-1);
                     82:
                     83:   checkRingIsR(f,g);
                     84:
                     85:   if (!isSameComponent_x(f,g)) return(-1);
                     86:   tf = f->m; tg = g->m; n = tf->ringp->n;
                     87:   for (i=1; i<n; i++) {
                     88:     if (tf->e[i].x < tg->e[i].x) return(-1);
                     89:     if (tf->e[i].D < tg->e[i].D) return(-1);
                     90:   }
                     91:   if (tf->e[0].D < tg->e[0].D) return(-1);  /*  h  */
1.2       takayama   92:   p = tf->e[0].x - tg->e[0].x;  /* H,  s */
1.1       takayama   93:   if (p >=0 ) return 0;
                     94:   else return(-p);
                     95: }
                     96:
                     97:
                     98: #define EP_SIZE 10
1.5       takayama   99: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf,POLY syz)
1.1       takayama  100: {
                    101:   struct ecartPolyArray *a;
                    102:   POLY *ep;
                    103:   int i;
                    104:   if (eparray == (struct ecartPolyArray *)NULL) {
                    105:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
                    106:     outofmemory(a);
                    107:     ep = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    108:     outofmemory(ep);
1.5       takayama  109:     a->cf = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    110:     outofmemory(a->cf);
                    111:     a->syz = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    112:     outofmemory(a->syz);
1.1       takayama  113:     a->limit = EP_SIZE;
                    114:     a->size = 0;
                    115:     a->pa = ep;
                    116:     eparray = a;
                    117:   }
                    118:   if (eparray->size >= eparray->limit) {
                    119:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
                    120:     outofmemory(a);
                    121:     ep = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    122:     outofmemory(ep);
1.5       takayama  123:     a->cf = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    124:     outofmemory(a->cf);
                    125:     a->syz = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    126:     outofmemory(a->syz);
1.1       takayama  127:     a->limit = (eparray->limit)*2;
                    128:     a->size = eparray->size;
                    129:     a->pa = ep;
                    130:     for (i=0; i<eparray->size; i++) {
                    131:       (a->pa)[i] = (eparray->pa)[i];
1.5       takayama  132:       (a->cf)[i] = (eparray->cf)[i];
                    133:       (a->syz)[i] = (eparray->syz)[i];
1.1       takayama  134:     }
                    135:     eparray  = a;
                    136:   }
                    137:   i = eparray->size;
                    138:   (eparray->pa)[i] = g;
1.5       takayama  139:   (eparray->cf)[i] = cf;
                    140:   (eparray->syz)[i] = syz;
1.1       takayama  141:   (eparray->size)++;
                    142:   return eparray;
                    143: }
                    144:
                    145: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,
                    146:                                             struct ecartPolyArray *epa)
                    147: {
                    148:   int grd;
                    149:   struct polySet *set;
                    150:   int minGrade = 0;
                    151:   int minGseti = 0;
                    152:   int minGgi   = 0;
                    153:   int ell1 = LARGE;
                    154:   int ell2 = LARGE;
                    155:   int ell;
                    156:   int i;
                    157:   struct ecartReducer er;
                    158:   /* Try to find a reducer in gset; */
                    159:   grd = 0;
                    160:   while (grd < gset->maxGrade) {
                    161:     set = gset->polys[grd];
                    162:     for (i=0; i<set->size; i++) {
1.2       takayama  163:       if (set->gh[i] == POLYNULL) {
                    164:         /* goHomogenize set->gh[i] */
1.4       takayama  165:           if (EcartAutomaticHomogenization) {
                    166:               set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
                    167:           }else{
                    168:               set->gh[i] = set->g[i];
                    169:           }
1.2       takayama  170:       }
                    171:       ell = ecartGetEll(r,set->gh[i]);
1.1       takayama  172:       if ((ell>=0) && (ell < ell1)) {
                    173:         ell1 = ell;
                    174:         minGrade = grd; minGseti=i;
                    175:       }
                    176:     }
                    177:     grd++;
                    178:   }
                    179:   if (epa != NULL) {
                    180:     /* Try to find in the second group. */
                    181:     for (i=0; i< epa->size; i++) {
                    182:       ell = ecartGetEll(r,(epa->pa)[i]);
                    183:       if ((ell>=0) && (ell < ell2)) {
                    184:         ell2 = ell;
                    185:         minGgi = i;
                    186:       }
                    187:     }
                    188:   }
                    189:
1.6     ! takayama  190:   if (DebugReductionRed || (DebugReductionEcart&1)) {
1.1       takayama  191:     printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi);
                    192:   }
                    193:   if (ell1 <= ell2) {
                    194:     if (ell1 == LARGE) {
                    195:       er.ell = -1;
                    196:       return er;
                    197:     }else{
                    198:       er.ell = ell1;
                    199:       er.first = 1;
                    200:       er.grade = minGrade;
                    201:       er.gseti = minGseti;
                    202:       return er;
                    203:     }
                    204:   }else{
                    205:       er.ell = ell2;
                    206:       er.first = 0;
                    207:       er.ggi = minGgi;
                    208:       return er;
                    209:   }
                    210: }
                    211:
1.5       takayama  212: static POLY  ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
                    213:                            struct gradedPolySet *gg,POLY r)
                    214: {
                    215:   POLY f;
                    216:   int grd,i;
                    217:   POLY q;
                    218:   struct coeff *c;
                    219:   f = ppMult(cf,r_0);
                    220:   while (syz != POLYNULL) {
                    221:     grd = syz->m->e[0].x;
                    222:     i = syz->m->e[0].D;
                    223:     c = syz->coeffp;
                    224:     if (c->tag == POLY_COEFF) {
                    225:       q = c->val.f;
                    226:     }else{
                    227:       q = POLYNULL;
                    228:     }
                    229:     f = ppAdd(f,ppMult(q,((gg->polys[grd])->g)[i]));
                    230:                                 /* not gh, works only for _ecart0 */
                    231:     syz = syz->next;
                    232:   }
                    233:   f = ppSub(f,r);
                    234:   return f;
                    235: }
                    236:
                    237:
                    238: POLY reduction_ecart(r,gset,needSyz,syzp)
                    239:      POLY r;
                    240:      struct gradedPolySet *gset;
                    241:      int needSyz;
                    242:      struct syz0 *syzp; /* set */
                    243: {
                    244:   if (EcartAutomaticHomogenization) {
                    245:     return reduction_ecart1(r,gset,needSyz,syzp);
                    246:   }else{
                    247:     return reduction_ecart0(r,gset,needSyz,syzp);
                    248:   }
                    249: }
                    250:
1.1       takayama  251: /*
                    252:   r and gset are assumed to be (0,1)-homogenized (h-homogenized)
1.5       takayama  253:   Polynomials r and gset are assumed
1.3       takayama  254:   to be double homogenized (h-homogenized and s(H)-homogenized)
1.1       takayama  255:  */
1.5       takayama  256: static POLY reduction_ecart0(r,gset,needSyz,syzp)
1.1       takayama  257:      POLY r;
                    258:      struct gradedPolySet *gset;
                    259:      int needSyz;
                    260:      struct syz0 *syzp; /* set */
                    261: {
                    262:   int reduced,reduced1,reduced2;
                    263:   int grd;
                    264:   struct polySet *set;
                    265:   POLY cf,syz;
                    266:   int i;
                    267:   POLY cc,cg;
                    268:   struct ecartReducer ells;
                    269:   struct ecartPolyArray *gg;
                    270:   POLY pp;
                    271:   int ell;
1.5       takayama  272:   POLY cf_o;
                    273:   POLY syz_o;
                    274:   POLY r_0;
1.1       takayama  275:
                    276:   extern struct ring *CurrentRingp;
                    277:   struct ring *rp;
                    278:
1.5       takayama  279:   r_0 = r;
1.1       takayama  280:   gg = NULL;
                    281:   if (needSyz) {
                    282:     if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
                    283:     cf = cxx(1,0,0,rp);
                    284:     syz = ZERO;
                    285:   }
                    286:
                    287:   if (r != POLYNULL) {
1.4       takayama  288:     rp = r->m->ringp;
                    289:     if (! rp->weightedHomogenization) {
                    290:       errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
                    291:     }
1.1       takayama  292:   }
                    293:
1.6     ! takayama  294:   if (DebugReductionEcart&1) printf("--------------------------------------\n");
1.5       takayama  295:   do {
                    296:     if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1));
1.6     ! takayama  297:     if (DebugReductionEcart & 1) printf("r=%s+...\n",POLYToString(head(r),'*',1));
1.5       takayama  298:     ells = ecartFindReducer(r,gset,gg);
                    299:     ell = ells.ell;
                    300:     if (ell > 0) {
1.6     ! takayama  301:       if (DebugReductionEcart & 2) printf("*");
1.5       takayama  302:       if (needSyz) {
                    303:         gg = ecartPutPolyInG(r,gg,cf,syz);
                    304:       }else{
                    305:         gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
                    306:       }
                    307:     }
                    308:     if (ell >= 0) {
                    309:       if (ells.first) {
                    310:         pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
                    311:       }else{
1.6     ! takayama  312:         if (DebugReductionEcart & 4) printf("#");
1.5       takayama  313:         pp = (gg->pa)[ells.ggi];
                    314:       }
                    315:       if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
                    316:       r = (*reduction1)(r,pp,needSyz,&cc,&cg);
                    317:       if (needSyz) {
                    318:         if (ells.first) {
                    319:           if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
                    320:           cf = ppMult(cc,cf);
                    321:           syz = cpMult(toSyzCoeff(cc),syz);
                    322:           syz = ppAdd(syz,toSyzPoly(cg,ells.grade,ells.gseti));
                    323:         }else{
                    324:           if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
                    325:           cf_o = (gg->cf)[ells.ggi];
                    326:           syz_o = (gg->syz)[ells.ggi];
                    327:           cf = ppMult(cc,cf);
                    328:           cf = ppAdd(cf,ppMult(cg,cf_o));
                    329:           syz = cpMult(toSyzCoeff(cc),syz);
                    330:           syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o));
1.6     ! takayama  331:           /* Note. 2003.07.19 */
1.5       takayama  332:         }
                    333:         if (DebugReductionRed) {
1.6     ! takayama  334:           POLY tp;
        !           335:           tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
        !           336:           if (tp != POLYNULL) {
        !           337:             fprintf(stderr,"reduction_ecart0(): sygyzy is broken. Return the Current values.\n");
        !           338:             fprintf(stderr,"%s\n",POLYToString(tp,'*',1));
        !           339:             syzp->cf = cf;
        !           340:             syzp->syz = syz;
        !           341:             return r;
        !           342:           }
1.5       takayama  343:         }
                    344:       }
                    345:       if (r ISZERO) goto ss;
                    346:       /*r = ecartDivideSv(r,&se); r = r/s^? Don't do this. */
                    347:     }
                    348:   }while (ell >= 0);
                    349:
                    350:  ss: ;
                    351:   if (needSyz) {
                    352:     syzp->cf = cf;   /* cf is in the CurrentRingp */
                    353:     syzp->syz = syz; /* syz is in the SyzRingp */
                    354:   }
                    355:
                    356:   return(r);
                    357: }
                    358:
                    359: /*
                    360:   r and gset are assumed to be (0,1)-homogenized (h-homogenized)
                    361:   r and gset are not assumed
                    362:   to be double homogenized (h-homogenized and s(H)-homogenized)
                    363:   They are automatically s(H)-homogenized, and s is set to 1
                    364:   when this function returns.
                    365:  */
                    366: static POLY reduction_ecart1(r,gset,needSyz,syzp)
                    367:      POLY r;
                    368:      struct gradedPolySet *gset;
                    369:      int needSyz;
                    370:      struct syz0 *syzp; /* set */
                    371: {
                    372:   int reduced,reduced1,reduced2;
                    373:   int grd;
                    374:   struct polySet *set;
                    375:   POLY cf,syz;
                    376:   int i;
                    377:   POLY cc,cg;
                    378:   struct ecartReducer ells;
                    379:   struct ecartPolyArray *gg;
                    380:   POLY pp;
                    381:   int ell;
                    382:   int se;
                    383:
                    384:   extern struct ring *CurrentRingp;
                    385:   struct ring *rp;
                    386:
                    387:   gg = NULL;
                    388:   if (needSyz) {
                    389:     if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
                    390:     cf = cxx(1,0,0,rp);
                    391:     syz = ZERO;
                    392:   }
                    393:
                    394:   if (r != POLYNULL) {
                    395:     rp = r->m->ringp;
                    396:     if (! rp->weightedHomogenization) {
                    397:       errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
                    398:     }
1.3       takayama  399:   }
1.5       takayama  400:
                    401:   r = goHomogenize11(r,DegreeShifto_vec,DegreeShifto_size,-1,1);
1.1       takayama  402:   /* 1 means homogenize only s */
                    403:
                    404:   do {
                    405:     if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1));
                    406:     ells = ecartFindReducer(r,gset,gg);
                    407:     ell = ells.ell;
                    408:     if (ell > 0) {
1.5       takayama  409:       gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
1.1       takayama  410:     }
                    411:     if (ell >= 0) {
                    412:       if (ells.first) {
                    413:         pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
                    414:       }else{
                    415:         pp = (gg->pa)[ells.ggi];
                    416:       }
                    417:       if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
                    418:       r = (*reduction1)(r,pp,needSyz,&cc,&cg);
                    419:       if (needSyz) {
                    420:         if (ells.first) {
                    421:           cf = ppMult(cc,cf);
                    422:           syz = cpMult(toSyzCoeff(cc),syz);
                    423:           syz = ppAddv(syz,toSyzPoly(cg,ells.grade,ells.gseti));
                    424:         }else{
1.6     ! takayama  425:           fprintf(stderr,"It has not yet implemented.\n");
        !           426:           exit(10);
1.1       takayama  427:           /* BUG: not yet */
                    428:         }
                    429:       }
1.5       takayama  430:       if (r ISZERO) goto ss1;
                    431:       r = ecartDivideSv(r,&se); /* r = r/s^? */
1.1       takayama  432:     }
                    433:   }while (ell >= 0);
                    434:
1.5       takayama  435:  ss1: ;
1.1       takayama  436:   if (needSyz) {
                    437:     syzp->cf = cf;   /* cf is in the CurrentRingp */
                    438:     syzp->syz = syz; /* syz is in the SyzRingp */
                    439:     /* BUG: dehomogenize the syzygy */
1.6     ! takayama  440:     fprintf(stderr,"It has not yet implemented.\n");
        !           441:     exit(10);
1.1       takayama  442:   }
1.5       takayama  443:
1.1       takayama  444:   r = goDeHomogenizeS(r);
1.5       takayama  445:
1.1       takayama  446:   return(r);
                    447: }

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