===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/E.c,v
retrieving revision 1.1.1.1
retrieving revision 1.8
diff -u -p -r1.1.1.1 -r1.8
--- OpenXM_contrib2/asir2000/engine/E.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/engine/E.c	2001/10/09 01:36:09	1.8
@@ -1,13 +1,55 @@
-/* $OpenXM: OpenXM/src/asir99/engine/E.c,v 1.1.1.1 1999/11/10 08:12:26 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/engine/E.c,v 1.7 2001/06/07 04:54:40 noro Exp $ 
+*/
 #include "ca.h"
 
-void henmv(vl,vn,f,g0,h0,a0,b0,lg,lh,lg0,lh0,q,k,gp,hp)
-VL vl;
-VN vn;
-P f,g0,h0,a0,b0,lg,lh,lg0,lh0;
-Q q;
-int k;
-P *gp,*hp;
+void henmv(VL vl,VN vn,P f,P g0,P h0,P a0,P b0,P lg,P lh,P lg0,P lh0,Q q,int k,P *gp,P *hp)
 {
 	P g1,h1,a1,b1;
 	N qn;
@@ -18,13 +60,7 @@ P *gp,*hp;
 	henmvmain(vl,vn,f,g1,h1,b1,a1,lg,lh,q,q2,k,gp,hp);
 }
 
-void henmvmain(vl,vn,f,fi0,fi1,gi0,gi1,l0,l1,mod,mod2,k,fr0,fr1)
-VL vl;
-VN vn;
-P f,fi0,fi1,gi0,gi1,l0,l1;
-Q mod,mod2;
-int k;
-P *fr0,*fr1;
+void henmvmain(VL vl,VN vn,P f,P fi0,P fi1,P gi0,P gi1,P l0,P l1,Q mod,Q mod2,int k,P *fr0,P *fr1)
 {
 	V v;
 	int n,i,j;
@@ -116,13 +152,7 @@ END:
 	output : gr1 * fr0 + gr0 * fr1 = 1 mod qr; ( qr = p^(2^(k-1)) )
 */
 
-void henzq(f,i0,fi0,i1,fi1,p,k,fr0p,fr1p,gr0p,gr1p,qrp)
-P f;
-UM fi0,fi1;
-int p,k;
-P i0,i1;
-P *fr0p,*fr1p,*gr0p,*gr1p;
-Q *qrp;
+void henzq(P f,P i0,UM fi0,P i1,UM fi1,int p,int k,P *fr0p,P *fr1p,P *gr0p,P *gr1p,Q *qrp)
 {
 	N qn;
 	Q q,qq,q2;
@@ -215,11 +245,7 @@ Q *qrp;
 	*fr0p = f0; *fr1p = f1; *gr0p = g0; *gr1p = g1; *qrp = q;
 }
 
-void henzq1(g,h,bound,gcp,hcp,qp)
-P g,h;
-Q bound;
-P *gcp,*hcp;
-Q *qp;
+void henzq1(P g,P h,Q bound,P *gcp,P *hcp,Q *qp)
 {
 	V v;
 	Q f,q,q1;
@@ -246,9 +272,7 @@ Q *qp;
 	wm = W_UMALLOC(m+n); 
 
 	for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
-		mod = lprime[index++];
-		if ( !mod )
-			error("henzq1 : lprime[] exhausted.");
+		mod = get_lprime(index++);
 		if ( !rem(NM((Q)LC(g)),mod) || !rem(NM((Q)LC(h)),mod) ) 
 			continue;
 		ptomp(mod,g,&tg); ptomp(mod,h,&th);
@@ -294,9 +318,7 @@ Q *qp;
 		mulp(CO,gc0,(P)q,&c); gc0 = c; mulp(CO,hc0,(P)q,&c); hc0 = c;
 	}
 	for ( index = 0; ; ) {
-		mod = lprime[index++];
-		if ( !mod )
-			error("henzq1 : lprime[] exhausted.");
+		mod = get_lprime(index++);
 		if ( !rem(NM((Q)zzz),mod) ||
 			!rem(NM((Q)LC(g)),mod) || 
 			!rem(NM((Q)LC(h)),mod) ) 
@@ -311,10 +333,7 @@ Q *qp;
 	}
 }
 
