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

1.21    ! takayama    1: /* $OpenXM: OpenXM/src/kan96xx/Kan/ecart.c,v 1.20 2003/09/21 02:19:43 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);
1.8       takayama   24: static struct ecartReducer ecartFindReducer_mod(POLY r,struct gradedPolySet *gset,struct ecartPolyArray *epa);
1.1       takayama   25: static int ecartCheckPoly(POLY f);  /* check if it does not contain s-variables */
                     26: static int ecartCheckEnv();         /* check if the environment is OK for ecart div*/
1.5       takayama   27: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf, POLY syz);
1.1       takayama   28: static int ecartGetEll(POLY r,POLY g);
1.21    ! takayama   29: static int ecartGetEllPartial(POLY r,POLY g);
1.5       takayama   30: static POLY ecartDivideSv(POLY r,int *d);
                     31: /* No automatic homogenization and s is used as a standart var. */
                     32: static POLY reduction_ecart0(POLY r,struct gradedPolySet *gset,
                     33:                              int needSyz, struct syz0 *syzp);
                     34: /* Automatic homogenization and s->1 */
                     35: static POLY reduction_ecart1(POLY r,struct gradedPolySet *gset,
                     36:                              int needSyz, struct syz0 *syzp);
1.8       takayama   37: static POLY reduction_ecart1_mod(POLY r,struct gradedPolySet *gset);
1.5       takayama   38: static POLY  ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
                     39:                             struct gradedPolySet *gg,POLY r);
1.16      takayama   40: static void  ecartCheckSyz0_printinfo(POLY cf,POLY r_0,POLY syz,
                     41:                             struct gradedPolySet *gg,POLY r);
1.1       takayama   42:
1.8       takayama   43: extern int DebugReductionRed;
                     44: extern int TraceLift;
                     45: struct ring *TraceLift_ringmod;
1.9       takayama   46: extern DoCancel;
1.6       takayama   47: int DebugReductionEcart = 0;
1.11      takayama   48: extern DebugContentReduction;
1.19      takayama   49: extern int Sugar;
1.1       takayama   50:
                     51: /* This is used for goHomogenization */
                     52: extern int DegreeShifto_size;
                     53: extern int *DegreeShifto_vec;
1.19      takayama   54: int Ecart_sugarGrade;
1.1       takayama   55:
1.3       takayama   56: /* It is used reduction_ecart() and ecartFindReducer()
                     57:    to determine if we homogenize in this function */
                     58: extern int EcartAutomaticHomogenization;
                     59:
1.1       takayama   60: #define LARGE 0x7fffffff
                     61:
                     62:
1.5       takayama   63: static POLY ecartDivideSv(POLY r,int *d) {
1.1       takayama   64:   POLY f;
                     65:   int k;
1.5       takayama   66:   *d = 0;
1.1       takayama   67:   if (r == POLYNULL) return r;
                     68:   f = r; k = LARGE;
                     69:   while (r != POLYNULL) {
                     70:     if (r->m->e[0].x < k) {
                     71:       k = r->m->e[0].x;
                     72:     }
                     73:     r = r->next;
                     74:   }
1.16      takayama   75:
                     76:   if (k > 0) {
                     77:     *d = k;
                     78:        return ppMult(cxx(1,0,-k,f->m->ringp),f);
                     79:   }else{
                     80:        return f;
                     81:   }
                     82:
                     83:   /* Do not do the below. It caused a bug. cf. misc-2003/07/ecart/b4.sm1 test2.
                     84:      Note. 2003.8.26
                     85:    */
                     86:   /*
1.1       takayama   87:   if (k > 0) {
1.5       takayama   88:     *d = k;
1.1       takayama   89:     f = r;
                     90:     while (r != POLYNULL) {
                     91:       r->m->e[0].x -= k;
                     92:       r = r->next;
                     93:     }
                     94:   }
                     95:   return f;
1.16      takayama   96:   */
1.1       takayama   97: }
                     98:
                     99: static int ecartGetEll(POLY f,POLY g) {
                    100:   int n,i,p;
                    101:   MONOMIAL tf;
                    102:   MONOMIAL tg;
                    103:
                    104:   if (f ISZERO) return(-1);
                    105:   if (g ISZERO) return(-1);
                    106:
                    107:   checkRingIsR(f,g);
                    108:
                    109:   if (!isSameComponent_x(f,g)) return(-1);
                    110:   tf = f->m; tg = g->m; n = tf->ringp->n;
                    111:   for (i=1; i<n; i++) {
                    112:     if (tf->e[i].x < tg->e[i].x) return(-1);
                    113:     if (tf->e[i].D < tg->e[i].D) return(-1);
                    114:   }
                    115:   if (tf->e[0].D < tg->e[0].D) return(-1);  /*  h  */
1.2       takayama  116:   p = tf->e[0].x - tg->e[0].x;  /* H,  s */
1.1       takayama  117:   if (p >=0 ) return 0;
                    118:   else return(-p);
                    119: }
                    120:
