=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Hgfs.c,v retrieving revision 1.27 retrieving revision 1.32 diff -u -p -r1.27 -r1.32 --- OpenXM_contrib2/asir2000/engine/Hgfs.c 2002/11/01 05:43:35 1.27 +++ OpenXM_contrib2/asir2000/engine/Hgfs.c 2015/08/08 14:19:41 1.32 @@ -1,8 +1,10 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.26 2002/10/25 02:43:40 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.31 2015/03/16 00:08:32 noro Exp $ */ #include "ca.h" #include "inline.h" +int debug_sfbfctr; + void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1); void extractcoefbm(BM f,int dx,UM r); @@ -79,6 +81,7 @@ void gensqfrsfum(UM p,struct oDUM *dc) { int n,i,j,d,mod; UM t,s,g,f,f1,b; + GFS u,v; if ( (n = DEG(p)) == 1 ) { dc[0].f = UMALLOC(DEG(p)); cpyum(p,dc[0].f); dc[0].n = 1; @@ -117,8 +120,11 @@ void gensqfrsfum(UM p,struct oDUM *dc) break; else { DEG(s) = DEG(t)/mod; - for ( j = 0; j <= DEG(t); j++ ) - COEF(s)[j] = COEF(t)[j*mod]; + for ( j = 0; j <= DEG(t); j++ ) { + iftogfs(COEF(t)[j*mod],&u); + pthrootgfs(u,&v); + COEF(s)[j] = v?FTOIF(CONT(v)):0; + } cpyum(s,b); d *= mod; } } @@ -506,21 +512,21 @@ void canzassf(UM f,int d,UM *r) /* Hensel related functions */ -int sfberle(VL,P,int,GFS *,DCP *); +int sfberle(V,V,P,int,GFS *,DCP *); void sfgcdgen(P,ML,ML *); void sfhenmain2(BM,UM,UM,int,BM *); void ptosfbm(int,P,BM); +void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp); /* f = f(x,y) */ -void sfhensel(int count,P f,V x,int degbound,GFS *evp,P *sfp,ML *listp) +void sfhensel(int count,P f,V x,V y,int degbound,GFS *evp,P *sfp,ML *listp) { int i; int fn; ML rlist; BM fl; VL vl,nvl; - V y; int dx,dy,bound; GFS ev; P f1,t,c,sf; @@ -534,14 +540,15 @@ void sfhensel(int count,P f,V x,int degbound,GFS *evp, reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&f1); vl = nvl; f = f1; } - y = vl->next->v; + if ( vl->next ) + y = vl->next->v; dx = getdeg(x,f); dy = getdeg(y,f); if ( dx == 1 ) { *listp = rlist = MLALLOC(1); rlist->n = 1; rlist->c[0] = 0; return; } - fn = sfberle(vl,f,count,&ev,&dc); + fn = sfberle(x,y,f,count,&ev,&dc); if ( fn <= 1 ) { /* fn == 0 => short of evaluation points */ *listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0; @@ -604,10 +611,12 @@ void sfhensel(int count,P f,V x,int degbound,GFS *evp, q = W_UMALLOC(dx); rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound; - fprintf(asir_out,"%d candidates\n",fn); + if ( debug_sfbfctr ) + fprintf(asir_out,"%d candidates\n",fn); init_eg(&eg_hensel); for ( i = 0; i < fn-1; i++ ) { - fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n", + if ( debug_sfbfctr ) + fprintf(asir_out,"deg(fm) = %d, deg(gm[%d]) = %d\n", DEG(fm),i,DEG(gm[i])); init_eg(&eg_hensel_t); get_eg(&tmp0); @@ -619,11 +628,15 @@ void sfhensel(int count,P f,V x,int degbound,GFS *evp, cpyum(hm,fm); get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1); add_eg(&eg_hensel,&tmp0,&tmp1); - print_eg("Hensel",&eg_hensel_t); + if ( debug_sfbfctr) { + print_eg("Hensel",&eg_hensel_t); + fprintf(asir_out,"\n"); + } + } + if ( debug_sfbfctr) { + print_eg("Hensel total",&eg_hensel); fprintf(asir_out,"\n"); } - print_eg("Hensel total",&eg_hensel); - fprintf(asir_out,"\n"); /* finally, fl must be the lift of gm[fn-1] */ rlist->c[i] = fl; @@ -640,20 +653,20 @@ void sfhensel(int count,P f,V x,int degbound,GFS *evp, /* main variable of f = x */ -int sfberle(VL vl,P f,int count,GFS *ev,DCP *dcp) +int sfberle(V x,V y,P f,int count,GFS *ev,DCP *dcp) { UM wf,wf1,wf2,wfs,gcd; int fn,n; GFS m,fm; DCP dc,dct,dc0; - VL nvl; - V x,y; + VL vl; P lc,lc0,f0; Obj obj; int j,q,index,i; - clctv(vl,f,&nvl); vl = nvl; - x = vl->v; y = vl->next->v; + NEWVL(vl); vl->v = x; + NEWVL(NEXT(vl)); NEXT(vl)->v = y; + NEXT(NEXT(vl)) =0; simp_ff((Obj)f,&obj); f = (P)obj; n = QTOS(DEG(DC(f))); wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n); @@ -756,9 +769,11 @@ void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp) wa = W_UMALLOC(2*dx); wb = W_UMALLOC(2*dx); /* XXX */ eucsfum(w1,w2,wa,wb); - fprintf(stderr,"dy=%d\n",dy); + if ( debug_sfbfctr) + fprintf(stderr,"dy=%d\n",dy); for ( k = 1; k <= dy; k++ ) { - fprintf(stderr,"."); + if ( debug_sfbfctr) + fprintf(stderr,"."); /* at this point, f = gk*hk mod y^k */ @@ -814,11 +829,15 @@ void sfhenmain2(BM f,UM g0,UM h0,int dy,BM *gp) cpyum(wg1,COEF(gk)[k]); cpyum(wh1,COEF(hk)[k]); } - fprintf(stderr,"\n"); + if ( debug_sfbfctr) + fprintf(stderr,"\n"); *gp = gk; DEG(f) = dy; for ( i = 0; i <= dy; i++ ) cpyum(COEF(hk)[i],COEF(f)[i]); +#if defined(__MINGW32__) || defined(__MINGW64__) + fflush(stderr); +#endif } /* a0*g+b0*h = 1 mod y -> a*g+b*h = 1 mod y^(dy+1) */ @@ -860,9 +879,11 @@ void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp) mulsfbm(a,g,c); mulsfbm(b,h,wz0); addtosfbm(wz0,c); COEF(COEF(c)[0])[0] = 0; - fprintf(stderr,"dy=%d\n",dy); + if ( debug_sfbfctr) + fprintf(stderr,"dy=%d\n",dy); for ( k = 1; k <= dy; k++ ) { - fprintf(stderr,"."); + if ( debug_sfbfctr) + fprintf(stderr,"."); /* at this point, a*g+b*h = 1 mod y^k, c = a*g+b*h-1 */ @@ -901,11 +922,15 @@ void sfexgcd_by_hensel(BM g,BM h,int dy,BM *ap,BM *bp) cpyum(wa1,COEF(a)[k]); cpyum(wb1,COEF(b)[k]); } - fprintf(stderr,"\n"); + if ( debug_sfbfctr) + fprintf(stderr,"\n"); DEG(a) = dy; DEG(b) = dy; *ap = a; *bp = b; +#if defined(__MINGW32__) || defined(__MINGW64__) + fflush(stderr); +#endif } /* fl->c[i] = coef_y(f,i) */ @@ -986,7 +1011,7 @@ void sfsqfr(P f,DCP *dcp) } else if ( !NEXT(vl) ) sfusqfr(f,dcp); else - sqfrsf(f,dcp); + sqfrsf(vl,f,dcp); } void sfusqfr(P f,DCP *dcp) @@ -1077,7 +1102,7 @@ void sfbfctr(P f,V x,V y,int degbound,DCP *dcp) int dx,dy; /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */ - sfhensel(5,f,x,degbound,&ev,&sf,&list); + sfhensel(5,f,x,y,degbound,&ev,&sf,&list); if ( list->n == 0 ) error("sfbfctr : short of evaluation points"); else if ( list->n == 1 ) { @@ -1111,7 +1136,7 @@ void sfbfctr_shift(P f,V x,V y,int degbound,GFS *evp,P int dx,dy; /* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */ - sfhensel(5,f,x,degbound,&ev,&sf,&list); + sfhensel(5,f,x,y,degbound,&ev,&sf,&list); if ( list->n == 0 ) error("sfbfctr_shift : short of evaluation points"); else if ( list->n == 1 ) { @@ -1199,10 +1224,11 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp) } } - fprintf(stderr,"np = %d\n",np); + if ( debug_sfbfctr) + fprintf(stderr,"np = %d\n",np); dtok = 0; for ( g = f, k = 1, dcf = dcf0 = 0, win[0] = 1, --np, z = 0; ; z++ ) { - if ( !(z % 1000) ) fprintf(stderr,"."); + if ( debug_sfbfctr && !(z % 1000) ) fprintf(stderr,"."); dt = sfdegtest(dy,bound,d1c,k,win); if ( dt ) dtok++; @@ -1267,9 +1293,13 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp) for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1; } - fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok); + if ( debug_sfbfctr ) + fprintf(stderr,"total %d, omitted by degtest %d\n",z,z-dtok); NEXTDC(dcf0,dcf); COEF(dcf) = g; DEG(dcf) = ONE; NEXT(dcf) = 0; *dcp = dcf0; +#if defined(__MINGW32__) || defined(__MINGW64__) + fflush(stderr); +#endif } void extractcoefbm(BM f,int dx,UM r)