===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.202
retrieving revision 1.207
diff -u -p -r1.202 -r1.207
--- OpenXM_contrib2/asir2000/engine/nd.c	2013/01/30 08:03:18	1.202
+++ OpenXM_contrib2/asir2000/engine/nd.c	2013/09/12 06:46:16	1.207
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.201 2013/01/30 07:44:48 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.206 2013/09/10 02:10:00 noro Exp $ */
 
 #include "nd.h"
 
@@ -74,6 +74,7 @@ void conv_ilist(int demand,int trace,NODE g,int **indp
 void parse_nd_option(NODE opt);
 void dltondl(int n,DL dl,UINT *r);
 DP ndvtodp(int mod,NDV p);
+DP ndtodp(int mod,ND p);
 
 extern int Denominator,DP_Multiple;
 
@@ -3953,20 +3954,31 @@ void nd_setup_parameters(int nvar,int max) {
         else if ( max < 65536 ) nd_bpe = 16;
         else nd_bpe = 32;
     }
-    if ( !do_weyl && weight_check && current_dl_weight_vector ) {
+    if ( !do_weyl && weight_check && (current_dl_weight_vector || nd_matrix) ) {
         UINT t;
-		/* t = max(weights) */
-        for ( i = 0, t = 0; i < nd_nvar; i++ )
-			if ( t < current_dl_weight_vector[i] )
-				t = current_dl_weight_vector[i];
-		/* i = bitsize of t */
-		for ( i = 0; t; t >>=1, i++ );
-		/* i += bitsize of nd_nvar */
-		for ( t = nd_nvar; t; t >>=1, i++);
-		/* nd_bpe+i = bitsize of max(weights)*max(exp)*nd_nvar */
-		if ( (nd_bpe+i) >= 31 )
-			error("nd_setup_parameters : too large weight");
-	}
+		int st;
+        int *v;
+	/* t = max(weights) */
+        t = 0;
+        if ( current_dl_weight_vector )
+            for ( i = 0, t = 0; i < nd_nvar; i++ ) {
+                if ( (st=current_dl_weight_vector[i]) < 0 ) st = -st;
+                if ( t < st ) t = st;
+            }
+        if ( nd_matrix )
+            for ( i = 0; i < nd_matrix_len; i++ )
+                for ( j = 0, v = nd_matrix[i]; j < nd_nvar; j++ ) {
+                    if ( (st=v[j]) < 0 ) st = -st;
+                    if ( t < st ) t = st;
+				}
+        /* i = bitsize of t */
+        for ( i = 0; t; t >>=1, i++ );
+        /* i += bitsize of nd_nvar */
+        for ( t = nd_nvar; t; t >>=1, i++);
+        /* nd_bpe+i = bitsize of max(weights)*max(exp)*nd_nvar */
+        if ( (nd_bpe+i) >= 31 )
+            error("nd_setup_parameters : too large weight");
+    }
     nd_epw = (sizeof(UINT)*8)/nd_bpe;
     elen = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
     nd_exporigin = nd_get_exporigin(nd_ord);
@@ -4878,6 +4890,27 @@ DP ndvtodp(int mod,NDV p)
     return d;
 }
 
