===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.64
retrieving revision 1.71
diff -u -p -r1.64 -r1.71
--- OpenXM_contrib2/asir2000/builtin/array.c	2013/11/05 02:55:02	1.64
+++ OpenXM_contrib2/asir2000/builtin/array.c	2017/02/21 09:20:23	1.71
@@ -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.63 2013/09/09 07:29:25 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.70 2017/01/08 03:05:39 noro Exp $
 */
 #include "ca.h"
 #include "base.h"
@@ -68,13 +68,15 @@
 extern int DP_Print; /* XXX */
 
 
-void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
+void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm(), Ptriangleq();
 void Pinvmat();
 void Pnewbytearray(),Pmemoryplot_to_coord();
 
 void Pgeneric_gauss_elim();
 void Pgeneric_gauss_elim_mod();
 
+void Pindep_rows_mod();
+
 void Pmat_to_gfmmat(),Plu_gfmmat(),Psolve_by_lu_gfmmat();
 void Pgeninvm_swap(), Premainder(), Psremainder(), Pvtol(), Pltov();
 void Pgeninv_sf_swap();
@@ -108,6 +110,7 @@ struct ftab array_tab[] = {
 	{"mat_to_gfmmat",Pmat_to_gfmmat,2},
 	{"generic_gauss_elim",Pgeneric_gauss_elim,1},
 	{"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
+	{"indep_rows_mod",Pindep_rows_mod,2},
 	{"newvect",Pnewvect,-2},
 	{"vect",Pvect,-99999999},
 	{"vector",Pnewvect,-2},
@@ -148,6 +151,7 @@ struct ftab array_tab[] = {
 	{"mat_col",Pmat_col,2},
 	{"lusolve_prep",Plusolve_prep,1},
 	{"lusolve_main",Plusolve_main,1},
+	{"triangleq",Ptriangleq,1},
 	{0,0,0},
 };
 
@@ -895,6 +899,10 @@ void Pvtol(NODE arg,LIST *rp)
 	pointer *a;
 	int len,i;
 
+	if ( OID(ARG0(arg)) == O_LIST ) {
+		*rp = ARG0(arg);
+		return;
+	}
 	asir_assert(ARG0(arg),O_VECT,"vtol");
 	v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
 	for ( i = len - 1, n = 0; i >= 0; i-- ) {
@@ -906,9 +914,18 @@ void Pvtol(NODE arg,LIST *rp)
 void Pltov(NODE arg,VECT *rp)
 {
 	NODE n;
-	VECT v;
+	VECT v,v0;
 	int len,i;
 
+	if ( OID(ARG0(arg)) == O_VECT ) {
+		v0 = (VECT)ARG0(arg); len = v0->len;
+		MKVECT(v,len);
+		for ( i = 0; i < len; i++ ) {
+			BDY(v)[i] = BDY(v0)[i];
+		}
+		*rp = v;
+		return;
+	}
 	asir_assert(ARG0(arg),O_LIST,"ltov");
 	n = (NODE)BDY((LIST)ARG0(arg));
 	len = length(n);
@@ -1159,6 +1176,46 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 	MKLIST(*rp,n0);
 }
 
+void Pindep_rows_mod(NODE arg,VECT *rp)
+{
+	MAT m,mat;
+	VECT rind;
+	Q **tmat;
+	int **wmat,**row0;
+	Q *rib;
+	int *rowstat,*p;
+	Q q;
+	int md,i,j,k,l,row,col,t,rank;
+
+	asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
+	asir_assert(ARG1(arg),O_N,"indep_rows_mod");
+	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];
+
+	rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
+	for ( i = 0; i < row; i++ )
+		for ( j = 0; j < col; j++ )
+			if ( q = (Q)tmat[i][j] ) {
+				t = rem(NM(q),md);
+				if ( t && SGN(q) < 0 )
+					t = (md - t) % md;
+				wmat[i][j] = t;
+			} else
+				wmat[i][j] = 0;
+	rank = indep_rows_mod(wmat,row,col,md,rowstat);
+
+	MKVECT(rind,rank);
+	rib = (Q *)rind->body;
+	for ( j = 0; j < rank; j++ ) {
+		STOQ(rowstat[j],rib[j]);
+	}
+    *rp = rind;
+}
+
 /*
 	input : a row x col matrix A
 		A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
@@ -2303,6 +2360,22 @@ void red_by_vect_sf(int m,unsigned int *p,unsigned int
 			*p = _addsf(_mulsf(*r,hc),*p);
 }
 
+extern GZ current_mod_lf;
+extern int current_mod_lf_size;
+
+void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
+{
+	mpz_set_ui(*p++,0); r++; len--;
+	for ( ; len; len--, r++, p++ ) {
+       mpz_addmul(*p,*r,hc);
+#if 0
+       if ( mpz_size(*p) > current_mod_lf_size )
+         mpz_mod(*p,*p,BDY(current_mod_lf));
+#endif
+    }
+}
+
+
 extern unsigned int **psca;
 
 void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
@@ -2391,6 +2464,103 @@ int generic_gauss_elim_mod(int **mat0,int row,int col,
 	return rank;
 }
 
+int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
+{
+	int i,j,k,l,inv,a,rank;
+	unsigned int *t,*pivot,*pk;
+	unsigned int **mat;
+
+	for ( i = 0; i < row; i++ ) rowstat[i] = i;
+	mat = (unsigned int **)mat0;
+	for ( rank = 0, j = 0; j < col; j++ ) {
+		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;
+		} else
+			colstat[j] = 1;
+		if ( i != rank ) {
+			t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
+			k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
+		}
+		pivot = mat[rank];
+		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] )
+				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 indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
+{
+	int i,j,k,l,inv,a,rank;
+	unsigned int *t,*pivot,*pk;
+	unsigned int **mat;
+
+	for ( i = 0; i < row; i++ ) rowstat[i] = i;
+	mat = (unsigned int **)mat0;
+	for ( rank = 0, j = 0; j < col; j++ ) {
+		for ( i = rank; i < row; i++ )
+			mat[i][j] %= md;
+		for ( i = rank; i < row; i++ )
+			if ( mat[i][j] )
+				break;
+		if ( i == row ) continue;
+		if ( i != rank ) {
+			t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
+			k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
+		}
+		pivot = mat[rank];
+		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] )
+				red_by_vect(md,t+j,pivot+j,md-a,col-j);
+		}
+		rank++;
+	}
+	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;
@@ -3677,4 +3847,45 @@ void Pmat_col(NODE arg,VECT *rp)
 		BDY(vect)[i] = BDY(mat)[i][j];
 	}
 	*rp = vect;
+}
+
+NODE triangleq(NODE e)
+{
+  int n,i,k;
+  V v;
+  VL vl;
+  P *p;
+  NODE r,r1;
+
+  n = length(e);
+  p = (P *)MALLOC(n*sizeof(P));
+  for ( i = 0; i < n; i++, e = NEXT(e) ) p[i] = (P)BDY(e);
+  i = 0;
+  while ( 1 ) {
+    for ( ; i < n && !p[i]; i++ );
+    if ( i == n ) break;
+    if ( OID(p[i]) == O_N ) return 0;
+    v = p[i]->v;
+    for ( k = i+1; k < n; k++ )
+      if ( p[k] ) {
+        if ( OID(p[k]) == O_N ) return 0;
+        if ( p[k]->v == v ) p[k] = 0;
+      }
+    i++;
+  }
+  for ( r = 0, i = 0; i < n; i++ ) {
+    if ( p[i] ) {
+      MKNODE(r1,p[i],r); r = r1;
+    }
+  }
+  return r;
+}
+
+void Ptriangleq(NODE arg,LIST *rp)
+{
+  NODE ret;
+
+  asir_assert(ARG0(arg),O_LIST,"sparseleq");
+  ret = triangleq(BDY((LIST)ARG0(arg)));
+  MKLIST(*rp,ret);
 }