===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.218
retrieving revision 1.227
diff -u -p -r1.218 -r1.227
--- OpenXM_contrib2/asir2000/engine/nd.c	2014/02/24 01:45:28	1.218
+++ OpenXM_contrib2/asir2000/engine/nd.c	2016/07/11 08:00:30	1.227
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.217 2014/02/03 02:43:05 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.226 2016/04/05 04:21:18 noro Exp $ */
 
 #include "nd.h"
 
@@ -8,11 +8,14 @@ int diag_period = 6;
 int weight_check = 1;
 int (*ndl_compare_function)(UINT *a1,UINT *a2);
 int nd_dcomp;
+int nd_rref2;
 NM _nm_free_list;
 ND _nd_free_list;
 ND_pairs _ndp_free_list;
 NODE nd_hcf;
 
+Obj nd_top_weight;
+
 static NODE nd_subst;
 static VL nd_vc;
 static int nd_ntrans;
@@ -53,6 +56,8 @@ static int nd_found,nd_create,nd_notfirst;
 static int nmv_adv;
 static int nd_demand;
 static int nd_module,nd_ispot,nd_mpos,nd_pot_nelim;
+static int nd_module_rank,nd_poly_weight_len;
+static int *nd_poly_weight,*nd_module_weight;
 static NODE nd_tracelist;
 static NODE nd_alltracelist;
 static int nd_gentrace,nd_gensyz,nd_nora,nd_newelim,nd_intersect;
@@ -85,8 +90,19 @@ NDV ndvgztondv(NDV p);
 ND ndtondgz(ND p);
 ND ndgztond(ND p);
 
+void Pdp_set_weight(NODE,VECT *);
+void Pox_cmo_rpc(NODE,Obj *);
+
 extern int Denominator,DP_Multiple;
 
+#define BLEN (8*sizeof(unsigned long))
+
+typedef struct matrix {
+  int row,col;
+  unsigned long **a;
+} *matrix;
+
+
 void nd_free_private_storage()
 {
     _nm_free_list = 0;
@@ -107,6 +123,20 @@ void _NM_alloc()
     }
 }
 
+matrix alloc_matrix(int row,int col)
+{
+  unsigned long **a;
+  int i,len,blen;
+  matrix mat;
+
+  mat = (matrix)MALLOC(sizeof(struct matrix));
+  mat->row = row;
+  mat->col = col;
+  mat->a = a = (unsigned long **)MALLOC(row*sizeof(unsigned long *));
+  return mat;
+}
+
+
 void _ND_alloc()
 {
     ND p;
@@ -497,11 +527,37 @@ int ndl_block_compare(UINT *d1,UINT *d2)
 
 int ndl_matrix_compare(UINT *d1,UINT *d2)
 {
-    int i,j,s;
+    int i,j,s,row;
     int *v;
+	Q **mat;
+    Q *w;
+    Q t,t1,t2;
 
     for ( j = 0; j < nd_nvar; j++ )
         nd_work_vector[j] = GET_EXP(d1,j)-GET_EXP(d2,j);
+	if ( nd_top_weight ) {
+		if ( OID(nd_top_weight) == O_VECT ) {
+		    mat = (Q **)&BDY((VECT)nd_top_weight);
+		    row = 1;
+		} else {
+			mat = (Q **)BDY((MAT)nd_top_weight);
+		    row = ((MAT)nd_top_weight)->row;
+		}
+		for ( i = 0; i < row; i++ ) {
+		    w = (Q *)mat[i];
+			for ( j = 0, t = 0; j < nd_nvar; j++ ) {
+				STOQ(nd_work_vector[j],t1);
+			    mulq(w[j],t1,&t2);
+				addq(t,t2,&t1);
+				t = t1;
+			}
+		    if ( t ) {
+		    	s = SGN(t); 
+		        if ( s > 0 ) return 1;
+		        else if ( s < 0 ) return -1;
+		    }
+		}
+	}
     for ( i = 0; i < nd_matrix_len; i++ ) {
         v = nd_matrix[i];
         for ( j = 0, s = 0; j < nd_nvar; j++ )
@@ -509,6 +565,8 @@ int ndl_matrix_compare(UINT *d1,UINT *d2)
         if ( s > 0 ) return 1;
         else if ( s < 0 ) return -1;
     }
+	if ( !ndl_equal(d1,d2) )
+		error("afo");
     return 0;
 }
 
@@ -591,10 +649,28 @@ int ndl_ww_lex_compare(UINT *d1,UINT *d2)
     return ndl_lex_compare(d1,d2);
 }
 
