===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.71
retrieving revision 1.73
diff -u -p -r1.71 -r1.73
--- OpenXM_contrib2/asir2000/engine/nd.c	2003/09/17 07:16:53	1.71
+++ OpenXM_contrib2/asir2000/engine/nd.c	2003/09/17 08:37:49	1.73
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.70 2003/09/15 10:51:45 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.72 2003/09/17 08:14:26 noro Exp $ */
 
 #include "ca.h"
 #include "inline.h"
@@ -2240,7 +2240,6 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe
 	int e,max,nvar;
 	NDV b;
 
-	nd_free_private_storage();
 	get_vars((Obj)f,&fv); pltovl(v,&vv);
 	nvar = length(vv);
 	nd_init_ord(ord);
@@ -2285,7 +2284,6 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru
 	nocheck = 0;
 	mindex = 0;
 
-	nd_free_private_storage();
 	/* setup modulus */
 	if ( trace < 0 ) {
 		trace = -trace;
@@ -2362,7 +2360,7 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,stru
 	}
 	/* dp->p */
 	nd_bpe = cbpe;
-	nd_setup_parameters(0,0);
+	nd_setup_parameters(nd_nvar,0);
 	for ( r = cand; r; r = NEXT(r) ) BDY(r) = (pointer)ndvtop(0,CO,vv,BDY(r));
 	MKLIST(*rp,cand);
 }
@@ -2744,23 +2742,23 @@ int nd_get_exporigin(struct order_spec *ord)
 void nd_setup_parameters(int nvar,int max) {
 	int i,j,n,elen,ord_o,ord_l,l,s;
 	struct order_pair *op;
+	int bpe;
 
-	/* if max == 0, don't touch nd_bpe */
-	if ( max > 0 ) {
-		if ( max < 2 ) nd_bpe = 1;
-		if ( max < 4 ) nd_bpe = 2;
-		else if ( max < 8 ) nd_bpe = 3;
-		else if ( max < 16 ) nd_bpe = 4;
-		else if ( max < 32 ) nd_bpe = 5;
-		else if ( max < 64 ) nd_bpe = 6;
-		else if ( max < 256 ) nd_bpe = 8;
-		else if ( max < 1024 ) nd_bpe = 10;
-		else if ( max < 65536 ) nd_bpe = 16;
-		else nd_bpe = 32;
-	}
-	/* nvar == 0, don't touch nd_nvar */
-	if ( nvar > 0 ) nd_nvar = nvar;
-
+	if ( !max ) bpe = nd_bpe;
+	else if ( max < 2 ) bpe = 1;
+	else if ( max < 4 ) bpe = 2;
+	else if ( max < 8 ) bpe = 3;
+	else if ( max < 16 ) bpe = 4;
+	else if ( max < 32 ) bpe = 5;
+	else if ( max < 64 ) bpe = 6;
+	else if ( max < 256 ) bpe = 8;
+	else if ( max < 1024 ) bpe = 10;
+	else if ( max < 65536 ) bpe = 16;
+	else bpe = 32;
+	if ( bpe != nd_bpe || nvar != nd_nvar )
+		nd_free_private_storage();
+	nd_bpe = bpe;
+	nd_nvar = nvar;
 	nd_epw = (sizeof(UINT)*8)/nd_bpe;
 	elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
 
@@ -2809,7 +2807,7 @@ ND_pairs nd_reconstruct(int mod,int trace,ND_pairs d)
 	else if ( obpe < 32 ) nd_bpe = 32;
 	else error("nd_reconstruct : exponent too large");
 
-	nd_setup_parameters(0,0);
+	nd_setup_parameters(nd_nvar,0);
 	prev_nm_free_list = _nm_free_list;
 	prev_ndp_free_list = _ndp_free_list;
 	_nm_free_list = 0;
@@ -3806,6 +3804,52 @@ void ndv_reduce_vect(int m,UINT *svect,int col,IndArra
 		if ( svect[i] >= (UINT)m ) svect[i] %= m;
 }
 
+void ndv_reduce_vect_sf(int m,UINT *svect,int col,IndArray *imat,NODE rp0)
+{
+	int i,j,k,len,pos,prev;
+	UINT c,c1,c2,c3,up,lo,dmy;
+	IndArray ivect;
+	unsigned char *ivc;
+	unsigned short *ivs;
+	unsigned int *ivi;
+	NDV redv;
+	NMV mr;
+	NODE rp;
+
+	for ( rp = rp0, i = 0; rp; i++, rp = NEXT(rp) ) {
+		ivect = imat[i];
+		k = ivect->head; svect[k] %= m;
+		if ( c = svect[k] ) {
+			c = _chsgnsf(c); redv = nd_ps[((NM_ind_pair)BDY(rp))->index];
+			len = LEN(redv); mr = BDY(redv);
+			svect[k] = 0; prev = k;
+			switch ( ivect->width ) {
+				case 1:
+					ivc = ivect->index.c;
+					for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
+						pos = prev+ivc[j]; prev = pos;
+						svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]);
+					}
+					break;
+				case 2:
+					ivs = ivect->index.s; 
+					for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
+						pos = prev+ivs[j]; prev = pos;
+						svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]);
+					}
+					break;
+				case 4:
+					ivi = ivect->index.i;
+					for ( j = 1, NMV_ADV(mr); j < len; j++, NMV_ADV(mr) ) {
+						pos = prev+ivi[j]; prev = pos;
+						svect[pos] = _addsf(_mulsf(CM(mr),c),svect[pos]);
+					}
+					break;
+			}
+		}
+	}
+}
+
 NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhead,UINT *s0vect)
 {
 	int j,k,len;
@@ -3881,7 +3925,8 @@ int nd_symbolic_preproc(PGeoBucket bucket,UINT **s0vec
 		}
 		col++;
 	}
-	NEXT(rp) = 0; NEXT(s) = 0;
+	if ( rp0 ) NEXT(rp) = 0;
+	NEXT(s) = 0;
 	s0v = (UINT *)MALLOC_ATOMIC(col*nd_wpd*sizeof(UINT));
 	for ( i = 0, p = s0v, s = s0; i < col;
 		i++, p += nd_wpd, s = NEXT(s) ) ndl_copy(DL(s),p);
@@ -4073,7 +4118,8 @@ NODE nd_f4(int m)
 		svect = (UINT *)MALLOC_ATOMIC(col*sizeof(UINT));
 		for ( a = sprow = 0, sp = sp0; a < nsp; a++, sp = NEXT(sp) ) {
 			nd_to_vect(m,s0vect,col,BDY(sp),svect);
-			ndv_reduce_vect(m,svect,col,imat,rp0);
+			if ( m == -1 ) ndv_reduce_vect_sf(m,svect,col,imat,rp0);
+			else ndv_reduce_vect(m,svect,col,imat,rp0);
 			for ( i = 0; i < col; i++ ) if ( svect[i] ) break;
 			if ( i < col ) {
 				spmat[sprow] = v = (UINT *)MALLOC_ATOMIC(spcol*sizeof(UINT));
@@ -4087,7 +4133,10 @@ NODE nd_f4(int m)
 
 		/* elimination (2nd step) */
 		colstat = (int *)ALLOCA(spcol*sizeof(int));
-		rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);
+		if ( m == -1 )
+			rank = generic_gauss_elim_sf(spmat,sprow,spcol,m,colstat);
+		else
+			rank = generic_gauss_elim_mod(spmat,sprow,spcol,m,colstat);
 
 		get_eg(&eg1); init_eg(&eg_f4); add_eg(&eg_f4,&eg0,&eg1);
 		if ( DP_Print ) {