===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.146
retrieving revision 1.149
diff -u -p -r1.146 -r1.149
--- OpenXM_contrib2/asir2000/engine/nd.c	2006/11/27 07:31:25	1.146
+++ OpenXM_contrib2/asir2000/engine/nd.c	2006/12/04 01:40:51	1.149
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.145 2006/08/26 05:38:06 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.148 2006/11/28 02:17:54 noro Exp $ */
 
 #include "nd.h"
 
@@ -52,6 +52,8 @@ void nd_det_reconstruct(NDV **dm,int n,int j,NDV d);
 int nd_monic(int m,ND *p);
 NDV plain_vect_to_ndv_q(Q *mat,int col,UINT *s0vect);
 
+extern int Denominator,DP_Multiple;
+
 void nd_free_private_storage()
 {
 	_nm_free_list = 0;
@@ -1685,6 +1687,7 @@ again:
 	else
 		for ( t = g; t; t = NEXT(t) )
 			BDY(t) = (pointer)nd_ps[(int)BDY(t)];
+	if ( !checkonly && DP_Print ) { printf("nd_gb done.\n"); fflush(stdout); }
 	return g;
 }
 
@@ -1733,6 +1736,17 @@ int do_diagonalize_trace(int sugar,int m)
 static struct oEGT eg_invdalg;
 struct oEGT eg_le;
 
+void nd_subst_vector(VL vl,P p,NODE subst,P *r)
+{
+	NODE tn;
+	P p1;
+
+	for ( tn = subst; tn; tn = NEXT(NEXT(tn)) ) {
+		substp(vl,p,BDY(tn),BDY(NEXT(tn)),&p1); p = p1;
+	}
+	*r = p;
+}
+
 NODE nd_gb_trace(int m,int ishomo)
 {
 	int i,nh,sugar,stat;
@@ -1742,6 +1756,7 @@ NODE nd_gb_trace(int m,int ishomo)
 	ND h,nf,nfq,s,head;
 	NDV nfv,nfqv;
 	Q q,den,num;
+	P hc;
 	union oNDC dn;
 	struct oEGT eg_monic,egm0,egm1;
 	int diag_count = 0;
@@ -1803,7 +1818,11 @@ again:
 			}
 			if ( nfq ) {
 				/* m|HC(nfq) => failure */
-				if ( !rem(NM(HCQ(nfq)),m) ) return 0;
+				if ( nd_vc ) {
+					nd_subst_vector(nd_vc,HCP(nfq),nd_subst,&hc); q = (Q)hc;
+				} else
+					q = HCQ(nfq);
+				if ( !rem(NM(q),m) ) return 0;
 
 				if ( DP_Print ) { printf("+"); fflush(stdout); }
 				if ( nd_nalg ) {
@@ -1848,6 +1867,7 @@ again:
 		print_eg("invdalg",&eg_invdalg);
 		print_eg("le",&eg_le);
 	}
+	if ( DP_Print ) { printf("nd_gb_trace done.\n"); fflush(stdout); }
 	return g;
 }
 
@@ -2395,6 +2415,8 @@ void nd_gr(LIST f,LIST v,int m,int f4,struct order_spe
 	if ( !m && Demand ) nd_demand = 1;
 	else nd_demand = 0;
 
+	if ( DP_Multiple )
+		nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1);
 #if 0
 	ndv_alloc = 0;
 #endif
@@ -2559,6 +2581,9 @@ void nd_gr_trace(LIST f,LIST v,int trace,int homo,int 
 	struct order_spec *ord1;
 	struct oEGT eg_check,eg0,eg1;
 
+	if ( DP_Multiple )
+		nd_scale = ((double)DP_Multiple)/(double)(Denominator?Denominator:1);
+
 	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 ) {
@@ -3298,10 +3323,12 @@ ND nd_copy(ND p)
 
 int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
 {
-	NM m;
+	NM m1,m2;
 	NDV p1,p2;
 	ND t1,t2;
 	UINT *lcm;
+	P gp,tp;
+	Q g,t;
 	int td;
 
 	if ( !mod && nd_demand ) {
@@ -3314,23 +3341,29 @@ int nd_sp(int mod,int trace,ND_pairs p,ND *rp)
 		}
 	}
 	lcm = LCM(p);
-	NEWNM(m);
-	CQ(m) = HCQ(p2);
-	ndl_sub(lcm,HDL(p1),DL(m));
-	if ( ndl_check_bound2(p->i1,DL(m)) )
-		return 0;
-	t1 = ndv_mul_nm(mod,m,p1);
-	if ( mod == -1 ) CM(m) = _chsgnsf(HCM(p1));
-	else if ( mod ) CM(m) = mod-HCM(p1); 
-	else chsgnp(HCP(p1),&CP(m));
-	ndl_sub(lcm,HDL(p2),DL(m));
-	if ( ndl_check_bound2(p->i2,DL(m)) ) {
-		nd_free(t1);
-		return 0;
+	NEWNM(m1); ndl_sub(lcm,HDL(p1),DL(m1));
+	if ( ndl_check_bound2(p->i1,DL(m1)) ) {
+		FREENM(m1); return 0;
 	}
-	t2 = ndv_mul_nm(mod,m,p2);
+	NEWNM(m2); ndl_sub(lcm,HDL(p2),DL(m2));
+	if ( ndl_check_bound2(p->i2,DL(m2)) ) {
+		FREENM(m1); FREENM(m2); return 0;
+	}
+
+	if ( mod == -1 ) {
+		CM(m1) = HCM(p2); CM(m2) = _chsgnsf(HCM(p1));
+	} else if ( mod ) {
+		CM(m1) = HCM(p2); CM(m2) = mod-HCM(p1); 
+	} else if ( nd_vc ) {
+		ezgcdpz(nd_vc,HCP(p1),HCP(p2),&gp);
+		divsp(nd_vc,HCP(p2),gp,&CP(m1));
+		divsp(nd_vc,HCP(p1),gp,&tp); chsgnp(tp,&CP(m2));
+	} else {
+		igcd_cofactor(HCQ(p1),HCQ(p2),&g,&t,&CQ(m1)); chsgnq(t,&CQ(m2));
+	}
+	t1 = ndv_mul_nm(mod,m1,p1); t2 = ndv_mul_nm(mod,m2,p2);
 	*rp = nd_add(mod,t1,t2);
-	FREENM(m);
+	FREENM(m1); FREENM(m2);
 	return 1;
 }
 
@@ -3710,9 +3743,8 @@ void ndv_mod(int mod,NDV p)
 	NMV t,d;
 	int r,s,u;
 	int i,len,dlen;
+	P cp;
 	Q c;
-	P cp,cp1;
-	NODE tn;
 	Obj gfs;
 
 	if ( !p ) return;
@@ -3730,9 +3762,7 @@ void ndv_mod(int mod,NDV p)
 	else
 		for ( t = d = BDY(p), i = 0; i < len; i++, NMV_ADV(t) ) {
 			if ( nd_vc ) {
-				for ( tn = nd_subst, cp = CP(t); tn; tn = NEXT(NEXT(tn)) ) {
-					substp(nd_vc,cp,BDY(tn),BDY(NEXT(tn)),&cp1); cp = cp1;
-				}
+				nd_subst_vector(nd_vc,CP(t),nd_subst,&cp);
 				c = (Q)cp;
 			} else
 				c = CQ(t);