+int ndl_module_weight_compare(UINT *d1,UINT *d2)
+{
+  int s,j;
+
+  if ( nd_nvar != nd_poly_weight_len )
+    error("invalid module weight : the length of polynomial weight != the number of variables");
+  s = 0;
+  for ( j = 0; j < nd_nvar; j++ )
+     s += (GET_EXP(d1,j)-GET_EXP(d2,j))*nd_poly_weight[j];
+  if ( MPOS(d1) >= 1 && MPOS(d2) >= 1 ) {
+    s += nd_module_weight[MPOS(d1)-1]-nd_module_weight[MPOS(d2)-1];
+  }
+  if ( s > 0 ) return 1;
+  else if ( s < 0 ) return -1;
+  else return 0;
+}
+
 int ndl_module_grlex_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
 		if ( nd_pot_nelim && MPOS(d1)>=nd_pot_nelim+1 && MPOS(d2) >= nd_pot_nelim+1 ) {
             if ( TD(d1) > TD(d2) ) return 1;
@@ -621,6 +697,7 @@ int ndl_module_glex_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;
@@ -639,6 +716,7 @@ int ndl_module_lex_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;
@@ -655,6 +733,7 @@ int ndl_module_block_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;
@@ -671,6 +750,7 @@ int ndl_module_matrix_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
         if ( MPOS(d1) < MPOS(d2) ) return 1;
         else if ( MPOS(d1) > MPOS(d2) ) return -1;
@@ -687,6 +767,7 @@ int ndl_module_composite_compare(UINT *d1,UINT *d2)
 {
     int i,c;
 
+    if ( nd_module_rank && (c = ndl_module_weight_compare(d1,d2)) ) return c;
     if ( nd_ispot ) {
         if ( MPOS(d1) > MPOS(d2) ) return 1;
         else if ( MPOS(d1) < MPOS(d2) ) return -1;
@@ -2279,9 +2360,11 @@ again:
         FREENDP(l);
     }
     if ( nd_nalg ) {
-        print_eg("monic",&eg_monic);
-        print_eg("invdalg",&eg_invdalg);
-        print_eg("le",&eg_le);
+        if ( DP_Print ) {
+          print_eg("monic",&eg_monic);
+          print_eg("invdalg",&eg_invdalg);
+          print_eg("le",&eg_le);
+        }
     }
 	conv_ilist(nd_demand,1,g,indp);
     if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); }
