===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.43
retrieving revision 1.48
diff -u -p -r1.43 -r1.48
--- OpenXM_contrib2/asir2000/builtin/array.c	2004/12/18 16:50:10	1.43
+++ OpenXM_contrib2/asir2000/builtin/array.c	2005/11/27 05:37:53	1.48
@@ -45,7 +45,7 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.42 2004/12/13 23:04:16 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.47 2005/11/27 00:07:05 noro Exp $
 */
 #include "ca.h"
 #include "base.h"
@@ -158,7 +158,7 @@ int generic_comp_obj(Obj *a,Obj *b)
 }
 
 
-void Pqsort(NODE arg,VECT *rp)
+void Pqsort(NODE arg,LIST *rp)
 {
 	VECT vect;
 	NODE n,n1;
@@ -205,9 +205,9 @@ void Pqsort(NODE arg,VECT *rp)
         for ( i = len - 1, n = 0; i >= 0; i-- ) {
             MKNODE(n1,a[i],n); n = n1;
         }
-        MKLIST((LIST)*rp,n);
+        MKLIST(*rp,n);
     }else {
-        *rp = vect;
+        *rp = (LIST)vect;
     }
 }
 
@@ -832,26 +832,44 @@ void Pinvmat(NODE arg,LIST *rp)
 	input : a row x col matrix A
 		A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
 
-	output : [B,R,C]
+	output : [B,D,R,C]
 		B : a rank(A) x col-rank(A) matrix
+		D : the denominator
 		R : a vector of length rank(A)
 		C : a vector of length col-rank(A)
-		B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
+		B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */
 
 void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 {
-	NODE n0;
+	NODE n0,opt,p;
 	MAT m,nm;
 	int *ri,*ci;
 	VECT rind,cind;
 	Q dn,q;
 	int i,j,k,l,row,col,t,rank;
+	int is_hensel = 0;
+	char *key;
+	Obj value;
 
+	if ( current_option ) {
+		for ( opt = current_option; opt; opt = NEXT(opt) ) {
+			p = BDY((LIST)BDY(opt));
+			key = BDY((STRING)BDY(p));
+			value = (Obj)BDY(NEXT(p));
+			if ( !strcmp(key,"hensel") && value ) {
+				is_hensel = value ? 1 : 0;
+				break;
+			}
+		}
+	}
 	asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
 	m = (MAT)ARG0(arg);
 	row = m->row; col = m->col;
-	rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
+	if ( is_hensel )
+		rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
+	else
+		rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
 	t = col-rank;
 	MKVECT(rind,rank);
 	MKVECT(cind,t);
@@ -875,6 +893,8 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 		B : a rank(A) x col-rank(A) matrix
 		R : a vector of length rank(A)
 		C : a vector of length col-rank(A)
+		RN : a vector of length rank(A) indicating useful rows
+
 		B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */
 
@@ -882,11 +902,11 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 {
 	NODE n0;
 	MAT m,mat;
-	VECT rind,cind;
+	VECT rind,cind,rnum;
 	Q **tmat;
-	int **wmat;
-	Q *rib,*cib;
-	int *colstat;
+	int **wmat,**row0;
+	Q *rib,*cib,*rnb;
+	int *colstat,*p;
 	Q q;
 	int md,i,j,k,l,row,col,t,rank;
 
@@ -895,6 +915,10 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 	m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
 	row = m->row; col = m->col; tmat = (Q **)m->body;
 	wmat = (int **)almat(row,col);
+
+	row0 = (int **)ALLOCA(row*sizeof(int *));
+	for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
+
 	colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
 	for ( i = 0; i < row; i++ )
 		for ( j = 0; j < col; j++ )
@@ -907,6 +931,13 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 				wmat[i][j] = 0;
 	rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
 
+	MKVECT(rnum,rank);
+	rnb = (Q *)rnum->body;
+	for ( i = 0; i < rank; i++ )
+		for ( j = 0, p = wmat[i]; j < row; j++ )
+			if ( p == row0[j] )
+				STOQ(j,rnb[i]);
+
 	MKMAT(mat,rank,col-rank);
 	tmat = (Q **)mat->body;
 	for ( i = 0; i < rank; i++ )
@@ -924,7 +955,7 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 		} else {
 			STOQ(j,cib[l]); l++;
 		}
-	n0 = mknode(3,mat,rind,cind);
+	n0 = mknode(4,mat,rind,cind,rnum);
 	MKLIST(*rp,n0);
 }
 
@@ -1178,6 +1209,10 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 	int ret;
 	struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
 	int period;
