===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Hgfs.c,v
retrieving revision 1.18
retrieving revision 1.20
diff -u -p -r1.18 -r1.20
--- OpenXM_contrib2/asir2000/engine/Hgfs.c	2001/10/09 01:36:10	1.18
+++ OpenXM_contrib2/asir2000/engine/Hgfs.c	2001/10/30 10:24:35	1.20
@@ -1,8 +1,9 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.17 2001/09/03 07:01:06 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/Hgfs.c,v 1.19 2001/10/30 07:25:58 noro Exp $ */
 
 #include "ca.h"
 
 void lnfsf(int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1);
+void extractcoefbm(BM f,int dx,UM r);
 
 int comp_dum(DUM a,DUM b)
 {
@@ -517,7 +518,7 @@ void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *li
 	BM fl;
 	VL vl,nvl;
 	V y;
-	int dx,dy;
+	int dx,dy,bound;
 	GFS ev;
 	P f1,t,c,sf;
 	DCP dc;
@@ -554,9 +555,14 @@ void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *li
 		ptosfum(dc->c,gm[i]);
 	}
 
+	/* set bound */
+	/* g | f, lc_y(g) = lc_y(f) => deg_y(g) <= deg_y(f) */
+	/* so, bound = dy is sufficient, but we use slightly large value */
+	bound = dy+2;
+
 	/* f(x,y) -> f(x,y+ev) */
-	fl = BMALLOC(dx,dy);
-	ptosfbm(dy,f,fl);
+	fl = BMALLOC(dx,bound);
+	ptosfbm(bound,f,fl);
 	if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
 
 	/* sf = f(x+ev) */
@@ -568,7 +574,7 @@ void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *li
 	hm = W_UMALLOC(dx);
 
 	q = W_UMALLOC(dx);
-	rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = dy;
+	rlist = MLALLOC(fn); rlist->n = fn; rlist->bound = bound;
 	fprintf(asir_out,"%d candidates\n",fn);
 	init_eg(&eg_hensel);
 	for ( i = 0; i < fn-1; i++ ) {
@@ -578,9 +584,9 @@ void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *li
 		get_eg(&tmp0);
 		/* fl = gm[i]*hm mod y */
 		divsfum(fm,gm[i],hm);
-		/* fl is replaced by the cofactor of gk mod y^dy */
+		/* fl is replaced by the cofactor of gk mod y^bound */
 		/* rlist->c[i] = gk */
-		sfhenmain2(fl,gm[i],hm,dy,(BM *)&rlist->c[i]);
+		sfhenmain2(fl,gm[i],hm,bound,(BM *)&rlist->c[i]);
 		cpyum(hm,fm);
 		get_eg(&tmp1); add_eg(&eg_hensel_t,&tmp0,&tmp1);
 		add_eg(&eg_hensel,&tmp0,&tmp1);
@@ -959,18 +965,19 @@ void sfbfctr(P f,V x,V y,DCP *dcp)
 void sfdtest(P f,ML list,V x,V y,DCP *dcp)
 {
 	int np,dx,dy;
-	int i,j,k;
+	int i,j,k,bound;
 	int *win;
 	P g,lcg,factor,cofactor,lcyx;
 	P csum;
 	DCP dcf,dcf0,dc;
 	BM *c;
 	BM lcy;
-	UM lcg0;
+	UM lcg0,lcy0,w;
+	UM *d1c;
 	ML wlist;
 	struct oVL vl1,vl0;
 	VL vl;
-	int z;
+	int z,dt,dtok;
 
 	/* vl = [x,y] */
 	vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;
@@ -982,7 +989,7 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp)
 	win = W_ALLOC(np+1);
 	wlist = W_MLALLOC(np);
 	wlist->n = list->n;
-	wlist->bound = dy;
+	bound = wlist->bound = list->bound;
 	c = (BM *)COEF(wlist);
 	bcopy((char *)COEF(list),(char *)c,(int)(sizeof(BM)*np));
 
@@ -1006,10 +1013,34 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp)
 	NEWP(lcyx); VR(lcyx) = x; DC(lcyx) = dc;
 	ptosfbm(dy,lcyx,lcy);
 
