===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/H.c,v
retrieving revision 1.2
retrieving revision 1.9
diff -u -p -r1.2 -r1.9
--- OpenXM_contrib2/asir2000/engine/H.c	2000/08/21 08:31:25	1.2
+++ OpenXM_contrib2/asir2000/engine/H.c	2015/08/06 10:01:52	1.9
@@ -23,7 +23,7 @@
  * 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@flab.fujitsu.co.jp of the detailed specification
+ * 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.
  * 
@@ -45,40 +45,13 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/engine/H.c,v 1.1.1.1 1999/12/03 07:39:08 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/H.c,v 1.8 2013/11/29 08:21:29 noro Exp $ 
 */
 #include "ca.h"
 #include "inline.h"
 #include "base.h"
 #include <math.h>
 
-void modfctrp(P,int,int,DCP *);
-void gensqfrum(int,UM,struct oDUM *);
-void srchum(int,UM,UM,UM);
-UM *resberle(int,UM,UM *);
-int substum(int,UM,int);
-void ddd(int,UM,UM *);
-void canzas(int,UM,int,UM *,UM *);
-int divum(int,UM,UM,UM);
-void randum(int,int,UM);
-void pwrmodum(int,UM,int,UM,UM);
-void spwrum(int,UM,UM *,UM,N,UM);
-void spwrum0(int,UM,UM,N,UM);
-void mulum(int,UM,UM,UM);
-void mulsum(int,UM,int,UM);
-void gcdum(int,UM,UM,UM);
-void mult_mod_tab(UM,int,UM *,UM,int);
-void make_qmat(UM,int,UM *,int ***);
-void null_mod(int **,int,int,int *);
-void null_to_sol(int **,int *,int,int,UM *);
-void newddd(int,UM,UM *);
-int irred_check(UM,int);
-int berlekamp(UM,int,int,UM *,UM *);
-int nfctr_mod(UM,int);
-void minipoly_mod(int,UM,UM,UM);
-int find_root(int,UM,int *);
-void showum(UM);
-void showumat(int **,int);
 #if 1
 #define Mulum mulum
 #define Divum divum  
@@ -93,18 +66,7 @@ void showumat(int **,int);
 			
 LUM LUMALLOC();
 
-struct p_pair {
-	UM p0;
-	UM p1;
-	struct p_pair *next;
-};
-
-void lnf_mod(int,int,UM,UM,struct p_pair *,UM,UM);
-
-void berle(index,count,f,listp)
-int index,count;
-P f;
-ML *listp;
+void berle(int index,int count,P f,ML *listp)
 {
 	UM wf,wf1,wf2,wfs,gcd;
 	ML flist;
@@ -146,9 +108,7 @@ ML *listp;
 		ptr[i] = ( ptr[i] * hd ) % m;
 }
 
-int berlecnt(mod,f)
-register int mod;
-UM f;
+int berlecnt(int mod,UM f)
 {
 	register int i,j,**c;
 	int d,dr,n;
@@ -176,10 +136,7 @@ UM f;
 
 /* XXX berlecntmain should not be used for large mod */
 
-int berlecntmain(mod,n,m,c)
-register int mod;
-int n,m;
-register int **c;
+int berlecntmain(int mod,int n,int m,int **c)
 {
 	register int *p1,*p2,i,j,k,l,a;
 	int *tmp,inv;
@@ -205,10 +162,7 @@ register int **c;
 	return ( cfs );
 }
 
-UM *berlemain(mod,f,fp)
-register int mod;
-UM f;
-UM *fp;
+UM *berlemain(int mod,UM f,UM *fp)
 {
 	UM wg,ws,wf,f0,gcd,q;
 	int n;
@@ -237,10 +191,7 @@ UM *fp;
 	return ( fp );
 }
 
-void hensel(index,count,f,listp)
-int index,count;
-P f;
-ML *listp;
+void hensel(int index,int count,P f,ML *listp)
 {
 	register int i,j;
 	int q,n,bound;
@@ -277,11 +228,248 @@ ML *listp;
 	}
 }
 
