===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.10
retrieving revision 1.44
diff -u -p -r1.10 -r1.44
--- OpenXM_contrib2/asir2000/builtin/array.c	2000/11/13 01:48:12	1.10
+++ OpenXM_contrib2/asir2000/builtin/array.c	2005/01/12 10:38:07	1.44
@@ -45,40 +45,33 @@
  * 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.9 2000/11/08 08:02:49 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);
 #endif
 
-extern int Print; /* XXX */
+extern int DP_Print; /* XXX */
 
-void inner_product_mat_int_mod(Q **,int **,int,int,int,Q *);
-void solve_by_lu_mod(int **,int,int,int **,int);
-void solve_by_lu_gfmmat(GFMMAT,unsigned int,unsigned int *,unsigned int *);
-int lu_gfmmat(GFMMAT,unsigned int,int *);
-void mat_to_gfmmat(MAT,unsigned int,GFMMAT *);
 
-int generic_gauss_elim_mod(int **,int,int,int,int *);
-int generic_gauss_elim(MAT ,MAT *,Q *,int **,int **);
-
-int gauss_elim_mod(int **,int,int,int);
-int gauss_elim_mod1(int **,int,int,int);
-int gauss_elim_geninv_mod(unsigned int **,int,int,int);
-int gauss_elim_geninv_mod_swap(unsigned int **,int,int,unsigned int,unsigned int ***,int **);
 void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
+void Pinvmat();
 void Pnewbytearray();
 
+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();
 void Pbconvmat_gf2n();
@@ -90,25 +83,44 @@ void Px962_irredpoly_up2();
 void Pirredpoly_up2();
 void Pnbpoly_up2();
 void Pqsort();
+void Pexponent_vector();
+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},
 	{"lu_gfmmat",Plu_gfmmat,2},
 	{"mat_to_gfmmat",Pmat_to_gfmmat,2},