1.21    ! takayama  121: static int ecartGetEllPartial(POLY f,POLY g) {
        !           122:   int n,i,p;
        !           123:   MONOMIAL tf;
        !           124:   MONOMIAL tg;
        !           125:   int nglob; int *glob; int k;
        !           126:
        !           127:   if (f ISZERO) return(-1);
        !           128:   if (g ISZERO) return(-1);
        !           129:
        !           130:   checkRingIsR(f,g);
        !           131:
        !           132:   if (!isSameComponent_x(f,g)) return(-1);
        !           133:   tf = f->m; tg = g->m; n = tf->ringp->n;
        !           134:   for (i=1; i<n; i++) {
        !           135:     if (tf->e[i].x < tg->e[i].x) return(-1);
        !           136:     if (tf->e[i].D < tg->e[i].D) return(-1);
        !           137:   }
        !           138:   if (tf->e[0].D < tg->e[0].D) return(-1);  /*  h  */
        !           139:
        !           140:   /* Find u and ell s.t. s^\ell f = u g.
        !           141:      When u contains variables in glob (index x var), then returns -1.
        !           142:      x' = glob, x'' = complement of glob.  Division in Q(x'')[x'].
        !           143:   */
        !           144:   nglob = tf->ringp->partialEcart;
        !           145:   glob = tf->ringp->partialEcartGlobalVarX;
        !           146:   for (i=0; i<nglob; i++) {
        !           147:     k = glob[i];
        !           148:     if (tf->e[k].x > tg->e[k].x) return(-1);
        !           149:   }
        !           150:
        !           151:   p = tf->e[0].x - tg->e[0].x;  /* H,  s */
        !           152:   if (p >=0 ) return 0;
        !           153:   else return(-p);
        !           154: }
        !           155:
1.1       takayama  156:
                    157: #define EP_SIZE 10
1.5       takayama  158: static struct ecartPolyArray *ecartPutPolyInG(POLY g,struct ecartPolyArray *eparray,POLY cf,POLY syz)
1.1       takayama  159: {
                    160:   struct ecartPolyArray *a;
                    161:   POLY *ep;
                    162:   int i;
                    163:   if (eparray == (struct ecartPolyArray *)NULL) {
                    164:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
                    165:     outofmemory(a);
                    166:     ep = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    167:     outofmemory(ep);
1.5       takayama  168:     a->cf = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    169:     outofmemory(a->cf);
                    170:     a->syz = (POLY *) sGC_malloc(sizeof(POLY)*EP_SIZE);
                    171:     outofmemory(a->syz);
1.1       takayama  172:     a->limit = EP_SIZE;
                    173:     a->size = 0;
                    174:     a->pa = ep;
                    175:     eparray = a;
                    176:   }
                    177:   if (eparray->size >= eparray->limit) {
                    178:     a = (struct ecartPolyArray *) sGC_malloc(sizeof(struct ecartPolyArray));
                    179:     outofmemory(a);
                    180:     ep = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    181:     outofmemory(ep);
1.5       takayama  182:     a->cf = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    183:     outofmemory(a->cf);
                    184:     a->syz = (POLY *) sGC_malloc(sizeof(POLY)*((eparray->limit)*2));
                    185:     outofmemory(a->syz);
1.1       takayama  186:     a->limit = (eparray->limit)*2;
                    187:     a->size = eparray->size;
                    188:     a->pa = ep;
                    189:     for (i=0; i<eparray->size; i++) {
                    190:       (a->pa)[i] = (eparray->pa)[i];
1.5       takayama  191:       (a->cf)[i] = (eparray->cf)[i];
                    192:       (a->syz)[i] = (eparray->syz)[i];
1.1       takayama  193:     }
                    194:     eparray  = a;
                    195:   }
                    196:   i = eparray->size;
                    197:   (eparray->pa)[i] = g;
1.5       takayama  198:   (eparray->cf)[i] = cf;
                    199:   (eparray->syz)[i] = syz;
1.1       takayama  200:   (eparray->size)++;
                    201:   return eparray;
                    202: }
                    203:
1.20      takayama  204: #ifndef EXPERIMENT
1.1       takayama  205: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,
                    206:                                             struct ecartPolyArray *epa)
                    207: {
                    208:   int grd;
                    209:   struct polySet *set;
                    210:   int minGrade = 0;
                    211:   int minGseti = 0;
                    212:   int minGgi   = 0;
                    213:   int ell1 = LARGE;
                    214:   int ell2 = LARGE;
                    215:   int ell;
                    216:   int i;
                    217:   struct ecartReducer er;
                    218:   /* Try to find a reducer in gset; */
                    219:   grd = 0;
                    220:   while (grd < gset->maxGrade) {
                    221:     set = gset->polys[grd];
                    222:     for (i=0; i<set->size; i++) {
1.2       takayama  223:       if (set->gh[i] == POLYNULL) {
                    224:         /* goHomogenize set->gh[i] */
1.4       takayama  225:           if (EcartAutomaticHomogenization) {
                    226:               set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
                    227:           }else{
                    228:               set->gh[i] = set->g[i];
                    229:           }
1.2       takayama  230:       }
                    231:       ell = ecartGetEll(r,set->gh[i]);
1.1       takayama  232:       if ((ell>=0) && (ell < ell1)) {
                    233:         ell1 = ell;
                    234:         minGrade = grd; minGseti=i;
                    235:       }
                    236:     }
                    237:     grd++;
                    238:   }
                    239:   if (epa != NULL) {
                    240:     /* Try to find in the second group. */
                    241:     for (i=0; i< epa->size; i++) {
1.21    ! takayama  242:       ell = ecartGetEllPartial(r,(epa->pa)[i]);
1.1       takayama  243:       if ((ell>=0) && (ell < ell2)) {
                    244:         ell2 = ell;
                    245:         minGgi = i;
                    246:       }
                    247:     }
                    248:   }
                    249:
1.18      takayama  250:   if ((DebugReductionRed&1) || (DebugReductionEcart&1)) {
1.1       takayama  251:     printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi);
                    252:   }
1.19      takayama  253:
1.1       takayama  254:   if (ell1 <= ell2) {
                    255:     if (ell1 == LARGE) {
                    256:       er.ell = -1;
                    257:       return er;
                    258:     }else{
                    259:       er.ell = ell1;
                    260:       er.first = 1;
                    261:       er.grade = minGrade;
                    262:       er.gseti = minGseti;
                    263:       return er;
                    264:     }
                    265:   }else{
                    266:       er.ell = ell2;
                    267:       er.first = 0;
                    268:       er.ggi = minGgi;
                    269:       return er;
                    270:   }
                    271: }
