===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.34
retrieving revision 1.44
diff -u -p -r1.34 -r1.44
--- OpenXM_contrib2/asir2000/builtin/array.c	2003/11/27 02:20:51	1.34
+++ OpenXM_contrib2/asir2000/builtin/array.c	2005/01/12 10:38:07	1.44
@@ -45,13 +45,15 @@
  * 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.33 2003/11/08 01:12:02 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.43 2004/12/18 16:50:10 saito Exp $
 */
 #include "ca.h"
 #include "base.h"
 #include "parse.h"
 #include "inline.h"
 
+#define F4_INTRAT_PERIOD 8 
+
 #if 0
 #undef DMAR
 #define DMAR(a1,a2,a3,d,r) (r)=dmar(a1,a2,a3,d);
@@ -87,6 +89,7 @@ void Pmat_swap_col_destructive();
 void Pvect();
 void Pmat();
 void Pmatc();
+void Pnd_det();
 
 struct ftab array_tab[] = {
 	{"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
@@ -111,6 +114,7 @@ struct ftab array_tab[] = {
 	{"ltov",Pltov,1},
 	{"size",Psize,1},
 	{"det",Pdet,-2},
+	{"nd_det",Pnd_det,-2},
 	{"invmat",Pinvmat,-2},
 	{"leqm",Pleqm,2},
 	{"leqm1",Pleqm1,2},
@@ -157,13 +161,28 @@ int generic_comp_obj(Obj *a,Obj *b)
 void Pqsort(NODE arg,VECT *rp)
 {
 	VECT vect;
-	NODE n;
+	NODE n,n1;
 	P p;
 	V v;
 	FUNC func;
+	int len,i;
+	pointer *a;
+	Obj t;
 
-	asir_assert(ARG0(arg),O_VECT,"qsort");
-	vect = (VECT)ARG0(arg);
+	t = ARG0(arg);
+    if (OID(t) == O_LIST) {
+        n = (NODE)BDY((LIST)t);
+        len = length(n);
+        MKVECT(vect,len);
+        for ( i = 0; i < len; i++, n = NEXT(n) ) {
+            BDY(vect)[i] = BDY(n);
+        }
+        
+    }else if (OID(t) != O_VECT) {
+        error("qsort : invalid argument");
+    }else {
+        vect = (VECT)t;
+    }
 	if ( argc(arg) == 1 )
 		qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
 	else {
@@ -181,7 +200,15 @@ void Pqsort(NODE arg,VECT *rp)
 		MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);	
 		qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
 	}
-	*rp = vect;
+    if (OID(t) == O_LIST) {
+        a = BDY(vect);
+        for ( i = len - 1, n = 0; i >= 0; i-- ) {
+            MKNODE(n1,a[i],n); n = n1;
+        }
+        MKLIST((LIST)*rp,n);
+    }else {
+        *rp = vect;
+    }
 }
 
 void PNBmul_gf2n(NODE arg,GF2N *rp)
@@ -733,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;
 		}
@@ -1093,7 +1124,7 @@ RESET:
 			add_eg(&eg_chrem_split,&tmp0,&tmp1);
 
 			get_eg(&tmp0);
-			if ( ind % 16 )
+			if ( ind % F4_INTRAT_PERIOD )
 				ret = 0;
 			else
 				ret = intmtoratm(crmat,m1,*nm,dn);
@@ -1144,7 +1175,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;
@@ -1192,9 +1229,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);
-			for ( q = ONE, count = 0; ; count++ ) {
-				fprintf(stderr,".");
+			init_eg(&eg_check); init_eg(&eg_intrat);
+			period = F4_INTRAT_PERIOD;
+			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] ) {
@@ -1204,17 +1249,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++ ) {
@@ -1232,18 +1279,50 @@ 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 % 16) && intmtoratm_q(xmat,NM(q),*nmmat,dn) ) {
-					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) ) {
-						fprintf(stderr,"\n");
-						print_eg("INV",&eg_inv);
-						print_eg("MUL",&eg_mul);
-						fflush(asir_out);
-						return rank;
+				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;
+						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 = 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;
 					}
 				}
 			}
@@ -1845,7 +1924,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;
@@ -1878,8 +1957,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];
 	}
 }
 
@@ -2890,4 +2973,12 @@ void printimat(int **mat,int row,int col)
 		}
 		printf("\n");
 	}
+}
+
+void Pnd_det(NODE arg,P *rp)
+{
+	if ( argc(arg) == 1 )
+		nd_det(0,ARG0(arg),rp);
+	else
+		nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
 }