===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.28
retrieving revision 1.40
diff -u -p -r1.28 -r1.40
--- OpenXM_contrib2/asir2000/builtin/array.c	2003/05/29 16:44:59	1.28
+++ OpenXM_contrib2/asir2000/builtin/array.c	2004/12/02 13:53:31	1.40
@@ -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.27 2003/01/06 01:16:37 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.39 2004/12/01 12:55:19 noro 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);
@@ -68,7 +70,7 @@ void Pgeneric_gauss_elim();
 void Pgeneric_gauss_elim_mod();
 
 void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
-void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol();
+void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
 void Pgeninv_sf_swap();
 void sepvect();
 void Pmulmat_gf2n();
@@ -86,6 +88,8 @@ void Pmat_swap_row_destructive();
 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},
@@ -100,13 +104,17 @@ struct ftab array_tab[] = {
 	{"newmat",Pnewmat,-3},
 	{"matrix",Pnewmat,-3},
 	{"mat",Pmat,-99999999},
+	{"matr",Pmat,-99999999},
+	{"matc",Pmatc,-99999999},
 	{"newbytearray",Pnewbytearray,-2},
 	{"sepmat_destructive",Psepmat_destructive,2},
 	{"sepvect",Psepvect,2},
 	{"qsort",Pqsort,-2},
 	{"vtol",Pvtol,1},
+	{"ltov",Pltov,1},
 	{"size",Psize,1},
 	{"det",Pdet,-2},
+	{"nd_det",Pnd_det,-2},
 	{"invmat",Pinvmat,-2},
 	{"leqm",Pleqm,2},
 	{"leqm1",Pleqm1,2},
@@ -153,12 +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 {
@@ -166,13 +190,25 @@ void Pqsort(NODE arg,VECT *rp)
 		if ( !p || OID(p)!=2 )
 			error("qsort : invalid argument");
 		v = VR(p);
-		if ( (int)v->attr != V_SR )
-			error("qsort : no such function");
-		generic_comp_obj_func = (FUNC)v->priv;
+		gen_searchf(NAME(v),&func);
+		if ( !func ) {
+			if ( (int)v->attr != V_SR )
+				error("qsort : no such function");
+			func = (FUNC)v->priv;
+		}
+		generic_comp_obj_func = func;
 		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)
@@ -371,6 +407,23 @@ void Pvect(NODE arg,VECT *rp) {
 	}
 
 	for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);
+	if ( len == 1 ) {
+		if ( ARG0(arg) != 0 ) {
+			switch ( OID(ARG0(arg)) ) {
+				case O_VECT:
+					*rp = ARG0(arg);
+					return;
+				case O_LIST:
+					for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );
+					MKVECT(vect,len-1);
+					for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);
+							tn; i++, tn = NEXT(tn) )
+						vb[i] = (pointer)BDY(tn);
+					*rp=vect;
+					return;
+			}
+		}
+	}
 	MKVECT(vect,len);
 	for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )
 		vb[i] = (pointer)BDY(tn);
@@ -462,9 +515,12 @@ void Pnewmat(NODE arg,MAT *rp)
 void Pmat(NODE arg, MAT *rp)
 {
 	int row,col;
+	int i;
 	MAT m;
 	pointer **mb;
+	pointer *ent;
 	NODE tn, sn;
+	VECT v;
 
 	if ( !arg ) {
 		*rp =0;
@@ -472,14 +528,92 @@ void Pmat(NODE arg, MAT *rp)
 	}
 
 	for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
-	for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
+	if ( row == 1 ) {
+		if ( OID(ARG0(arg)) == O_MAT ) {
+			*rp=ARG0(arg);
+			return;
+		} else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
+			error("mat : invalid argument");
+		}
+	}
+	if ( OID(ARG0(arg)) == O_VECT ) {
+		v = ARG0(arg);
+		col = v->len;
+	} else if ( OID(ARG0(arg)) == O_LIST ) {
+		for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
+	} else {
+		error("mat : invalid argument");
+	}
+
 	MKMAT(m,row,col);
-	for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++)
-		for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
-			mb[row][col] = (pointer)BDY(sn);
+	for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {
+		if ( BDY(tn) == 0 ) {
+			error("mat : invalid argument");
+		} else if ( OID(BDY(tn)) == O_VECT ) {
+			v = tn->body;
+			ent = BDY(v);
+			for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];
+		} else if ( OID(BDY(tn)) == O_LIST ) {
+			for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
+				mb[row][col] = (pointer)BDY(sn);
+		} else {
+			error("mat : invalid argument");
+		}
+	}
 	*rp = m;
 }
 
