===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/distm.c,v
retrieving revision 1.11
retrieving revision 1.13
diff -u -p -r1.11 -r1.13
--- OpenXM_contrib2/asir2000/engine/distm.c	2003/07/18 10:13:13	1.11
+++ OpenXM_contrib2/asir2000/engine/distm.c	2003/07/21 01:57:51	1.13
@@ -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/engine/distm.c,v 1.10 2002/01/28 00:54:43 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.12 2003/07/20 08:55:23 noro Exp $ 
 */
 #include "ca.h"
 #include "inline.h"
@@ -925,6 +925,8 @@ static ND *nps;
 int nd_mod,nd_nvar;
 int is_rlex;
 int nd_epw,nd_bpe,nd_wpd;
+unsigned int nd_mask[32];
+unsigned int nd_mask0;
 NM _nm_free_list;
 ND _nd_free_list;
 ND_pairs _ndp_free_list;
@@ -972,14 +974,16 @@ void _ND_alloc();
 int ndl_td(unsigned int *d);
 ND nd_add(ND p1,ND p2);
 ND nd_mul_nm(ND p,NM m0);
+ND nd_mul_term(ND p,int td,unsigned int *d);
 ND nd_sp(ND_pairs p);
-ND nd_reducer(ND p1,ND p2);
-ND nd_nf(NODE b,ND g,ND *ps,int full);
+ND nd_find_reducer(ND g);
+ND nd_nf(ND g,int full);
 void ndl_print(unsigned int *dl);
 void nd_print(ND p);
 void ndp_print(ND_pairs d);
 int nd_length(ND p);
 void nd_monic(ND p);
+void nd_mul_c(ND p,int mul);
 
 void nd_free_private_storage()
 {
@@ -993,7 +997,7 @@ void _NM_alloc()
 	NM p;
 	int i;
 
-	for ( i = 0; i < 1024; i++ ) {
+	for ( i = 0; i < 10240; i++ ) {
 		p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int));
 		p->next = _nm_free_list; _nm_free_list = p;
 	}
@@ -1015,7 +1019,7 @@ void _NDP_alloc()
 	ND_pairs p;
 	int i;
 
-	for ( i = 0; i < 1024; i++ ) {
+	for ( i = 0; i < 10240; i++ ) {
 		p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs)
 			+(nd_wpd-1)*sizeof(unsigned int));
 		p->next = _ndp_free_list; _ndp_free_list = p;
@@ -1038,7 +1042,7 @@ INLINE nd_length(ND p)
 INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2)
 {
 	unsigned int u1,u2;
-	int i;
+	int i,j;
 
 	switch ( nd_bpe ) {
 		case 4:
@@ -1055,6 +1059,17 @@ INLINE int ndl_reducible(unsigned int *d1,unsigned int
 			}
 			return 1;
 			break;
+		case 6:
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i];
+				if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0;
+				if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0;
+				if ( (u1&0x3f000) < (u2&0x3f000) ) return 0;
+				if ( (u1&0xfc0) < (u2&0xfc0) ) return 0;
+				if ( (u1&0x3f) < (u2&0x3f) ) return 0;
+			}
+			return 1;
+			break;
 		case 8:
 			for ( i = 0; i < nd_wpd; i++ ) {
 				u1 = d1[i]; u2 = d2[i];
@@ -1079,14 +1094,19 @@ INLINE int ndl_reducible(unsigned int *d1,unsigned int
 			return 1;
 			break;
 		default:
-			error("ndl_reducible : not implemented yet");
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i];
+				for ( j = 0; j < nd_epw; j++ )
+					if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0;
+			}
+			return 1;
 	}
 }
 
-INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
+void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d)
 {
 	unsigned int t1,t2,u,u1,u2;
-	int i;
+	int i,j;
 
 	switch ( nd_bpe ) {
 		case 4:
@@ -1103,6 +1123,17 @@ INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsi
 				d[i] = u;
 			}
 			break;
+		case 6:
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i];
+				t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2;
+				t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2;
+				t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2;
+				t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2;
+				t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2;
+				d[i] = u;
+			}
+			break;
 		case 8:
 			for ( i = 0; i < nd_wpd; i++ ) {
 				u1 = d1[i]; u2 = d2[i];
@@ -1128,43 +1159,26 @@ INLINE void ndl_lcm(unsigned int *d1,unsigned *d2,unsi
 			}
 			break;
 		default:
-			error("ndl_lcm : not implemented yet");
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i]; 
+				for ( j = 0, u = 0; j < nd_epw; j++ ) {
+					t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2;
+				}
+				d[i] = u;
+			}
+			break;
 	}
 }
 