-void addm2p(vl,mod,mod2,n1,n2,nr)
-VL vl;
-Q mod,mod2;
-P n1,n2,*nr;
+void addm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 {
 	P t;
 
@@ -325,10 +344,7 @@ P n1,n2,*nr;
 		cm2p(mod,mod2,t,nr); 
 }
 
-void subm2p(vl,mod,mod2,n1,n2,nr)
-VL vl;
-Q mod,mod2;
-P n1,n2,*nr;
+void subm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 {
 	P t;
 
@@ -339,10 +355,7 @@ P n1,n2,*nr;
 		cm2p(mod,mod2,t,nr); 
 }
 
-void mulm2p(vl,mod,mod2,n1,n2,nr)
-VL vl;
-Q mod,mod2;
-P n1,n2,*nr;
+void mulm2p(VL vl,Q mod,Q mod2,P n1,P n2,P *nr)
 {
 	P t;
 
@@ -353,9 +366,7 @@ P n1,n2,*nr;
 		cm2p(mod,mod2,t,nr); 
 }
 
-void cmp(mod,p,pr)
-Q mod;
-P p,*pr;
+void cmp(Q mod,P p,P *pr)
 {
 	P t;
 	DCP dc,dcr,dcr0;
@@ -382,9 +393,7 @@ P p,*pr;
 	}
 }
 
-void cm2p(mod,m,p,pr)
-Q mod,m;
-P p,*pr;
+void cm2p(Q mod,Q m,P p,P *pr)
 {
 	P t;
 	DCP dc,dcr,dcr0;
@@ -411,9 +420,7 @@ P p,*pr;
 	}
 }
 
-void addm2q(mod,mod2,n1,n2,nr)
-Q mod,mod2;
-Q n1,n2,*nr;
+void addm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 {
 	Q t;
 
@@ -424,9 +431,7 @@ Q n1,n2,*nr;
 		rem2q(t,mod,mod2,nr); 
 }
 
-void subm2q(mod,mod2,n1,n2,nr)
-Q mod,mod2;
-Q n1,n2,*nr;
+void subm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 {
 	Q t;
 
@@ -437,9 +442,7 @@ Q n1,n2,*nr;
 		rem2q(t,mod,mod2,nr); 
 }
 
-void mulm2q(mod,mod2,n1,n2,nr)
-Q mod,mod2;
-Q n1,n2,*nr;
+void mulm2q(Q mod,Q mod2,Q n1,Q n2,Q *nr)
 {
 	Q t;
 
@@ -450,8 +453,7 @@ Q n1,n2,*nr;
 		rem2q(t,mod,mod2,nr); 
 }
 
-void rem2q(n,m,m2,nr)
-Q n,m,m2,*nr;
+void rem2q(Q n,Q m,Q m2,Q *nr)
 {
 	N q,r,s;
 	int sgn;
@@ -468,12 +470,38 @@ Q n,m,m2,*nr;
 	}
 }
 
-void exthp(vl,p,d,pr)
-VL vl;
-P p;
-int d;
-P *pr;
+/* 
+	extract d-homogeneous part with respect to vl - {v}
+*/
+
+void exthpc_generic(VL vl,P p,int d,V v,P *pr)
 {
+	P w,x,t,t1,a,xd;
+	V v0;
+	DCP dc;
+
+	if ( d < 0 || !p )
+		*pr = 0;
+	else if ( NUM(p) )
+		if ( d == 0 )
+			*pr = p;
+		else
+			*pr = 0;
+	else if ( v == VR(p) )
+		exthpc(vl,v,p,d,pr);
+	else {
+		v0 = VR(p);
+		for ( MKV(v0,x), dc = DC(p), w = 0; dc; dc = NEXT(dc) ) {
+			exthpc_generic(vl,COEF(dc),d-QTOS(DEG(dc)),v,&t);
+			pwrp(vl,x,DEG(dc),&xd);
+			mulp(vl,xd,t,&t1); addp(vl,w,t1,&a); w = a;
+		}
+		*pr = w;
+	}
+}
+
+void exthp(VL vl,P p,int d,P *pr)
+{
 	P t,t1,a,w,x,xd;
 	DCP dc;
 
@@ -494,12 +522,7 @@ P *pr;
 	}
 }
 