1.20      takayama  272: #endif
1.1       takayama  273:
1.5       takayama  274: static POLY  ecartCheckSyz0(POLY cf,POLY r_0,POLY syz,
                    275:                            struct gradedPolySet *gg,POLY r)
                    276: {
                    277:   POLY f;
                    278:   int grd,i;
                    279:   POLY q;
                    280:   struct coeff *c;
                    281:   f = ppMult(cf,r_0);
                    282:   while (syz != POLYNULL) {
                    283:     grd = syz->m->e[0].x;
                    284:     i = syz->m->e[0].D;
                    285:     c = syz->coeffp;
                    286:     if (c->tag == POLY_COEFF) {
                    287:       q = c->val.f;
                    288:     }else{
                    289:       q = POLYNULL;
                    290:     }
                    291:     f = ppAdd(f,ppMult(q,((gg->polys[grd])->g)[i]));
                    292:                                 /* not gh, works only for _ecart0 */
                    293:     syz = syz->next;
                    294:   }
                    295:   f = ppSub(f,r);
                    296:   return f;
                    297: }
                    298:
1.16      takayama  299: static void  ecartCheckSyz0_printinfo(POLY cf,POLY r_0,POLY syz,
                    300:                                       struct gradedPolySet *gg,POLY r)
                    301: {
                    302:   POLY f;
                    303:   int grd,i;
                    304:   POLY q;
                    305:   struct coeff *c;
                    306:   fprintf(stderr,"cf=%s\n",POLYToString(cf,'*',1));
                    307:   fprintf(stderr,"r_0=%s\n",POLYToString(r_0,'*',1));
                    308:   fprintf(stderr,"syz=%s\n",POLYToString(syz,'*',1));
                    309:   fprintf(stderr,"r=%s\n",POLYToString(r,'*',1));
                    310:   f = ppMult(cf,r_0);
                    311:   while (syz != POLYNULL) {
                    312:     grd = syz->m->e[0].x;
                    313:     i = syz->m->e[0].D;
                    314:     c = syz->coeffp;
                    315:     if (c->tag == POLY_COEFF) {
                    316:       q = c->val.f;
                    317:     }else{
                    318:       q = POLYNULL;
                    319:     }
                    320:        fprintf(stderr,"[grd,idx]=[%d,%d], %s\n",grd,i,
                    321:                        POLYToString(((gg->polys[grd])->g)[i],'*',1));
                    322:     /* f = ppAdd(f,ppMult(q,((gg->polys[grd])->g)[i])); */
                    323:     syz = syz->next;
                    324:   }
                    325:   /* f = ppSub(f,r); */
                    326: }
                    327:
1.5       takayama  328:
                    329: POLY reduction_ecart(r,gset,needSyz,syzp)
                    330:      POLY r;
                    331:      struct gradedPolySet *gset;
                    332:      int needSyz;
                    333:      struct syz0 *syzp; /* set */
                    334: {
1.8       takayama  335:   POLY rn;
1.12      takayama  336:   if (TraceLift && needSyz) {
                    337:     warningGradedSet("TraceLift cannot be used to get syzygy. TraceLift is turned off.\n");
                    338:     TraceLift = 0;
                    339:   }
1.8       takayama  340:   if (TraceLift) {
                    341:        if (EcartAutomaticHomogenization) {
                    342:          if (TraceLift_ringmod == NULL) {
                    343:                warningPoly("reduction_ecart: TraceLift_ringmod is not set.\n");
                    344:                return reduction_ecart1(r,gset,needSyz,syzp);
                    345:          }
1.12      takayama  346:          rn = reduction_ecart1_mod(r,gset);
1.8       takayama  347:          if (rn == POLYNULL) return rn;
                    348:          else return reduction_ecart1(r,gset,needSyz,syzp);
                    349:        }else{
                    350:          return reduction_ecart0(r,gset,needSyz,syzp);
                    351:        }
1.5       takayama  352:   }else{
1.8       takayama  353:        if (EcartAutomaticHomogenization) {
                    354:          return reduction_ecart1(r,gset,needSyz,syzp);
                    355:        }else{
                    356:          return reduction_ecart0(r,gset,needSyz,syzp);
                    357:        }
1.5       takayama  358:   }
                    359: }
                    360:
1.1       takayama  361: /*
                    362:   r and gset are assumed to be (0,1)-homogenized (h-homogenized)
1.5       takayama  363:   Polynomials r and gset are assumed
1.3       takayama  364:   to be double homogenized (h-homogenized and s(H)-homogenized)
1.1       takayama  365:  */
1.5       takayama  366: static POLY reduction_ecart0(r,gset,needSyz,syzp)
1.1       takayama  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;
1.5       takayama  382:   POLY cf_o;
                    383:   POLY syz_o;
                    384:   POLY r_0;
1.15      takayama  385:   int se;
1.14      takayama  386:   struct coeff *cont;
1.1       takayama  387:
                    388:   extern struct ring *CurrentRingp;
                    389:   struct ring *rp;
                    390:
1.5       takayama  391:   r_0 = r;
1.1       takayama  392:   gg = NULL;
                    393:   if (needSyz) {
                    394:     if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
                    395:     cf = cxx(1,0,0,rp);
                    396:     syz = ZERO;
                    397:   }
                    398:
                    399:   if (r != POLYNULL) {
1.4       takayama  400:     rp = r->m->ringp;
                    401:     if (! rp->weightedHomogenization) {
                    402:       errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
                    403:     }
1.1       takayama  404:   }
                    405:
1.14      takayama  406:   if (DoCancel && (r != POLYNULL)) shouldReduceContent(r,1);
                    407:
1.20      takayama  408:   if (DebugReductionEcart&1) printf("reduction_ecart0--------------------------------------\n");
1.5       takayama  409:   do {
1.18      takayama  410:     if (DebugReductionRed&1) printf("r=%s\n",POLYToString(r,'*',1));
1.6       takayama  411:     if (DebugReductionEcart & 1) printf("r=%s+...\n",POLYToString(head(r),'*',1));
1.5       takayama  412:     ells = ecartFindReducer(r,gset,gg);
                    413:     ell = ells.ell;
                    414:     if (ell > 0) {
1.6       takayama  415:       if (DebugReductionEcart & 2) printf("*");
1.5       takayama  416:       if (needSyz) {
                    417:         gg = ecartPutPolyInG(r,gg,cf,syz);
                    418:       }else{
                    419:         gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
                    420:       }
                    421:     }
                    422:     if (ell >= 0) {
                    423:       if (ells.first) {
                    424:         pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
                    425:       }else{
1.6       takayama  426:         if (DebugReductionEcart & 4) printf("#");
1.5       takayama  427:         pp = (gg->pa)[ells.ggi];
                    428:       }
1.16      takayama  429:       if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
1.5       takayama  430:       r = (*reduction1)(r,pp,needSyz,&cc,&cg);
1.14      takayama  431:
                    432:       if (DoCancel && (r != POLYNULL)) {
                    433:         if (shouldReduceContent(r,0)) {
                    434:           r = reduceContentOfPoly(r,&cont);
                    435:           shouldReduceContent(r,1);
                    436:           if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("CONT=%s ",coeffToString(cont));
                    437:         }
                    438:       }
                    439:
                    440:
1.5       takayama  441:       if (needSyz) {
                    442:         if (ells.first) {
                    443:           if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
                    444:           cf = ppMult(cc,cf);
                    445:           syz = cpMult(toSyzCoeff(cc),syz);
                    446:           syz = ppAdd(syz,toSyzPoly(cg,ells.grade,ells.gseti));
                    447:         }else{
                    448:           if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
                    449:           cf_o = (gg->cf)[ells.ggi];
                    450:           syz_o = (gg->syz)[ells.ggi];
                    451:           cf = ppMult(cc,cf);
                    452:           cf = ppAdd(cf,ppMult(cg,cf_o));
                    453:           syz = cpMult(toSyzCoeff(cc),syz);
                    454:           syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o));
1.6       takayama  455:           /* Note. 2003.07.19 */
1.5       takayama  456:         }
1.18      takayama  457:         if (DebugReductionRed && needSyz) {
1.6       takayama  458:           POLY tp;
                    459:           tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
                    460:           if (tp != POLYNULL) {
                    461:             fprintf(stderr,"reduction_ecart0(): sygyzy is broken. Return the Current values.\n");
1.16      takayama  462:             fprintf(stderr,"tp=%s\n",POLYToString(tp,'*',1));
                    463:                        ecartCheckSyz0_printinfo(cf,r_0,syz,gset,r);
1.6       takayama  464:             syzp->cf = cf;
                    465:             syzp->syz = syz;
                    466:             return r;
                    467:           }
1.5       takayama  468:         }
                    469:       }
                    470:       if (r ISZERO) goto ss;
1.15      takayama  471:
                    472:       /* r = r/s^? Don't do this?? */
                    473:       r = ecartDivideSv(r,&se);
                    474:          if (needSyz && (se > 0)) {
                    475:                POLY tt;
                    476:                tt = cxx(1,0,-se,rp);
1.16      takayama  477:                cf = ppMult(tt,cf);
1.15      takayama  478:                syz = cpMult(toSyzCoeff(tt),syz);
                    479:          }
                    480:
1.5       takayama  481:     }
                    482:   }while (ell >= 0);
                    483:
                    484:  ss: ;
                    485:   if (needSyz) {
1.16      takayama  486:     syzp->cf = cf;   /* cf is in the ring of r */
                    487:     syzp->syz = syz; /* syz is in the syzRing of r */
1.14      takayama  488:   }
                    489:
                    490:   if (DoCancel && (r != POLYNULL)) {
                    491:     if (r->m->ringp->p == 0) {
                    492:          r = reduceContentOfPoly(r,&cont);
                    493:          if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("cont=%s ",coeffToString(cont));
                    494:     }
1.5       takayama  495:   }
                    496:
                    497:   return(r);
                    498: }
                    499:
                    500: /*
                    501:   r and gset are assumed to be (0,1)-homogenized (h-homogenized)
                    502:   r and gset are not assumed
                    503:   to be double homogenized (h-homogenized and s(H)-homogenized)
                    504:   They are automatically s(H)-homogenized, and s is set to 1
                    505:   when this function returns.
                    506:  */
                    507: static POLY reduction_ecart1(r,gset,needSyz,syzp)
                    508:      POLY r;
                    509:      struct gradedPolySet *gset;
                    510:      int needSyz;
                    511:      struct syz0 *syzp; /* set */
                    512: {
                    513:   int reduced,reduced1,reduced2;
                    514:   int grd;
                    515:   struct polySet *set;
                    516:   POLY cf,syz;
                    517:   int i;
                    518:   POLY cc,cg;
                    519:   struct ecartReducer ells;
                    520:   struct ecartPolyArray *gg;
                    521:   POLY pp;
                    522:   int ell;
1.12      takayama  523:   POLY cf_o;
                    524:   POLY syz_o;
                    525:   POLY r_0;
1.18      takayama  526:   POLY r0;
1.5       takayama  527:   int se;
1.9       takayama  528:   struct coeff *cont;
1.5       takayama  529:
                    530:   extern struct ring *CurrentRingp;
                    531:   struct ring *rp;
1.9       takayama  532:   extern struct ring *SmallRingp;
1.5       takayama  533:
1.12      takayama  534:   r_0 = r;
1.5       takayama  535:   gg = NULL;
                    536:   if (needSyz) {
                    537:     if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; }
                    538:     cf = cxx(1,0,0,rp);
                    539:     syz = ZERO;
                    540:   }
                    541:
                    542:   if (r != POLYNULL) {
                    543:     rp = r->m->ringp;
                    544:     if (! rp->weightedHomogenization) {
                    545:       errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
                    546:     }
1.3       takayama  547:   }
1.5       takayama  548:
                    549:   r = goHomogenize11(r,DegreeShifto_vec,DegreeShifto_size,-1,1);