-void hsq(index,count,f,nindex,dcp)
-int index,count,*nindex;
-P f;
-DCP *dcp;
+void hensel2(int index,int count,P f,ML *listp)
 {
+	register int i,j;
+	int mod,q,n,bound,dx;
+	ML blist,clist,bqlist,cqlist,rlist;
+	UM fm,qfm,gm,qgm,hm,qhm,qam,qbm,w;
+	UM *b;
+	LUM fl,tl;
+	int k;
+
+	dx = UDEG(f);
+	if ( dx == 1 ) {
+		*listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
+		return;
+	}
+	berle(index,count,f,&blist);
+	n = blist->n;
+	mod = blist->mod;
+
+	if ( n == 1 ) {
+		*listp = blist = MLALLOC(1); blist->n = 1; blist->c[0] = 0;
+		return;
+	}
+	
+	/* find k s.t. mod^k <= 2^27 < mod^(k+1); set q = mod^k */
+	for ( k = 1, q = mod; q <= ((1<<27)/mod); q *= mod, k++ );
+
+	/* mignotte bound */
+	bound = mignotte(q,f);
+
+	*listp = rlist = MLALLOC(n); 
+	rlist->n = n;
+	rlist->mod = q;
+	rlist->bound = bound;
+
+	if ( bound == 1 ) {
+		gcdgen(f,blist,&clist);
+		henprep(f,blist,clist,&bqlist,&cqlist);
+
+		for ( i = 0, b = (UM *)bqlist->c; i < n; i++ ) {
+			COEF(rlist)[i] = tl = LUMALLOC(DEG(b[i]),1);
+			for ( j = 0; j <= DEG(tl); j++ )
+					COEF(tl)[j][0] = COEF(b[i])[j];
+			COEF(rlist)[i] = tl;
+		}
+	} else {
+		/* fl = f mod q */
+		fl = LUMALLOC(dx,bound);
+		ptolum(q,bound,f,fl);
+		/* fm = f mod mod */
+		fm = W_UMALLOC(dx);
+		ptoum(mod,f,fm);
+		/* fm = f mod q */
+		qfm = W_UMALLOC(dx);
+		ptoum(q,f,qfm);
+
+		gm = W_UMALLOC(dx); qgm = W_UMALLOC(dx);
+		hm = W_UMALLOC(dx); qhm = W_UMALLOC(dx);
+		qam = W_UMALLOC(dx); qbm = W_UMALLOC(dx);
+		w = W_UMALLOC(dx);
+		for ( i = 0; i < n-1; i++ ) {
+			cpyum(COEF(blist)[i],gm);
+			cpyum(fm,w);
+			divum(mod,w,gm,hm);
+
+			/* find am,bm s.t. qam*qgm+qbm*qhm=1 mod q, qgm=gm mod mod, qhm=hm mod mod */
+			henprep2(mod,q,k,qfm,gm,hm,qgm,qhm,qam,qbm);
+
+			henmain2(fl,qgm,qhm,qam,qbm,q,bound,&tl);
+			rlist->c[i] = (pointer)tl;
+			cpyum(hm,fm);
+			cpyum(qhm,qfm);
+		}
+		rlist->c[i] = fl;
+	}
+}
+
+/*
+	f = g0*h0 mod m -> f = gk*hk mod m^(bound), f is replaced by hk
+*/
+
+void henmain2(LUM f,UM g0,UM h0,UM a0,UM b0,int m,int bound,LUM *gp)
+{
+	int n,dg,dh,i,k,j,dg1,dh1;
+	UM wu,wr,ws,wt,q,wh1,wg1,wc,wd,we,wz;
+	LUM wb0,wb1,wb2,fk,gk,hk;
+
+	n = DEG(f); dg = DEG(g0); dh = DEG(h0);
+
+	W_LUMALLOC(n,bound,wb0);
+	W_LUMALLOC(n,bound,wb1);
+	W_LUMALLOC(n,bound,wb2);
+
+	wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n);
+	wr = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
+	q = W_UMALLOC(2*n);
+	wh1 = W_UMALLOC(2*n); wg1 = W_UMALLOC(2*n);
+	
+	/* gk = g0 */
+	gk = LUMALLOC(n,bound);
+	DEG(gk) = dg;
+	for ( i = 0; i <= dg; i++ )
+		COEF(gk)[i][0] = COEF(g0)[i];
+
+	/* hk = h0 */
+	W_LUMALLOC(n,bound,hk);
+	DEG(hk) = dh;
+	for ( i = 0; i <= dh; i++ )
+		COEF(hk)[i][0] = COEF(h0)[i];
+
+	/* fk = gk*hk */
+	W_LUMALLOC(n,bound,fk);
+	mullum(m,bound,gk,hk,fk);
+
+	wc = W_UMALLOC(2*n); wd = W_UMALLOC(2*n);
+	we = W_UMALLOC(2*n); wz = W_UMALLOC(2*n);
+
+#if 0
+	mulum(m,a0,g0,wc);
+	mulum(m,b0,h0,wd);
+	addum(m,wc,wd,wz);
+	if ( DEG(wz) != 0 || COEF(wz)[0] != 1 )
+		error("henmain2 : cannot happen(extgcd)");
+#endif
+
+#if 1
+	fprintf(stderr,"bound=%d\n",bound);
+#endif
+	for ( k = 1; k < bound; k++ ) {
+#if 1
+		fprintf(stderr,".");
+#endif
+
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
+		check_intr();
+#endif
+		/* at this point, f = gk*hk mod y^k */
+
+#if 0
+		for ( j = 0; j < k; j++ )
+			for ( i = 0; i <= n; i++ )
+				if ( COEF(f)[i][j] != COEF(f)[i][j] )
+					error("henmain2 : cannot happen(f=fk)");
+#endif
+
+		/* wt = (f-gk*hk)/y^k */
+		for ( i = 0; i <= n; i++ )
+			COEF(ws)[i] = COEF(f)[i][k];
+		degum(ws,n);
+		for ( i = 0; i <= n; i++ )
+			COEF(wu)[i] = COEF(fk)[i][k];
+		degum(wu,n);
+		subum(m,ws,wu,wt);
+		
+		/* compute wf1,wg1 s.t. wh1*g0+wg1*h0 = wt */
+		mulum(m,a0,wt,wh1); DEG(wh1) = divum(m,wh1,h0,q);
+		mulum(m,wh1,g0,wc); subum(m,wt,wc,wd); DEG(wd) = divum(m,wd,h0,wg1);
+
+		/* check */
+#if 0
+		if ( DEG(wd) >= 0 || DEG(wg1) > dg )
+			error("henmain2 : cannot happen(adj)");
+
+		mulum(m,wg1,h0,wc); mulum(m,wh1,g0,wd); addum(m,wc,wd,we);
+		subum(m,we,wt,wz);
+		if ( DEG(wz) >= 0 )
+			error("henmain2 : cannot happen(coef)");
+#endif
+
+		/* fk += ((wg1*hk+wh1*gk)*y^k+wg1*wh1*y^(2*k) mod m^bound */
+
+		/* wb0 = wh1*y^k */
+		clearlum(n,bound,wb0);
+		DEG(wb0) = dh1 = DEG(wh1);
+		for ( i = 0; i <= dh1; i++ )
+			COEF(wb0)[i][k] = COEF(wh1)[i];
+
+		/* wb2 = gk*wb0 mod y^bound */
+		clearlum(n,bound,wb2);
+		mullum(m,bound,gk,wb0,wb2);
+
+		/* fk += wb2 */
+		addtolum(m,bound,wb2,fk);
+
+		/* wb1 = wg1*y^k */
+		clearlum(n,bound,wb1);
+		DEG(wb1) = dg1 = DEG(wg1);
+		for ( i = 0; i <= n; i++ )
+			COEF(wb1)[i][k] = COEF(wg1)[i];
+
+		/* wb2 = hk*wb1 mod y^bound */
+		clearlum(n,bound,wb2);
+		mullum(m,bound,hk,wb1,wb2);
+
+		/* fk += wb2 */
+		addtolum(m,bound,wb2,fk);
+
+		/* fk += wg1*wh1*y^(2*k) mod y^bound) */
+		if ( 2*k < bound ) {
+			clearlum(n,bound,wb2);
+			mullum(m,bound,wb0,wb1,wb2);
+			addtolum(m,bound,wb2,fk);
+		}
+
+		/* gk += wg1*y^k, hk += wh1*y^k */
+		for ( i = 0; i <= DEG(wg1); i++ )
+			COEF(gk)[i][k] = COEF(wg1)[i];
+		for ( i = 0; i <= DEG(wh1); i++ )
+			COEF(hk)[i][k] = COEF(wh1)[i];
+	}
+#if 1
+	fprintf(stderr,"\n");
+#endif
+	*gp = gk;
+	clearlum(n,bound,f);
+	DEG(f) = dh;
+	for ( i = 0; i <= dh; i++ )
+		for ( j = 0; j < bound; j++ )
+			COEF(f)[i][j] = COEF(hk)[i][j];
+}
+
+void clearlum(int n,int bound,LUM f)
+{
+	int i;
+
+	for ( i = 0; i <= n; i++ )
+		bzero(COEF(f)[i],bound*sizeof(int));
+}
+
+/* g += f */
+
+void addtolum(int m,int bound,LUM f,LUM g)
+{
+	int n,i;
+
+	n = DEG(f);
+	for ( i = 0; i <= n; i++ )
+		addpadic(m,bound,COEF(f)[i],COEF(g)[i]);
+}
+
+void hsq(int index,int count,P f,int *nindex,DCP *dcp)
+{
 	register int i,j,k;
 	register int **pp,**fpp;
 	register int *px,*py;
@@ -406,9 +594,7 @@ DCP *dcp;
 	*dcp = 0;
 }
 