-void exthpc(vl,v,p,d,pr)
-VL vl;
-V v;
-P p;
-int d;
-P *pr;
+void exthpc(VL vl,V v,P p,int d,P *pr)
 {
 	P t,t1,a,w,x,xd;
 	DCP dc;
@@ -518,10 +541,7 @@ P *pr;
 	}
 }
 
-void cbound(vl,p,b)
-VL vl;
-P p;
-Q *b;
+void cbound(VL vl,P p,Q *b)
 {
 	Q n,e,t,m;
 	int k;
@@ -536,9 +556,7 @@ Q *b;
 	mulq(m,n,b); 
 }
 
-int geldb(vl,p)
-VL vl;
-P p;
+int geldb(VL vl,P p)
 {
 	int m;
 
@@ -547,11 +565,9 @@ P p;
 	return ( m );
 }
 
-int getdeg(v,p)
-V v;
-P p;
+int getdeg(V v,P p)
 {
-	int m;
+	int m,t;
 	DCP dc;
 	
 	if ( !p || NUM(p) ) 
@@ -559,15 +575,15 @@ P p;
 	else if ( v == VR(p) ) 
 		return ( deg(v,p) );
 	else {
-		for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) 
-			m = MAX(m,getdeg(v,COEF(dc)));
+		for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
+			t = getdeg(v,COEF(dc));
+			m = MAX(m,t);
+		}
 		return ( m );
 	}
 }
 
-void cmax(p,b)
-P p;
-Q *b;
+void cmax(P p,Q *b)
 {
 	DCP dc;
 	Q m,m1;
@@ -586,9 +602,7 @@ Q *b;
 	}
 }
 
-int nextbin(vn,n)
-VN vn;
-int n;
+int nextbin(VN vn,int n)
 {
 	int tmp,i,carry;
 
@@ -603,9 +617,7 @@ int n;
 	return ( carry );
 }
 
-void mulsgn(vn,vnt,n,vn1)
-VN vn,vnt,vn1;
-int n;
+void mulsgn(VN vn,VN vnt,int n,VN vn1)
 {
 	int i;
 
@@ -616,8 +628,7 @@ int n;
 			vn1[(int)vnt[i].v].n *= -1;
 }
 
-void next(vn)
-VN vn;
+void next(VN vn)
 {
 	int i,m,n,tmp,carry;
 
@@ -644,10 +655,7 @@ VN vn;
 	}
 }
 		
-void clctv(vl,p,nvlp)
-VL vl;
-P p;
-VL *nvlp;
+void clctv(VL vl,P p,VL *nvlp)
 {
 	int i,n;
 	VL tvl;
@@ -669,10 +677,7 @@ VL *nvlp;
 	vntovl(tvn,n,nvlp);
 }
 
-void markv(vn,n,p)
-VN vn;
-int n;
-P p;
+void markv(VN vn,int n,P p)
 {
 	V v;
 	DCP dc;
@@ -690,10 +695,7 @@ P p;
 		markv(vn,n,COEF(dc));
 }
 			
-void vntovl(vn,n,vlp)
-VN vn;
-int n;
-VL *vlp;
+void vntovl(VN vn,int n,VL *vlp)
 {
 	int i;
 	VL tvl,tvl0;
@@ -718,11 +720,9 @@ VL *vlp;
 	*vlp = tvl0;
 }
 
-int dbound(v,f)
-V v;
-P f;
+int dbound(V v,P f)
 {
-	int m;
+	int m,t;
 	DCP dc;
 
 	if ( !f ) 
@@ -730,16 +730,17 @@ P f;
 	else if ( v != VR(f) ) 
 		return homdeg(f);
 	else {
-		for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) 
-			m = MAX(m,homdeg(COEF(dc)));
+		for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
+			t = homdeg(COEF(dc));
+			m = MAX(m,t);
+		}
 		return ( m );
 	}
 }
 
