===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.97
retrieving revision 1.102
diff -u -p -r1.97 -r1.102
--- OpenXM_contrib2/asir2000/engine/nd.c	2004/03/25 01:31:03	1.97
+++ OpenXM_contrib2/asir2000/engine/nd.c	2004/09/15 01:43:33	1.102
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.96 2004/03/17 08:16:24 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.101 2004/09/14 10:00:26 noro Exp $ */
 
 #include "nd.h"
 
@@ -2702,6 +2702,35 @@ UINT *ndv_compute_bound(NDV p)
 	return t;
 }
 
+UINT *nd_compute_bound(ND p)
+{
+	UINT *d1,*d2,*t;
+	UINT u;
+	int i,j,k,l,len,ind;
+	NM m;
+
+	if ( !p )
+		return 0;
+	d1 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	d2 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	len = LEN(p);
+	m = BDY(p); ndl_copy(DL(m),d1); m = NEXT(m);
+	for ( m = NEXT(m); m; m = NEXT(m) ) {
+		ndl_lcm(DL(m),d1,d2);
+		t = d1; d1 = d2; d2 = t;
+	}
+	l = nd_nvar+31;
+	t = (UINT *)MALLOC_ATOMIC(l*sizeof(UINT));
+	for ( i = nd_exporigin, ind = 0; i < nd_wpd; i++ ) {
+		u = d1[i];
+		k = (nd_epw-1)*nd_bpe;
+		for ( j = 0; j < nd_epw; j++, k -= nd_bpe, ind++ )
+			t[ind] = (u>>k)&nd_mask0;
+	}
+	for ( ; ind < l; ind++ ) t[ind] = 0;
+	return t;
+}
+
 int nd_get_exporigin(struct order_spec *ord)
 {
 	switch ( ord->id ) {
@@ -3167,6 +3196,44 @@ ND ndv_mul_nm(int mod,NM m0,NDV p)
 	}
 }
 