+	{"generic_gauss_elim",Pgeneric_gauss_elim,1},
 	{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
 	{"newvect",Pnewvect,-2},
+	{"vect",Pvect,-99999999},
+	{"vector",Pnewvect,-2},
+	{"exponent_vector",Pexponent_vector,-99999999},
 	{"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},
 	{"geninvm",Pgeninvm,2},
 	{"geninvm_swap",Pgeninvm_swap,2},
+	{"geninv_sf_swap",Pgeninv_sf_swap,1},
 	{"remainder",Premainder,2},
 	{"sremainder",Psremainder,2},
 	{"mulmat_gf2n",Pmulmat_gf2n,1},
@@ -119,11 +131,12 @@ struct ftab array_tab[] = {
 	{"x962_irredpoly_up2",Px962_irredpoly_up2,2},
 	{"irredpoly_up2",Pirredpoly_up2,2},
 	{"nbpoly_up2",Pnbpoly_up2,2},
+	{"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
+	{"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
 	{0,0,0},
 };
 
-int comp_obj(a,b)
-Obj *a,*b;
+int comp_obj(Obj *a,Obj *b)
 {
 	return arf_comp(CO,*a,*b);
 }
@@ -131,8 +144,7 @@ Obj *a,*b;
 static FUNC generic_comp_obj_func;
 static NODE generic_comp_obj_arg;
 
-int generic_comp_obj(a,b)
-Obj *a,*b;
+int generic_comp_obj(Obj *a,Obj *b)
 {
 	Q r;
 	
@@ -146,19 +158,31 @@ Obj *a,*b;
 }
 
 
-void Pqsort(arg,rp)
-NODE arg;
-VECT *rp;
+void Pqsort(NODE arg,VECT *rp)
 {
 	VECT vect;
-	char buf[BUFSIZ];
-	char *fname;
-	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,18 +190,28 @@ 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(arg,rp)
-NODE arg;
-GF2N *rp;
+void PNBmul_gf2n(NODE arg,GF2N *rp)
 {
 	GF2N a,b;
 	GF2MAT mat;
@@ -214,9 +248,7 @@ GF2N *rp;
 	}
 }
 
-void Pmul_vect_mat_gf2n(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
 {
 	GF2N a;
 	GF2MAT mat;
@@ -247,9 +279,7 @@ GF2N *rp;
 	}
 }
 
-void Pbconvmat_gf2n(arg,rp)
-NODE arg;
-LIST *rp;
+void Pbconvmat_gf2n(NODE arg,LIST *rp)
 {
 	P p0,p1;
 	int to;
@@ -269,9 +299,7 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Pmulmat_gf2n(arg,rp)
-NODE arg;
-GF2MAT *rp;
+void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
 {
 	GF2MAT m;
 
@@ -280,9 +308,7 @@ GF2MAT *rp;
 	*rp = m;
 }
 
-void Psepmat_destructive(arg,rp)
-NODE arg;
-LIST *rp;
+void Psepmat_destructive(NODE arg,LIST *rp)
 {
 	MAT mat,mat1;
 	int i,j,row,col;
@@ -313,17 +339,12 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Psepvect(arg,rp)
-NODE arg;
-VECT *rp;
+void Psepvect(NODE arg,VECT *rp)
 {
 	sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
 }
 
-void sepvect(v,d,rp)
-VECT v;
-int d;
-VECT *rp;
+void sepvect(VECT v,int d,VECT *rp)
 {
 	int i,j,k,n,q,q1,r;
 	pointer *pv,*pw,*pu;
@@ -347,9 +368,7 @@ VECT *rp;
 	}
 }
 
-void Pnewvect(arg,rp)
-NODE arg;
-VECT *rp;
+void Pnewvect(NODE arg,VECT *rp)
 {
 	int len,i,r;
 	VECT vect;
@@ -376,10 +395,48 @@ VECT *rp;
 	*rp = vect;
 }
 
-void Pnewbytearray(arg,rp)
-NODE arg;
-BYTEARRAY *rp;
+void Pvect(NODE arg,VECT *rp) {
+	int len,i,r;
+	VECT vect;
+	pointer *vb;
+	NODE tn;
+
+	if ( !arg ) {
+		*rp =0;
+		return;
+	}
+
+	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);
+	*rp = vect;
+}
+
+void Pexponent_vector(NODE arg,DP *rp)
 {
+	nodetod(arg,rp);
+}
+
+void Pnewbytearray(NODE arg,BYTEARRAY *rp)
+{
 	int len,i,r;
 	BYTEARRAY array;
 	unsigned char *vb;
@@ -420,9 +477,7 @@ BYTEARRAY *rp;
 	*rp = array;
 }
 
-void Pnewmat(arg,rp)
-NODE arg;
-MAT *rp;
+void Pnewmat(NODE arg,MAT *rp)
 {
 	int row,col;
 	int i,j,r,c;
@@ -457,10 +512,110 @@ MAT *rp;
 	*rp = m;
 }
 
-void Pvtol(arg,rp)
-NODE arg;
-LIST *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;
+		return;
+	}
+
+	for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
+	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++) {
+		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;
 	VECT v;
 	pointer *a;
@@ -474,10 +629,23 @@ LIST *rp;
 	MKLIST(*rp,n);
 }
 
-void Premainder(arg,rp)
-NODE arg;
-Obj *rp;
+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;
 	VECT v,w;
 	MAT m,l;
@@ -524,9 +692,7 @@ Obj *rp;
 	}
 }
 
-void Psremainder(arg,rp)
-NODE arg;
-Obj *rp;
+void Psremainder(NODE arg,Obj *rp)
 {
 	Obj a;
 	VECT v,w;
@@ -575,9 +741,7 @@ Obj *rp;
 	}
 }
 
-void Psize(arg,rp)
-NODE arg;
-LIST *rp;
+void Psize(NODE arg,LIST *rp)
 {
 
 	int n,m;
@@ -596,6 +760,10 @@ 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;
 		}
@@ -603,9 +771,7 @@ LIST *rp;
 	MKLIST(*rp,t);
 }
 
-void Pdet(arg,rp)
-NODE arg;
-P *rp;
+void Pdet(NODE arg,P *rp)
 {
 	MAT m;
 	int n,i,j,mod;
@@ -629,6 +795,39 @@ P *rp;
 	}
 }
 
+void Pinvmat(NODE arg,LIST *rp)
+{
+	MAT m,r;
+	int n,i,j,mod;
+	P dn;
+	P **mat,**imat,**w;
+	NODE nd;
+
+	m = (MAT)ARG0(arg);
+	asir_assert(m,O_MAT,"invmat");
+	if ( m->row != m->col )
+		error("invmat : non-square matrix");
+	else if ( argc(arg) == 1 ) {
+		n = m->row;
+		invmatp(CO,(P **)BDY(m),n,&imat,&dn);
+		NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
+		nd = mknode(2,r,dn);
+		MKLIST(*rp,nd);
+	} else {
+		n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
+		w = (P **)almat_pointer(n,n);
+		for ( i = 0; i < n; i++ )
+			for ( j = 0; j < n; j++ )
+				ptomp(mod,mat[i][j],&w[i][j]);
+#if 0
+		detmp(CO,mod,w,n,&d);
+		mptop(d,rp);
+#else
+		error("not implemented yet");
+#endif
+	}
+}
+
 /*
 	input : a row x col matrix A
 		A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
@@ -640,11 +839,48 @@ P *rp;
 		B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */
 