1.1       takayama  550:   /* 1 means homogenize only s */
1.10      takayama  551:   if (DoCancel && (r != POLYNULL)) shouldReduceContent(r,1);
1.1       takayama  552:
1.20      takayama  553:   if (DebugReductionEcart&1) printf("reduction_ecart1=======================================\n");
1.1       takayama  554:   do {
1.18      takayama  555:     if (DebugReductionRed & 1) printf("(ecart1(d)) r=%s\n",POLYToString(r,'*',1));
1.7       takayama  556:     if (DebugReductionEcart & 1) printf("r=%s+,,,\n",POLYToString(head(r),'*',1));
                    557:
1.1       takayama  558:     ells = ecartFindReducer(r,gset,gg);
                    559:     ell = ells.ell;
                    560:     if (ell > 0) {
1.7       takayama  561:       if (DebugReductionEcart & 2) printf("%");
1.16      takayama  562:       if (needSyz) {
                    563:         gg = ecartPutPolyInG(r,gg,cf,syz);
                    564:       }else{
                    565:         gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
                    566:       }
1.1       takayama  567:     }
                    568:     if (ell >= 0) {
                    569:       if (ells.first) {
                    570:         pp = ((gset->polys[ells.grade])->gh)[ells.gseti];
                    571:       }else{
1.9       takayama  572:         if (DebugReductionEcart & 4) {printf("+"); fflush(NULL);}
1.1       takayama  573:         pp = (gg->pa)[ells.ggi];
                    574:       }
1.16      takayama  575:       if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
1.18      takayama  576:          r0 = r;
1.1       takayama  577:       r = (*reduction1)(r,pp,needSyz,&cc,&cg);
1.18      takayama  578:       if (DebugReductionEcart & 8) {
                    579:         if (ell > 0) {printf("ell+ "); fflush(NULL);
                    580:         }else {
                    581:                  printf("ell0 "); fflush(NULL);
                    582:           if ((*mmLarger)(r,r0) >= 1) {
                    583:             printf("error in reduction.");
                    584:             printf("   r0=%s\n",POLYToString(r0,'*',1));
                    585:             printf("==>r=%s\n",POLYToString(r,'*',1));
                    586:             getchar(); getchar();
                    587:           }
                    588:         }
                    589:       }
1.10      takayama  590:
1.12      takayama  591:       if (DoCancel && (r != POLYNULL)) {
1.10      takayama  592:         if (shouldReduceContent(r,0)) {
                    593:           r = reduceContentOfPoly(r,&cont);
                    594:           shouldReduceContent(r,1);
1.11      takayama  595:           if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("CONT=%s ",coeffToString(cont));
1.10      takayama  596:         }
                    597:       }
                    598:
1.1       takayama  599:       if (needSyz) {
                    600:         if (ells.first) {
1.16      takayama  601:           if (ell > 0) cc = ppMult(cc,cxx(1,0,ell,rp));
1.1       takayama  602:           cf = ppMult(cc,cf);
                    603:           syz = cpMult(toSyzCoeff(cc),syz);
1.16      takayama  604:           syz = ppAdd(syz,toSyzPoly(cg,ells.grade,ells.gseti));
1.1       takayama  605:         }else{
1.12      takayama  606:           if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp));
                    607:           cf_o = (gg->cf)[ells.ggi];
                    608:           syz_o = (gg->syz)[ells.ggi];
                    609:           cf = ppMult(cc,cf);
                    610:           cf = ppAdd(cf,ppMult(cg,cf_o));
                    611:           syz = cpMult(toSyzCoeff(cc),syz);
                    612:           syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o));
                    613:           /* Note. 2003.07.19 */
1.1       takayama  614:         }
                    615:       }
1.12      takayama  616:
1.18      takayama  617:       if (DebugReductionRed && needSyz) {
1.12      takayama  618:         POLY tp;
                    619:         tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
1.16      takayama  620:         tp = goDeHomogenizeS(tp);
1.12      takayama  621:         if (tp != POLYNULL) {
1.16      takayama  622:           fprintf(stderr,"Error: reduction_ecart1(): sygyzy is broken. Return the Current values.\n");
1.12      takayama  623:           fprintf(stderr,"%s\n",POLYToString(tp,'*',1));
                    624:           syzp->cf = cf;
                    625:           syzp->syz = syz;
                    626:           return r;
                    627:         }
                    628:       }
                    629:
1.5       takayama  630:       if (r ISZERO) goto ss1;
                    631:       r = ecartDivideSv(r,&se); /* r = r/s^? */