-void gcdgen(f,blist,clistp)
-P f;
-ML blist,*clistp;
+void gcdgen(P f,ML blist,ML *clistp)
 {
 	register int i;
 	int n,d,mod,np;
@@ -435,13 +621,57 @@ ML blist,*clistp;
 	}
 }
 
+/* find a,b s.t. qa*qg+qb*qh=1 mod q, qg=g mod mod, qh=h mod mod */
+/* q = mod^k */
+
+void henprep2(int mod,int q,int k,UM f,UM g,UM h,UM qg,UM qh,UM qa,UM qb)
+{
+	int n;
+	UM wg,wh,wa,wb;
+	ML bl,cl,bql,cql;
+	P ff;
+
+	n = DEG(f);
+	wg = W_UMALLOC(2*n); wh = W_UMALLOC(2*n);
+	wa = W_UMALLOC(2*n); wb = W_UMALLOC(2*n);
+	cpyum(g,wg); cpyum(h,wh);
+
+	/* wa*g+wb*h = 1 mod mod */
+	eucum(mod,wg,wh,wa,wb);
+
+#if 0
+	/* check */
+	wt = W_UMALLOC(2*n); ws = W_UMALLOC(2*n); wu = W_UMALLOC(2*n);
+	mulum(mod,wa,g,wt);
+	mulum(mod,wb,h,ws);
+	addum(mod,wt,ws,wu);
+ 	if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
+		error("henprep 1");
+#endif
+
+	bl = MLALLOC(2); bl->n = 2; bl->mod = mod; bl->c[0] = g; bl->c[1] = h;
+	cl = MLALLOC(2); cl->n = 2; cl->mod = mod; cl->c[0] = wb; cl->c[1] = wa;
+	umtop(CO->v,f,&ff); /* XXX */
+	henprep(ff,bl,cl,&bql,&cql); /* XXX */
+
+	cpyum(bql->c[0],qg); cpyum(bql->c[1],qh);
+	cpyum(cql->c[0],qb); cpyum(cql->c[1],qa);
+
+#if 0
+	/* check */
+	mulum(q,qa,qg,wt);
+	mulum(q,qb,qh,ws);
+	addum(q,wt,ws,wu);
+ 	if ( DEG(wu) != 0 || COEF(wu)[0] != 1 )
+		error("henprep 2");
+#endif
+}
+
 /*
-	henprep(f,&blist,&clist,&bqlist,&cqlist);
+	henprep(f,blist,clist,&bqlist,&cqlist);
  */
 