+DP ndtodp(int mod,ND p)
+{
+    MP m,m0;
+	DP d;
+    NM t;
+    int i,len;
+
+    if ( !p ) return 0;
+    m0 = 0;
+    len = p->len;
+    for ( t = BDY(p); t; t = NEXT(t) ) {
+        NEXTMP(m0,m);
+    	m->dl = ndltodl(nd_nvar,DL(t));
+    	m->c = ndctop(mod,t->c);
+    }
+    NEXT(m) = 0;
+	MKDP(nd_nvar,m0,d);
+    SG(d) = SG(p);
+    return d;
+}
+
 void ndv_print(NDV p)
 {
     NMV m;
@@ -4946,6 +4979,8 @@ NODE ndv_reducebase(NODE x,int *perm)
 void nd_init_ord(struct order_spec *ord)
 {
     nd_module = (ord->id >= 256);
+	nd_matrix = 0;
+	nd_matrix_len = 0;
     switch ( ord->id ) {
         case 0:
             switch ( ord->ord.simple ) {
@@ -6644,6 +6679,32 @@ void ndv_save(NDV p,int index)
     fclose(s);
 }
 
+void nd_save_mod(ND p,int index)
+{
+    FILE *s;
+    char name[BUFSIZ];
+    int nv,sugar,len,c;
+    NM m;
+
+    sprintf(name,"%s/%d",Demand,index);
+    s = fopen(name,"w");
+    if ( !p ) { 
+		len = 0;
+    	write_int(s,&len);
+		fclose(s);
+        return;
+    }
+    nv = NV(p);
+    sugar = SG(p);
+    len = LEN(p);
+    write_int(s,&nv); write_int(s,&sugar); write_int(s,&len);
+	for ( m = BDY(p); m; m = NEXT(m) ) {
+	  c = CM(m); write_int(s,&c);
+	  write_intarray(s,DL(m),nd_wpd);
+	}
+    fclose(s);
+}
+
 NDV ndv_load(int index)
 {
     FILE *s;
@@ -6688,6 +6749,36 @@ NDV ndv_load(int index)
     return d;
 }
 
+ND nd_load_mod(int index)
+{
+    FILE *s;
+    char name[BUFSIZ];
+    int nv,sugar,len,i,c;
+	ND d;
+    NM m0,m;
+
+    sprintf(name,"%s/%d",Demand,index);
+    s = fopen(name,"r");
+	/* if the file does not exist, it means p[index]=0 */
+    if ( !s ) return 0;
+
+    read_int(s,&nv);
+	if ( !nv ) { fclose(s); return 0; }
+
+    read_int(s,&sugar);
+    read_int(s,&len);
+	for ( m0 = 0, i = 0; i < len; i++ ) {
+		NEXTNM(m0,m);
+		read_int(s,&c); CM(m) = c;
+		read_intarray(s,DL(m),nd_wpd);
+	}
+	NEXT(m) = 0;
+    MKND(nv,m0,len,d);
+    SG(d) = sugar;
+	fclose(s);
+    return d;
+}
+
 void nd_det(int mod,MAT f,P *rp)
 {
     VL fv,tv;
@@ -7182,4 +7273,281 @@ void parse_nd_option(NODE opt)
 		else if ( !strcmp(key,"intersect") )
             nd_intersect = value?1:0;
     }
+}
+
+ND mdptond(DP d);
+ND nd_mul_nm(int mod,NM m0,ND p);
+ND *btog(NODE ti,ND **p,int nb,int mod);
+ND btog_one(NODE ti,ND *p,int nb,int mod);
+MAT nd_btog(LIST f,LIST v,int m,struct order_spec *ord,LIST tlist,MAT *rp);
+VECT nd_btog_one(LIST f,LIST v,int m,struct order_spec *ord,LIST tlist,int pos,MAT *rp);
+
+/* d:monomial */
+ND mdptond(DP d)
+{
+  NM m;
+  ND r;
+
+  if ( OID(d) == 1 )
+	r = ptond(CO,CO,(P)d);
+  else {
+    NEWNM(m);
+    dltondl(NV(d),BDY(d)->dl,DL(m));
+    CQ(m) = (Q)BDY(d)->c;
+    NEXT(m) = 0;
+    MKND(NV(d),m,1,r);
+  }
+    return r;
+}
+
+ND nd_mul_nm(int mod,NM m0,ND p)
+{
+  UINT *d0;
+  int c0,c1,c;
+  NM tm,mr,mr0;
+  ND r;
+
+  if ( !p ) return 0;
+  d0 = DL(m0);
+  c0 = CM(m0);
+  mr0 = 0;
+  for ( tm = BDY(p); tm; tm = NEXT(tm) ) {
+    NEXTNM(mr0,mr);
+	c = CM(tm); DMAR(c0,c,0,mod,c1); CM(mr) = c1;
+	ndl_add(d0,DL(tm),DL(mr));
+  }
+  NEXT(mr) = 0;
+  MKND(NV(p),mr0,LEN(p),r);
+  return r;
+}
+
+ND *btog(NODE ti,ND **p,int nb,int mod)
+{
+  PGeoBucket *r;
+  int i,ci;
+  NODE t,s;
+  ND m,tp;
+  ND *pi,*rd;
+  P c;
+
+  r = (PGeoBucket *)MALLOC(nb*sizeof(PGeoBucket));
+  for ( i = 0; i < nb; i++ )
+  	r[i] = create_pbucket();
+  for ( t = ti; t; t = NEXT(t) ) {
+	s = BDY((LIST)BDY(t));
+    if ( ARG0(s) ) {
+	  m = mdptond((DP)ARG2(s));
+	  ptomp(mod,(P)HCQ(m),&c);
+	  if ( ci = ((MQ)c)->cont ) {
+	    HCM(m) = ci;
+	    pi = p[QTOS((Q)ARG1(s))];
+	    for ( i = 0; i < nb; i++ ) {
+		  tp = nd_mul_nm(mod,BDY(m),pi[i]);
+	      add_pbucket(mod,r[i],tp);
+	    }
+	  }
+	  ci = 1;
+    } else {
+	  ptomp(mod,(P)ARG3(s),&c); ci = ((MQ)c)->cont;
+	  ci = invm(ci,mod);
+	}
+  }
+  rd = (ND *)MALLOC(nb*sizeof(ND));
+  for ( i = 0; i < nb; i++ )
+	rd[i] = normalize_pbucket(mod,r[i]);
+  if ( ci != 1 )
+    for ( i = 0; i < nb; i++ ) nd_mul_c(mod,rd[i],ci);
+   return rd;
+}
+
+ND btog_one(NODE ti,ND *p,int nb,int mod)
+{
+  PGeoBucket r;
+  int i,ci,j;
+  NODE t,s;
+  ND m,tp;
+  ND pi,rd;
+  P c;
+
+  r = create_pbucket();
+  for ( t = ti; t; t = NEXT(t) ) {
+	s = BDY((LIST)BDY(t));
+    if ( ARG0(s) ) {
+	  m = mdptond((DP)ARG2(s));
+	  ptomp(mod,(P)HCQ(m),&c);
+	  if ( ci = ((MQ)c)->cont ) {
+	    HCM(m) = ci;
+	    pi = p[j=QTOS((Q)ARG1(s))];
+		if ( !pi ) {
+		  pi = nd_load_mod(j);
+		  tp = nd_mul_nm(mod,BDY(m),pi);
+		  nd_free(pi);
+	      add_pbucket(mod,r,tp);
+		} else {
+		  tp = nd_mul_nm(mod,BDY(m),pi);
+	      add_pbucket(mod,r,tp);
+	    }
+	  }
+	  ci = 1;
+    } else {
+	  ptomp(mod,(P)ARG3(s),&c); ci = ((MQ)c)->cont;
+	  ci = invm(ci,mod);
+	}
+  }
+  rd = normalize_pbucket(mod,r);
+  free_pbucket(r);
+  if ( ci != 1 ) nd_mul_c(mod,rd,ci);
+  return rd;
+}
+
+MAT nd_btog(LIST f,LIST v,int mod,struct order_spec *ord,LIST tlist,MAT *rp)
+{
+  int i,j,n,m,nb,pi0,pi1,nvar;
+  VL fv,tv,vv;
+  NODE permtrace,perm,trace,intred,ind,t,pi,ti;
+  ND **p;
+  ND *c;
+  ND u;
+  P inv;
+  MAT mat;
+
+  parse_nd_option(current_option);
+  get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);
+  for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
+  switch ( ord->id ) {
+    case 1:
+      if ( ord->nv != nvar )
+        error("nd_check : invalid order specification");
+      break;
+    default:
+      break;
+  }
+  nd_init_ord(ord);
+#if 0
+  nd_bpe = QTOS((Q)ARG7(BDY(tlist)));
+#else
+  nd_bpe = 32;
+#endif
+  nd_setup_parameters(nvar,0);
+  permtrace = BDY((LIST)ARG2(BDY(tlist))); 
+  intred = BDY((LIST)ARG3(BDY(tlist))); 
+  ind = BDY((LIST)ARG4(BDY(tlist)));
+  perm = BDY((LIST)BDY(permtrace)); trace =NEXT(permtrace);
+  for ( i = length(perm)-1, t = trace; t; t = NEXT(t) ) {
+  	j = QTOS((Q)BDY(BDY((LIST)BDY(t))));
+	if ( j > i ) i = j;
+  }
+  n = i+1;
+  nb = length(BDY(f));
+  p = (ND **)MALLOC(n*sizeof(ND *));
+  for ( t = perm, i = 0; t; t = NEXT(t), i++ ) {
+    pi = BDY((LIST)BDY(t)); 
+	pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));
+	p[pi0] = c = (ND *)MALLOC(nb*sizeof(ND));
+	ptomp(mod,(P)ARG2(pi),&inv);
+	u = ptond(CO,vv,(P)ONE);
+	HCM(u) = ((MQ)inv)->cont;
+	c[pi1] = u;
+  }
+  for ( t = trace,i=0; t; t = NEXT(t), i++ ) {
+	printf("%d ",i); fflush(stdout);
+    ti = BDY((LIST)BDY(t));
+    p[j=QTOS((Q)ARG0(ti))] = btog(BDY((LIST)ARG1(ti)),p,nb,mod);
+    if ( j == 441 )
+		printf("afo");
+  }
+  for ( t = intred, i=0; t; t = NEXT(t), i++ ) {
+	printf("%d ",i); fflush(stdout);
+    ti = BDY((LIST)BDY(t));
+    p[j=QTOS((Q)ARG0(ti))] = btog(BDY((LIST)ARG1(ti)),p,nb,mod);
+    if ( j == 441 )
+		printf("afo");
+  }
+  m = length(ind);
+  MKMAT(mat,nb,m);
+  for ( j = 0, t = ind; j < m; j++, t = NEXT(t) ) 
+    for ( i = 0, c = p[QTOS((Q)BDY(t))]; i < nb; i++ ) 
+		BDY(mat)[i][j] = ndtodp(mod,c[i]);
+  return mat;
+}
+
+VECT nd_btog_one(LIST f,LIST v,int mod,struct order_spec *ord,
+  LIST tlist,int pos,MAT *rp)
+{
+  int i,j,n,m,nb,pi0,pi1,nvar;
+  VL fv,tv,vv;
+  NODE permtrace,perm,trace,intred,ind,t,pi,ti;
+  ND *p;
+  ND *c;
+  ND u;
+  P inv;
+  VECT vect;
+
+  parse_nd_option(current_option);
+  get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&nd_vc);
+  for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
+  switch ( ord->id ) {
+    case 1:
+      if ( ord->nv != nvar )
+        error("nd_check : invalid order specification");
+      break;
+    default:
+      break;
+  }
+  nd_init_ord(ord);
+#if 0
+  nd_bpe = QTOS((Q)ARG7(BDY(tlist)));
+#else
+  nd_bpe = 32;
+#endif
+  nd_setup_parameters(nvar,0);
+  permtrace = BDY((LIST)ARG2(BDY(tlist))); 
+  intred = BDY((LIST)ARG3(BDY(tlist))); 
+  ind = BDY((LIST)ARG4(BDY(tlist)));
+  perm = BDY((LIST)BDY(permtrace)); trace =NEXT(permtrace);
+  for ( i = length(perm)-1, t = trace; t; t = NEXT(t) ) {
+  	j = QTOS((Q)BDY(BDY((LIST)BDY(t))));
+	if ( j > i ) i = j;
+  }
+  n = i+1;
+  nb = length(BDY(f));
+  p = (ND *)MALLOC(n*sizeof(ND *));
+  for ( t = perm, i = 0; t; t = NEXT(t), i++ ) {
+    pi = BDY((LIST)BDY(t)); 
+	pi0 = QTOS((Q)ARG0(pi)); pi1 = QTOS((Q)ARG1(pi));
+	if ( pi1 == pos ) {
+	  ptomp(mod,(P)ARG2(pi),&inv);
+	  u = ptond(CO,vv,(P)ONE);
+	  HCM(u) = ((MQ)inv)->cont;
+	  p[pi0] = u;
+    }
+  }
+  for ( t = trace,i=0; t; t = NEXT(t), i++ ) {
+	printf("%d ",i); fflush(stdout);
+    ti = BDY((LIST)BDY(t));
+    p[j=QTOS((Q)ARG0(ti))] = btog_one(BDY((LIST)ARG1(ti)),p,nb,mod);
+    if ( Demand ) {
+        nd_save_mod(p[j],j); nd_free(p[j]); p[j] = 0;
+	}
+  }
+  for ( t = intred, i=0; t; t = NEXT(t), i++ ) {
+	printf("%d ",i); fflush(stdout);
+    ti = BDY((LIST)BDY(t));
+    p[j=QTOS((Q)ARG0(ti))] = btog_one(BDY((LIST)ARG1(ti)),p,nb,mod);
+    if ( Demand ) {
+        nd_save_mod(p[j],j); nd_free(p[j]); p[j] = 0;
+	}
+  }
+  m = length(ind);
+  MKVECT(vect,m);
+  for ( j = 0, t = ind; j < m; j++, t = NEXT(t) ) {
+	u = p[QTOS((Q)BDY(t))];
+	if ( !u ) {
+	  u = nd_load_mod(QTOS((Q)BDY(t)));
+	  BDY(vect)[j] = ndtodp(mod,u);
+	  nd_free(u);
+	} else
+	  BDY(vect)[j] = ndtodp(mod,u);
+  }
+  return vect;
 }