1.15      takayama  632:
1.16      takayama  633:       if (needSyz && (se > 0)) { /* It may not necessary because of dehomo*/
                    634:         POLY tt;
                    635:         /*printf("!1/H!"); fflush(NULL);*/ /* misc-2003/ecart/t1.sm1 foo4 */
                    636:         tt = cxx(1,0,-se,rp);
                    637:         cf = ppMult(tt,cf);
                    638:         syz = cpMult(toSyzCoeff(tt),syz);
                    639:       }
                    640:
                    641:       /* For debug */
                    642:       if (DebugReductionRed && needSyz) {
                    643:         POLY tp;
                    644:         tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
                    645:         tp = goDeHomogenizeS(tp);
                    646:         if (tp != POLYNULL) {
                    647:           fprintf(stderr,"Error: reduction_ecart1() after divide: sygyzy is broken. Return the Current values.\n");
                    648:           fprintf(stderr,"%s\n",POLYToString(tp,'*',1));
                    649:           syzp->cf = cf;
                    650:           syzp->syz = syz;
                    651:           return r;
                    652:         }
                    653:       }
1.15      takayama  654:
1.1       takayama  655:     }
                    656:   }while (ell >= 0);
                    657:
1.5       takayama  658:  ss1: ;
1.1       takayama  659:   if (needSyz) {
1.12      takayama  660:     /* dehomogenize the syzygy. BUG, this may be inefficient.  */
1.16      takayama  661:     cf = goDeHomogenizeS(cf);
                    662:     syz = goDeHomogenizeS(syz);
                    663:     /*printf("cf=%s\n",POLYToString(cf,'*',1));
                    664:       printf("syz=%s\n",POLYToString(syz,'*',1));*/
1.1       takayama  665:     syzp->cf = cf;   /* cf is in the CurrentRingp */
                    666:     syzp->syz = syz; /* syz is in the SyzRingp */
                    667:   }
1.8       takayama  668:
                    669:   r = goDeHomogenizeS(r);
1.12      takayama  670:   if (DoCancel && (r != POLYNULL)) {
1.10      takayama  671:     if (r->m->ringp->p == 0) {
1.16      takayama  672:       r = reduceContentOfPoly(r,&cont);
                    673:       if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("cont=%s ",coeffToString(cont));
                    674:     }
                    675:   }
                    676:
                    677:   /* For debug */
                    678:   if (DebugReductionRed && needSyz) {
                    679:     POLY tp;
                    680:     tp = ecartCheckSyz0(cf,r_0,syz,gset,r);
                    681:     tp = goDeHomogenizeS(tp);
                    682:     if (tp != POLYNULL) {
                    683:       fprintf(stderr,"Error: reduction_ecart1() last step: sygyzy is broken. Return the Current values.\n");
                    684:       fprintf(stderr,"%s\n",POLYToString(tp,'*',1));
                    685:       syzp->cf = cf;
                    686:       syzp->syz = syz;
                    687:       return r;
1.10      takayama  688:     }
1.9       takayama  689:   }
1.8       takayama  690:
                    691:   return(r);
                    692: }
                    693:
                    694: /* Functions for trace lift  */
                    695: static struct ecartReducer ecartFindReducer_mod(POLY r,
                    696:                                                 struct gradedPolySet *gset,
                    697:                                                 struct ecartPolyArray *epa)
                    698: {
                    699:   int grd;
                    700:   struct polySet *set;
                    701:   int minGrade = 0;
                    702:   int minGseti = 0;
                    703:   int minGgi   = 0;
                    704:   int ell1 = LARGE;
                    705:   int ell2 = LARGE;
                    706:   int ell;
                    707:   int i;
                    708:   struct ecartReducer er;
                    709:   /* Try to find a reducer in gset; */
                    710:   grd = 0;
                    711:   while (grd < gset->maxGrade) {
                    712:     set = gset->polys[grd];
                    713:     for (i=0; i<set->size; i++) {
                    714:       if (set->gh[i] == POLYNULL) {
                    715:         /* goHomogenize set->gh[i] */
                    716:           if (EcartAutomaticHomogenization) {
                    717:               set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
                    718:           }else{
                    719:               set->gh[i] = set->g[i];
                    720:           }
                    721:       }
                    722:          if (TraceLift && (set->gmod[i] == POLYNULL)) {
                    723:                set->gmod[i] = modulop(set->gh[i],TraceLift_ringmod);
                    724:          }
                    725:          if (TraceLift) {
                    726:                ell = ecartGetEll(r,set->gmod[i]);
                    727:          }else{
                    728:                ell = ecartGetEll(r,set->gh[i]);
                    729:          }
                    730:       if ((ell>=0) && (ell < ell1)) {
                    731:         ell1 = ell;
                    732:         minGrade = grd; minGseti=i;
                    733:       }
                    734:     }
                    735:     grd++;
                    736:   }
                    737:   if (epa != NULL) {
                    738:     /* Try to find in the second group. */
                    739:     for (i=0; i< epa->size; i++) {
1.21    ! takayama  740:       ell = ecartGetEllPartial(r,(epa->pa)[i]);
1.8       takayama  741:       if ((ell>=0) && (ell < ell2)) {
                    742:         ell2 = ell;
                    743:         minGgi = i;
                    744:       }
                    745:     }
                    746:   }
                    747:
1.18      takayama  748:   if ((DebugReductionRed & 1)|| (DebugReductionEcart&1)) {
1.8       takayama  749:     printf("ecartFindReducer_mod(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d, p=%d\n",ell1,ell2,minGrade,minGseti,minGgi,TraceLift_ringmod->p);
                    750:   }
                    751:   if (ell1 <= ell2) {
                    752:     if (ell1 == LARGE) {
                    753:       er.ell = -1;
                    754:       return er;
                    755:     }else{
                    756:       er.ell = ell1;
                    757:       er.first = 1;
                    758:       er.grade = minGrade;
                    759:       er.gseti = minGseti;
                    760:       return er;
                    761:     }
                    762:   }else{
                    763:       er.ell = ell2;
                    764:       er.first = 0;
                    765:       er.ggi = minGgi;
                    766:       return er;
                    767:   }
                    768: }
                    769:
                    770: static POLY reduction_ecart1_mod(r,gset)
                    771:      POLY r;
                    772:      struct gradedPolySet *gset;
                    773: {
                    774:   int reduced,reduced1,reduced2;
                    775:   int grd;
                    776:   struct polySet *set;
                    777:   int i;
                    778:   POLY cc,cg;
                    779:   struct ecartReducer ells;
                    780:   struct ecartPolyArray *gg;
                    781:   POLY pp;
1.18      takayama  782:   POLY r0;
1.8       takayama  783:   int ell;
                    784:   int se;
                    785:
                    786:   extern struct ring *CurrentRingp;
                    787:   struct ring *rp;
                    788:
                    789:   gg = NULL;
                    790:
                    791:   if (r != POLYNULL) {
                    792:     rp = r->m->ringp;
                    793:     if (! rp->weightedHomogenization) {
                    794:       errorKan1("%s\n","ecart.c: the given ring must be declared with [(weightedHomogenization) 1]");
                    795:     }
                    796:   }
                    797:
                    798:   r = goHomogenize11(r,DegreeShifto_vec,DegreeShifto_size,-1,1);
                    799:   /* 1 means homogenize only s */
                    800:   /*printf("r=%s (mod 0)\n",POLYToString(head(r),'*',1));
                    801:        KshowRing(TraceLift_ringmod); **/
                    802:
                    803:   r = modulop(r,TraceLift_ringmod);
1.17      takayama  804:   if (r != POLYNULL) rp = r->m->ringp; /* reset rp */
1.8       takayama  805:
                    806:   /* printf("r=%s (mod p)\n",POLYToString(head(r),'*',1)); **/
                    807:
                    808:   if (DebugReductionEcart&1) printf("=====================================mod\n");
                    809:   do {
1.18      takayama  810:     if (DebugReductionRed & 1) printf("(ecart1_mod(d)) r=%s\n",POLYToString(r,'*',1));
1.8       takayama  811:     if (DebugReductionEcart & 1) printf("r=%s+,,,\n",POLYToString(head(r),'*',1));
                    812:
                    813:     ells = ecartFindReducer_mod(r,gset,gg);
                    814:     ell = ells.ell;
                    815:     if (ell > 0) {
                    816:       if (DebugReductionEcart & 2) printf("%");
                    817:       gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL);
                    818:     }
                    819:     if (ell >= 0) {
                    820:       if (ells.first) {
                    821:         pp = ((gset->polys[ells.grade])->gmod)[ells.gseti];
                    822:       }else{
1.9       takayama  823:         if (DebugReductionEcart & 4) {printf("+"); fflush(NULL);}
1.8       takayama  824:         pp = (gg->pa)[ells.ggi];
                    825:       }
1.16      takayama  826:       if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */
1.18      takayama  827:
                    828:          r0 = r;
1.8       takayama  829:       r = (*reduction1)(r,pp,0,&cc,&cg);
1.18      takayama  830:
                    831:       if (DebugReductionEcart & 8) {
                    832:         if (ell > 0) {printf("ell+ "); fflush(NULL);
                    833:         }else {
                    834:                  printf("ell0 "); fflush(NULL);
                    835:           if ((*mmLarger)(r,r0) >= 1) {
                    836:             printf("error in reduction.");
                    837:             printf("   r0=%s\n",POLYToString(r0,'*',1));
                    838:             printf("==>r=%s\n",POLYToString(r,'*',1));
                    839:             getchar(); getchar();
                    840:           }
                    841:         }
                    842:       }
                    843:
1.8       takayama  844:       if (r ISZERO) goto ss1;
                    845:       r = ecartDivideSv(r,&se); /* r = r/s^? */
                    846:     }
                    847:   }while (ell >= 0);
                    848:
                    849:  ss1: ;