-void Pgeneric_gauss_elim_mod(arg,rp)
-NODE arg;
-LIST *rp;
+void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 {
 	NODE n0;
+	MAT m,nm;
+	int *ri,*ci;
+	VECT rind,cind;
+	Q dn,q;
+	int i,j,k,l,row,col,t,rank;
+
+	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);
+	t = col-rank;
+	MKVECT(rind,rank);
+	MKVECT(cind,t);
+	for ( i = 0; i < rank; i++ ) {
+		STOQ(ri[i],q);
+		BDY(rind)[i] = (pointer)q;
+	}
+	for ( i = 0; i < t; i++ ) {
+		STOQ(ci[i],q);
+		BDY(cind)[i] = (pointer)q;
+	}
+	n0 = mknode(4,nm,dn,rind,cind);
+	MKLIST(*rp,n0);
+}
+
+/*
+	input : a row x col matrix A
+		A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
+
+	output : [B,R,C]
+		B : a rank(A) x col-rank(A) matrix
+		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]}+...
+*/
+
+void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
+{
+	NODE n0;
 	MAT m,mat;
 	VECT rind,cind;
 	Q **tmat;
@@ -652,7 +888,7 @@ LIST *rp;
 	Q *rib,*cib;
 	int *colstat;
 	Q q;
-	int md,i,j,k,l,row,col,t,n,rank;
+	int md,i,j,k,l,row,col,t,rank;
 
 	asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
 	asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
@@ -692,9 +928,7 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Pleqm(arg,rp)
-NODE arg;
-VECT *rp;
+void Pleqm(NODE arg,VECT *rp)
 {
 	MAT m;
 	VECT vect;
@@ -733,9 +967,7 @@ VECT *rp;
 	}
 }
 