-void henprep(f,blist,clist,bqlistp,cqlistp)
-P f;
-ML blist,clist,*bqlistp,*cqlistp;
+void henprep(P f,ML blist,ML clist,ML *bqlistp,ML *cqlistp)
 {
 	register int i,j,k,*px,*py,*pz;
 	int n,pmax,dr,tmp,p,p1,mod,np,b,q;
@@ -522,9 +752,7 @@ ML blist,clist,*bqlistp,*cqlistp;
 	henmain(fl,bqlist,cqlist,listp)
 */
 
-void henmain(f,bqlist,cqlist,listp)
-LUM f;
-ML bqlist,cqlist,*listp;
+void henmain(LUM f,ML bqlist,ML cqlist,ML *listp)
 {
 	register int i,j,k;
 	int *px,*py;
@@ -549,7 +777,16 @@ ML bqlist,cqlist,*listp;
 		for ( j = DEG(b[i]), pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- ) 
 			pp[j][0] = px[j];
 	}
+#if 0
+	fprintf(stderr,"bound=%d\n",bound);
+#endif
 	for ( i = 1; i < bound; i++ ) {
+#if 0
+		fprintf(stderr,".");
+#endif
+#if defined(VISUAL) || defined(__MINGW32__) || defined(__MINGW64__)
+		check_intr();
+#endif
 		mullum(mod,i+1,l[0],l[1],wb0);
 		for ( j = 2; j < np; j++ ) {
 			mullum(mod,i+1,l[j],wb0,wb1);
@@ -578,14 +815,87 @@ ML bqlist,cqlist,*listp;
 		for ( j = n, px = COEF(wq0); j >= 0; j-- ) 
 			px[j] = 0;
 	}
+#if 0
+	fprintf(stderr,"\n");
+#endif
 }
 
+/*
+	henmain_incremental(fl,bqlist,cqlist,start)
+    fl = bqlist[0]*... mod q^start
+*/
+
+void henmain_incremental(LUM f,LUM *bqlist,ML cqlist,
+	int np, int mod, int start, int bound)
+{
+	register int i,j,k;
+	int *px,*py;
+	int **pp,**pp1;
+	int n,dr,tmp;
+	UM wt,wq0,wq,wr,wm,wm0,wa,q;
+	LUM wb0,wb1,tlum;
+	UM *b,*c;
+	LUM *l;
+	ML list;
+
+	n = DEG(f);
+	W_LUMALLOC(n,bound,wb0); W_LUMALLOC(n,bound,wb1);
+	wt = W_UMALLOC(n); wq0 = W_UMALLOC(n); wq = W_UMALLOC(n);
+	wr = W_UMALLOC(n); wm = W_UMALLOC(2*n); wm0 = W_UMALLOC(2*n);
+	wa = W_UMALLOC(2*n); q = W_UMALLOC(2*n);
+	c = (UM *)cqlist->c; l = bqlist;
+	b = (UM *)ALLOCA(n*sizeof(UM));
+	for ( i = 0; i < np; i++ ) {
+		j = DEG(l[i]);
+		b[i] = W_UMALLOC(j);
+		DEG(b[i]) = j;
+		for ( pp = COEF(l[i]), px = COEF(b[i]); j >= 0; j-- )
+			px[j] = pp[j][0];
+	}
+#if 0
+	fprintf(stderr,"bound=%d\n",bound);
+#endif
+	for ( i = start; i < bound; i++ ) {
+#if 0
+		fprintf(stderr,".");
+#endif
+		mullum(mod,i+1,l[0],l[1],wb0);
+		for ( j = 2; j < np; j++ ) {
+			mullum(mod,i+1,l[j],wb0,wb1);
+			tlum = wb0; wb0 = wb1; wb1 = tlum;
+		}
+		for ( j = n, px = COEF(wt); j >= 0; j-- )
+			px[j] = 0;
+		for ( j = n, pp = COEF(f), pp1 = COEF(wb0); j >= 0; j-- ) {
+			tmp = ( pp[j][i] - pp1[j][i] ) % mod;
+			COEF(wt)[j] = ( tmp < 0 ? tmp + mod : tmp );
+		}
+		degum(wt,n);
+		for ( j = n, px = COEF(wq0); j >= 0; j-- )
+			px[j] = 0;
+		for ( j = 1; j < np; j++ ) {
+			mulum(mod,wt,c[j],wm); dr = divum(mod,wm,b[j],q);
+			for ( k = DEG(q), px = COEF(wq0), py = COEF(q); k >= 0; k-- ) 
+				px[k] = ( px[k] + py[k] ) % mod;
+			for ( k = dr, pp = COEF(l[j]), px = COEF(wm); k >= 0; k-- ) 
+				pp[k][i] = px[k];
+		}
+		degum(wq0,n); mulum(mod,wq0,b[0],wm);
+		mulum(mod,wt,c[0],wm0); addum(mod,wm,wm0,wa);
+		for ( j = DEG(wa), pp = COEF(l[0]), px = COEF(wa); j >= 0; j-- ) 
+			pp[j][i] = px[j];
+		for ( j = n, px = COEF(wq0); j >= 0; j-- ) 
+			px[j] = 0;
+	}
+#if 0
+	fprintf(stderr,"\n");
+#endif
+}
+
 static double M;
 static int E;
 
-int mignotte(q,f)
-int q;
-P f;
+int mignotte(int q,P f)
 {
 	int p;
 	unsigned int *b;
@@ -601,9 +911,7 @@ P f;
 	return (int)ceil( (0.31*(E+UDEG(f)+1)+log10((double)M)) / log10((double)q) );
 }
 
-int mig(q,d,f)
-int q,d;
-P f;
+int mig(int q,int d,P f)
 {
 	int p;
 	unsigned int *b;
@@ -620,9 +928,7 @@ P f;
 	return (int)ceil( (0.31*(E+d+1)+log10((double)M)) / log10((double)q) );
 }
 
-void sqad(man,exp)
-unsigned int man;
-int exp;
+void sqad(unsigned int man,int exp)
 {
 	int e,sqe;
 	unsigned int t;
@@ -658,10 +964,7 @@ int exp;
 	}
 }
 
