===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/Q.c,v
retrieving revision 1.4
retrieving revision 1.9
diff -u -p -r1.4 -r1.9
--- OpenXM_contrib2/asir2000/engine/Q.c	2000/08/21 08:31:27	1.4
+++ OpenXM_contrib2/asir2000/engine/Q.c	2002/08/14 04:49: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,14 +45,13 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.3 2000/05/29 08:54:46 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/Q.c,v 1.8 2002/08/08 10:57:01 noro Exp $ 
 */
 #include "ca.h"
 #include "base.h"
 #include "inline.h"
 
-void addq(n1,n2,nr)
-Q n1,n2,*nr;
+void addq(Q n1,Q n2,Q *nr)
 {
 	N nm,dn,nm1,nm2,nm3,dn0,dn1,dn2,g,g1,g0,m;
 	int sgn;
@@ -135,8 +134,7 @@ Q n1,n2,*nr;
 	}
 }
 
-void subq(n1,n2,nr)
-Q n1,n2,*nr;
+void subq(Q n1,Q n2,Q *nr)
 {
 	Q m;
 
@@ -155,8 +153,7 @@ Q n1,n2,*nr;
 	}
 }
 
-void mulq(n1,n2,nr)
-Q n1,n2,*nr;
+void mulq(Q n1,Q n2,Q *nr)
 {
 	N nm,nm1,nm2,dn,dn1,dn2,g;
 	int sgn;
@@ -202,8 +199,7 @@ Q n1,n2,*nr;
 	}
 }
 
-void divq(n1,n2,nq)
-Q n1,n2,*nq;
+void divq(Q n1,Q n2,Q *nq)
 {
 	Q m;
 
@@ -220,9 +216,28 @@ Q n1,n2,*nq;
 	}
 }
 
-void invq(n,nr)
-Q n,*nr;
+void divsq(Q n1,Q n2,Q *nq)
 {
+	int sgn;
+	N tn;
+
+	if ( !n2 ) {
+		error("division by 0");
+		*nq = 0;
+		return;
+	} else if ( !n1 ) 
+		*nq = 0;
+	else if ( n1 == n2 ) 
+		*nq = ONE;
+	else {
+		divsn(NM(n1),NM(n2),&tn);
+		sgn = SGN(n1)*SGN(n2);
+		NTOQ(tn,sgn,*nq);
+	}
+}
+
+void invq(Q n,Q *nr)
+{
 	N nm,dn;
 
 	if ( INT(n) )
@@ -242,8 +257,7 @@ Q n,*nr;
 	}
 }
 
-void chsgnq(n,nr)
-Q n,*nr;
+void chsgnq(Q n,Q *nr)
 {
 	Q t;
 
@@ -254,8 +268,7 @@ Q n,*nr;
 	}
 }
 
-void pwrq(n1,n,nr)
-Q n1,n,*nr;
+void pwrq(Q n1,Q n,Q *nr)
 {
 	N nm,dn;
 	int sgn;
@@ -295,8 +308,7 @@ Q n1,n,*nr;
 	}
 }
 
-int cmpq(q1,q2)
-Q q1,q2;
+int cmpq(Q q1,Q q2)
 {
 	int sgn;
 	N t,s;
@@ -330,9 +342,7 @@ Q q1,q2;
 	}
 }
 