-int gauss_elim_mod(mat,row,col,md)
-int **mat;
-int row,col,md;
+int gauss_elim_mod(int **mat,int row,int col,int md)
 {
 	int i,j,k,inv,a,n;
 	int *t,*pivot;
@@ -773,12 +1005,9 @@ int row,col,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,nm,dn,rindp,cindp)
-MAT mat;
-MAT *nm;
-Q *dn;
-int **rindp,**cindp;
+int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
 {
 	int **wmat;
 	Q **bmat;
@@ -805,10 +1034,10 @@ int **rindp,**cindp;
 	colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
 	wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
 	for ( ind = 0; ; ind++ ) {
-		if ( Print ) {
+		if ( DP_Print ) {
 			fprintf(asir_out,"."); fflush(asir_out);
 		}
-		md = lprime[ind];
+		md = get_lprime(ind);
 		get_eg(&tmp0);
 		for ( i = 0; i < row; i++ )
 			for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
@@ -842,13 +1071,13 @@ RESET:
 					}
 		} else {
 			if ( rank < rank0 ) {
-				if ( Print ) {
+				if ( DP_Print ) {
 					fprintf(asir_out,"lower rank matrix; continuing...\n");
 					fflush(asir_out);
 				}
 				continue;
 			} else if ( rank > rank0 ) {
-				if ( Print ) {
+				if ( DP_Print ) {
 					fprintf(asir_out,"higher rank matrix; resetting...\n");
 					fflush(asir_out);
 				}
@@ -856,7 +1085,7 @@ RESET:
 			} else {
 				for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
 				if ( j < col ) {
-					if ( Print ) {
+					if ( DP_Print ) {
 						fprintf(asir_out,"inconsitent colstat; resetting...\n");
 						fflush(asir_out);
 					}
@@ -895,7 +1124,10 @@ RESET:
 			add_eg(&eg_chrem_split,&tmp0,&tmp1);
 
 			get_eg(&tmp0);
-			ret = intmtoratm(crmat,m1,*nm,dn);
+			if ( ind % F4_INTRAT_PERIOD )
+				ret = 0;
+			else
+				ret = intmtoratm(crmat,m1,*nm,dn);
 			get_eg(&tmp1);
 			add_eg(&eg_intrat,&tmp0,&tmp1);
 			add_eg(&eg_intrat_split,&tmp0,&tmp1);
@@ -912,7 +1144,7 @@ RESET:
 					get_eg(&tmp1);
 					add_eg(&eg_gschk,&tmp0,&tmp1);
 					add_eg(&eg_gschk_split,&tmp0,&tmp1);
-					if ( Print ) {
+					if ( DP_Print ) {
 						print_eg("Mod",&eg_mod_split);
 						print_eg("Elim",&eg_elim_split);
 						print_eg("ChRem",&eg_chrem_split);
@@ -927,11 +1159,7 @@ RESET:
 	}
 }
 
-int generic_gauss_elim_hensel(mat,nmmat,dn,rindp,cindp)
-MAT mat;
-MAT *nmmat;
-Q *dn;
-int **rindp,**cindp;
+int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
 {
 	MAT bmat,xmat;
 	Q **a0,**a,**b,**x,**nm;
@@ -947,13 +1175,19 @@ int **rindp,**cindp;
 	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;
 	w = (int **)almat(row,col);
 	for ( ind = 0; ; ind++ ) {
-		md = lprime[ind];
+		md = get_lprime(ind);
 		STOQ(md,mdq);
 		for ( i = 0; i < row; i++ )
 			for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
@@ -965,7 +1199,7 @@ int **rindp,**cindp;
 				} else
 					wi[j] = 0;
 
-		rank = find_lhs_and_lu_mod(w,row,col,md,&rinfo,&cinfo);
+		rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
 		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++ )
@@ -995,9 +1229,17 @@ int **rindp,**cindp;
 			*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] ) {
@@ -1007,17 +1249,19 @@ int **rindp,**cindp;
 							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++ ) {
@@ -1035,18 +1279,50 @@ int **rindp,**cindp;
 				add_eg(&eg_mul,&tmp0,&tmp1);
 				/* q = q*md */
 				mulq(q,mdq,&u); q = u;
-				if ( !(count % 2) && 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;
 					}
 				}
 			}
@@ -1055,13 +1331,10 @@ int **rindp,**cindp;
 
 int f4_nocheck;
 
-int gensolve_check(mat,nm,dn,rind,cind)
-MAT mat,nm;
-Q dn;
-int *rind,*cind;
+int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
 {
 	int row,col,rank,clen,i,j,k,l;
-	Q s,t,u;
+	Q s,t;
 	Q *w;
 	Q *mati,*nmk;
 
@@ -1106,13 +1379,10 @@ int *rind,*cind;
 
 /* assuming 0 < c < m */
 
-int inttorat(c,m,b,sgnp,nmp,dnp)
-N c,m,b;
-int *sgnp;
-N *nmp,*dnp;
+int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
 {
-	Q qq,t,u1,v1,r1,nm;
-	N q,r,u2,v2,r2;
+	Q qq,t,u1,v1,r1;
+	N q,u2,v2,r2;
 
 	u1 = 0; v1 = ONE; u2 = m; v2 = c;
 	while ( cmpn(v2,b) >= 0 ) {
@@ -1131,14 +1401,10 @@ N *nmp,*dnp;
 
 /* mat->body = N ** */
 
-int intmtoratm(mat,md,nm,dn)
-MAT mat;
-N md;
-MAT nm;
-Q *dn;
+int intmtoratm(MAT mat,N md,MAT nm,Q *dn)
 {
 	N t,s,b;
-	Q bound,dn0,dn1,nm1,q,tq;
+	Q dn0,dn1,nm1,q;
 	int i,j,k,l,row,col;
 	Q **rmat;
 	N **tmat;
@@ -1188,14 +1454,10 @@ Q *dn;
 
 /* mat->body = Q ** */
 
-int intmtoratm_q(mat,md,nm,dn)
-MAT mat;
-N md;
-MAT nm;
-Q *dn;
+int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)
 {
 	N t,s,b;
-	Q bound,dn0,dn1,nm1,q,tq;
+	Q dn0,dn1,nm1,q;
 	int i,j,k,l,row,col;
 	Q **rmat;
 	Q **tmat;
@@ -1247,10 +1509,7 @@ Q *dn;
 
 #define ONE_STEP1  if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
 
-void reduce_reducers_mod(mat,row,col,md)
-int **mat;
-int row,col;
-int md;
+void reduce_reducers_mod(int **mat,int row,int col,int md)
 {
 	int i,j,k,l,hc,zzz;
 	int *t,*s,*tj,*ind;
@@ -1289,7 +1548,7 @@ int md;
 					ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
 					ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
 				}
-				for ( ; k >= 0; k-- ) {
+				for ( ; k > 0; k-- ) {
 					if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
 				}
 			}
@@ -1305,10 +1564,7 @@ int md;
 	2. reduce spolys by the reduced reducers
 */
 
-void pre_reduce_mod(mat,row,col,nred,md)
-int **mat;
-int row,col,nred;
-int md;
+void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
 {
 	int i,j,k,l,hc,inv;
 	int *t,*s,*tk,*ind;
@@ -1358,14 +1614,10 @@ int md;
 	mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
 */
 
-void reduce_sp_by_red_mod(sp,redmat,ind,nred,col,md)
-int *sp,**redmat;
-int *ind;
-int nred,col;
-int md;
+void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
 {
 	int i,j,k,hc,zzz;
-	int *t,*s,*tj;
+	int *s,*tj;
 
 	/* reduce the spolys by redmat */
 	for ( i = nred-1; i >= 0; i-- ) {
@@ -1376,25 +1628,105 @@ int md;
 			hc = md-hc;
 			s = redmat[i]+j;
 			tj = sp+j;
-			for ( k = col-j; k >= 0; k-- ) {
+			for ( k = col-j; k > 0; k-- ) {
 				if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
 			}
 		}
 	}
 }
 
+/*
+	mat[i] : compressed reducers (i=0,...,nred-1)
+	mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
+*/
+
+void red_by_compress(int m,unsigned int *p,unsigned int *r,
+	unsigned int *ri,unsigned int hc,int len)
+{
+	unsigned int up,lo;
+	unsigned int dmy;
+	unsigned int *pj;
+
+	p[*ri] = 0; r++; ri++;
+	for ( len--; len; len--, r++, ri++ ) {
+		pj = p+ *ri;
+		DMA(*r,hc,*pj,up,lo);
+		if ( up ) {
+			DSAB(m,up,lo,dmy,*pj);
+		} else
+			*pj = lo;
+	}
+}
+
+/* p -= hc*r */
+
+void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
+{
+	register unsigned int up,lo;
+	unsigned int dmy;
+
+	*p++ = 0; r++; len--;
+	for ( ; len; len--, r++, p++ )
+		if ( *r ) {
+			DMA(*r,hc,*p,up,lo);
+			if ( up ) {
+				DSAB(m,up,lo,dmy,*p);
+			} else
+				*p = lo;
+		}
+}
+
+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,
+	int nred,int col,int md)
+{
+	int i,len;
+	CDP ri;
+	unsigned int hc;
+	unsigned int *usp;
+
+	usp = (unsigned int *)sp;
+	/* reduce the spolys by redmat */
+	for ( i = nred-1; i >= 0; i-- ) {
+		/* reduce sp by redmat[i] */
+		usp[ind[i]] %= md;
+		if ( hc = usp[ind[i]] ) {
+			/* sp = sp-hc*redmat[i] */
+			hc = md-hc;
+			ri = redmat[i];
+			len = ri->len;
+			red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
+		}
+	}
+	for ( i = 0; i < col; i++ )
+		if ( usp[i] >= (unsigned int)md )
+			usp[i] %= md;
+}
+
 #define ONE_STEP2  if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
 
-int generic_gauss_elim_mod(mat,row,col,md,colstat)
-int **mat;
-int row,col,md;
-int *colstat;
+int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
 {
-	int i,j,k,l,inv,a,rank,zzz;
-	int *t,*pivot,*pk,*tk;
+	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 && !mat[i][j]; i++ );
+		for ( i = rank; i < row; i++ )
+			mat[i][j] %= md;
+		for ( i = rank; i < row; i++ )
+			if ( mat[i][j] )
+				break;
 		if ( i == row ) {
 			colstat[j] = 0;
 			continue;
@@ -1407,36 +1739,68 @@ int *colstat;
 		inv = invm(pivot[j],md);
 		for ( k = j, pk = pivot+k; k < col; k++, pk++ )
 			if ( *pk ) {
+				if ( *pk >= (unsigned int)md )
+					*pk %= md;
 				DMAR(*pk,inv,0,md,*pk)
 			}
 		for ( i = rank+1; i < row; i++ ) {
 			t = mat[i];
-			if ( a = t[j] ) {
-				a = md - a; pk = pivot+j; tk = t+j;
-				k = col-j;
-				for ( ; k >= 64; k -= 64 ) {
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-				}
-				for ( ; k >= 0; k -- ) {
-					if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
-				}
+			if ( a = t[j] )
+				red_by_vect(md,t+j,pivot+j,md-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];
+				t[j] %= md;
+				if ( a = t[j] )
+					red_by_vect(md,t+j,pivot+j,md-a,col-j);
 			}
+			l--;
 		}
+	for ( j = 0, l = 0; l < rank; j++ )
+		if ( colstat[j] ) {
+			t = mat[l];
+			for ( k = j; k < col; k++ )
+				if ( t[k] >= (unsigned int)md )
+					t[k] %= md;
+			l++;
+		}
+	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-- )
@@ -1444,31 +1808,8 @@ int *colstat;
 			pivot = mat[l];
 			for ( i = 0; i < l; i++ ) {
 				t = mat[i];
-				if ( a = t[j] ) {
-					a = md-a; pk = pivot+j; tk = t+j;
-					k = col-j;
-					for ( ; k >= 64; k -= 64 ) {
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					ONE_STEP2 ONE_STEP2 ONE_STEP2 ONE_STEP2
-					}
-					for ( ; k >= 0; k -- ) {
-						if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
-					}
-				}
+				if ( a = t[j] )
+					red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
 			}
 			l--;
 		}
@@ -1477,13 +1818,10 @@ int *colstat;
 
 /* LU decomposition; a[i][i] = 1/U[i][i] */
 
-int lu_gfmmat(mat,md,perm)
-GFMMAT mat;
-unsigned int md;
-int *perm;
+int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
 {
 	int row,col;
-	int i,j,k,l;
+	int i,j,k;
 	unsigned int *t,*pivot;
 	unsigned int **a;
 	unsigned int inv,m;
@@ -1535,12 +1873,10 @@ int *perm;
 	cinfo[j]=1 <=> j-th column is contained in the LU decomp.
 */
 
-int find_lhs_and_lu_mod(a,row,col,md,rinfo,cinfo)
-unsigned int **a;
-unsigned int md;
-int **rinfo,**cinfo;
+int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
+	unsigned int md,int **rinfo,int **cinfo)
 {
-	int i,j,k,l,d;
+	int i,j,k,d;
 	int *rp,*cp;
 	unsigned int *t,*pivot;
 	unsigned int inv,m;
@@ -1588,12 +1924,7 @@ int **rinfo,**cinfo;
 	b = a^(-1)b
  */
 
-void solve_by_lu_mod(a,n,md,b,l)
-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;
@@ -1626,14 +1957,16 @@ int l;
 			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];
 	}
 }
 
-void Pleqm1(arg,rp)
-NODE arg;
-VECT *rp;
+void Pleqm1(NODE arg,VECT *rp)
 {
 	MAT m;
 	VECT vect;
@@ -1672,9 +2005,7 @@ VECT *rp;
 	}
 }
 
-gauss_elim_mod1(mat,row,col,md)
-int **mat;
-int row,col,md;
+int gauss_elim_mod1(int **mat,int row,int col,int md)
 {
 	int i,j,k,inv,a,n;
 	int *t,*pivot;
@@ -1711,9 +2042,7 @@ int row,col,md;
 		return -1;
 }
 
-void Pgeninvm(arg,rp)
-NODE arg;
-LIST *rp;
+void Pgeninvm(NODE arg,LIST *rp)
 {
 	MAT m;
 	pointer **mat;
@@ -1747,17 +2076,15 @@ LIST *rp;
 		MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
 		for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
 			for ( j = 0; j < row; j++ )
-				STOQ(wmat[i][j+col],tmat[i][j]);
+				UTOQ(wmat[i][j+col],tmat[i][j]);
 		for ( tmat = (Q **)mat2->body; i < row; i++ )
 			for ( j = 0; j < row; j++ )
-				STOQ(wmat[i][j+col],tmat[i-col][j]);
+				UTOQ(wmat[i][j+col],tmat[i-col][j]);
 	 	MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
 	}
 }
 
-int gauss_elim_geninv_mod(mat,row,col,md)
-unsigned int **mat;
-int row,col,md;
+int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
 {
 	int i,j,k,inv,a,n,m;
 	unsigned int *t,*pivot;
@@ -1793,9 +2120,7 @@ int row,col,md;
 	return 0;
 }
 
-void Psolve_by_lu_gfmmat(arg,rp)
-NODE arg;
-VECT *rp;
+void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
 {
 	GFMMAT lu;
 	Q *perm,*rhs,*v;
@@ -1816,15 +2141,12 @@ VECT *rp;
 	solve_by_lu_gfmmat(lu,md,b,sol);
 	MKVECT(r,n);
 	for ( i = 0, v = (Q *)r->body; i < n; i++ )
-			STOQ(sol[i],v[i]);
+			UTOQ(sol[i],v[i]);
 	*rp = r;
 }
 
-void solve_by_lu_gfmmat(lu,md,b,x)
-GFMMAT lu;
-unsigned int md;
-unsigned int *b;
-unsigned int *x;
+void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
+	unsigned int *b,unsigned int *x)
 {
 	int n;
 	unsigned int **a;
@@ -1856,9 +2178,7 @@ unsigned int *x;
 	}
 }
 
-void Plu_gfmmat(arg,rp)
-NODE arg;
-LIST *rp;
+void Plu_gfmmat(NODE arg,LIST *rp)
 {
 	MAT m;
 	GFMMAT mm;
@@ -1888,9 +2208,7 @@ LIST *rp;
 	MKLIST(*rp,n0);
 }
 
-void Pmat_to_gfmmat(arg,rp)
-NODE arg;
-GFMMAT *rp;
+void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
 {
 	MAT m;
 	unsigned int md;
@@ -1901,10 +2219,7 @@ GFMMAT *rp;
 	mat_to_gfmmat(m,md,rp);	
 }
 
-void mat_to_gfmmat(m,md,rp)
-MAT m;
-unsigned int md;
-GFMMAT *rp;
+void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
 {
 	unsigned int **wmat;
 	unsigned int t;
@@ -2030,14 +2345,108 @@ int **indexp;
 	return 0;
 }
 
+void Pgeninv_sf_swap(NODE arg,LIST *rp)
+{
+	MAT m;
+	GFS **mat,**tmat;
+	Q *tvect;
+	GFS q;
+	int **wmat,**invmat;
+	int *index;
+	unsigned int t;
+	int i,j,row,col,status;
+	MAT mat1;
+	VECT vect1;
+	NODE node1,node2;
+
+	asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
+	m = (MAT)ARG0(arg);
+	row = m->row; col = m->col; mat = (GFS **)m->body;
+	wmat = (int **)almat(row,col+row);
+	for ( i = 0; i < row; i++ ) {
+		bzero((char *)wmat[i],(col+row)*sizeof(int));
+		for ( j = 0; j < col; j++ )
+			if ( q = (GFS)mat[i][j] )
+				wmat[i][j] = FTOIF(CONT(q));
+		wmat[i][col+i] = _onesf();
+	}
+	status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
+	if ( status > 0 )
+		*rp = 0;
+	else {
+		MKMAT(mat1,col,col);
+		for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
+			for ( j = 0; j < col; j++ )
+				if ( t = invmat[i][j] ) {
+					MKGFS(IFTOF(t),tmat[i][j]);
+				}
+		MKVECT(vect1,row);
+		for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
+			STOQ(index[i],tvect[i]);
+	 	MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
+	}
+}
+
+int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
+	int ***invmatp,int **indexp)
+{
+	int i,j,k,inv,a,n,m,u;
+	int *t,*pivot,*s;
+	int *index;
+	int **invmat;
+
+	n = col; m = row+col;
+	*indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
+	for ( i = 0; i < row; i++ )
+		index[i] = i;
+	for ( j = 0; j < n; j++ ) {
+		for ( i = j; i < row && !mat[i][j]; i++ );
+		if ( i == row ) {
+			*indexp = 0; *invmatp = 0; return 1;
+		}
+		if ( i != j ) {
+			t = mat[i]; mat[i] = mat[j]; mat[j] = t;
+			k = index[i]; index[i] = index[j]; index[j] = k;	
+		}
+		pivot = mat[j];
+		inv = _invsf(pivot[j]);
+		for ( k = j; k < m; k++ )
+			if ( pivot[k] )
+				pivot[k] = _mulsf(pivot[k],inv);
+		for ( i = j+1; i < row; i++ ) {
+			t = mat[i];
+			if ( a = t[j] )
+				for ( k = j, a = _chsgnsf(a); k < m; k++ )
+					if ( pivot[k] ) {
+						u = _mulsf(pivot[k],a);
+						t[k] = _addsf(u,t[k]);
+					}
+		}
+	}
+	for ( j = n-1; j >= 0; j-- ) {
+		pivot = mat[j];
+		for ( i = j-1; i >= 0; i-- ) {
+			t = mat[i];
+			if ( a = t[j] )
+				for ( k = j, a = _chsgnsf(a); k < m; k++ )
+					if ( pivot[k] ) {
+						u = _mulsf(pivot[k],a);
+						t[k] = _addsf(u,t[k]);
+					}
+		}
+	}
+	*invmatp = invmat = (int **)almat(col,col);
+	for ( i = 0; i < col; i++ )
+		for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
+			s[j] = t[col+index[j]];
+	return 0;
+}
+
 void _addn(N,N,N);
 int _subn(N,N,N);
 void _muln(N,N,N);
 
