===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Hgfs.c,v
retrieving revision 1.18
retrieving revision 1.23
diff -u -p -r1.18 -r1.23
--- OpenXM_contrib2/asir2000/engine/Hgfs.c	2001/10/09 01:36:10	1.18
+++ OpenXM_contrib2/asir2000/engine/Hgfs.c	2002/09/27 08:40:48	1.23
@@ -1,8 +1,10 @@
-/* $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.22 2002/09/27 04:24:04 noro Exp $ */
 
 #include "ca.h"
+#include "inline.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)
 {
@@ -26,7 +28,9 @@ void fctrsf(P p,DCP *dcp)
 
 	simp_ff((Obj)p,&obj); p = (P)obj;
 	if ( !p ) {
-		*dcp = 0; return;
+		NEWDC(dc); COEF(dc) = 0; DEG(dc) = ONE; 
+		NEXT(dc) = 0; *dcp = dc;
+		return;
 	}
 	mp = W_UMALLOC(UDEG(p));
 	ptosfum(p,mp);
@@ -509,7 +513,7 @@ void ptosfbm(int,P,BM);
 
 /* f = f(x,y) */
 
-void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *listp)
+void sfhensel(int count,P f,V x,int degbound,GFS *evp,P *sfp,ML *listp)
 {
 	int i;
 	int fn;
@@ -517,10 +521,10 @@ 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;
+	DCP dc,dct,dc0;
 	UM q,fm,hm;
 	UM *gm;
 	struct oEGT tmp0,tmp1,eg_hensel,eg_hensel_t;
@@ -543,9 +547,35 @@ void sfhensel(int count,P f,V x,GFS *evp,P *sfp,ML *li
 		*listp = rlist = MLALLOC(1); rlist->n = fn; rlist->c[0] = 0;
 		return;
 	}
-	/* pass the the leading coeff. to the first element */
-	c = dc->c; dc = NEXT(dc);
-	mulp(vl,dc->c,c,&t); dc->c = t;
+	if ( degbound >= 0 ) {
+		/*
+		 * reconstruct dc so that
+		 * dc[1],... : factors satisfy degree bound
+		 * dc[0]     : product of others
+		 */
+		c = dc->c; dc = NEXT(dc);
+		dc0 = 0;
+		fn = 0;
+		while ( dc ) {
+			if ( getdeg(x,COEF(dc)) <= degbound ) {
+				dct = NEXT(dc); NEXT(dc) = dc0; dc0 = dc; dc = dct;
+				fn++;
+			} else {
+				mulp(vl,COEF(dc),c,&t); c = t;
+				dc = NEXT(dc);
+			}
+		}
+		if ( OID(c) == O_P ) {
+			NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = dc0;
+			fn++;
+		} else {
+			mulp(vl,dc0->c,c,&t); dc0->c = t; dc = dc0;
+		}
+	} else {
+		/* pass the the leading coeff. to the first element */
+		c = dc->c; dc = NEXT(dc);
+		mulp(vl,dc->c,c,&t); dc->c = t;
+	}
 
 	/* convert mod y-a factors into UM */
 	gm = (UM *)ALLOCA(fn*sizeof(UM));
@@ -554,9 +584,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 +603,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 +613,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);
@@ -615,7 +650,7 @@ int sfberle(VL vl,P f,int count,GFS *ev,DCP *dcp)
 	V x,y;
 	P lc,lc0,f0;
 	Obj obj;
-	int j,q1,index,i;
+	int j,q,index,i;
 
 	clctv(vl,f,&nvl); vl = nvl;
 	x = vl->v; y = vl->next->v;
@@ -623,11 +658,11 @@ int sfberle(VL vl,P f,int count,GFS *ev,DCP *dcp)
 	n = QTOS(DEG(DC(f)));
 	wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); wf2 = W_UMALLOC(n);
 	wfs = W_UMALLOC(n); gcd = W_UMALLOC(n);
-	q1 = field_order_sf()-1;
+	q = field_order_sf();
 	lc = DC(f)->c;
 	for ( j = 0, fn = n + 1, index = 0; 
-		index < q1 && j < count && fn > 1; index++ ) {
-		MKGFS(index,m);
+		index < q && j < count && fn > 1; index++ ) {
+		indextogfs(index,&m);
 		substp(vl,lc,y,(P)m,&lc0);
 		if ( lc0 ) {
 			substp(vl,f,y,(P)m,&f0);
@@ -643,7 +678,7 @@ int sfberle(VL vl,P f,int count,GFS *ev,DCP *dcp)
 			}
 		}
 	}
-	if ( index == q1 )
+	if ( index == q )
 		return 0;
 	else if ( fn == 1 )
 		return 1;