+ND nd_quo(int mod,ND p,NDV d)
+{
+	NM mq0,mq;
+	NMV tm;
+	Q q;
+	int i,nv,sg,c,c1,c2;
+	ND t,r;
+	
+	if ( !p ) return 0;
+	else {
+		nv = NV(p);
+		sg = SG(p);
+		mq0 = 0;
+		tm = (NMV)ALLOCA(nmv_adv);
+		while ( p ) {
+			NEXTNM(mq0,mq);
+			ndl_sub(HDL(p),HDL(d),DL(mq));
+			ndl_copy(DL(mq),DL(tm));
+			if ( mod ) {
+				c1 = invm(HCM(d),mod); c2 = HCM(p);
+				DMAR(c1,c2,0,mod,c); CM(mq) = c;
+				CM(tm) = mod-c;
+			} else {
+				divq(HCQ(p),HCQ(d),&CQ(mq));
+				chsgnq(CQ(mq),&CQ(tm));
+			}
+			t = ndv_mul_nmv_trunc(mod,tm,d,HDL(d));
+			p = nd_add(mod,p,t);
+		}
+		NEXT(mq) = 0;
+		for ( i = 0, mq = mq0; mq; mq = NEXT(mq), i++ );
+		MKND(nv,mq0,i,r);
+		/* XXX */
+		SG(r) = sg-SG(d);
+		return r;
+	}
+}
+
 void ndv_realloc(NDV p,int obpe,int oadv,EPOS oepos)
 {
 	NMV m,mr,mr0,t;
@@ -3316,7 +3383,7 @@ ND ptond(VL vl,VL dvl,P p)
 		w = (DCP *)ALLOCA(k*sizeof(DCP));
 		for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ ) w[j] = dc;
 		for ( i = 0, tvl = dvl, v = VR(p);
-			vl && tvl->v != v; tvl = NEXT(tvl), i++ );
+			tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
 		if ( !tvl ) {
 			for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
 				t = ptond(vl,dvl,COEF(w[j]));
@@ -4586,4 +4653,150 @@ NDV ndv_load(int index)
 	MKNDV(nv,m0,len,d);
 	SG(d) = sugar;
 	return d;
+}
+
+void nd_det(int mod,MAT f,P *rp)
+{
+	VL fv,tv;
+	int n,i,j,max,e,nvar,sgn,k0,l0,len0,len,k,l,a;
+	pointer **m;
+	Q mone;
+	NDV **dm;
+	NDV *t,*mi,*mj;
+	NDV d,s,mij,mjj;
+	ND u;
+	NMV nmv;
+	PGeoBucket bucket;
+	struct order_spec *ord;
+
+	create_order_spec(0,0,&ord);
+	nd_init_ord(ord);
+	get_vars((Obj)f,&fv);
+	if ( f->row != f->col )
+		error("nd_det : non-square matrix");
+	n = f->row;
+	for ( nvar = 0, tv = fv; tv; tv = NEXT(tv), nvar++ );
+	m = f->body;
+	for ( i = 0, max = 0; i < n; i++ )
+		for ( j = 0; j < n; j++ )
+			for ( tv = fv; tv; tv = NEXT(tv) ) {
+				e = getdeg(tv->v,(P)m[i][j]);
+				max = MAX(e,max);
+			}
+	nd_setup_parameters(nvar,1024);
+	dm = (NDV **)almat_pointer(n,n);
+	for ( i = 0, max = 0; i < n; i++ )
+		for ( j = 0; j < n; j++ ) {
+			dm[i][j] = ptondv(CO,fv,m[i][j]);
+			if ( mod ) ndv_mod(mod,dm[i][j]);
+			if ( dm[i][j] && !LEN(dm[i][j]) ) dm[i][j] = 0;
+		}
+	d = ptondv(CO,fv,(P)ONE);
+	if ( mod ) ndv_mod(mod,d);
+	chsgnq(ONE,&mone);
+	for ( j = 0, sgn = 1; j < n; j++ ) {
+		for ( i = j; i < n && !dm[i][j]; i++ );
+		if ( i == n ) {
+			*rp = 0;
+			return;
+		}
+		k0 = i; l0 = j; len0 = LEN(dm[k0][l0]);
+		for ( k = j; k < n; k++ )
+			for ( l = j; l < n; l++ )
+				if ( dm[k][l] && LEN(dm[k][l]) < len0 ) {
+					k0 = k; l0 = l; len0 = LEN(dm[k][l]);
+				}
+		if ( k0 != j ) {
+			t = dm[j]; dm[j] = dm[k0]; dm[k0] = t;
+			sgn = -sgn;
+		}
+		if ( l0 != j ) {
+			for ( k = j; k < n; k++ ) {
+				s = dm[k][j]; dm[k][j] = dm[k][l0]; dm[k][l0] = s;
+			}
+			sgn = -sgn;
+		}
+		for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) {
+			mi = dm[i]; mij = mi[j];
+			if ( mod )
+				ndv_mul_c(mod,mij,mod-1);
+			else
+				ndv_mul_c_q(mij,mone);
+			for ( k = j+1; k < n; k++ ) {
+				bucket = create_pbucket();
+				if ( mi[k] )
+					nmv = BDY(mjj); len = LEN(mjj);
+					for ( a = 0; a < len; a++, NMV_ADV(nmv) ) {
+						u = ndv_mul_nmv_trunc(mod,nmv,mi[k],DL(BDY(d)));
+						add_pbucket(mod,bucket,u);
+					}
+				if ( mj[k] && mij ) {
+					nmv = BDY(mij); len = LEN(mij);
+					for ( a = 0; a < len; a++, NMV_ADV(nmv) ) {
+						u = ndv_mul_nmv_trunc(mod,nmv,mj[k],DL(BDY(d)));
+						add_pbucket(mod,bucket,u);
+					}
+				}
+				u = normalize_pbucket(mod,bucket);
+				u = nd_quo(mod,u,d);
+				mi[k] = ndtondv(mod,u);
+			}
+		}
+		d = mjj;
+	}
+	if ( sgn < 0 )
+		if ( mod )
+			ndv_mul_c(mod,d,mod-1);
+		else
+			ndv_mul_c_q(d,mone);
+	*rp = ndvtop(mod,CO,fv,d);
+}
+
+ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d)
+{
+	NM mr,mr0;
+	NM tnm;
+	NMV m;
+	UINT *d0,*dt,*dm;
+	int c,n,td,i,c1,c2,len;
+	Q q;
+	ND r;
+
+	if ( !p ) return 0;
+	else {
+		n = NV(p); m = BDY(p); len = LEN(p);
+		d0 = DL(m0);
+		td = TD(d);
+		mr0 = 0;
+		NEWNM(tnm);
+		if ( mod ) {
+			c = CM(m0);
+			for ( i = 0; i < len; i++, NMV_ADV(m) ) {
+				ndl_add(DL(m),d0,DL(tnm));
+				if ( ndl_reducible(DL(tnm),d) ) {
+					NEXTNM(mr0,mr);
+					c1 = CM(m); DMAR(c1,c,0,mod,c2); CM(mr) = c2;
+					ndl_add(DL(m),d0,DL(mr));
+				}
+			}
+		} else {
+			q = CQ(m0);
+			for ( i = 0; i < len; i++, NMV_ADV(m) ) {
+				ndl_add(DL(m),d0,DL(tnm));
+				if ( ndl_reducible(DL(tnm),d) ) {
+					NEXTNM(mr0,mr);
+					mulq(CQ(m),q,&CQ(mr));
+					ndl_add(DL(m),d0,DL(mr));
+				}
+			}
+		}
+		if ( !mr0 )
+			return 0;
+		else {
+			NEXT(mr) = 0; 
+			MKND(NV(p),mr0,len,r);
+			SG(r) = SG(p) + TD(d0);
+			return r;
+		}
+	}
 }