@@ -5193,7 +5276,16 @@ NODE ndv_reducebase(NODE x,int *perm)
 
 void nd_init_ord(struct order_spec *ord)
 {
-    nd_module = (ord->id >= 256);
+  nd_module = (ord->id >= 256);
+  if ( nd_module ) {
+    nd_dcomp = -1;
+    nd_ispot = ord->ispot;
+    nd_pot_nelim = ord->pot_nelim;
+    nd_poly_weight_len = ord->nv;
+    nd_poly_weight = ord->top_weight;
+    nd_module_rank = ord->module_rank;
+    nd_module_weight = ord->module_top_weight;
+  }
 	nd_matrix = 0;
 	nd_matrix_len = 0;
     switch ( ord->id ) {
@@ -5249,9 +5341,6 @@ void nd_init_ord(struct order_spec *ord)
 
         /* module order */
         case 256:
-            nd_ispot = ord->ispot;
-            nd_pot_nelim = ord->pot_nelim;
-            nd_dcomp = -1;
             switch ( ord->ord.simple ) {
                 case 0:
                     nd_isrlex = 1;
@@ -5271,17 +5360,11 @@ void nd_init_ord(struct order_spec *ord)
             break;
         case 257:
             /* block order */
-            nd_ispot = ord->ispot;
-            nd_pot_nelim = ord->pot_nelim;
-            nd_dcomp = -1;
             nd_isrlex = 0;
             ndl_compare_function = ndl_module_block_compare;
             break;
         case 258:
             /* matrix order */
-            nd_ispot = ord->ispot;
-            nd_pot_nelim = ord->pot_nelim;
-            nd_dcomp = -1;
             nd_isrlex = 0;
             nd_matrix_len = ord->ord.matrix.row;
             nd_matrix = ord->ord.matrix.matrix;
@@ -5289,9 +5372,6 @@ void nd_init_ord(struct order_spec *ord)
             break;
         case 259:
             /* composite order */
-            nd_ispot = ord->ispot;
-            nd_pot_nelim = ord->pot_nelim;
-            nd_dcomp = -1;
             nd_isrlex = 0;
             nd_worb_len = ord->ord.composite.length;
             nd_worb = ord->ord.composite.w_or_b;
@@ -5476,6 +5556,44 @@ int nd_to_vect_q(UINT *s0,int n,ND d,Q *r)
     return i;
 }
 
+unsigned long *nd_to_vect_2(UINT *s0,int n,int *s0hash,ND p)
+{
+    NM m;
+    unsigned long *v;
+    int i,j,h,size;
+	UINT *s,*t;
+
+	size = sizeof(unsigned long)*(n+BLEN-1)/BLEN;
+    v = (unsigned long *)MALLOC_ATOMIC_IGNORE_OFF_PAGE(size);
+    bzero(v,size);
+    for ( i = j = 0, s = s0, m = BDY(p); m; j++, m = NEXT(m) ) {
+		t = DL(m);
+		h = ndl_hash_value(t);
+        for ( ; h != s0hash[i] || !ndl_equal(t,s); s += nd_wpd, i++ );
+	    v[i/BLEN] |= 1L <<(i%BLEN);
+    }
+    return v;
+}
+
+int nd_nm_to_vect_2(UINT *s0,int n,int *s0hash,NDV p,NM m,unsigned long *v)
+{
+    NMV mr;
+    UINT *d,*t,*s;
+    int i,j,len,h,head;
+
+    d = DL(m);
+    len = LEN(p);
+    t = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+    for ( i = j = 0, s = s0, mr = BDY(p); j < len; j++, NMV_ADV(mr) ) {
+        ndl_add(d,DL(mr),t);    
+		h = ndl_hash_value(t);
+        for ( ; h != s0hash[i] || !ndl_equal(t,s); s += nd_wpd, i++ );
+		if ( j == 0 ) head = i;
+	    v[i/BLEN] |= 1L <<(i%BLEN);
+    }
+    return head;
+}
+
 Q *nm_ind_pair_to_vect(int mod,UINT *s0,int n,NM_ind_pair pair)
 {
     NM m;
@@ -5875,6 +5993,28 @@ NDV vect_to_ndv(UINT *vect,int spcol,int col,int *rhea
     }
 }
 
+NDV vect_to_ndv_2(unsigned long *vect,int col,UINT *s0vect)
+{
+    int j,k,len;
+    UINT *p;
+    NDV r;
+    NMV mr0,mr;
+
+    for ( j = 0, len = 0; j < col; j++ ) if ( vect[j/BLEN] & (1L<<(j%BLEN)) ) len++;
+    if ( !len ) return 0;
+    else {
+        mr0 = (NMV)MALLOC_ATOMIC_IGNORE_OFF_PAGE(nmv_adv*len);
+        mr = mr0; 
+        p = s0vect;
+        for ( j = 0; j < col; j++, p += nd_wpd )
+		  if ( vect[j/BLEN] & (1L<<(j%BLEN)) ) {
+            ndl_copy(p,DL(mr)); CM(mr) = 1; NMV_ADV(mr);
+          }
+        MKNDV(nd_nvar,mr0,len,r);
+        return r;
+    }
+}
+
 /* for preprocessed vector */
 
 NDV vect_to_ndv_q(Q *vect,int spcol,int col,int *rhead,UINT *s0vect)
@@ -6348,6 +6488,176 @@ NODE nd_f4_pseudo_trace(int m,int **indp)
     return g;
 }
 
+int rref(matrix mat,int *sugar)
+{
+  int row,col,i,j,k,l,s,wcol,wj;
+  unsigned long bj;
+  unsigned long **a;
+  unsigned long *ai,*ak,*as,*t;
+  int *pivot;
+
+  row = mat->row;
+  col = mat->col;
+  a = mat->a;
+  wcol = (col+BLEN-1)/BLEN; 
+  pivot = (int *)MALLOC_ATOMIC(row*sizeof(int));
+  i = 0;
+  for ( j = 0; j < col; j++ ) {
+	wj = j/BLEN; bj = 1L<<(j%BLEN);
+    for ( k = i; k < row; k++ )
+	  if ( a[k][wj] & bj ) break;
+    if ( k == row ) continue;
+	pivot[i] = j;
+    if ( k != i ) {
+	 t = a[i]; a[i] = a[k]; a[k] = t;
+	 s = sugar[i]; sugar[i] = sugar[k]; sugar[k] = s;
+	}
+	ai = a[i];
+    for ( k = i+1; k < row; k++ ) {
+	  ak = a[k];
+	  if ( ak[wj] & bj ) {
+	    for ( l = wj; l < wcol; l++ )
+		  ak[l] ^= ai[l];
+	    sugar[k] = MAX(sugar[k],sugar[i]);
+	  }
+	}
+	i++;
+  }
+  for ( k = i-1; k >= 0; k-- ) {
+    j = pivot[k]; wj = j/BLEN; bj = 1L<<(j%BLEN);
+	ak = a[k];
+    for ( s = 0; s < k; s++ ) {
+	  as = a[s];
+      if ( as[wj] & bj ) {
+        for ( l = wj; l < wcol; l++ )
+		  as[l] ^= ak[l];
+	    sugar[s] = MAX(sugar[s],sugar[k]);
+	  }
+	}
+  }
+  return i;
+}
+
+void print_matrix(matrix mat)
+{
+  int row,col,i,j;
+  unsigned long *ai;
+
+  row = mat->row;
+  col = mat->col;
+  printf("%d x %d\n",row,col);
+  for ( i = 0; i < row; i++ ) {
+	ai = mat->a[i];
+    for ( j = 0; j < col; j++ ) {
+	  if ( ai[j/BLEN] & (1L<<(j%BLEN)) ) putchar('1');
+	  else putchar('0');
+	}
+	putchar('\n');
+  }
+}
+
+NDV vect_to_ndv_2(unsigned long *vect,int col,UINT *s0vect);
+
+void red_by_vect_2(matrix mat,int *sugar,unsigned long *v,int rhead,int rsugar)
+{
+  int row,col,wcol,wj,i,j;
+  unsigned long bj;
+  unsigned long *ai;
+  unsigned long **a;
+  int len;
+  int *pos;
+
+  row = mat->row;
+  col = mat->col;
+  wcol = (col+BLEN-1)/BLEN; 
+  pos = (int *)ALLOCA(wcol*sizeof(int));
+  bzero(pos,wcol*sizeof(int));
+  for ( i = j = 0; i < wcol; i++ )
+    if ( v[i] ) pos[j++] = i;;
+  len = j;
+  wj = rhead/BLEN;
+  bj = 1L<<rhead%BLEN;
+  a = mat->a;
+  for ( i = 0; i < row; i++ ) {
+	ai = a[i];
+    if ( ai[wj]&bj ) {
+	  for ( j = 0; j < len; j++ )
+	    ai[pos[j]] ^= v[pos[j]];
+	  sugar[i] = MAX(sugar[i],rsugar); 
+	}
+  }
+}
+
+NODE nd_f4_red_2(ND_pairs sp0,UINT *s0vect,int col,NODE rp0,ND_pairs *nz)
+{
+    int nsp,nred,i,i0,k,rank,row;
+    NODE r0,rp;
+    ND_pairs sp;
+	ND spol;
+	NM_ind_pair rt;
+    int *s0hash;
+	UINT *s;
+	int *pivot,*sugar,*head;
+	matrix mat;
+    NM m;
+    NODE r;
+	struct oEGT eg0,eg1,eg2,eg_elim1,eg_elim2;
+	int rhead,rsugar,size;
+    unsigned long *v;
+
+    get_eg(&eg0);
+init_eg(&eg_search);
+    for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
+    nred = length(rp0);
+    mat = alloc_matrix(nsp,col);
+    s0hash = (int *)ALLOCA(col*sizeof(int));
+    for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd )
+        s0hash[i] = ndl_hash_value(s);
+
+	sugar = (int *)ALLOCA(nsp*sizeof(int));
+	for ( i = 0, sp = sp0; sp; sp = NEXT(sp) ) {
+		nd_sp(2,0,sp,&spol);
+		if ( spol ) {
+	      mat->a[i] = nd_to_vect_2(s0vect,col,s0hash,spol);
+		  sugar[i] = SG(spol);
+		  i++;
+		}
+	}
+	mat->row = i;
+    if ( DP_Print ) {
+      fprintf(asir_out,"%dx%d,",mat->row,mat->col); fflush(asir_out);
+    }
+	size = ((col+BLEN-1)/BLEN)*sizeof(unsigned long);
+	v = CALLOC((col+BLEN-1)/BLEN,sizeof(unsigned long));
+    for ( rp = rp0, i = 0; rp; rp = NEXT(rp), i++ ) {
+		rt = (NM_ind_pair)BDY(rp);
+		bzero(v,size);
+        rhead = nd_nm_to_vect_2(s0vect,col,s0hash,nd_ps[rt->index],rt->mul,v);
+		rsugar = SG(nd_ps[rt->index])+TD(DL(rt->mul));
+	    red_by_vect_2(mat,sugar,v,rhead,rsugar);
+	}
+
+    get_eg(&eg1);
+    init_eg(&eg_elim1); add_eg(&eg_elim1,&eg0,&eg1);
+	rank = rref(mat,sugar);
+
+    for ( i = 0, r0 = 0; i < rank; i++ ) {
+      NEXTNODE(r0,r); 
+	  BDY(r) = (pointer)vect_to_ndv_2(mat->a[i],col,s0vect);
+      SG((NDV)BDY(r)) = sugar[i];
+    }
+    if ( r0 ) NEXT(r) = 0;
+    get_eg(&eg2);
+    init_eg(&eg_elim2); add_eg(&eg_elim2,&eg1,&eg2);
+    if ( DP_Print ) {
+        fprintf(asir_out,"elim1=%fsec,elim2=%fsec\n",
+		  eg_elim1.exectime+eg_elim1.gctime,eg_elim2.exectime+eg_elim2.gctime);
+        fflush(asir_out);
+	}
+    return r0;
+}
+
+
 NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0vect,int col,NODE rp0,ND_pairs *nz)
 {
     IndArray *imat;
@@ -6359,6 +6669,9 @@ NODE nd_f4_red(int m,ND_pairs sp0,int trace,UINT *s0ve
     UINT *s;
     int *s0hash;
 
+    if ( m == 2 && nd_rref2 )
+	   return nd_f4_red_2(sp0,s0vect,col,rp0,nz);
+
 init_eg(&eg_search);
     for ( sp = sp0, nsp = 0; sp; sp = NEXT(sp), nsp++ );
     nred = length(rp0);
@@ -6367,7 +6680,9 @@ init_eg(&eg_search);
     for ( i = 0; i < col; i++ ) rhead[i] = 0;
 
     /* construction of index arrays */
-	fprintf(stderr,"%dx%d,",nsp+nred,col);
+    if ( DP_Print ) {
+	  fprintf(stderr,"%dx%d,",nsp+nred,col);
+    }
     rvect = (NM_ind_pair *)ALLOCA(nred*sizeof(NM_ind_pair));
     s0hash = (int *)ALLOCA(col*sizeof(int));
     for ( i = 0, s = s0vect; i < col; i++, s += nd_wpd )
@@ -6381,7 +6696,7 @@ init_eg(&eg_search);
         r0 = nd_f4_red_main(m,sp0,nsp,s0vect,col,rvect,rhead,imat,nred,nz);
     else
         r0 = nd_f4_red_gz_main(sp0,nsp,trace,s0vect,col,rvect,rhead,imat,nred);
-print_eg("search",&eg_search);
+    if ( DP_Print ) print_eg("search",&eg_search);
     return r0;
 }
 
@@ -7354,7 +7669,9 @@ void nd_det(int mod,MAT f,P *rp)
     if ( mod ) ndv_mod(mod,d);
     chsgnq(ONE,&mone);
     for ( j = 0, sgn = 1; j < n; j++ ) {
-        if ( DP_Print ) fprintf(stderr,".",j);
+        if ( DP_Print ) {
+	  fprintf(stderr,".",j);
+	}
         for ( i = j; i < n && !dm[i][j]; i++ );
         if ( i == n ) {
             *rp = 0;
@@ -7413,7 +7730,9 @@ void nd_det(int mod,MAT f,P *rp)
         }
         d = mjj;
     }
-    if ( DP_Print ) fprintf(stderr,"\n",k);
+    if ( DP_Print ) {
+      fprintf(stderr,"\n",k);
+    }
     if ( sgn < 0 )
         if ( mod )
             ndv_mul_c(mod,d,mod-1);