+	int *wx,*ptr;
+	int wxsize,nsize;
+	N wn;
+	Q wq;
 
 	a0 = (Q **)mat->body;
 	row = mat->row; col = mat->col;
@@ -1195,7 +1230,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 				} else
 					wi[j] = 0;
 
+		if ( DP_Print ) {
+			fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
+		}
 		rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
+		if ( DP_Print ) {
+			fprintf(asir_out,"done.\n"); fflush(asir_out);
+		}
 		a = (Q **)almat_pointer(rank,rank); /* lhs mat */
 		MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
 		for ( j = li = ri = 0; j < col; j++ )
@@ -1227,10 +1268,15 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 			init_eg(&eg_mul); init_eg(&eg_inv);
 			init_eg(&eg_check); init_eg(&eg_intrat);
 			period = F4_INTRAT_PERIOD;
-			for ( q = ONE, count = 0; ; count++ ) {
-				if ( DP_Print > 3 )
+			nsize = period;
+			wxsize = rank*ri*nsize;
+			wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
+			for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
+			for ( q = ONE, count = 0; ; ) {
+				if ( DP_Print )
 					fprintf(stderr,"o");
 				/* wc = -b mod md */
+				get_eg(&tmp0);
 				for ( i = 0; i < rank; i++ )
 					for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
 						if ( u = (Q)bi[j] ) {
@@ -1240,17 +1286,19 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 							wi[j] = t;
 						} else
 							wi[j] = 0;
-				/* wc = A^(-1)wc; wc is normalized */
-				get_eg(&tmp0);
-				solve_by_lu_mod(w,rank,md,wc,ri);
-				get_eg(&tmp1);
-				add_eg(&eg_inv,&tmp0,&tmp1);
-				/* x = x-q*wc */
+				/* wc = A^(-1)wc; wc is not normalized */
+				solve_by_lu_mod(w,rank,md,wc,ri,0);
+				/* wx += q*wc */
+				ptr = wx;
 				for ( i = 0; i < rank; i++ )
-					for ( j = 0, xi = x[i], wi = wc[i]; j < ri; j++ ) {
-						STOQ(wi[j],u); mulq(q,u,&s);
-						subq(xi[j],s,&u); xi[j] = u;
+					for ( j = 0, wi = wc[i]; j < ri; j++ ) {
+						if ( wi[j] )
+							muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
+						ptr += nsize;
 					}
+				count++;
+				get_eg(&tmp1);
+				add_eg(&eg_inv,&tmp0,&tmp1);
 				get_eg(&tmp0);
 				for ( i = 0; i < rank; i++ )
 					for ( j = 0; j < ri; j++ ) {
@@ -1268,8 +1316,21 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 				add_eg(&eg_mul,&tmp0,&tmp1);
 				/* q = q*md */
 				mulq(q,mdq,&u); q = u;
-				if ( !(count % period) ) {
+				if ( count == period ) {
 					get_eg(&tmp0);
+					ptr = wx;
+					for ( i = 0; i < rank; i++ )
+						for ( j = 0, xi = x[i]; j < ri;
+							j++, ptr += nsize ) {
+							for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
+							if ( k >= 0 ) {
+								wn = NALLOC(k+1);
+								PL(wn) = k+1;
+								for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
+								NTOQ(wn,1,wq);
+								subq(xi[j],wq,&u); xi[j] = u;
+							}
+						}
 					ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
 					get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
 					if ( ret ) {
@@ -1292,8 +1353,14 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 							}
 							return rank;
 						}
-					} else
-						period *=2;
+					} else {
+						period = period*3/2;
+						count = 0;
+						nsize += period;
+						wxsize += rank*ri*nsize;
+						wx = (int *)REALLOC(wx,wxsize*sizeof(int));
+						for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
+					}
 				}
 			}
 	}
@@ -1894,7 +1961,7 @@ int find_lhs_and_lu_mod(unsigned int **a,int row,int c
 	b = a^(-1)b
  */
 
-void solve_by_lu_mod(int **a,int n,int md,int **b,int l)
+void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
 {
 	unsigned int *y,*c;
 	int i,j,k;
@@ -1927,8 +1994,12 @@ void solve_by_lu_mod(int **a,int n,int md,int **b,int 
 			DMAR(t,a[i][i],0,md,c[i])
 		}
 		/* copy c to b[.][k] with normalization */
-		for ( i = 0; i < n; i++ )
-			b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
+		if ( normalize )
+			for ( i = 0; i < n; i++ )
+				b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
+		else
+			for ( i = 0; i < n; i++ )
+				b[i][k] = c[i];
 	}
 }