===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/gf.c,v
retrieving revision 1.1.1.1
retrieving revision 1.15
diff -u -p -r1.1.1.1 -r1.15
--- OpenXM_contrib2/asir2000/builtin/gf.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/gf.c	2002/03/15 02:52:09	1.15
@@ -1,4 +1,52 @@
-/* $OpenXM: OpenXM/src/asir99/builtin/gf.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
+/*
+ * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
+ * All rights reserved.
+ * 
+ * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
+ * non-exclusive and royalty-free license to use, copy, modify and
+ * redistribute, solely for non-commercial and non-profit purposes, the
+ * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
+ * conditions of this Agreement. For the avoidance of doubt, you acquire
+ * only a limited right to use the SOFTWARE hereunder, and FLL or any
+ * third party developer retains all rights, including but not limited to
+ * copyrights, in and to the SOFTWARE.
+ * 
+ * (1) FLL does not grant you a license in any way for commercial
+ * purposes. You may use the SOFTWARE only for non-commercial and
+ * non-profit purposes only, such as academic, research and internal
+ * business use.
+ * (2) The SOFTWARE is protected by the Copyright Law of Japan and
+ * international copyright treaties. If you make copies of the SOFTWARE,
+ * with or without modification, as permitted hereunder, you shall affix
+ * to all such copies of the SOFTWARE the above copyright notice.
+ * (3) An explicit reference to this SOFTWARE and its copyright owner
+ * shall be made on your publication or presentation in any form of the
+ * results obtained by use of the SOFTWARE.
+ * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
+ * for such modification or the source code of the modified part of the
+ * SOFTWARE.
+ * 
+ * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
+ * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
+ * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
+ * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
+ * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
+ * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
+ * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
+ * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
+ * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
+ * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
+ * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
+ * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
+ * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
+ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
+ * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
+ *
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/gf.c,v 1.14 2001/10/09 01:36:06 noro Exp $ 
+*/
 #include "ca.h"
 #include "parse.h"
 
@@ -8,6 +56,7 @@ struct resf_dlist {
 
 int resf_degtest(int,int *,int,struct resf_dlist *);
 void uhensel(P,NODE,int,int,NODE *);
+void uhensel_incremental(P,NODE,int,int,int,NODE *);
 void resf_hensel(int,P,int,P *,ML *);
 void resf_dtest(P,ML,int,int *,int *,DCP *);
 void resf_dtest_special(P,ML,int,int *,int *,DCP *);
@@ -18,6 +67,8 @@ void nullspace(UM **,UM,int,int,int *);
 void nullspace_lm(LM **,int,int *);
 void nullspace_gf2n(GF2N **,int,int *);
 void nullspace_gfpn(GFPN **,int,int *);
+void nullspace_gfs(GFS **,int,int *);
+void nullspace_gfsn(GFSN **,int,int *);
 void null_to_sol(int **,int *,int,int,UM *);
 
 void showgfmat(UM **,int);
@@ -25,6 +76,8 @@ void pwr_mod(P,P,V,P,int,N,P *);
 void rem_mod(P,P,V,P,int,P *);
 
 void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
+void Puhensel_incremental();
+void Psfuhensel();
 
 void Pnullspace_ff();
 
@@ -34,6 +87,7 @@ void Plinear_form_to_vect(),Pvect_to_linear_form();
 void solve_linear_equation_gf2n(GF2N **,int,int,int *);
 void linear_form_to_array(P,VL,int,Num *);
 void array_to_linear_form(Num *,VL,int,P *);
+void sfuhensel(P,NODE,GFS,int,NODE *);
 
 extern int current_ff;
 
@@ -46,6 +100,8 @@ struct ftab gf_tab[] = {
 	{"resfmain",Presfmain,4},
 	{"pwr_mod",Ppwr_mod,6},
 	{"uhensel",Puhensel,4},
+	{"uhensel_incremental",Puhensel_incremental,5},
+	{"sfuhensel",Psfuhensel,4},
 	{0,0,0},
 };
 
@@ -98,6 +154,23 @@ LIST *rp;
 	MKLIST(*rp,r);
 }
 
+void Puhensel_incremental(arg,rp)
+NODE arg;
+LIST *rp;
+{
+	P f;
+	NODE mfl,r;
+	int mod,bound,start;
+
+	f = (P)ARG0(arg);
+	mfl = BDY((LIST)ARG1(arg));
+	mod = QTOS((Q)ARG2(arg));
+	start = QTOS((Q)ARG3(arg));
+	bound = QTOS((Q)ARG4(arg));
+	uhensel_incremental(f,mfl,mod,start,bound,&r);
+	MKLIST(*rp,r);
+}
+
 void uhensel(f,mfl,mod,bound,rp)
 P f;
 NODE mfl;
@@ -130,6 +203,141 @@ NODE *rp;
 	*rp = top;
 }
 