+void Pmatc(NODE arg, MAT *rp)
+{
+	int row,col;
+	int i;
+	MAT m;
+	pointer **mb;
+	pointer *ent;
+	NODE tn, sn;
+	VECT v;
+
+	if ( !arg ) {
+		*rp =0;
+		return;
+	}
+
+	for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);
+	if ( col == 1 ) {
+		if ( OID(ARG0(arg)) == O_MAT ) {
+			*rp=ARG0(arg);
+			return;
+		} else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
+			error("matc : invalid argument");
+		}
+	}
+	if ( OID(ARG0(arg)) == O_VECT ) {
+		v = ARG0(arg);
+		row = v->len;
+	} else if ( OID(ARG0(arg)) == O_LIST ) {
+		for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);
+	} else {
+		error("matc : invalid argument");
+	}
+
+	MKMAT(m,row,col);
+	for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {
+		if ( BDY(tn) == 0 ) {
+			error("matc : invalid argument");
+		} else if ( OID(BDY(tn)) == O_VECT ) {
+			v = tn->body;
+			ent = BDY(v);
+			for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];
+		} else if ( OID(BDY(tn)) == O_LIST ) {
+			for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )
+				mb[row][col] = (pointer)BDY(sn);
+		} else {
+			error("matc : invalid argument");
+		}
+	}
+	*rp = m;
+}
+
 void Pvtol(NODE arg,LIST *rp)
 {
 	NODE n,n1;
@@ -495,6 +629,21 @@ void Pvtol(NODE arg,LIST *rp)
 	MKLIST(*rp,n);
 }
 
+void Pltov(NODE arg,VECT *rp)
+{
+	NODE n;
+	VECT v;
+	int len,i;
+
+	asir_assert(ARG0(arg),O_LIST,"ltov");
+	n = (NODE)BDY((LIST)ARG0(arg));
+	len = length(n);
+	MKVECT(v,len);
+	for ( i = 0; i < len; i++, n = NEXT(n) )
+		BDY(v)[i] = BDY(n);
+	*rp = v;
+}
+
 void Premainder(NODE arg,Obj *rp)
 {
 	Obj a;
@@ -852,6 +1001,7 @@ int gauss_elim_mod(int **mat,int row,int col,int md)
 }
 
 struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
+struct oEGT eg_conv;
 
 int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
 {
@@ -970,7 +1120,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);
@@ -1022,6 +1172,7 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 	int *rind,*cind;
 	int count;
 	struct oEGT eg_mul,eg_inv,tmp0,tmp1;
+	int period;
 
 	a0 = (Q **)mat->body;
 	row = mat->row; col = mat->col;
@@ -1069,8 +1220,10 @@ 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);
+			period = F4_INTRAT_PERIOD;
 			for ( q = ONE, count = 0; ; count++ ) {
-				fprintf(stderr,".");
+				if ( DP_Print )
+					fprintf(stderr,".");
 				/* wc = -b mod md */
 				for ( i = 0; i < rank; i++ )
 					for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
@@ -1109,20 +1262,24 @@ 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) )
+					if ( 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) ) {
+							if ( DP_Print ) {
+								fprintf(stderr,"\n");
+								print_eg("INV",&eg_inv);
+								print_eg("MUL",&eg_mul);
+								fflush(asir_out);
+							}
+							return rank;
+						}
+					} else
+						period *=2;
 			}
 	}
 }
@@ -1474,6 +1631,14 @@ void red_by_vect(int m,unsigned int *p,unsigned int *r
 		}
 }
 
+void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
+{
+	*p++ = 0; r++; len--;
+	for ( ; len; len--, r++, p++ )
+		if ( *r )
+			*p = _addsf(_mulsf(*r,hc),*p);
+}
+
 extern unsigned int **psca;
 
 void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
@@ -1562,6 +1727,50 @@ int generic_gauss_elim_mod(int **mat0,int row,int col,
 	return rank;
 }
 
+int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
+{
+	int i,j,k,l,inv,a,rank;
+	unsigned int *t,*pivot,*pk;
+	unsigned int **mat;
+
+	mat = (unsigned int **)mat0;
+	for ( rank = 0, j = 0; j < col; j++ ) {
+		for ( i = rank; i < row; i++ )
+			if ( mat[i][j] )
+				break;
+		if ( i == row ) {
+			colstat[j] = 0;
+			continue;
+		} else
+			colstat[j] = 1;
+		if ( i != rank ) {
+			t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
+		}
+		pivot = mat[rank];
+		inv = _invsf(pivot[j]);
+		for ( k = j, pk = pivot+k; k < col; k++, pk++ )
+			if ( *pk )
+				*pk = _mulsf(*pk,inv);
+		for ( i = rank+1; i < row; i++ ) {
+			t = mat[i];
+			if ( a = t[j] )
+				red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
+		}
+		rank++;
+	}
+	for ( j = col-1, l = rank-1; j >= 0; j-- )
+		if ( colstat[j] ) {
+			pivot = mat[l];
+			for ( i = 0; i < l; i++ ) {
+				t = mat[i];
+				if ( a = t[j] )
+					red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
+			}
+			l--;
+		}
+	return rank;
+}
+
 /* LU decomposition; a[i][i] = 1/U[i][i] */
 
 int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
@@ -2715,4 +2924,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);
 }