===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.40
retrieving revision 1.45
diff -u -p -r1.40 -r1.45
--- OpenXM_contrib2/asir2000/builtin/array.c	2004/12/02 13:53:31	1.40
+++ OpenXM_contrib2/asir2000/builtin/array.c	2005/01/23 14:03:47	1.45
@@ -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.39 2004/12/01 12:55:19 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.44 2005/01/12 10:38:07 noro Exp $
 */
 #include "ca.h"
 #include "base.h"
@@ -760,6 +760,10 @@ void Psize(NODE arg,LIST *rp)
 				n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
 				STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
 				break;
+			case O_IMAT:
+				n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
+				STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
+				break;
 			default:
 				error("size : invalid argument"); break;
 		}
@@ -828,11 +832,12 @@ 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)
@@ -1171,8 +1176,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 	int *cinfo,*rinfo;
 	int *rind,*cind;
 	int count;
-	struct oEGT eg_mul,eg_inv,tmp0,tmp1;
+	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;
@@ -1220,11 +1230,17 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 			*cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
 
 			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 )
-					fprintf(stderr,".");
+			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 > 3 )
+					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] ) {
@@ -1234,17 +1250,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++ ) {
@@ -1262,24 +1280,52 @@ 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 ( intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
+				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 ) {
 						for ( j = k = l = 0; j < col; j++ )
 							if ( cinfo[j] )
 								rind[k++] = j;	
 							else
 								cind[l++] = j;
-						if ( gensolve_check(mat,*nmmat,*dn,rind,cind) ) {
-							if ( DP_Print ) {
+						get_eg(&tmp0);
+						ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
+						get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
+						if ( ret ) {
+							if ( DP_Print > 3 ) {
 								fprintf(stderr,"\n");
 								print_eg("INV",&eg_inv);
 								print_eg("MUL",&eg_mul);
+								print_eg("INTRAT",&eg_intrat);
+								print_eg("CHECK",&eg_check);
 								fflush(asir_out);
 							}
 							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;
+					}
+				}
 			}
 	}
 }
@@ -1879,7 +1925,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;
@@ -1912,8 +1958,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];
 	}
 }