-void inner_product_int(a,b,n,r)
-Q *a,*b;
-int n;
-Q *r;
+void inner_product_int(Q *a,Q *b,int n,Q *r)
 {
 	int la,lb,i;
 	int sgn,sgn1;
@@ -2090,11 +2499,7 @@ Q *r;
 
 /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
 
-void inner_product_mat_int_mod(a,b,n,k,l,r)
-Q **a;
-int **b;
-int n,k,l;
-Q *r;
+void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
 {
 	int la,lb,i;
 	int sgn,sgn1;
@@ -2150,9 +2555,7 @@ Q *r;
 		NTOQ(sum,sgn,*r);
 }
 
-void Pmul_mat_vect_int(arg,rp)
-NODE arg;
-VECT *rp;
+void Pmul_mat_vect_int(NODE arg,VECT *rp)
 {
 	MAT mat;
 	VECT vect,r;
@@ -2163,14 +2566,13 @@ VECT *rp;
 	row = mat->row;
 	col = mat->col;
 	MKVECT(r,row);
-	for ( i = 0; i < row; i++ )
-		inner_product_int(mat->body[i],vect->body,col,&r->body[i]);
+	for ( i = 0; i < row; i++ ) {
+		inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
+	}
 	*rp = r;
 }
 
-void Pnbpoly_up2(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pnbpoly_up2(NODE arg,GF2N *rp)
 {
 	int m,type,ret;
 	UP2 r;
@@ -2184,11 +2586,9 @@ GF2N *rp;
 		*rp = 0;
 }
 
-void Px962_irredpoly_up2(arg,rp)
-NODE arg;
-GF2N *rp;
+void Px962_irredpoly_up2(NODE arg,GF2N *rp)
 {
-	int m,type,ret,w;
+	int m,ret,w;
 	GF2N prev;
 	UP2 r;
 
@@ -2204,18 +2604,16 @@ GF2N *rp;
 			bzero((char *)r->b,w*sizeof(unsigned int));
 		}
 	}
-	ret = _generate_irreducible_polynomial(r,m,type);
+	ret = _generate_irreducible_polynomial(r,m);
 	if ( ret == 0 )
 		MKGF2N(r,*rp);
 	else
 		*rp = 0;
 }
 
