=================================================================== RCS file: /home/cvs/OpenXM/src/kan96xx/Kan/ecart.c,v retrieving revision 1.10 retrieving revision 1.19 diff -u -p -r1.10 -r1.19 --- OpenXM/src/kan96xx/Kan/ecart.c 2003/08/20 05:18:35 1.10 +++ OpenXM/src/kan96xx/Kan/ecart.c 2003/09/20 09:57:29 1.19 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM/src/kan96xx/Kan/ecart.c,v 1.9 2003/08/20 01:39:17 takayama Exp $ */ +/* $OpenXM: OpenXM/src/kan96xx/Kan/ecart.c,v 1.18 2003/09/12 02:52:50 takayama Exp $ */ #include #include "datatype.h" #include "extern2.h" @@ -36,17 +36,21 @@ static POLY reduction_ecart1(POLY r,struct gradedPolyS static POLY reduction_ecart1_mod(POLY r,struct gradedPolySet *gset); static POLY ecartCheckSyz0(POLY cf,POLY r_0,POLY syz, struct gradedPolySet *gg,POLY r); -static int shouldReduceContent(POLY f,int ss); +static void ecartCheckSyz0_printinfo(POLY cf,POLY r_0,POLY syz, + struct gradedPolySet *gg,POLY r); extern int DebugReductionRed; extern int TraceLift; struct ring *TraceLift_ringmod; extern DoCancel; int DebugReductionEcart = 0; +extern DebugContentReduction; +extern int Sugar; /* This is used for goHomogenization */ extern int DegreeShifto_size; extern int *DegreeShifto_vec; +int Ecart_sugarGrade; /* It is used reduction_ecart() and ecartFindReducer() to determine if we homogenize in this function */ @@ -67,8 +71,20 @@ static POLY ecartDivideSv(POLY r,int *d) { } r = r->next; } + if (k > 0) { *d = k; + return ppMult(cxx(1,0,-k,f->m->ringp),f); + }else{ + return f; + } + + /* Do not do the below. It caused a bug. cf. misc-2003/07/ecart/b4.sm1 test2. + Note. 2003.8.26 + */ + /* + if (k > 0) { + *d = k; f = r; while (r != POLYNULL) { r->m->e[0].x -= k; @@ -76,6 +92,7 @@ static POLY ecartDivideSv(POLY r,int *d) { } } return f; + */ } static int ecartGetEll(POLY f,POLY g) { @@ -193,9 +210,39 @@ static struct ecartReducer ecartFindReducer(POLY r,str } } - if (DebugReductionRed || (DebugReductionEcart&1)) { + if ((DebugReductionRed&1) || (DebugReductionEcart&1)) { printf("ecartFindReducer(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d\n",ell1,ell2,minGrade,minGseti,minGgi); } + if (Sugar) { /* experimental */ + if (ell1 <= ell2) { + if (ell1 == LARGE) { + er.ell = -1; + return er; + }else{ + int new_s; + er.ell = ell1; + er.first = 1; + er.grade = minGrade; + er.gseti = minGseti; + /* reduce if and only if Ecart_sugarGrade does not increase. */ + new_s = grade_gen(r)-grade_gen(gset->polys[minGrade]->gh[minGseti]); + if (new_s + minGrade <= Ecart_sugarGrade) { + return er; + }else{ + printf("new_s=%d, minGrade=%d, sugarGrade=%d\n",new_s,minGrade, + Ecart_sugarGrade); + er.ell = -1; + return er; + } + } + }else{ + er.ell = ell2; + er.first = 0; + er.ggi = minGgi; + return er; + } + } + if (ell1 <= ell2) { if (ell1 == LARGE) { er.ell = -1; @@ -240,7 +287,36 @@ static POLY ecartCheckSyz0(POLY cf,POLY r_0,POLY syz, return f; } +static void ecartCheckSyz0_printinfo(POLY cf,POLY r_0,POLY syz, + struct gradedPolySet *gg,POLY r) +{ + POLY f; + int grd,i; + POLY q; + struct coeff *c; + fprintf(stderr,"cf=%s\n",POLYToString(cf,'*',1)); + fprintf(stderr,"r_0=%s\n",POLYToString(r_0,'*',1)); + fprintf(stderr,"syz=%s\n",POLYToString(syz,'*',1)); + fprintf(stderr,"r=%s\n",POLYToString(r,'*',1)); + f = ppMult(cf,r_0); + while (syz != POLYNULL) { + grd = syz->m->e[0].x; + i = syz->m->e[0].D; + c = syz->coeffp; + if (c->tag == POLY_COEFF) { + q = c->val.f; + }else{ + q = POLYNULL; + } + fprintf(stderr,"[grd,idx]=[%d,%d], %s\n",grd,i, + POLYToString(((gg->polys[grd])->g)[i],'*',1)); + /* f = ppAdd(f,ppMult(q,((gg->polys[grd])->g)[i])); */ + syz = syz->next; + } + /* f = ppSub(f,r); */ +} + POLY reduction_ecart(r,gset,needSyz,syzp) POLY r; struct gradedPolySet *gset; @@ -248,13 +324,17 @@ POLY reduction_ecart(r,gset,needSyz,syzp) struct syz0 *syzp; /* set */ { POLY rn; + if (TraceLift && needSyz) { + warningGradedSet("TraceLift cannot be used to get syzygy. TraceLift is turned off.\n"); + TraceLift = 0; + } if (TraceLift) { if (EcartAutomaticHomogenization) { if (TraceLift_ringmod == NULL) { warningPoly("reduction_ecart: TraceLift_ringmod is not set.\n"); return reduction_ecart1(r,gset,needSyz,syzp); } - rn = reduction_ecart1_mod(r,gset); /* BUG: syzygy is not obtained. */ + rn = reduction_ecart1_mod(r,gset); if (rn == POLYNULL) return rn; else return reduction_ecart1(r,gset,needSyz,syzp); }else{ @@ -293,6 +373,8 @@ static POLY reduction_ecart0(r,gset,needSyz,syzp) POLY cf_o; POLY syz_o; POLY r_0; + int se; + struct coeff *cont; extern struct ring *CurrentRingp; struct ring *rp; @@ -312,9 +394,11 @@ static POLY reduction_ecart0(r,gset,needSyz,syzp) } } + if (DoCancel && (r != POLYNULL)) shouldReduceContent(r,1); + if (DebugReductionEcart&1) printf("--------------------------------------\n"); do { - if (DebugReductionRed) printf("r=%s\n",POLYToString(r,'*',1)); + if (DebugReductionRed&1) printf("r=%s\n",POLYToString(r,'*',1)); if (DebugReductionEcart & 1) printf("r=%s+...\n",POLYToString(head(r),'*',1)); ells = ecartFindReducer(r,gset,gg); ell = ells.ell; @@ -333,8 +417,18 @@ static POLY reduction_ecart0(r,gset,needSyz,syzp) if (DebugReductionEcart & 4) printf("#"); pp = (gg->pa)[ells.ggi]; } - if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ + if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ r = (*reduction1)(r,pp,needSyz,&cc,&cg); + + if (DoCancel && (r != POLYNULL)) { + if (shouldReduceContent(r,0)) { + r = reduceContentOfPoly(r,&cont); + shouldReduceContent(r,1); + if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("CONT=%s ",coeffToString(cont)); + } + } + + if (needSyz) { if (ells.first) { if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp)); @@ -351,12 +445,13 @@ static POLY reduction_ecart0(r,gset,needSyz,syzp) syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o)); /* Note. 2003.07.19 */ } - if (DebugReductionRed) { + if (DebugReductionRed && needSyz) { POLY tp; tp = ecartCheckSyz0(cf,r_0,syz,gset,r); if (tp != POLYNULL) { fprintf(stderr,"reduction_ecart0(): sygyzy is broken. Return the Current values.\n"); - fprintf(stderr,"%s\n",POLYToString(tp,'*',1)); + fprintf(stderr,"tp=%s\n",POLYToString(tp,'*',1)); + ecartCheckSyz0_printinfo(cf,r_0,syz,gset,r); syzp->cf = cf; syzp->syz = syz; return r; @@ -364,16 +459,32 @@ static POLY reduction_ecart0(r,gset,needSyz,syzp) } } if (r ISZERO) goto ss; - /*r = ecartDivideSv(r,&se); r = r/s^? Don't do this. */ + + /* r = r/s^? Don't do this?? */ + r = ecartDivideSv(r,&se); + if (needSyz && (se > 0)) { + POLY tt; + tt = cxx(1,0,-se,rp); + cf = ppMult(tt,cf); + syz = cpMult(toSyzCoeff(tt),syz); + } + } }while (ell >= 0); ss: ; if (needSyz) { - syzp->cf = cf; /* cf is in the CurrentRingp */ - syzp->syz = syz; /* syz is in the SyzRingp */ + syzp->cf = cf; /* cf is in the ring of r */ + syzp->syz = syz; /* syz is in the syzRing of r */ } + if (DoCancel && (r != POLYNULL)) { + if (r->m->ringp->p == 0) { + r = reduceContentOfPoly(r,&cont); + if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("cont=%s ",coeffToString(cont)); + } + } + return(r); } @@ -400,6 +511,10 @@ static POLY reduction_ecart1(r,gset,needSyz,syzp) struct ecartPolyArray *gg; POLY pp; int ell; + POLY cf_o; + POLY syz_o; + POLY r_0; + POLY r0; int se; struct coeff *cont; @@ -407,6 +522,7 @@ static POLY reduction_ecart1(r,gset,needSyz,syzp) struct ring *rp; extern struct ring *SmallRingp; + r_0 = r; gg = NULL; if (needSyz) { if (r ISZERO) { rp = CurrentRingp; } else { rp = r->m->ringp; } @@ -427,14 +543,18 @@ static POLY reduction_ecart1(r,gset,needSyz,syzp) if (DebugReductionEcart&1) printf("=======================================\n"); do { - if (DebugReductionRed) printf("(ecart1(d)) r=%s\n",POLYToString(r,'*',1)); + if (DebugReductionRed & 1) printf("(ecart1(d)) r=%s\n",POLYToString(r,'*',1)); if (DebugReductionEcart & 1) printf("r=%s+,,,\n",POLYToString(head(r),'*',1)); ells = ecartFindReducer(r,gset,gg); ell = ells.ell; if (ell > 0) { if (DebugReductionEcart & 2) printf("%"); - gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL); + if (needSyz) { + gg = ecartPutPolyInG(r,gg,cf,syz); + }else{ + gg = ecartPutPolyInG(r,gg,POLYNULL,POLYNULL); + } } if (ell >= 0) { if (ells.first) { @@ -443,52 +563,122 @@ static POLY reduction_ecart1(r,gset,needSyz,syzp) if (DebugReductionEcart & 4) {printf("+"); fflush(NULL);} pp = (gg->pa)[ells.ggi]; } - if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ + if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ + r0 = r; r = (*reduction1)(r,pp,needSyz,&cc,&cg); + if (DebugReductionEcart & 8) { + if (ell > 0) {printf("ell+ "); fflush(NULL); + }else { + printf("ell0 "); fflush(NULL); + if ((*mmLarger)(r,r0) >= 1) { + printf("error in reduction."); + printf(" r0=%s\n",POLYToString(r0,'*',1)); + printf("==>r=%s\n",POLYToString(r,'*',1)); + getchar(); getchar(); + } + } + } - if (DoCancel && (r != POLYNULL)) { /* BUG: syzygy should be corrected. */ + if (DoCancel && (r != POLYNULL)) { if (shouldReduceContent(r,0)) { r = reduceContentOfPoly(r,&cont); shouldReduceContent(r,1); - if (DebugReductionEcart || DebugReductionRed) printf("CONT=%d ",coeffToString(cont)); + if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("CONT=%s ",coeffToString(cont)); } } if (needSyz) { if (ells.first) { + if (ell > 0) cc = ppMult(cc,cxx(1,0,ell,rp)); cf = ppMult(cc,cf); syz = cpMult(toSyzCoeff(cc),syz); - syz = ppAddv(syz,toSyzPoly(cg,ells.grade,ells.gseti)); + syz = ppAdd(syz,toSyzPoly(cg,ells.grade,ells.gseti)); }else{ - fprintf(stderr,"It has not yet implemented.\n"); - exit(10); - /* BUG: not yet */ + if (ell >0) cc = ppMult(cc,cxx(1,0,ell,rp)); + cf_o = (gg->cf)[ells.ggi]; + syz_o = (gg->syz)[ells.ggi]; + cf = ppMult(cc,cf); + cf = ppAdd(cf,ppMult(cg,cf_o)); + syz = cpMult(toSyzCoeff(cc),syz); + syz = ppAdd(syz,cpMult(toSyzCoeff(cg),syz_o)); + /* Note. 2003.07.19 */ } } + + if (DebugReductionRed && needSyz) { + POLY tp; + tp = ecartCheckSyz0(cf,r_0,syz,gset,r); + tp = goDeHomogenizeS(tp); + if (tp != POLYNULL) { + fprintf(stderr,"Error: reduction_ecart1(): sygyzy is broken. Return the Current values.\n"); + fprintf(stderr,"%s\n",POLYToString(tp,'*',1)); + syzp->cf = cf; + syzp->syz = syz; + return r; + } + } + if (r ISZERO) goto ss1; r = ecartDivideSv(r,&se); /* r = r/s^? */ + + if (needSyz && (se > 0)) { /* It may not necessary because of dehomo*/ + POLY tt; + /*printf("!1/H!"); fflush(NULL);*/ /* misc-2003/ecart/t1.sm1 foo4 */ + tt = cxx(1,0,-se,rp); + cf = ppMult(tt,cf); + syz = cpMult(toSyzCoeff(tt),syz); + } + + /* For debug */ + if (DebugReductionRed && needSyz) { + POLY tp; + tp = ecartCheckSyz0(cf,r_0,syz,gset,r); + tp = goDeHomogenizeS(tp); + if (tp != POLYNULL) { + fprintf(stderr,"Error: reduction_ecart1() after divide: sygyzy is broken. Return the Current values.\n"); + fprintf(stderr,"%s\n",POLYToString(tp,'*',1)); + syzp->cf = cf; + syzp->syz = syz; + return r; + } + } + } }while (ell >= 0); ss1: ; if (needSyz) { + /* dehomogenize the syzygy. BUG, this may be inefficient. */ + cf = goDeHomogenizeS(cf); + syz = goDeHomogenizeS(syz); + /*printf("cf=%s\n",POLYToString(cf,'*',1)); + printf("syz=%s\n",POLYToString(syz,'*',1));*/ syzp->cf = cf; /* cf is in the CurrentRingp */ syzp->syz = syz; /* syz is in the SyzRingp */ - /* BUG: dehomogenize the syzygy */ - fprintf(stderr,"It has not yet implemented.\n"); - exit(10); } r = goDeHomogenizeS(r); - if (DoCancel && (r != POLYNULL)) { /* BUG: syzygy should be corrected. */ + if (DoCancel && (r != POLYNULL)) { if (r->m->ringp->p == 0) { - if (coeffSizeMin(r) >= DoCancel) { - r = reduceContentOfPoly(r,&cont); - if (DebugReductionEcart || DebugReductionRed) printf("cont=%d ",coeffToString(cont)); - } + r = reduceContentOfPoly(r,&cont); + if (DebugReductionEcart || DebugReductionRed || DebugContentReduction) printf("cont=%s ",coeffToString(cont)); } } + /* For debug */ + if (DebugReductionRed && needSyz) { + POLY tp; + tp = ecartCheckSyz0(cf,r_0,syz,gset,r); + tp = goDeHomogenizeS(tp); + if (tp != POLYNULL) { + fprintf(stderr,"Error: reduction_ecart1() last step: sygyzy is broken. Return the Current values.\n"); + fprintf(stderr,"%s\n",POLYToString(tp,'*',1)); + syzp->cf = cf; + syzp->syz = syz; + return r; + } + } + return(r); } @@ -546,7 +736,7 @@ static struct ecartReducer ecartFindReducer_mod(POLY r } } - if (DebugReductionRed || (DebugReductionEcart&1)) { + if ((DebugReductionRed & 1)|| (DebugReductionEcart&1)) { printf("ecartFindReducer_mod(): ell1=%d, ell2=%d, minGrade=%d, minGseti=%d, minGgi=%d, p=%d\n",ell1,ell2,minGrade,minGseti,minGgi,TraceLift_ringmod->p); } if (ell1 <= ell2) { @@ -580,6 +770,7 @@ static POLY reduction_ecart1_mod(r,gset) struct ecartReducer ells; struct ecartPolyArray *gg; POLY pp; + POLY r0; int ell; int se; @@ -601,13 +792,13 @@ static POLY reduction_ecart1_mod(r,gset) KshowRing(TraceLift_ringmod); **/ r = modulop(r,TraceLift_ringmod); - rp = r->m->ringp; /* reset rp */ + if (r != POLYNULL) rp = r->m->ringp; /* reset rp */ /* printf("r=%s (mod p)\n",POLYToString(head(r),'*',1)); **/ if (DebugReductionEcart&1) printf("=====================================mod\n"); do { - if (DebugReductionRed) printf("(ecart1_mod(d)) r=%s\n",POLYToString(r,'*',1)); + if (DebugReductionRed & 1) printf("(ecart1_mod(d)) r=%s\n",POLYToString(r,'*',1)); if (DebugReductionEcart & 1) printf("r=%s+,,,\n",POLYToString(head(r),'*',1)); ells = ecartFindReducer_mod(r,gset,gg); @@ -623,8 +814,24 @@ static POLY reduction_ecart1_mod(r,gset) if (DebugReductionEcart & 4) {printf("+"); fflush(NULL);} pp = (gg->pa)[ells.ggi]; } - if (ell > 0) r = mpMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ + if (ell > 0) r = ppMult(cxx(1,0,ell,rp),r); /* r = s^ell r */ + + r0 = r; r = (*reduction1)(r,pp,0,&cc,&cg); + + if (DebugReductionEcart & 8) { + if (ell > 0) {printf("ell+ "); fflush(NULL); + }else { + printf("ell0 "); fflush(NULL); + if ((*mmLarger)(r,r0) >= 1) { + printf("error in reduction."); + printf(" r0=%s\n",POLYToString(r0,'*',1)); + printf("==>r=%s\n",POLYToString(r,'*',1)); + getchar(); getchar(); + } + } + } + if (r ISZERO) goto ss1; r = ecartDivideSv(r,&se); /* r = r/s^? */ } @@ -637,20 +844,3 @@ static POLY reduction_ecart1_mod(r,gset) return(r); } -static int shouldReduceContent(POLY f,int ss) { - static int prevSize = 1; - int size; - if (f == POLYNULL) return 0; - if (f->m->ringp->p != 0) return 0; - if (f->coeffp->tag != MP_INTEGER) return 0; - size = mpz_size(f->coeffp->val.bigp); - if (ss > 0) { - prevSize = size; - return 0; - } - if (size > 2*prevSize) { - return 1; - }else{ - return 0; - } -}