-int homdeg(f)
-P f;
+int homdeg(P f)
 {
-	int m;
+	int m,t;
 	DCP dc;
 
 	if ( !f ) 
@@ -747,16 +748,17 @@ P f;
 	else if ( NUM(f) ) 
 		return( 0 );
 	else {
-		for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) 
-			m = MAX(m,QTOS(DEG(dc))+homdeg(COEF(dc)));
+		for ( dc = DC(f), m = 0; dc; dc = NEXT(dc) ) {
+			t = QTOS(DEG(dc))+homdeg(COEF(dc));
+			m = MAX(m,t);
+		}
 		return ( m );
 	}
 }
 
-int minhomdeg(f)
-P f;
+int minhomdeg(P f)
 {
-	int m;
+	int m,t;
 	DCP dc;
 
 	if ( !f ) 
@@ -764,17 +766,15 @@ P f;
 	else if ( NUM(f) ) 
 		return( 0 );
 	else {
-		for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) 
-			m = MIN(m,QTOS(DEG(dc))+minhomdeg(COEF(dc)));
+		for ( dc = DC(f), m = homdeg(f); dc; dc = NEXT(dc) ) {
+			t = QTOS(DEG(dc))+minhomdeg(COEF(dc));
+			m = MIN(m,t);
+		}
 		return ( m );
 	}
 }
 
-void adjc(vl,f,a,lc0,q,fr,ar)
-VL vl;
-P f,a,lc0;
-Q q;
-P *fr,*ar;
+void adjc(VL vl,P f,P a,P lc0,Q q,P *fr,P *ar)
 {
 	P m,m1;
 	Q t;
@@ -788,13 +788,7 @@ P *fr,*ar;
 }
 
 #if 1
-void affinemain(vl,p,v0,n,pl,pr)
-VL vl;
-V v0;
-int n;
-P *pl;
-P p;
-P *pr;
+void affinemain(VL vl,P p,V v0,int n,P *pl,P *pr)
 {
 	P x,t,m,c,s,a;
 	DCP dc;
@@ -833,11 +827,7 @@ P *pr;
 #endif
 
 #if 0
-affine(vl,p,vn,r)
-VL vl;
-P p;
-VN vn;
-P *r;
+affine(VL vl,P p,VN vn,P *r)
 {
 	int n,d,d1,i;
 	Q *q;
@@ -855,12 +845,7 @@ P *r;
 	}
 }
 
-afmain(vl,bc,p,vn,r)
-VL vl;
-Q **bc;
-P p;
-VN vn;
-P *r;
+afmain(VL vl,Q **bc,P p,VN vn,P *r)
 {
 	P t,s,u;
 	P *c,*rc;
@@ -901,11 +886,7 @@ P *r;
 }
 #endif
 
-void restore(vl,f,vn,fr)
-VL vl;
-P f;
-VN vn;
-P *fr;
+void restore(VL vl,P f,VN vn,P *fr)
 {
 	int i;
 	P vv,g,g1,t;
@@ -925,8 +906,7 @@ P *fr;
 	*fr = g;
 }
 
-void mergev(vl,vl1,vl2,nvlp)
-VL vl,vl1,vl2,*nvlp;
+void mergev(VL vl,VL vl1,VL vl2,VL *nvlp)
 {
 	int i,n;
 	VL tvl;
@@ -962,11 +942,7 @@ VL vl,vl1,vl2,*nvlp;
 }
 
 #if 0
-void substvp(vl,f,vn,g)
-VL vl;
-P f;
-VN vn;
-P *g;
+void substvp(VL vl,P f,VN vn,P *g)
 {
 	V v;
 	int i;
@@ -981,11 +957,7 @@ P *g;
 	*g = h;
 }
 
-void affine(vl,f,vn,fr)
-VL vl;
-P f;
-VN vn;
-P *fr;
+void affine(VL vl,P f,VN vn,P *fr)
 {
 	int i,j,n;
 	P vv,g,g1,t,u;