@@ -830,7 +865,7 @@ void sfbmtop(BM f,V x,V y,P *fp)
 			if ( DEG(c[j]) >= i && (a = COEF(c[j])[i]) ) {
 				NEWDC(dct);
 				STOQ(j,DEG(dct));
-				MKGFS(IFTOF(a),b);
+				iftogfs(a,&b);
 				COEF(dct) = (P)b;
 				NEXT(dct) = dc;
 				dc = dct;
@@ -863,8 +898,10 @@ void sfsqfr(P f,DCP *dcp)
 		NEWDC(dc); DEG(dc) = ONE; COEF(dc) = f; NEXT(dc) = 0; *dcp = dc;
 	} else if ( !NEXT(vl) )
 		sfusqfr(f,dcp);
+#if 0
 	else if ( !NEXT(NEXT(vl)) )
 		sfbsqfr(f,vl->v,NEXT(vl)->v,dcp);
+#endif
 	else
 		error("sfsqfr : not implemented yet");
 }
@@ -897,16 +934,21 @@ void sfusqfr(P f,DCP *dcp)
 	*dcp = dct;
 }
 
+void sfbsqfrmain(P f,V x,V y,DCP *dcp)
+{
+	/* XXX*/
+}
+
+/* f is bivariate */
+
 void sfbsqfr(P f,V x,V y,DCP *dcp)
 {
 	P t,rf,cx,cy;
 	VL vl,rvl;
-	DCP dcx,dcy;
+	DCP dcx,dcy,dct,dc;
 	struct oVL vl0,vl1;
 
-	/* vl = [x,y] */
-	vl0.v = x; vl0.next = &vl1; vl1.v = y; vl1.next = 0; vl = &vl0;
-	/* cy(y) = cont(f,x), f /= cx */
+	/* cy(y) = cont(f,x), f /= cy */
 	cont_pp_sfp(vl,f,&cy,&t); f = t;
 	/* rvl = [y,x] */
 	reordvar(vl,y,&rvl); reorderp(rvl,vl,f,&rf);
@@ -915,13 +957,32 @@ void sfbsqfr(P f,V x,V y,DCP *dcp)
 	reorderp(vl,rvl,rf,&f);
 
 	/* f -> cx*cy*f */
-	sfusqfr(cx,&dcx);
-	sfusqfr(cy,&dcy);
+	sfsqfr(cx,&dcx); dcx = NEXT(dcx);
+	sfsqfr(cy,&dcy); dcy = NEXT(dcy);
+	if ( dcx ) {
+		for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
+		NEXT(dct) = dcy;
+	} else
+		dcx = dcy;
+	if ( OID(f) == O_N )
+		*dcp = dcx;
+	else {
+		/* f must be bivariate */
+		sfbsqfrmain(f,x,y,&dc);
+		if ( dcx ) {
+			for ( dct = dcx; NEXT(dct); dct = NEXT(dct) );
+			NEXT(dct) = dc;
+		} else
+			dcx = dc;
+		*dcp = dcx;
+	}
 }
 
 void sfdtest(P,ML,V,V,DCP *);
 
-void sfbfctr(P f,V x,V y,DCP *dcp)
+/* if degbound >= 0 find factor s.t. deg_x(factor) <= degbound */
+
+void sfbfctr(P f,V x,V y,int degbound,DCP *dcp)
 {
 	ML list;
 	P sf;
@@ -931,7 +992,7 @@ void sfbfctr(P f,V x,V y,DCP *dcp)
 	int dx,dy;
 
 	/* sf(x) = f(x+ev) = list->c[0]*list->c[1]*... */
-	sfhensel(5,f,x,&ev,&sf,&list);
+	sfhensel(5,f,x,degbound,&ev,&sf,&list);
 	if ( list->n == 0 )
 		error("sfbfctr : short of evaluation points");
 	else if ( list->n == 1 ) {
@@ -959,18 +1020,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 +1044,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 +1068,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 +1108,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 +1132,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 +1155,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)
 {
@@ -1184,7 +1322,7 @@ void cont_pp_sfp(VL vl,P f,P *cp,P *fp)
 	y = vl->next->v;
 	d = getdeg(y,f);
 	if ( d == 0 ) {
-		MKGFS(0,g);
+		itogfs(1,&g);
 		*cp = (P)g;
 		*fp = f;  /* XXX */
 	} else {
@@ -1250,6 +1388,7 @@ int divtp_by_sfbm(VL vl,P f,P g,P *qp)
 /* XXX generate an irreducible poly of degree n */
 
 extern int current_gfs_q1;
+extern int *current_gfs_ntoi;
 
 void generate_defpoly_sfum(int n,UM *dp)
 {
@@ -1277,8 +1416,12 @@ void generate_defpoly_sfum(int n,UM *dp)
 		for ( j = 0; j < i; j++ )
 			w[j] = 0;
 		w[i]++;
-		for ( i = 0; i < n; i++ )
-			c[i] = w[i]?FTOIF(w[i]-1):0;
+		if ( !current_gfs_ntoi )
+			for ( i = 0; i < n; i++ )
+				c[i] = w[i]?FTOIF(w[i]):0;
+		else
+			for ( i = 0; i < n; i++ )
+				c[i] = w[i]?FTOIF(w[i]-1):0;
 		if ( !c[0] )
 			continue;
 		diffsfum(r,dr); cpyum(r,t); gcdsfum(t,dr,g);