1.5       takayama  850:
1.1       takayama  851:   r = goDeHomogenizeS(r);
1.5       takayama  852:
1.1       takayama  853:   return(r);
1.10      takayama  854: }
1.20      takayama  855:
                    856: #ifdef EXPERIMENT
                    857: static struct ecartReducer ecartFindReducer(POLY r,struct gradedPolySet *gset,
                    858:                                             struct ecartPolyArray *epa)
                    859: {
                    860:   int grd;
                    861:   struct polySet *set;
                    862:   int minGrade = 0;
                    863:   int minGseti = 0;
                    864:   int minGgi   = 0;
                    865:   int ell1 = LARGE;
                    866:   int ell2 = LARGE;
                    867:   int ell;
                    868:   int i;
                    869:   struct ecartReducer er;
                    870:   int ellmin ;
                    871:   ell1 = ell2 = ellmin = LARGE;
                    872:   /* Try to find a reducer in gset; */
                    873:   grd = 0;
                    874:   while (grd < gset->maxGrade) {
                    875:     set = gset->polys[grd];
                    876:     for (i=0; i<set->size; i++) {
                    877:       if (set->gh[i] == POLYNULL) {
                    878:         /* goHomogenize set->gh[i] */
                    879:           if (EcartAutomaticHomogenization) {
                    880:               set->gh[i] = goHomogenize11(set->g[i],DegreeShifto_vec,DegreeShifto_size,-1,1);
                    881:           }else{
                    882:               set->gh[i] = set->g[i];
                    883:           }
                    884:       }
                    885:       ell = ecartGetEll(r,set->gh[i]);
                    886:       if ((ell>=0) && (ell < ell1)) {
                    887:         ell1 = ell;
                    888:         minGrade = grd; minGseti=i;
                    889:       }
                    890:     }
                    891:     grd++;
                    892:   }
                    893:   if (epa != NULL) {
                    894:     /* Try to find in the second group. */
                    895:     for (i=0; i< epa->size; i++) {
                    896:       ell = ecartGetEll(r,(epa->pa)[i]);
                    897:       if ((ell>=0) && (ell < ell2)) {
                    898:         ell2 = ell;
                    899:         minGgi = i;
                    900:       }
                    901:     }
                    902:   }
                    903:
                    904:   if ((DebugReductionRed&1) || (DebugReductionEcart&1)) {
                    905:     printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi);
                    906:   }
                    907:
                    908: #ifdef EXPERIMENTAL1
                    909:   if (Sugar) {  /* experimental */
                    910:        if (ell1 <= ell2) {
                    911:          if (ell1 == LARGE) {
                    912:                er.ell = -1;
                    913:                return er;
                    914:          }else{
                    915:                int new_s;
                    916:                er.ell = ell1;
                    917:                er.first = 1;
                    918:                er.grade = minGrade;
                    919:                er.gseti = minGseti;
                    920:                /* reduce if and only if Ecart_sugarGrade does not increase. */
                    921:                new_s = grade_gen(r)-grade_gen(gset->polys[minGrade]->gh[minGseti]);
                    922:                if (new_s + minGrade <= Ecart_sugarGrade+1) {
                    923:                  return er;
                    924:                }else{
                    925:                  printf("new_s=%d, minGrade=%d, sugarGrade=%d\n",new_s,minGrade,
                    926:                                 Ecart_sugarGrade);
                    927:                  er.ell = -1;
                    928:                  return er;
                    929:                }
                    930:          }
                    931:        }else{
                    932:       er.ell = ell2;
                    933:       er.first = 0;
                    934:       er.ggi = minGgi;
                    935:          return er;
                    936:        }
                    937:   }
                    938: #endif
                    939:
                    940:   if ((ell1 == LARGE) &&  (ell2 == LARGE)) {
                    941:        er.ell = -1;
                    942:        return er;
                    943:   }
                    944:
                    945:   if (ell1 < ell2) {
                    946:        ellmin = ell1;
                    947:   }else{
                    948:        ellmin = ell2;
                    949:   }
                    950:   {
                    951: #define M 100
                    952:        int aminGrade[M];
                    953:        int aminGseti[M];
                    954:        int aminGgi[M];
                    955:        int sp1,sp2;
                    956:        int i;
                    957:        int gmin;
                    958:        int gtmp;
                    959:        sp1 = sp2 = 0;
                    960:
                    961:        if (ell1 == ellmin) {
                    962:          grd = 0;
                    963:          while (grd < gset->maxGrade) {
                    964:                set = gset->polys[grd];
                    965:                for (i=0; i<set->size; i++) {
                    966:                  ell = ecartGetEll(r,set->gh[i]);
                    967:                  if (ell == ellmin) {
                    968:                        aminGrade[sp1] = grd; aminGseti[sp1]=i;
                    969:                        sp1++;
                    970:                        if (sp1 >= M) {
                    971:                          fprintf(stderr,"aminGrade, buffer overflow. sp1 is set to 1.\n");
                    972:                          sp1 = 1;
                    973:                        }
                    974:                  }
                    975:                }
                    976:                grd++;
                    977:          }
                    978:        }
                    979:        if (ell2 == ellmin) {
                    980:          if (epa != NULL) {
                    981:                for (i=0; i< epa->size; i++) {
                    982:                  ell = ecartGetEll(r,(epa->pa)[i]);
                    983:                  if (ell == ellmin) {
                    984:                        aminGgi[sp2] = i;
                    985:                        sp2++;
                    986:                        if (sp2 >= M) {
                    987:                          fprintf(stderr,"aminGgi, buffer overflow. sp2 is set to 1.\n");
                    988:                          sp2 = 1;
                    989:                        }
                    990:                  }
                    991:                }
                    992:          }
                    993:        }
                    994:        /* print summary */
                    995:        printf("summary -----------------------\n");
                    996:        for (i=0; i<sp1; i++) {
                    997:          printf("i=%d: minGrade=%d, minGseti=%d,",i,aminGrade[i],aminGseti[i]);
                    998:          printf("grade=%d\n",grade_gen(gset->polys[aminGrade[i]]->gh[aminGseti[i]]));
                    999:        }
                   1000:        gmin = LARGE;
                   1001:        for (i=0; i<sp2; i++) {
                   1002:          printf("i=%d: minGgi=%d,",i,aminGgi[i]);
                   1003:          gtmp = grade_gen((epa->pa)[aminGgi[i]]);
                   1004:          printf("grade=%d\n",gtmp);
                   1005:          if (gtmp < gmin) {
                   1006:                gmin = gtmp;
                   1007:                minGgi = aminGgi[i];
                   1008:          }
                   1009:        }
                   1010:
                   1011:        if (ell1 <= ell2) {
                   1012:       er.ell = ell1;
                   1013:       er.first = 1;
                   1014:       er.grade = minGrade;
                   1015:       er.gseti = minGseti;
                   1016:       return er;
                   1017:        }else{
                   1018:       er.ell = ell2;
                   1019:       er.first = 0;
                   1020:       er.ggi = minGgi;
                   1021:       return er;
                   1022:        }
                   1023:
                   1024:   }
                   1025:
                   1026:   if (ell1 <= ell2) {
                   1027:       er.ell = ell1;
                   1028:       er.first = 1;
                   1029:       er.grade = minGrade;
                   1030:       er.gseti = minGseti;
                   1031:       return er;
                   1032:   }else{
                   1033:       er.ell = ell2;
                   1034:       er.first = 0;
                   1035:       er.ggi = minGgi;
                   1036:       return er;
                   1037:   }
                   1038: }
                   1039: #endif
1.10      takayama 1040:

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