-void ptolum(q,bound,f,fl)
-int q,bound;
-P f;
-LUM fl;
+void ptolum(int q,int bound,P f,LUM fl)
 {
 	DCP dc;
 	int i,j;
@@ -676,7 +979,7 @@ LUM fl;
 		c = pp[QTOS(DEG(dc))]; w = (unsigned int *)W_ALLOC(d);
 		for ( i = 0; i < d; i++ )
 			w[i] = m[i];
-		for ( i = 0; d >= 1; ) {
+		for ( i = 0; i < bound && d >= 1; ) {
 			for ( j = d - 1, r = 0; j >= 0; j-- ) {
 				DSAB(q,r,w[j],w[j],r)
 			}
@@ -694,10 +997,7 @@ LUM fl;
 	}
 }
 
-void modfctrp(p,mod,flag,dcp)
-P p;
-int mod,flag;
-DCP *dcp;
+void modfctrp(P p,int mod,int flag,DCP *dcp)
 {
 	int cm,n,i,j,k;
 	DCP dc,dc0;
@@ -789,10 +1089,7 @@ DCP *dcp;
 	NEXT(dc) = 0; *dcp = dc0;
 }
 
-void gensqfrum(mod,p,dc)
-int mod;
-UM p;
-struct oDUM *dc;
+void gensqfrum(int mod,UM p,struct oDUM *dc)
 {
 	int n,i,j,d;
 	UM t,s,g,f,f1,b;
@@ -845,9 +1142,7 @@ struct oDUM *dc;
 }
 
 #if 0
-void srchum(mod,p1,p2,gr)
-int mod;
-UM p1,p2,gr;
+void srchum(int mod,UM p1,UM p2,UM gr)
 {
 	UM m,m1,m2,q,r,t,g1,g2;
 	int lc,d,d1,d2,i,j,k,l,l1,l2,l3,tmp,adj;
@@ -917,10 +1212,7 @@ UM p1,p2,gr;
 	}
 }
 