-void Pirredpoly_up2(arg,rp)
-NODE arg;
-GF2N *rp;
+void Pirredpoly_up2(NODE arg,GF2N *rp)
 {
-	int m,type,ret,w;
+	int m,ret,w;
 	GF2N prev;
 	UP2 r;
 
@@ -2231,13 +2629,55 @@ GF2N *rp;
 			bzero((char *)r->b,w*sizeof(unsigned int));
 		}
 	}
-	ret = _generate_good_irreducible_polynomial(r,m,type);
+	ret = _generate_good_irreducible_polynomial(r,m);
 	if ( ret == 0 )
 		MKGF2N(r,*rp);
 	else
 		*rp = 0;
 }
 
+void Pmat_swap_row_destructive(NODE arg, MAT *m)
+{
+	int i1,i2;
+	pointer *t;
+	MAT mat;
+
+	asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
+	asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
+	asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
+	mat = (MAT)ARG0(arg);
+	i1 = QTOS((Q)ARG1(arg));
+	i2 = QTOS((Q)ARG2(arg));
+	if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
+		error("mat_swap_row_destructive : Out of range");
+	t = mat->body[i1];
+	mat->body[i1] = mat->body[i2];
+	mat->body[i2] = t;
+	*m = mat;
+}
+
+void Pmat_swap_col_destructive(NODE arg, MAT *m)
+{
+	int j1,j2,i,n;
+	pointer *mi;
+	pointer t;
+	MAT mat;
+
+	asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
+	asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
+	asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
+	mat = (MAT)ARG0(arg);
+	j1 = QTOS((Q)ARG1(arg));
+	j2 = QTOS((Q)ARG2(arg));
+	if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
+		error("mat_swap_col_destructive : Out of range");
+	n = mat->row;
+	for ( i = 0; i < n; i++ ) {
+		mi = mat->body[i];
+		t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
+	}
+	*m = mat;
+}
 /*
  * f = type 'type' normal polynomial of degree m if exists
  * IEEE P1363 A.7.2
@@ -2511,9 +2951,7 @@ PENTA:
 	return 1;
 }
 
-printqmat(mat,row,col)
-Q **mat;
-int row,col;
+void printqmat(Q **mat,int row,int col)
 {
 	int i,j;
 
@@ -2525,9 +2963,7 @@ int row,col;
 	}
 }
 
-printimat(mat,row,col)
-int **mat;
-int row,col;
+void printimat(int **mat,int row,int col)
 {
 	int i,j;
 
@@ -2537,4 +2973,12 @@ int row,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);
 }