+void uhensel_incremental(f,mfl,mod,start,bound,rp)
+P f;
+NODE mfl;
+int mod,start,bound;
+NODE *rp;
+{
+	ML blist,clist,rlist;
+	LUM fl;
+	LUM *lblist;
+	int nf,i,j,k;
+	int **p;
+	P s;
+	V v;
+	NODE t,top;
+
+	nf = length(mfl);
+	blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
+	lblist = (LUM *)MALLOC(nf*sizeof(LUM));
+	for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
+		blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
+		ptoum(mod,(P)BDY(t),blist->c[i]);
+		W_LUMALLOC((int)UDEG((P)BDY(t)),bound,lblist[i]);
+		ptolum(mod,start,(P)BDY(t),lblist[i]);
+		p = lblist[i]->c;
+		for ( j = DEG(lblist[i]); j >= 0; j-- )
+			for ( k = start; k < bound; k++ )
+				p[j][k] = 0;
+	}
+	gcdgen(f,blist,&clist);	
+	clist->bound = bound;
+	W_LUMALLOC((int)UDEG(f),bound,fl);
+	ptolum(mod,bound,f,fl);
+	henmain_incremental(fl,lblist,clist,nf,mod,start,bound);
+	v = VR(f);
+	for ( i = nf-1, top = 0; i >= 0; i-- ) {
+		lumtop_unsigned(v,mod,bound,lblist[i],&s);
+		MKNODE(t,s,top); top = t;
+	}
+	*rp = top;
+}
+
+void Psfuhensel(arg,rp)
+NODE arg;
+LIST *rp;
+{
+	P f;
+	int bound;
+	NODE r,mfl;
+	GFS ev;
+
+	f = (P)ARG0(arg);
+	mfl = BDY((LIST)ARG1(arg));
+	ev = (GFS)ARG2(arg);
+	bound = QTOS((Q)ARG3(arg));
+	sfuhensel(f,mfl,ev,bound,&r);
+	MKLIST(*rp,r);
+}
+
+void sfuhensel(f,mfl,ev,bound,rp)
+P f;
+NODE mfl;
+GFS ev;
+int bound;
+NODE *rp;
+{
+	BM fl;
+	BM *r;
+	VL vl,nvl;
+	int i,fn,dx,dy,d;
+	NODE t,top;
+	UM fm,hm,q;
+	UM *gm;
+	V x,y;
+	P g,s,u;
+
+	clctv(CO,f,&vl);
+	if ( !vl || !vl->next || vl->next->next )
+		error("sfuhensel : f must be a bivariate poly");
+
+	for ( i = 0, t = mfl; t; i++, t = NEXT(t) );
+	fn = i;
+
+	gm = (UM *)MALLOC(fn*sizeof(UM));
+
+	/* XXX : more severe check is necessary */
+	x = VR((P)BDY(mfl));
+	y = vl->v == x ? vl->next->v : vl->v;
+
+	for ( i = 0, t = mfl, d = 0; i < fn; i++, t = NEXT(t) ) {
+		gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));
+		ptosfum((P)BDY(t),gm[i]);
+		d += DEG(gm[i]);
+	}
+
+	/* reorder f if necessary */
+	if ( vl->v != x ) {
+		reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);
+		vl = nvl; f = g;
+	}
+	dx = getdeg(x,f);
+	if ( dx != d )
+		error("sfuhensel : product of factors has incompatible degree");
+
+	dy = getdeg(y,f);
+	dy = MAX(dy,bound);
+	fl = BMALLOC(dx,dy);
+	ptosfbm(dy,f,fl);
+	if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
+
+	/* fm = fl mod y */
+	fm = W_UMALLOC(dx);
+	cpyum(COEF(fl)[0],fm);
+	hm = W_UMALLOC(dx);
+
+	q = W_UMALLOC(dx);
+	r = (BM *)MLALLOC(fn*sizeof(BM));
+	for ( i = 0; i < fn-1; i++ ) {
+		/* fl = gm[i]*hm mod y */
+		divsfum(fm,gm[i],hm);
+		/* fl is replaced by the cofactor of gk mod y^bound */
+		/* r[i] = gk */
+		sfhenmain2(fl,gm[i],hm,bound,r+i);
+		cpyum(hm,fm);
+	}
+	/* finally, fl must be the lift of gm[fn-1] */
+	r[i] = fl;
+
+	for ( i = fn-1, top = 0; i >= 0; i-- ) {
+		sfbmtop(r[i],x,y,&s);
+		reorderp(CO,vl,s,&u);
+		MKNODE(t,u,top); top = t;
+	}
+	*rp = top;
+}
+
 void Presfmain(arg,rp)
 NODE arg;
 LIST *rp;