-UM *resberle(mod,f,fp)
-register int mod;
-UM f;
-UM *fp;
+UM *resberle(int mod,UM f,UM *fp)
 {
 	UM w,wg,ws,wf,f0,gcd,q,res;
 	int n;
@@ -953,10 +1245,7 @@ UM *fp;
 	return ( fp );
 }
 
-int substum(mod,p,a)
-int mod;
-UM p;
-int a;
+int substum(int mod,UM p,int a)
 {
 	int i,j,s;
 	int *c;
@@ -976,9 +1265,7 @@ int a;
 }
 #endif
 
-void ddd(mod,f,r)
-int mod;
-UM f,*r;
+void ddd(int mod,UM f,UM *r)
 {
 	register int i,j;
 	int d,n;
@@ -1024,9 +1311,7 @@ UM f,*r;
 }
 
 #if 0
-void canzas(mod,f,d,base,r)
-int mod;
-UM f,*base,*r;
+void canzas(int mod,UM f,int d,UM *base,UM *r)
 {
 	UM t,s,u,w,g,o,q;
 	N n1,n2,n3,n4,n5;
@@ -1064,9 +1349,7 @@ UM f,*base,*r;
 	}
 }
 #else
-void canzas(mod,f,d,base,r)
-int mod;
-UM f,*base,*r;
+void canzas(int mod,UM f,int d,UM *base,UM *r)
 {
 	UM t,s,u,w,g,o,q;
 	N n1,n2,n3,n4,n5;
@@ -1105,9 +1388,7 @@ UM f,*base,*r;
 }
 #endif
 
