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

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

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