-void remq(n,m,nr)
-Q n,m;
-Q *nr;
+void remq(Q n,Q m,Q *nr)
 {
 	N q,r,s;
 
@@ -352,9 +362,7 @@ Q *nr;
 
 /* t = [nC0 nC1 ... nCn] */
 
-void mkbc(n,t)
-int n;
-Q *t;
+void mkbc(int n,Q *t)
 {
 	int i;
 	N c,d;
@@ -374,9 +382,7 @@ Q *t;
  *  where W(k,l,i) = i! * kCi * lCi
  */
 
-void mkwc(k,l,t)
-int k,l;
-Q *t;
+void mkwc(int k,int l,Q *t)
 {
 	int i,n,up,low;
 	N nm,d,c;
@@ -396,9 +402,7 @@ Q *t;
 /* mod m table */
 /* XXX : should be optimized */
 
-void mkwcm(k,l,m,t)
-int k,l,m;
-int *t;
+void mkwcm(int k,int l,int m,int *t)
 {
 	int i,n;
 	Q *s;
@@ -411,9 +415,8 @@ int *t;
 	}
 }
 
-void factorial(n,r)
-int n;
-Q *r;
+#if 0
+void factorial(int n,Q *r)
 {
 	Q t,iq,s;
 	unsigned int i,m,m0;
@@ -436,10 +439,67 @@ Q *r;
 		*r = t;
 	}
 }
+#endif
 
-void invl(a,mod,ar)
-Q a,mod,*ar;
+void partial_factorial(int s,int e,N *r);
+
+void factorial(int n,Q *r)
 {
+	N nm;
+
+	if ( !n )
+		*r = ONE;
+	else if ( n < 0 ) 
+		*r = 0;
+	else {
+		partial_factorial(1,n,&nm);
+		NTOQ(nm,1,*r);
+	}
+}
+
+/* s*(s+1)*...*e */
+void partial_factorial(int s,int e,N *r)
+{
+	unsigned int len,i,m,m0,c;
+	unsigned int *p,*p1,*p2;
+	N nm,r1,r2;
+
+	/* XXX */
+	if ( e-s > 2000 ) {
+		m = (e+s)/2;
+		partial_factorial(s,m,&r1);
+		partial_factorial(m+1,e,&r2);
+		kmuln(r1,r2,r);	
+		return;
+	}
+	/* find i s.t. 2^(i-1) < m <= 2^i */
+	for ( i = 0, m = 1; m < e; m <<=1, i++ );
+	 /* XXX estimate of word length of the result */
+	len = (i*(e-s+1)+1+31)/32;
+	p = ALLOCA((len+1)*sizeof(int));
+	p1 = ALLOCA((len+1)*sizeof(int));
+	p[0] = s;
+	len = 1;
+	for ( i = s+1; (int)i <= e; ) {
+		for ( m0 = m = 1, c = 0; ((int)i <= e) && !c; i++ ) {
+			m0 = m; DM(m0,i,c,m)
+		}
+		if ( c ) {
+			m = m0; i--;
+		}
+		bzero(p1,(len+1)*sizeof(int));
+		muln_1(p,len,m,p1);
+		if ( p1[len] )
+			len++;
+		p2 = p; p = p1; p1 = p2;
+	}
+	*r = nm = NALLOC(len);
+	bcopy(p,BD(nm),sizeof(int)*len);
+	PL(nm) = len;
+}
+
+void invl(Q a,Q mod,Q *ar)
+{
 	Q f1,f2,a1,a2,q,m,s;
 	N qn,rn;
 
@@ -478,8 +538,7 @@ Q a,mod,*ar;
 
 int kara_mag=100;
 
-void kmuln(n1,n2,nr)
-N n1,n2,*nr;
+void kmuln(N n1,N n2,N *nr)
 {
 	N n,t,s,m,carry;
 	int d,d1,d2,len,i,l;
@@ -522,10 +581,7 @@ N n1,n2,*nr;
 	bcopy((char *)r0,(char *)BD(t),l*sizeof(int));
 }
 
-void extractn(n,index,len,nr)
-N n;
-int index,len;
-N *nr;
+void extractn(N n,int index,int len,N *nr)
 {
 	unsigned int *m;
 	int l;
@@ -547,25 +603,19 @@ N *nr;
 	}
 }
 
-void copyn(n,len,p)
-N n;
-int len;
-int *p;
+void copyn(N n,int len,int *p)
 {
 	if ( n )
 		bcopy((char *)BD(n),(char *)p,MIN(PL(n),len)*sizeof(int));
 }
 
-void dupn(n,p)
-N n;
-N p;
+void dupn(N n,N p)
 {
 	if ( n )
 		bcopy((char *)n,(char *)p,(1+PL(n))*sizeof(int));
 }
 
-void kmulnmain(n1,n2,nr)
-N n1,n2,*nr;
+void kmulnmain(N n1,N n2,N *nr)
 {
 	int d1,d2,h,sgn,sgn1,sgn2,len;
 	N n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;