+	/* initialize lcy0 by LC(f) */
+	lcy0 = W_UMALLOC(bound);
+	ptosfum(COEF(DC(g)),lcy0);
+
+	/* ((d-1 coefs)*lcy0 */
+	d1c = (UM *)W_ALLOC(np*sizeof(UM));
+	w = W_UMALLOC(2*bound);
+	for ( i = 1; i < np; i++ ) {
+		extractcoefbm(c[i],degbm(c[i])-1,w);
+		d1c[i] = W_UMALLOC(2*bound);
+		mulsfum(w,lcy0,d1c[i]);
+		/* d1c[i] = d1c[i] mod y^(bound+1) */
+		if ( DEG(d1c[i]) > bound ) {
+			for ( j = DEG(d1c[i]); j > bound; j-- )
+				COEF(d1c[i])[j] = 0;
+			degum(d1c[i],bound);
+		}
+	}
+
 	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 ( sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,k,win,&factor,&cofactor) ) {
+		dt = sfdegtest(dy,bound,d1c,k,win);
+		if ( dt )
+			dtok++;
+		if ( dt && sfdtestmain(vl,lcg,lcg0,lcy,csum,wlist,
+				k,win,&factor,&cofactor) ) {
 			NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
 			g = cofactor;
 
@@ -1022,6 +1053,9 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp)
 			/* update csum */
 			sfcsump(vl,lcg,&csum);
 
+			/* update dy */
+			dy = getdeg(y,g);
+
 			/* update lcy */
 			clearbm(0,lcy);
 			COEF(dc) = COEF(DC(g));
@@ -1043,6 +1077,22 @@ void sfdtest(P f,ML list,V x,V y,DCP *dcp)
 			else 
 				for ( i = 1; i < k; i++ ) 
 					win[i] = win[0] + i;
+
+
+			/* update lcy0  */
+			ptosfum(COEF(DC(g)),lcy0);
+		
+			/* update d-1 coeffs */
+			for ( i = 1; i <= np; i++ ) {
+				extractcoefbm(c[i],degbm(c[i])-1,w);
+				mulsfum(w,lcy0,d1c[i]);
+				/* d1c[i] = d1c[1] mod y^(bound+1) */
+				if ( DEG(d1c[i]) > bound ) {
+					for ( j = DEG(d1c[i]); j > bound; j-- )
+						COEF(d1c[i])[j] = 0;
+					degum(d1c[i],bound);
+				}
+			}
 		} else if ( !ncombi(1,np,k,win) ) 
 			if ( k == np ) 
 				break;
@@ -1050,12 +1100,45 @@ 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,"\n");
+	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;
 }
 
+void extractcoefbm(BM f,int dx,UM r)
+{
+	int j;
+	UM fj;
+
+	for ( j = DEG(f); j >= 0; j-- ) {
+		fj = COEF(f)[j];
+		if ( fj && DEG(fj) >= dx ) {
+			COEF(r)[j] = COEF(fj)[dx];
+		} else
+			COEF(r)[j] = 0;
+	}
+	degum(r,DEG(f));
+}
+
+/* deg_y(prod mod y^(bound+1)) <= dy ? */
+
+int sfdegtest(int dy,int bound,UM *d1c,int k,int *in)
+{
+	int i,j;
+	UM w,w1,wt;
+	BM t;
+
+	w = W_UMALLOC(bound);
+	w1 = W_UMALLOC(bound);
+	clearum(w,bound);
+	for ( i = 0; i < k; i++ ) {
+		addsfum(w,d1c[in[i]],w1); wt = w; w = w1; w1 = wt;
+	}
+	return DEG(w) <= dy ? 1 : 0;
+}
+
 /* lcy = LC(g), lcg = lcy*g, lcg0 = const part of lcg */
+
 int sfdtestmain(VL vl,P lcg,UM lcg0,BM lcy,P csum,ML list,
 	int k,int *in,P *fp,P *cofp)
 {