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

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

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