===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.110
retrieving revision 1.116
diff -u -p -r1.110 -r1.116
--- OpenXM_contrib2/asir2000/engine/nd.c	2004/09/21 05:23:14	1.110
+++ OpenXM_contrib2/asir2000/engine/nd.c	2004/12/01 12:36:17	1.116
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.109 2004/09/21 04:50:15 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.115 2004/11/18 08:29:11 noro Exp $ */
 
 #include "nd.h"
 
@@ -41,6 +41,9 @@ static int nd_found,nd_create,nd_notfirst;
 static int nmv_adv;
 static int nd_demand;
 
+UINT *nd_det_compute_bound(NDV **dm,int n,int j);
+void nd_det_reconstruct(NDV **dm,int n,int j,NDV d);
+
 void nd_free_private_storage()
 {
 	_nm_free_list = 0;
@@ -724,13 +727,11 @@ int ndl_disjoint(UINT *d1,UINT *d2)
 #endif
 }
 
-int ndl_check_bound2(int index,UINT *d2)
+int ndl_check_bound(UINT *d1,UINT *d2)
 {
 	UINT u2;
-	UINT *d1;
 	int i,j,ind,k;
 
-	d1 = nd_bound[index];
 	ind = 0;
 #if USE_UNROLL
 	switch ( nd_bpe ) {
@@ -819,6 +820,11 @@ int ndl_check_bound2(int index,UINT *d2)
 #endif
 }
 
+int ndl_check_bound2(int index,UINT *d2)
+{
+	return ndl_check_bound(nd_bound[index],d2);
+}
+
 INLINE int ndl_hash_value(UINT *d)
 {
 	int i;
@@ -2581,14 +2587,8 @@ void removecont_array(Q *c,int n)
 	Q d0,d1,a,u,u1,gcd;
 	int i,j;
 	N qn,rn,gn;
-	Q *t,*q,*r;
+	Q *q,*r;
 
-	t = (Q *)ALLOCA(n*sizeof(Q));
-	for ( i = j = 0; i < n; i++ )
-		if ( c[i] )
-			t[j++] = c[i];
-	n = j;
-	c = t;
 	q = (Q *)ALLOCA(n*sizeof(Q));
 	r = (Q *)ALLOCA(n*sizeof(Q));
 	v.id = O_VECT; v.len = n; v.body = (pointer *)c;
@@ -2622,6 +2622,7 @@ void nd_mul_c(int mod,ND p,int mul)
 	int c,c1;
 
 	if ( !p ) return;
+	if ( mul == 1 ) return;
 	if ( mod == -1 )
 		for ( m = BDY(p); m; m = NEXT(m) )
 			CM(m) = _mulsf(CM(m),mul);
@@ -2637,6 +2638,7 @@ void nd_mul_c_q(ND p,Q mul)
 	Q c;
 
 	if ( !p ) return;
+	if ( UNIQ(mul) ) return;
 	for ( m = BDY(p); m; m = NEXT(m) ) {
 		mulq(CQ(m),mul,&c); CQ(m) = c;
 	}
@@ -3715,6 +3717,10 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec
 	int stat,nvar,max,e;
 	union oNDC dn;
 
+	if ( !f ) {
+		*rp = 0;
+		return;
+	}
 	pltovl(v,&vv);
 	for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
 
@@ -4893,6 +4899,7 @@ void nd_det(int mod,MAT f,P *rp)
 	NDV d,s,mij,mjj;
 	ND u;
 	NMV nmv;
+	UINT *bound;
 	PGeoBucket bucket;
 	struct order_spec *ord;
 
@@ -4910,7 +4917,7 @@ void nd_det(int mod,MAT f,P *rp)
 				e = getdeg(tv->v,(P)m[i][j]);
 				max = MAX(e,max);
 			}
-	nd_setup_parameters(nvar,1024);
+	nd_setup_parameters(nvar,max);
 	dm = (NDV **)almat_pointer(n,n);
 	for ( i = 0, max = 0; i < n; i++ )
 		for ( j = 0; j < n; j++ ) {
@@ -4944,15 +4951,19 @@ void nd_det(int mod,MAT f,P *rp)
 			}
 			sgn = -sgn;
 		}
+		bound = nd_det_compute_bound(dm,n,j);
+		if ( ndl_check_bound(bound,bound) )
+			nd_det_reconstruct(dm,n,j,d);
+
 		for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) {
-			if ( DP_Print ) fprintf(stderr,"	i=%d\n		",i);
+/*			if ( DP_Print ) fprintf(stderr,"	i=%d\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++ ) {
-				if ( DP_Print ) fprintf(stderr,"k=%d ",k);
+/*				if ( DP_Print ) fprintf(stderr,"k=%d ",k); */
 				bucket = create_pbucket();
 				if ( mi[k] ) {
 					nmv = BDY(mjj); len = LEN(mjj);
@@ -4971,7 +4982,7 @@ void nd_det(int mod,MAT f,P *rp)
 				u = nd_quo(mod,bucket,d);
 				mi[k] = ndtondv(mod,u);
 			}
-			if ( DP_Print ) fprintf(stderr,"\n",k);
+/*			if ( DP_Print ) fprintf(stderr,"\n",k); */
 		}
 		d = mjj;
 	}
@@ -5031,4 +5042,57 @@ ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d)
 			return r;
 		}
 	}
+}
+
+void nd_det_reconstruct(NDV **dm,int n,int j,NDV d)
+{
+	int i,obpe,oadv,h,k,l;
+	static NM prev_nm_free_list;
+	EPOS oepos;
+
+	obpe = nd_bpe;
+	oadv = nmv_adv;
+	oepos = nd_epos;
+	if ( obpe < 2 ) nd_bpe = 2;
+	else if ( obpe < 3 ) nd_bpe = 3;
+	else if ( obpe < 4 ) nd_bpe = 4;
+	else if ( obpe < 5 ) nd_bpe = 5;
+	else if ( obpe < 6 ) nd_bpe = 6;
+	else if ( obpe < 8 ) nd_bpe = 8;
+	else if ( obpe < 10 ) nd_bpe = 10;
+	else if ( obpe < 16 ) nd_bpe = 16;
+	else if ( obpe < 32 ) nd_bpe = 32;
+	else error("nd_det_reconstruct : exponent too large");
+
+	nd_setup_parameters(nd_nvar,0);
+	prev_nm_free_list = _nm_free_list;
+	_nm_free_list = 0;
+	for ( k = j; k < n; k++ )
+		for (l = j; l < n; l++ )
+			ndv_realloc(dm[k][l],obpe,oadv,oepos);
+	ndv_realloc(d,obpe,oadv,oepos);
+	prev_nm_free_list = 0;
+#if 0
+	GC_gcollect();
+#endif
+}
+
+UINT *nd_det_compute_bound(NDV **dm,int n,int j)
+{
+	UINT *d0,*d1,*d,*t,*r;
+	int k,l;
+
+	d0 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	d1 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	for ( k = 0; k < nd_wpd; k++ ) d0[k] = 0;
+	for ( k = j; k < n; k++ )
+		for ( l = j; l < n; l++ )
+			if ( dm[k][l] ) {
+				d = ndv_compute_bound(dm[k][l]);
+				ndl_lcm(d,d0,d1);
+				t = d1; d1 = d0; d0 = t;
+			}
+	r = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	for ( k = 0; k < nd_wpd; k++ ) r[k] = d0[k];
+	return r;
 }