-void randum(mod,d,p)
-int mod,d;
-UM p;
+void randum(int mod,int d,UM p)
 {
 	unsigned int n;
 	int i;
@@ -1117,9 +1398,7 @@ UM p;
 		COEF(p)[i] = ((unsigned int)random()) % mod;
 }
 
-void pwrmodum(mod,p,e,f,pr)
-int mod,e;
-UM p,f,pr;
+void pwrmodum(int mod,UM p,int e,UM f,UM pr)
 {
 	UM wt,ws,q;
 
@@ -1145,11 +1424,7 @@ UM p,f,pr;
 	}
 }
 
-void spwrum(mod,m,base,f,e,r)
-int mod;
-UM f,m,r;
-UM *base;
-N e;
+void spwrum(int mod,UM m,UM *base,UM f,N e,UM r)
 {
 	int a,n,i;
 	N e1,an;
@@ -1177,10 +1452,7 @@ N e;
 	}
 }
 
-void spwrum0(mod,m,f,e,r)
-int mod;
-UM f,m,r;
-N e;
+void spwrum0(int mod,UM m,UM f,N e,UM r)
 {
 	UM t,s,q;
 	N e1;
@@ -1203,9 +1475,7 @@ N e;
 }
 
 #if 0
-void Mulum(mod,p1,p2,pr)
-register int mod;
-UM p1,p2,pr;
+void Mulum(int mod,UM p1,UM p2,UM pr)
 {
 	register int *pc1,*pcr;
 	register int mul,i,j,d1,d2;
@@ -1224,9 +1494,7 @@ UM p1,p2,pr;
 	DEG(pr) = d1 + d2;
 }
 
-void Mulsum(mod,p,n,pr)
-register int mod,n;
-UM p,pr;
+void Mulsum(int mod,UM p,int n,UM pr)
 {
 	register int *sp,*dp;
 	register int i;
@@ -1236,9 +1504,7 @@ UM p,pr;
 		*dp = (*sp * n) % mod;
 }
 	
-int Divum(mod,p1,p2,pq)
-register int mod;
-UM p1,p2,pq;
+int Divum(int mod,UM p1,UM p2,UM pq)
 {
 	register int *pc1,*pct;
 	register int tmp,i,j,inv;
@@ -1274,9 +1540,7 @@ UM p1,p2,pq;
 	return( i );
 }
 
-void Gcdum(mod,p1,p2,pr)
-register int mod;
-UM p1,p2,pr;
+void Gcdum(int mod,UM p1,UM p2,UM pr)
 {
 	register int *sp,*dp;
 	register int i,inv;
@@ -1303,10 +1567,7 @@ UM p1,p2,pr;
 }
 #endif
 