-INLINE int ndl_td(unsigned int *d)
+int ndl_td(unsigned int *d)
 {
 	unsigned int t,u;
-	int i;
+	int i,j;
 
-	switch ( nd_bpe ) {
-		case 4:
-			for ( t = 0, i = 0; i < nd_wpd; i++ ) {
-				u = d[i];
-				t += ((u&0xf0000000)>>28)+((u&0xf000000)>>24)
-					+((u&0xf00000)>>20)+((u&0xf0000)>>16)
-					+((u&0xf000)>>12)+((u&0xf00)>>8)+((u&0xf0)>>4)+(u&0xf);
-			}
-			break;
-		case 8:
-			for ( t = 0, i = 0; i < nd_wpd; i++ ) {
-				u = d[i];
-				t += ((u&0xff000000)>>24)+((u&0xff0000)>>16)
-					+((u&0xff00)>>8)+(u&0xff);
-			}
-			break;
-		case 16:
-			for ( t = 0, i = 0; i < nd_wpd; i++ ) {
-				u = d[i];
-				t += ((u&0xffff0000)>>16)+(u&0xffff);
-			}
-			break;
-		case 32:
-			for ( t = 0, i = 0; i < nd_wpd; i++ )
-				t += d[i];
-			break;
-		default:
-			error("ndl_td : not implemented yet");
+	for ( t = 0, i = 0; i < nd_wpd; i++ ) {
+		u = d[i];
+		for ( j = 0; j < nd_epw; j++, u>>=nd_bpe )
+			t += (u&nd_mask0); 
 	}
 	return t;
 }
@@ -1173,10 +1187,10 @@ INLINE int ndl_compare(unsigned int *d1,unsigned int *
 {
 	int i;
 
-	for ( i = 0; i < nd_wpd; i++ )
-		if ( d1[i] > d2[i] )
+	for ( i = 0; i < nd_wpd; i++, d1++, d2++ )
+		if ( *d1 > *d2 )
 			return is_rlex ? -1 : 1;
-		else if ( d1[i] < d2[i] )
+		else if ( *d1 < *d2 )
 			return is_rlex ? 1 : -1;
 	return 0;
 }
@@ -1199,7 +1213,7 @@ INLINE void ndl_add(unsigned int *d1,unsigned int *d2,
 		d[i] = d1[i]+d2[i];
 }
 
-INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
+void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d)
 {
 	int i;
 
@@ -1207,10 +1221,10 @@ INLINE void ndl_sub(unsigned int *d1,unsigned int *d2,
 		d[i] = d1[i]-d2[i];
 }
 
-INLINE int ndl_disjoint(unsigned int *d1,unsigned int *d2)
+int ndl_disjoint(unsigned int *d1,unsigned int *d2)
 {
 	unsigned int t1,t2,u,u1,u2;
-	int i;
+	int i,j;
 
 	switch ( nd_bpe ) {
 		case 4:
@@ -1227,6 +1241,17 @@ INLINE int ndl_disjoint(unsigned int *d1,unsigned int 
 			}
 			return 1;
 			break;
+		case 6:
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i];
+				t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0;
+				t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0;
+				t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0;
+				t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0;
+				t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0;
+			}
+			return 1;
+			break;
 		case 8:
 			for ( i = 0; i < nd_wpd; i++ ) {
 				u1 = d1[i]; u2 = d2[i];
@@ -1251,107 +1276,30 @@ INLINE int ndl_disjoint(unsigned int *d1,unsigned int 
 			return 1;
 			break;
 		default:
-			error("ndl_disjoint : not implemented yet");
-	}
-}
-
-ND nd_add(ND p1,ND p2)
-{
-	int n,c;
-	int t;
-	ND r;
-	NM m1,m2,mr0,mr,s;
-
-	if ( !p1 )
-		return p2;
-	else if ( !p2 )
-		return p1;
-	else {
-		for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
-			if ( m1->td > m2->td )
-				c = 1;
-			else if ( m1->td < m2->td )
-				c = -1;
-			else
-				c = ndl_compare(m1->dl,m2->dl);
-			switch ( c ) {
-				case 0:
-					t = ((C(m1))+(C(m2))) - nd_mod;
-					if ( t < 0 )
-						t += nd_mod;
-					s = m1; m1 = NEXT(m1);
-					if ( t ) {
-						NEXTNM2(mr0,mr,s); C(mr) = (t);
-					} else {
-						FREENM(s);
-					}
-					s = m2; m2 = NEXT(m2); FREENM(s);
-					break;
-				case 1:
-					s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
-					break;
-				case -1:
-					s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
-					break;
+			for ( i = 0; i < nd_wpd; i++ ) {
+				u1 = d1[i]; u2 = d2[i];
+				for ( j = 0; j < nd_epw; j++ ) {
+					if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0;
+					u1 >>= nd_bpe; u2 >>= nd_bpe;
+				}
 			}
-		}
-		if ( !mr0 )
-			if ( m1 )
-				mr0 = m1;
-			else if ( m2 )
-				mr0 = m2;
-			else
-				return 0;
-		else if ( m1 )
-			NEXT(mr) = m1;
-		else if ( m2 )
-			NEXT(mr) = m2;
-		else
-			NEXT(mr) = 0;
-		BDY(p1) = mr0;
-		p1->sugar = MAX(p1->sugar,p2->sugar);
-		FREEND(p2);
-		return p1;
+			return 1;
+			break;
 	}
 }
 
-INLINE ND nd_mul_nm(ND p,NM m0)
-{
-	NM m,mr,mr0;
-	unsigned int *d,*dt,*dm;
-	int c,n,td;
-	int *pt,*p1,*p2;
-	ND r;
-
-	if ( !p )
-		return 0;
-	else {
-		n = NV(p); m = BDY(p);
-		d = m0->dl; td = m0->td; c = C(m0);
-		mr0 = 0;
-		for ( ; m; m = NEXT(m) ) {
-			NEXTNM(mr0,mr);
-			C(mr) = (C(m)*c)%nd_mod;
-			mr->td = m->td+td;
-			ndl_add(m->dl,d,mr->dl);
-		}
-		NEXT(mr) = 0; 
-		MKND(NV(p),mr0,r);
-		r->sugar = p->sugar + td;
-		return r;
-	}
-}
-
 ND nd_reduce(ND p1,ND p2)
 {
-	int c,t,td,td2,mul;
+	int c,c1,c2,t,td,td2,mul;
 	NM m2,prev,head,cur,new;
 	unsigned int *d;
 
 	if ( !p1 )
 		return 0;
 	else {
-		mul = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod;
+		c2 = invm(HC(p2),nd_mod);
+		c1 = nd_mod-HC(p1);
+		DMAR(c1,c2,0,nd_mod,mul);
 		td = HTD(p1)-HTD(p2);
 		d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int));
 		ndl_sub(HDL(p1),HDL(p2),d);
@@ -1361,7 +1309,9 @@ ND nd_reduce(ND p1,ND p2)
 			td2 = new->td = m2->td+td;
 			ndl_add(m2->dl,d,new->dl);
 			if ( !cur ) {
-				C(new) = (C(m2)*mul)%nd_mod;
+				c1 = C(m2);
+				DMAR(c1,mul,0,nd_mod,c2);
+				C(new) = c2;
 				if ( !prev ) {
 					prev = new;
 					NEXT(prev) = 0;
@@ -1383,7 +1333,9 @@ ND nd_reduce(ND p1,ND p2)
 				c = ndl_compare(cur->dl,new->dl);
 			switch ( c ) {
 				case 0:
-					t = (C(cur)+C(m2)*mul)%nd_mod;
+					c2 = C(m2);
+					c1 = C(cur);
+					DMAR(c2,mul,c1,nd_mod,t);
 					if ( t )
 						C(cur) = t;
 					else if ( !prev ) {
@@ -1405,11 +1357,15 @@ ND nd_reduce(ND p1,ND p2)
 					if ( !prev ) {
 						/* cur = head */
 						prev = new;
-						C(prev) = (C(m2)*mul)%nd_mod;
+						c2 = C(m2);
+						DMAR(c2,mul,0,nd_mod,c1);
+						C(prev) = c1;
 						NEXT(prev) = head;
 						head = prev;
 					} else {
-						C(new) = (C(m2)*mul)%nd_mod;
+						c2 = C(m2);
+						DMAR(c2,mul,0,nd_mod,c1);
+						C(new) = c1;
 						NEXT(prev) = new;
 						NEXT(new) = cur;
 						prev = new;
@@ -1452,74 +1408,204 @@ ND nd_sp(ND_pairs p)
 	return nd_add(t1,t2);
 }
 
-ND nd_reducer(ND p1,ND p2)
+ND nd_find_reducer(ND g)
 {
 	NM m;
+	ND r,p;
+	int i,c1,c2,c;
+
+	for ( i = 0; i < nd_psn; i++ ) {
+		p = nps[i];
+		if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
+#if 1
+			NEWNM(m);
+			c1 = invm(HC(p),nd_mod);
+			c2 = nd_mod-HC(g);
+			DMAR(c1,c2,0,nd_mod,c);
+			C(m) = c;
+			m->td = HTD(g)-HTD(p);
+			ndl_sub(HDL(g),HDL(p),m->dl);
+			NEXT(m) = 0;
+			r = nd_mul_nm(p,m);
+			FREENM(m);
+			r->sugar = m->td + p->sugar;
+			return r;
+#else
+			return p;
+#endif
+		}
+	}
+	return 0;
+}
+
+ND nd_find_monic_reducer(ND g)
+{
+	int *d;
+	ND p,r;
+	int i;
+
+	for ( i = 0; i < nd_psn; i++ ) {
+		p = nps[i];
+		if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
+			d = (int *)ALLOCA(nd_wpd*sizeof(int));
+			ndl_sub(HDL(g),HDL(p),d);
+			r = nd_mul_term(p,HTD(g)-HTD(p),d);
+			return r;
+		}
+	}
+	return 0;
+}
+
+ND nd_add(ND p1,ND p2)
+{
+	int n,c;
+	int t;
 	ND r;
+	NM m1,m2,mr0,mr,s;
 
-	NEWNM(m);
-	C(m) = ((nd_mod-HC(p1))*invm(HC(p2),nd_mod))%nd_mod;
-	m->td = HTD(p1)-HTD(p2);
-	ndl_sub(HDL(p1),HDL(p2),m->dl);
-	NEXT(m) = 0;
-	r = nd_mul_nm(p2,m);
-	FREENM(m);
-	return r;
+	if ( !p1 )
+		return p2;
+	else if ( !p2 )
+		return p1;
+	else {
+		for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
+			if ( m1->td > m2->td )
+				c = 1;
+			else if ( m1->td < m2->td )
+				c = -1;
+			else
+				c = ndl_compare(m1->dl,m2->dl);
+			switch ( c ) {
+				case 0:
+					t = ((C(m1))+(C(m2))) - nd_mod;
+					if ( t < 0 )
+						t += nd_mod;
+					s = m1; m1 = NEXT(m1);
+					if ( t ) {
+						NEXTNM2(mr0,mr,s); C(mr) = (t);
+					} else {
+						FREENM(s);
+					}
+					s = m2; m2 = NEXT(m2); FREENM(s);
+					break;
+				case 1:
+					s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s);
+					break;
+				case -1:
+					s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s);
+					break;
+			}
+		}
+		if ( !mr0 )
+			if ( m1 )
+				mr0 = m1;
+			else if ( m2 )
+				mr0 = m2;
+			else
+				return 0;
+		else if ( m1 )
+			NEXT(mr) = m1;
+		else if ( m2 )
+			NEXT(mr) = m2;
+		else
+			NEXT(mr) = 0;
+		BDY(p1) = mr0;
+		p1->sugar = MAX(p1->sugar,p2->sugar);
+		FREEND(p2);
+		return p1;
+	}
 }
 
-#if 1
-ND nd_nf(NODE b,ND g,ND *ps,int full)
+ND nd_mul_nm(ND p,NM m0)
 {
-	ND u,p,d,red;
-	NODE l;
-	NM m,mrd;
-	int sugar,psugar,n,h_reducible;
+	NM m,mr,mr0;
+	unsigned int *d,*dt,*dm;
+	int c,n,td,i,c1,c2;
+	int *pt,*p1,*p2;
+	ND r;
 
-	if ( !g ) {
+	if ( !p )
 		return 0;
+	else {
+		n = NV(p); m = BDY(p);
+		d = m0->dl; td = m0->td; c = C(m0);
+		mr0 = 0;
+		for ( ; m; m = NEXT(m) ) {
+			NEXTNM(mr0,mr);
+			c1 = C(m);
+			DMAR(c1,c,0,nd_mod,c2);
+			C(mr) = c2;
+			mr->td = m->td+td;
+			ndl_add(m->dl,d,mr->dl);
+		}
+		NEXT(mr) = 0; 
+		MKND(NV(p),mr0,r);
+		r->sugar = p->sugar + td;
+		return r;
 	}
+}
+
+ND nd_mul_term(ND p,int td,unsigned int *d)
+{
+	NM m,mr,mr0;
+	int c,n;
+	ND r;
+
+	if ( !p )
+		return 0;
+	else {
+		n = NV(p); m = BDY(p);
+		mr0 = 0;
+		for ( ; m; m = NEXT(m) ) {
+			NEXTNM(mr0,mr);
+			C(mr) = C(m);
+			mr->td = m->td+td;
+			ndl_add(m->dl,d,mr->dl);
+		}
+		NEXT(mr) = 0; 
+		MKND(NV(p),mr0,r);
+		r->sugar = p->sugar + td;
+		return r;
+	}
+}
+
+#if 1
+ND nd_nf(ND g,int full)
+{
+	ND p,d,red;
+	NM m,mrd,tail;
+	int n,sugar,psugar;
+
+	if ( !g )
+		return 0;
 	sugar = g->sugar;
-	n = g->nv;
+	n = NV(g);
 	for ( d = 0; g; ) {
-		for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
-			p = ps[(int)BDY(l)];
-			if ( HTD(g)>=HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) {
-				h_reducible = 1;
-				psugar = HTD(g)-HTD(p) + p->sugar;
-#if 0
-				red = nd_reducer(g,p);
-				g = nd_add(g,red);
+		red = nd_find_reducer(g);
+		if ( red ) {
+#if 1
+			g = nd_add(g,red);
+			sugar = MAX(sugar,red->sugar);
 #else
-				g = nd_reduce(g,p);
+			psugar = (HTD(g)-HTD(red))+red->sugar;
+			g = nd_reduce(g,red);
+			sugar = MAX(sugar,psugar);
 #endif
-				sugar = MAX(sugar,psugar);
-				if ( !g ) {
-					if ( d )
-						d->sugar = sugar;
-					return d;
-				}
-				break;
+		} else if ( !full )
+			return g;
+		else {
+			m = BDY(g); 
+			if ( NEXT(m) ) {
+				BDY(g) = NEXT(m); NEXT(m) = 0;
+			} else {
+				FREEND(g); g = 0;
 			}
-		}
-		if ( !h_reducible ) {
-			/* head term is not reducible */
-			if ( !full ) {
-				if ( g )
-					g->sugar = sugar;
-				return g;
+			if ( d ) {
+				NEXT(tail)=m;
+				tail=m;
 			} else {
-				m = BDY(g); 
-				if ( NEXT(m) ) {
-					BDY(g) = NEXT(m); NEXT(m) = 0;
-				} else {
-					FREEND(g); g = 0;
-				}
-				if ( d ) {
-					for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
-					NEXT(mrd) = m;
-				} else {
-					MKND(n,m,d);
-				}
+				MKND(n,m,d);
+				tail = BDY(d);
 			}
 		}
 	}
@@ -1632,7 +1718,7 @@ ND normalize_pbucket(PGeoBucket g)
 	return r;
 }
 
-ND nd_nf(NODE b,ND g,ND *ps,int full)
+ND nd_nf(ND g,int full)
 {
 	ND u,p,d,red;
 	NODE l;
@@ -1656,41 +1742,31 @@ ND nd_nf(NODE b,ND g,ND *ps,int full)
 			return d;
 		}
 		g = bucket->body[h];
-		for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
-			p = ps[(int)BDY(l)];
-			if ( ndl_reducible(HDL(g),HDL(p)) ) {
-				h_reducible = 1;
-				psugar = HTD(g)-HTD(p) + p->sugar;
-				red = nd_reducer(g,p);
-				bucket->body[h] = nd_remove_head(g);
-				red = nd_remove_head(red);
-				add_pbucket(bucket,red);
-				sugar = MAX(sugar,psugar);
-				break;
+		red = nd_find_reducer(g);
+		if ( red ) {
+			bucket->body[h] = nd_remove_head(g);
+			red = nd_remove_head(red);
+			add_pbucket(bucket,red);
+			sugar = MAX(sugar,red->sugar);
+		} else if ( !full ) {
+			g = normalize_pbucket(bucket);
+			if ( g )
+				g->sugar = sugar;
+			return g;
+		} else {
+			m = BDY(g); 
+			if ( NEXT(m) ) {
+				BDY(g) = NEXT(m); NEXT(m) = 0;
+			} else {
+				FREEND(g); g = 0;
 			}
-		}
-		if ( !h_reducible ) {
-			/* head term is not reducible */
-			if ( !full ) {
-				g = normalize_pbucket(bucket);
-				if ( g )
-					g->sugar = sugar;
-				return g;
+			bucket->body[h] = g;
+			NEXT(m) = 0;
+			if ( d ) {
+				for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
+				NEXT(mrd) = m;
 			} else {
-				m = BDY(g); 
-				if ( NEXT(m) ) {
-					BDY(g) = NEXT(m); NEXT(m) = 0;
-				} else {
-					FREEND(g); g = 0;
-				}
-				bucket->body[h] = g;
-				NEXT(m) = 0;
-				if ( d ) {
-					for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
-					NEXT(mrd) = m;
-				} else {
-					MKND(n,m,d);
-				}
+				MKND(n,m,d);
 			}
 		}
 	}
@@ -1717,7 +1793,7 @@ NODE nd_gb(NODE f)
 #endif
 		l = nd_minp(d,&d);
 		h = nd_sp(l);
-		nf = nd_nf(gall,h,nps,!Top);
+		nf = nd_nf(h,!Top);
 		if ( nf ) {
 			printf("+"); fflush(stdout);
 #if 0
@@ -2022,6 +2098,14 @@ NODE nd_setup(NODE f)
 	nd_bpe = 4;
 	nd_epw = (sizeof(unsigned int)*8)/nd_bpe;
 	nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0);
+	if ( nd_bpe < 32 ) {
+		nd_mask0 = (1<<nd_bpe)-1;
+	} else {
+		nd_mask0 = 0xffffffff;
+	}
+	bzero(nd_mask,sizeof(nd_mask));
+	for ( i = 0; i < nd_epw; i++ )
+		nd_mask[nd_epw-i-1] = (nd_mask0<<(i*nd_bpe));
 	nd_free_private_storage();
 	for ( i = 0; i < nd_psn; i++, f = NEXT(f) ) {
 		nps[i] = dptond((DP)BDY(f));
@@ -2085,6 +2169,7 @@ void dltondl(int n,DL dl,unsigned int *r)
 	int i;
 
 	d = dl->d;
+	bzero(r,nd_wpd*sizeof(unsigned int));
 	if ( is_rlex )
 		for ( i = 0; i < n; i++ )
 			r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe));
@@ -2206,12 +2291,22 @@ void ndp_print(ND_pairs d)
 
 void nd_monic(ND p)
 {
-	int mul;
+	if ( !p )
+		return;
+	else
+		nd_mul_c(p,invm(HC(p),nd_mod));
+}
+
+void nd_mul_c(ND p,int mul)
+{
 	NM m;
+	int c,c1;
 
 	if ( !p )
 		return;
-	mul = invm(HC(p),nd_mod);
-	for ( m = BDY(p); m; m = NEXT(m) )
-		C(m) = (C(m)*mul)%nd_mod;
+	for ( m = BDY(p); m; m = NEXT(m) ) {
+		c1 = C(m);
+		DMAR(c1,mul,0,nd_mod,c);
+		C(m) = c;
+	}
 }