@@ -180,19 +388,42 @@ P *mfl;
 ML *listp;
 {
 	register int i,j;
-	int q,n,bound;
+	int q,n,bound,inv,lc;
 	int *p;
 	int **pp;
 	ML blist,clist,bqlist,cqlist,rlist;
 	UM *b;
 	LUM fl,tl;
 	LUM *l;
+	UM w;
 
+	w = W_UMALLOC(UDEG(f));
 	blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
-	for ( i = 0; i < nf; i++ ) {
+
+	/* c[0] must have lc(f) */
+	blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
+	ptoum(mod,mfl[0],w);
+	inv = invm(w->c[UDEG(mfl[0])],mod);
+	lc = rem(NM((Q)LC(f)),mod);
+	if ( SGN((Q)LC(f)) < 0 )
+		lc = (mod-lc)%mod;
+	lc = dmar(inv,lc,0,mod);
+	if ( lc == 1 )
+		copyum(w,blist->c[0]);
+	else
+		mulsum(mod,w,lc,blist->c[0]);
+
+	/* c[i] (i=1,...,nf-1) must be monic */
+	for ( i = 1; i < nf; i++ ) {
 		blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
-		ptoum(mod,mfl[i],blist->c[i]);
+		ptoum(mod,mfl[i],w);
+		inv = invm(w->c[UDEG(mfl[i])],mod);
+		if ( inv == 1 )
+			copyum(w,blist->c[i]);
+		else
+			mulsum(mod,w,inv,blist->c[i]);
 	}
+
 	gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
 	n = bqlist->n; q = bqlist->mod;
 	bqlist->bound = cqlist->bound = bound = mignotte(q,f);
@@ -714,7 +945,6 @@ NODE arg;
 LIST *rp;
 {
 	int i,j,n;
-	Q mod;
 	MAT mat,r;
 	VECT u;
 	Q q;
@@ -738,6 +968,10 @@ LIST *rp;
 			nullspace_gf2n((GF2N **)w,n,ind); break;
 		case FF_GFPN:
 			nullspace_gfpn((GFPN **)w,n,ind); break;
+		case FF_GFS:
+			nullspace_gfs((GFS **)w,n,ind); break;
+		case FF_GFSN:
+			nullspace_gfsn((GFSN **)w,n,ind); break;
 		default:
 			error("nullspace_ff : current_ff is not set");
 	}
@@ -876,6 +1110,89 @@ int *ind;
 		}
 	}
 }
+
+void nullspace_gfs(mat,n,ind)
+GFS **mat;
+int n;
+int *ind;
+{
+	int i,j,l,s;
+	GFS w,w1,h,inv;
+	GFS *t,*u;
+	GFS one;
+
+	bzero(ind,n*sizeof(int));
+	ind[0] = 0;
+	mqtogfs(ONEM,&one);
+
+	for ( i = j = 0; j < n; i++, j++ ) {
+		for ( ; j < n; j++ ) {
+			for ( l = i; l < n; l++ )
+				if ( mat[l][j] )
+					break;
+			if ( l < n ) {
+				t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
+			} else
+				ind[j] = 1;
+		}
+		if ( j == n )
+			break;
+		divgfs(one,mat[i][j],&inv);
+		for ( s = j, t = mat[i]; s < n; s++ ) {
+			mulgfs(t[s],inv,&w); t[s] = w;
+		}
+		for ( l = 0; l < n; l++ ) {
+			if ( l == i )
+				continue;
+			u = mat[l]; 
+			chsgngfs(u[j],&h);
+			for ( s = j; s < n; s++ ) {
+				mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
+			}
+		}
+	}
+}
+
+void nullspace_gfsn(mat,n,ind)
+GFSN **mat;
+int n;
+int *ind;
+{
+	int i,j,l,s;
+	GFSN w,w1,h,inv;
+	GFSN *t,*u;
+
+	bzero(ind,n*sizeof(int));
+	ind[0] = 0;
+
+	for ( i = j = 0; j < n; i++, j++ ) {
+		for ( ; j < n; j++ ) {
+			for ( l = i; l < n; l++ )
+				if ( mat[l][j] )
+					break;
+			if ( l < n ) {
+				t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
+			} else
+				ind[j] = 1;
+		}
+		if ( j == n )
+			break;
+		invgfsn(mat[i][j],&inv);
+		for ( s = j, t = mat[i]; s < n; s++ ) {
+			mulgfsn(t[s],inv,&w); t[s] = w;
+		}
+		for ( l = 0; l < n; l++ ) {
+			if ( l == i )
+				continue;
+			u = mat[l]; 
+			chsgngfsn(u[j],&h);
+			for ( s = j; s < n; s++ ) {
+				mulgfsn(h,t[s],&w); addgfsn(w,u[s],&w1); u[s] = w1;
+			}
+		}
+	}
+}
+
 /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
 
 void linear_form_to_array(p,vl,m,array)