-void mult_mod_tab(p,mod,tab,r,d)
-UM p,r;
-UM *tab;
-int mod,d;
+void mult_mod_tab(UM p,int mod,UM *tab,UM r,int d)
 {
 	UM w,w1,c;
 	int n,i;
@@ -1324,11 +1585,7 @@ int mod,d;
 		}
 }
 
-void make_qmat(p,mod,tab,mp)
-UM p;
-int mod;
-UM *tab;
-int ***mp;
+void make_qmat(UM p,int mod,UM *tab,int ***mp)
 {
 	int n,i,j;
 	int *c;
@@ -1347,10 +1604,7 @@ int ***mp;
 		mat[i][i] = (mat[i][i]+mod-1) % mod;
 }
 
-void null_mod(mat,mod,n,ind)
-int **mat;
-int *ind;
-int mod,n;
+void null_mod(int **mat,int mod,int n,int *ind)
 {
 	int i,j,l,s,h,inv;
 	int *t,*u;
@@ -1381,11 +1635,7 @@ int mod,n;
 	}
 }
 
-void null_to_sol(mat,ind,mod,n,r)
-int **mat;
-int *ind;
-int mod,n;
-UM *r;
+void null_to_sol(int **mat,int *ind,int mod,int n,UM *r)
 {
 	int i,j,k,l;
 	int *c;
@@ -1414,9 +1664,7 @@ null_mod(mat,mod,n,ind)
 null_to_sol(mat,ind,mod,n,r)
 */
 
-void newddd(mod,f,r)
-int mod;
-UM f,*r;
+void newddd(int mod,UM f,UM *r)
 {
 	register int i,j;
 	int d,n;
@@ -1459,9 +1707,7 @@ UM f,*r;
 	r[j] = 0;
 }
 
-int nfctr_mod(f,mod)
-int mod;
-UM f;
+int nfctr_mod(UM f,int mod)
 {
 	register int i,j;
 	int d,n;
@@ -1502,9 +1748,7 @@ UM f;
 	return j;
 }
 
-int irred_check(f,mod)
-UM f;
-int mod;
+int irred_check(UM f,int mod)
 {
 	register int i,j;
 	int d,n;
@@ -1547,10 +1791,7 @@ int mod;
 	return 1;
 }
 
-int berlekamp(p,mod,df,tab,r)
-UM p;
-int mod,df;
-UM *tab,*r;
+int berlekamp(UM p,int mod,int df,UM *tab,UM *r)
 {
 	int n,i,j,k,nf,d,nr;
 	int **mat;
@@ -1596,11 +1837,12 @@ UM *tab,*r;
 			}
 		}
 	}
+	/* NOTREACHED */
+	error("berlekamp : cannot happen");
+	return -1;
 }
 
-void minipoly_mod(mod,f,p,mp)
-int mod;
-UM f,p,mp;
+void minipoly_mod(int mod,UM f,UM p,UM mp)
 {
 	struct p_pair *list,*l,*l1,*lprev;	
 	int n,d;
@@ -1635,11 +1877,7 @@ UM f,p,mp;
 	}
 }
 
-void lnf_mod(mod,n,p0,p1,list,np0,np1)
-int mod,n;
-UM p0,p1;
-struct p_pair *list;
-UM np0,np1;
+void lnf_mod(int mod,int n,UM p0,UM p1,struct p_pair *list,UM np0,UM np1)
 {
 	int inv,h,d1;
 	UM t0,t1,s0,s1;
@@ -1659,10 +1897,7 @@ UM np0,np1;
 	}
 }
 
-int find_root(mod,p,root)
-int mod;
-UM p;
-int *root;
+int find_root(int mod,UM p,int *root)
 {
 	UM *r;
 	int i,j;
@@ -1675,8 +1910,7 @@ int *root;
 	return j;
 }
 
-void showum(p)
-UM p;
+void showum(UM p)
 {
 	int i;
 	int *c;
@@ -1687,9 +1921,7 @@ UM p;
 	printf("\n");
 }
 
-void showumat(mat,n)
-int **mat;
-int n;
+void showumat(int **mat,int n)
 {
 	int i,j;