===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/gr.c,v
retrieving revision 1.38
retrieving revision 1.69
diff -u -p -r1.38 -r1.69
--- OpenXM_contrib2/asir2000/builtin/gr.c	2001/11/16 10:35:07	1.38
+++ OpenXM_contrib2/asir2000/builtin/gr.c	2013/11/05 11:36:58	1.69
@@ -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/gr.c,v 1.37 2001/10/09 01:36:06 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/gr.c,v 1.68 2012/12/17 07:20:44 noro Exp $
 */
 #include "ca.h"
 #include "parse.h"
@@ -53,7 +53,7 @@
 #include "ox.h"
 
 #if defined(__GNUC__)
-#define INLINE inline
+#define INLINE static inline
 #elif defined(VISUAL)
 #define INLINE __inline
 #else
@@ -69,7 +69,7 @@ double get_rtime();	
 
 struct oEGT eg_nf,eg_nfm;
 struct oEGT eg_znfm,eg_pz,eg_np,eg_ra,eg_mc,eg_gc;
-int TP,NBP,NMP,NFP,NDP,ZR,NZR;
+int TP,N_BP,NMP,NFP,NDP,ZR,NZR;
 
 extern int (*cmpdl)();
 extern int do_weyl;
@@ -89,26 +89,29 @@ static VL VC;
 
 int PCoeffs;
 int DP_Print = 0;
+int DP_PrintShort = 0;
 int DP_Multiple = 0;
 int DP_NFStat = 0;
 LIST Dist = 0;
 int NoGCD = 0;
 int GenTrace = 0;
+int GenSyz = 0;
 int OXCheck = -1;
+int OneZeroHomo = 0;
 
-static int NoSugar = 0;
+int NoSugar = 0;
 static int NoCriB = 0;
 static int NoGC = 0;
 static int NoMC = 0;
 static int NoRA = 0;
-static int DP_PrintShort = 0;
 static int ShowMag = 0;
 static int Stat = 0;
-static int Denominator = 1;
-static int Top = 0;
-static int Reverse = 0;
+int Denominator = 1;
+int Top = 0;
+int Reverse = 0;
 static int Max_mag = 0;
-static char *Demand = 0;
+static int Max_coef = 0;
+char *Demand = 0;
 static int PtozpRA = 0;
 
 int doing_f4;
@@ -122,17 +125,70 @@ void Pox_pop_local(NODE,Obj *);
 INLINE int eqdl(int nv,DL dl1,DL dl2)
 {
 	int i;
-	int *b1,*b2;
+	int *p1,*p2;
 
 	if ( dl1->td != dl2->td )
 		return 0;
-	for ( i = 0, b1 = dl1->d, b2 = dl2->d; i < nv; i++ )
-		if ( b1[i] != b2[i] )
-			break;
-	if ( i == nv )
-		return 1;
-	else
-		return 0;
+	i = nv-1;
+	p1 = dl1->d;
+	p2 = dl2->d;
+	while ( i >= 7 ) {
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		if ( *p1++ != *p2++ ) return 0;
+		i -= 8;
+	}
+	switch ( i ) {
+		case 6:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 5:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 4:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 3:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 2:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 1:
+			if ( *p1++ != *p2++ ) return 0;
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		case 0:
+			if ( *p1++ != *p2++ ) return 0;
+			return 1;
+		default:
+			return 1;
+	}
 }
 
 /* b[] should be cleared */
@@ -244,12 +300,13 @@ void pdl(NODE f)
 void dp_gr_main(LIST f,LIST v,Num homo,int modular,int field,struct order_spec *ord,LIST *rp)
 {
 	int i,mindex,m,nochk;
-	struct order_spec ord1;
+	struct order_spec *ord1;
 	Q q;
 	VL fv,vv,vc;
 	NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
 	NODE ind,ind0;
 	LIST trace,gbindex;
+	int input_is_dp = 0;
 
 	mindex = 0; nochk = 0; dp_fcoeffs = field;
 	get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
@@ -262,14 +319,22 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int
 		homogenize_order(ord,NVars,&ord1);
 		for ( fd0 = fi0 = 0, t = BDY(f); t; t = NEXT(t) ) {
 			NEXTNODE(fd0,fd); NEXTNODE(fi0,fi);
-			ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi)); dp_homo((DP)BDY(fi),(DP *)&BDY(fd));
+			if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+				dp_sort((DP)BDY(t),(DP *)&BDY(fi)); input_is_dp = 1;
+			} else 
+				ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fi));
+			dp_homo((DP)BDY(fi),(DP *)&BDY(fd));
 		}
 		if ( fd0 ) NEXT(fd) = 0;
 		if ( fi0 ) NEXT(fi) = 0;
-		initd(&ord1);
+		initd(ord1);
 	} else {
 		for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
-			NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
+			NEXTNODE(fd0,fd);
+			if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+				dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1;
+			} else
+				ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
 		}
 		if ( fd0 ) NEXT(fd) = 0;
 		fi0 = fd0;
@@ -286,7 +351,7 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int
 	init_stat();
 	while ( 1 ) {
 		if ( homo ) {
-			initd(&ord1); CNVars = NVars+1;
+			initd(ord1); CNVars = NVars+1;
 		}
 		if ( DP_Print && modular ) {
 			fprintf(asir_out,"mod= %d, eval = ",m); printsubst(subst);
@@ -317,7 +382,10 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int
 	}
 	for ( r0 = 0, ind0 = 0; x; x = NEXT(x) ) {
 		NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
-		dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+		if ( input_is_dp )
+			BDY(r) = (pointer)ps[(int)BDY(x)];
+		else
+			dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
 		NEXTNODE(ind0,ind);
 		STOQ((int)BDY(x),q); BDY(ind) = q;
 	}
@@ -338,16 +406,63 @@ void dp_gr_main(LIST f,LIST v,Num homo,int modular,int
 	}
 	print_stat();
 	if ( ShowMag )
-		fprintf(asir_out,"\nMax_mag=%d\n",Max_mag);
+		fprintf(asir_out,"\nMax_mag=%d, Max_coef=%d\n",Max_mag, Max_coef);
 }
 
+void dp_interreduce(LIST f,LIST v,int field,struct order_spec *ord,LIST *rp)
+{
+	int i,mindex,m,nochk;
+	struct order_spec *ord1;
+	Q q;
+	VL fv,vv,vc;
+	NODE fd,fd0,fi,fi0,r,r0,t,subst,x,s,xx;
+	NODE ind,ind0;
+	LIST trace,gbindex;
+	int input_is_dp = 0;
+
+	mindex = 0; nochk = 0; dp_fcoeffs = field;
+	get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
+	NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
+	CNVars = NVars;
+	if ( ord->id && NVars != ord->nv )
+		error("dp_interreduce : invalid order specification");
+	initd(ord);
+	for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
+		NEXTNODE(fd0,fd);
+		if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+			dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1;
+		} else
+			ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
+	}
+	if ( fd0 ) NEXT(fd) = 0;
+	fi0 = fd0;
+	setup_arrays(fd0,0,&s);
+	init_stat();
+	x = s;
+	reduceall(x,&xx); x = xx;
+	for ( r0 = 0, ind0 = 0; x; x = NEXT(x) ) {
+		NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
+		if ( input_is_dp )
+			BDY(r) = (pointer)ps[(int)BDY(x)];
+		else
+			dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+		NEXTNODE(ind0,ind);
+		STOQ((int)BDY(x),q); BDY(ind) = q;
+	}
+	if ( r0 ) NEXT(r) = 0;
+	if ( ind0 ) NEXT(ind) = 0;
+	MKLIST(*rp,r0);
+	MKLIST(gbindex,ind0);
+}
+
 void dp_gr_mod_main(LIST f,LIST v,Num homo,int m,struct order_spec *ord,LIST *rp)
 {
-	struct order_spec ord1;
+	struct order_spec *ord1;
 	VL fv,vv,vc;
 	NODE fd,fd0,r,r0,t,x,s,xx;
 	DP a,b,c;
-extern struct oEGT eg_red_mod;
+	extern struct oEGT eg_red_mod;
+	int input_is_dp = 0;
 
 	get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
 	NVars = length((NODE)vv); PCoeffs = vc ? 1 : 0; VC = vc;
@@ -357,7 +472,11 @@ extern struct oEGT eg_red_mod;
 	initd(ord);
 	if ( homo ) {
 		for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
-			ptod(CO,vv,(P)BDY(t),&a); dp_homo(a,&b);
+			if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+				dp_sort((DP)BDY(t),&a); input_is_dp = 1;
+			} else 
+				ptod(CO,vv,(P)BDY(t),&a);
+			dp_homo(a,&b);
 			if ( PCoeffs )
 				dp_mod(b,m,0,&c);
 			else
@@ -366,10 +485,13 @@ extern struct oEGT eg_red_mod;
 				NEXTNODE(fd0,fd); BDY(fd) = (pointer)c;
 			}
 		}
-		homogenize_order(ord,NVars,&ord1); initd(&ord1);
+		homogenize_order(ord,NVars,&ord1); initd(ord1);
 	} else {
 		for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) {
-			ptod(CO,vv,(P)BDY(t),&b);
+			if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+				dp_sort((DP)BDY(t),&b); input_is_dp = 1;
+			} else 
+				ptod(CO,vv,(P)BDY(t),&b);
 			if ( PCoeffs )
 				dp_mod(b,m,0,&c);
 			else
@@ -383,7 +505,7 @@ extern struct oEGT eg_red_mod;
 	setup_arrays(fd0,m,&s);
 	init_stat();
 	if ( homo ) {
-		initd(&ord1); CNVars = NVars+1;
+		initd(ord1); CNVars = NVars+1;
 	}
 /* init_eg(&eg_red_mod); */
 	x = gb_mod(s,m);
@@ -395,11 +517,19 @@ extern struct oEGT eg_red_mod;
 	reduceall_mod(x,m,&xx); x = xx;
 	if ( PCoeffs )
 		for ( r0 = 0; x; x = NEXT(x) ) {
-			NEXTNODE(r0,r); mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+			NEXTNODE(r0,r); 
+			if ( input_is_dp )
+				mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r));
+			else
+				mdtop(CO,m,vv,ps[(int)BDY(x)],(P *)&BDY(r));
 		}
 	else
 		for ( r0 = 0; x; x = NEXT(x) ) {
-			NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+			NEXTNODE(r0,r); 
+			if ( input_is_dp )
+				_mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r));
+			else
+				_dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
 		}
 	print_stat();
 	if ( r0 ) NEXT(r) = 0;
@@ -411,6 +541,7 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L
 	int homogen;
 	VL fv,vv,vc;
 	NODE fd,fd0,r,r0,t,x,s,xx;
+	int input_is_dp = 0;
 
 	dp_fcoeffs = 0;
 	get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
@@ -420,7 +551,11 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L
 		error("dp_f4_main : invalid order specification");
 	initd(ord);
 	for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
-		NEXTNODE(fd0,fd); ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
+		NEXTNODE(fd0,fd); 
+		if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+			dp_sort((DP)BDY(t),(DP *)&BDY(fd)); input_is_dp = 1;
+		} else 
+			ptod(CO,vv,(P)BDY(t),(DP *)&BDY(fd));
 		if ( homogen )
 			homogen = dp_homogeneous(BDY(fd));
 	}
@@ -432,7 +567,10 @@ void dp_f4_main(LIST f,LIST v,struct order_spec *ord,L
 	}
 	for ( r0 = 0; x; x = NEXT(x) ) {
 		NEXTNODE(r0,r); dp_load((int)BDY(x),&ps[(int)BDY(x)]);
-		dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+		if ( input_is_dp )
+			BDY(r) = (pointer)ps[(int)BDY(x)];
+		else
+			dtop(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
 	}
 	if ( r0 ) NEXT(r) = 0;
 	MKLIST(*rp,r0);
@@ -444,6 +582,7 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s
 	VL fv,vv,vc;
 	DP b,c,c1;
 	NODE fd,fd0,r,r0,t,x,s,xx;
+	int input_is_dp = 0;
 
 	dp_fcoeffs = 0;
 	get_vars((Obj)f,&fv); pltovl(v,&vv); vlminus(fv,vv,&vc);
@@ -453,7 +592,10 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s
 		error("dp_f4_mod_main : invalid order specification");
 	initd(ord);
 	for ( fd0 = 0, t = BDY(f), homogen = 1; t; t = NEXT(t) ) {
-		ptod(CO,vv,(P)BDY(t),&b);
+		if ( BDY(t) && OID(BDY(t)) == O_DP ) {
+			dp_sort((DP)BDY(t),&b); input_is_dp = 1;
+		} else 
+			ptod(CO,vv,(P)BDY(t),&b);
 		if ( homogen )
 			homogen = dp_homogeneous(b);
 		_dp_mod(b,m,0,&c);
@@ -473,7 +615,11 @@ void dp_f4_mod_main(LIST f,LIST v,int m,struct order_s
 		reduceall_mod(x,m,&xx); x = xx;
 	}
 	for ( r0 = 0; x; x = NEXT(x) ) {
-		NEXTNODE(r0,r); _dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
+		NEXTNODE(r0,r);
+		if ( input_is_dp )
+			_mdtodp(ps[(int)BDY(x)],(DP *)&BDY(r));
+		else
+			_dtop_mod(CO,vv,ps[(int)BDY(x)],(P *)&BDY(r));
 	}
 	if ( r0 ) NEXT(r) = 0;
 	MKLIST(*rp,r0);
@@ -600,6 +746,8 @@ NODE gb_f4(NODE f)
 /* initial bases are monic */
 
 unsigned int **psca;
+GeoBucket create_bucket();
+DL remove_head_bucket(GeoBucket,int);
 
 NODE gb_f4_mod(NODE f,int m)
 {
@@ -618,14 +766,16 @@ NODE gb_f4_mod(NODE f,int m)
 	int *indred,*isred;
 	CDP ri;
 	int pscalen;
-	struct oEGT tmp0,tmp1,eg_split_symb,eg_split_elim1,eg_split_elim2;
-	extern struct oEGT eg_symb,eg_elim1,eg_elim2;
+	GeoBucket bucket;
+	DL head;
+	struct oEGT tmp0,tmp1,eg_split_symb,eg_split_conv,eg_split_elim1,eg_split_elim2;
+	extern struct oEGT eg_symb,eg_conv,eg_elim1,eg_elim2;
 
 	/* initialize coeffcient array list of ps[] */
 	pscalen = pslen;
 	psca = (unsigned int **)MALLOC(pscalen*sizeof(unsigned int *));
 
-	init_eg(&eg_symb); init_eg(&eg_elim1); init_eg(&eg_elim2);
+	init_eg(&eg_symb); init_eg(&eg_conv); init_eg(&eg_elim1); init_eg(&eg_elim2);
 	for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
 		i = (int)BDY(r);
 		d = updpairs(d,g,i);
@@ -640,41 +790,57 @@ NODE gb_f4_mod(NODE f,int m)
 		minsugar(d,&dm,&dr); d = dr;
 		if ( DP_Print )
 			fprintf(asir_out,"sugar=%d\n",dm->sugar);
-		blist = 0; s0 = 0;
+		blist = 0;
+		bucket = create_bucket();
 		/* asph : sum of all head terms of spoly */
 		for ( t = dm; t; t = NEXT(t) ) {
 			_dp_sp_mod(ps[t->dp1],ps[t->dp2],m,&sp);
 /*			fprintf(stderr,"splen=%d-",dp_nt(sp)); */
 			if ( sp ) {
 				MKNODE(bt,sp,blist); blist = bt;
-				s0 = symb_merge(s0,dp_dllist(sp),nv);
+				add_bucket(bucket,dp_dllist(sp),nv);
 /*				fprintf(stderr,"%d-",length(s0)); */
 			}
 		}
+#if 0
 		if ( DP_Print )
 			fprintf(asir_out,"initial spmat : %d x %d ",length(blist),length(s0));
+#endif
 		/* s0 : all the terms appeared in symbolic reduction */
-		for ( s = s0, nred = 0; s; s = NEXT(s) ) {
+		nred = 0;
+		s0 = 0;
+		while ( 1 ) {
+			head = remove_head_bucket(bucket,nv);
+			if ( !head ) break;
+			else {
+				NEXTNODE(s0,s);
+				BDY(s) = (pointer)head;
+			}
 			for ( r = gall;	r; r = NEXT(r) )
-				if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,BDY(s),nv) )
+				if ( _dl_redble(BDY(ps[(int)BDY(r)])->dl,head,nv) )
 					break;
 			if ( r ) {
-				dltod(BDY(s),nv,&tdp);
+				dltod(head,nv,&tdp);
 				dp_subd(tdp,ps[(int)BDY(r)],&sd);
 				dt = mul_dllist(BDY(sd)->dl,ps[(int)BDY(r)]);
+				add_bucket(bucket,NEXT(dt),nv);
 /*				fprintf(stderr,"[%d]",length(dt)); */
 				/* list of [t,f] */
 				bt1 = mknode(2,BDY(sd)->dl,BDY(r));
 				MKNODE(bt,bt1,blist); blist = bt;
-				symb_merge(s,dt,nv);
-/*				fprintf(stderr,"%d-",length(s)); */
+/*				fprintf(stderr,"%d-",length(s0));  */
 				nred++;
 			}
 		}
+		if ( s0 ) NEXT(s) = 0;
 /*		fprintf(stderr,"\n"); */
+		get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
+		init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
+
 		if ( DP_Print )
 			fprintf(asir_out,"number of reducers : %d\n",nred);
 
+		get_eg(&tmp0);
 		/* the first nred polys in blist are reducers */
 		/* row = the number of all the polys */
 		for ( r = blist, row = 0; r; r = NEXT(r), row++ );
@@ -713,8 +879,8 @@ NODE gb_f4_mod(NODE f,int m)
 		for ( j = 0, k = 0; j < col; j++ )
 			if ( !isred[j] )
 				st[k++] = at[j];
-		get_eg(&tmp1); add_eg(&eg_symb,&tmp0,&tmp1);
-		init_eg(&eg_split_symb); add_eg(&eg_split_symb,&tmp0,&tmp1);
+		get_eg(&tmp1); add_eg(&eg_conv,&tmp0,&tmp1);
+		init_eg(&eg_split_conv); add_eg(&eg_split_conv,&tmp0,&tmp1);
 
 		get_eg(&tmp1);
 		/* spoly matrix; stored in reduced form; terms in ht[] are omitted */
@@ -745,8 +911,8 @@ NODE gb_f4_mod(NODE f,int m)
 
 		/* XXX free redmat explicitly */
 		for ( k = 0; k < nred; k++ ) {
-			GC_free(BDY(redmat[k]));
-			GC_free(redmat[k]);
+			GCFREE(BDY(redmat[k]));
+			GCFREE(redmat[k]);
 		}
 
 		get_eg(&tmp0); add_eg(&eg_elim1,&tmp1,&tmp0);
@@ -771,6 +937,7 @@ NODE gb_f4_mod(NODE f,int m)
 		if ( DP_Print ) {
 			fprintf(asir_out,"done rank = %d\n",rank,row,col);
 			print_eg("Symb",&eg_split_symb);
+			print_eg("Conv",&eg_split_conv);
 			print_eg("Elim1",&eg_split_elim1);
 			print_eg("Elim2",&eg_split_elim2);
 			fprintf(asir_out,"\n");
@@ -808,11 +975,12 @@ NODE gb_f4_mod(NODE f,int m)
 
 		/* XXX free spmat[] explicitly */
 		for ( j = 0; j < nsp; j++ ) {
-			GC_free(spmat[j]);
+			GCFREE(spmat[j]);
 		}
 	}
 	if ( DP_Print ) {
 		print_eg("Symb",&eg_symb);
+		print_eg("Conv",&eg_conv);
 		print_eg("Elim1",&eg_elim1);
 		print_eg("Elim2",&eg_elim2);
 		fflush(asir_out);
@@ -1019,6 +1187,20 @@ void pltovl(LIST l,VL *vl)
 	*vl = r0;
 }
 
+void vltopl(VL vl,LIST *l)
+{
+	VL n;
+	NODE r,r0;
+	P p;
+
+	n = vl;
+	for ( r0 = 0; n; n = NEXT(n) ) {
+		NEXTNODE(r0,r); MKV(n->v,p); BDY(r) = (pointer)p;
+	}
+	if ( r0 ) NEXT(r) = 0;
+	MKLIST(*l,r0);
+}
+
 void makesubst(VL v,NODE *s)
 {
 	NODE r,r0;
@@ -1158,11 +1340,9 @@ void prim_part(DP f,int m,DP *r)
 		if ( PCoeffs )
 			dp_prim_mod(f,m,r);
 		else
-			*r = f;
+			_dp_monic(f,m,r);
 	} else {
-		if ( dp_fcoeffs )
-			*r = f;
-		else if ( PCoeffs )
+		if ( dp_fcoeffs || PCoeffs )
 			dp_prim(f,r);
 		else
 			dp_ptozp(f,r);
@@ -1251,9 +1431,10 @@ void reduceall(NODE in,NODE *h)
 	w = (int *)ALLOCA(n*sizeof(int));
 	for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
 		w[i] = (int)BDY(t);
+	/* w[i] < 0 : reduced to 0 */
 	for ( i = 0; i < n; i++ ) {
 		for ( top = 0, j = n-1; j >= 0; j-- )
-			if ( j != i ) {
+			if ( w[j] >= 0 && j != i ) {
 				MKNODE(t,(pointer)w[j],top); top = t;
 			}
 		get_eg(&tmp0);
@@ -1275,10 +1456,16 @@ void reduceall(NODE in,NODE *h)
 		if ( DP_Print || DP_PrintShort ) {
 			fprintf(asir_out,"."); fflush(asir_out);
 		}
-		w[i] = newps(g1,0,(NODE)0);
+		if ( g1 ) {
+			w[i] = newps(g1,0,(NODE)0);
+		} else {
+			w[i] = -1;
+		}
 	}
 	for ( top = 0, j = n-1; j >= 0; j-- ) {
-		MKNODE(t,(pointer)w[j],top); top = t;
+		if ( w[j] >= 0 ) {
+			MKNODE(t,(pointer)w[j],top); top = t;
+		}
 	}
 	*h = top;
 	if ( DP_Print || DP_PrintShort )
@@ -1304,9 +1491,10 @@ void reduceall_mod(NODE in,int m,NODE *h)
 	w = (int *)ALLOCA(n*sizeof(int));
 	for ( i = 0, t = r; i < n; i++, t = NEXT(t) )
 		w[i] = (int)BDY(t);
+	/* w[i] < 0 : reduced to 0 */
 	for ( i = 0; i < n; i++ ) {
 		for ( top = 0, j = n-1; j >= 0; j-- )
-			if ( j != i ) {
+			if ( w[j] >= 0 && j != i ) {
 				MKNODE(t,(pointer)w[j],top); top = t;
 			}
 		get_eg(&tmp0);
@@ -1320,10 +1508,16 @@ void reduceall_mod(NODE in,int m,NODE *h)
 		if ( DP_Print || DP_PrintShort ) {
 			fprintf(asir_out,"."); fflush(asir_out);
 		}
-		w[i] = newps_mod(g,m);
+		if ( g ) {
+			w[i] = newps_mod(g,m);
+		} else {
+			w[i] = -1;
+		}
 	}
 	for ( top = 0, j = n-1; j >= 0; j-- ) {
-		MKNODE(t,(pointer)w[j],top); top = t;
+		if ( w[j] >= 0 ) {
+			MKNODE(t,(pointer)w[j],top); top = t;
+		}
 	}
 	*h = top;
 	if ( DP_Print || DP_PrintShort )
@@ -1462,7 +1656,7 @@ void reducebase_dehomo(NODE f,NODE *g)
 				NODE node;
 
 				STOQ(r[i],q);
-				node = mknode(4,0,q,0,0);
+				node = mknode(4,NULLP,q,NULLP,NULLP);
 				MKLIST(hist,node);
 				MKNODE(TraceList,hist,0);
 			}
@@ -1537,7 +1731,7 @@ void minsugar(DP_pairs d,DP_pairs *dm,DP_pairs *dr)
 
 NODE gb(NODE f,int m,NODE subst)
 {
-	int i,nh,prev,mag;
+	int i,nh,prev,mag,mag0,magt;
 	NODE r,g,gall;
 	DP_pairs d;
 	DP_pairs l;
@@ -1551,8 +1745,10 @@ NODE gb(NODE f,int m,NODE subst)
 	static prev_sugar = -1;
 
 	Max_mag = 0;
+	Max_coef = 0;
 	prev = 1;
 	doing_f4 = 0;
+	init_denomlist();
 	if ( m ) {
 		psm = (DP *)MALLOC(pslen*sizeof(DP));
 		for ( i = 0; i < psn; i++ )
@@ -1602,7 +1798,7 @@ NODE gb(NODE f,int m,NODE subst)
 				_dp_nf(gall,h,ps,!Top,&nf);
 			else
 				_dp_nf_z(gall,h,ps,!Top,DP_Multiple,&nf);
-			if ( DP_Print )
+			if ( DP_Print && nf )
 				fprintf(asir_out,"(%.3g)",get_rtime()-t_0);
 			get_eg(&tnf1); add_eg(&eg_nf,&tnf0,&tnf1);
 		} else
@@ -1613,6 +1809,7 @@ skip_nf:
 			get_eg(&tpz0);
 			prim_part(nf,0,&h);
 			get_eg(&tpz1); add_eg(&eg_pz,&tpz0,&tpz1);
+			add_denomlist(BDY(h)->c);
 			get_eg(&tnp0);
 			if ( Demand && skip_nf_flag )
 				nh = newps_nosave(h,m,subst);
@@ -1625,8 +1822,12 @@ skip_nf:
 			g = updbase(g,nh);
 			gall = append_one(gall,nh);
 			if ( !dp_fcoeffs && ShowMag ) {
-				for ( mag = 0, mp = BDY(h); mp; mp = NEXT(mp) )
-					mag += p_mag((P)mp->c);
+				for ( mag = 0, mag0 = 0, mp = BDY(h); mp; mp = NEXT(mp) ) {
+					magt = p_mag((P)mp->c);
+					mag0 = MAX(mag0,magt);
+					mag += magt;
+				}
+				Max_coef = MAX(Max_coef,mag0);
 				Max_mag = MAX(Max_mag,mag);
 			}
 			if ( DP_Print ) {
@@ -1638,7 +1839,7 @@ skip_nf:
 					l->dp1,l->dp2,length(g),length(gall),DPPlength(d),
 					pss[nh]);
 				if ( ShowMag )
-					fprintf(asir_out,",mag=%d",mag);
+					fprintf(asir_out,",mag=(%d,%d)",mag,mag0);
 				fprintf(asir_out,"\n"); fflush(asir_out);
 			} else if ( DP_PrintShort ) {
 				fprintf(asir_out,"+"); fflush(asir_out);
@@ -1733,7 +1934,7 @@ DP_pairs updpairs( DP_pairs d, NODE /* of index */ g, 
 	if ( !NoCriB && d ) {
 		dl = DPPlength(d);
 		d = criterion_B( d, t );
-		dl -= DPPlength(d); NBP += dl;
+		dl -= DPPlength(d); N_BP += dl;
 	}
 	d1 = newpairs( g, t );
 	if ( NEXT(d1) ) {
@@ -1779,7 +1980,7 @@ DP_pairs newpairs( NODE /* of index */ g, int t )
 		p->lcm = lcm_of_DL(CNVars, dl = psh[dp], tdl, (DL)0 );
 #if 0
 		if ( do_weyl )
-			p->sugar = dl_weight(p->lcm);
+			p->sugar = dl_weyl_weight(p->lcm);
 		else
 #endif
 			p->sugar = (ts > (s = pss[dp] - dl->td) ? ts : s) + p->lcm->td;
@@ -1895,7 +2096,7 @@ DP_pairs criterion_F( DP_pairs d1 )
 	for ( head = last = 0, p = d1; NEXT(p); ) {
 		s = (r = w = collect_pairs_of_hdlcm( p, &rest ))->sugar;
 		while ( w = NEXT(w) )
-			if ( criterion_2( w->dp1, w->dp2 ) ) {
+			if ( !do_weyl && criterion_2( w->dp1, w->dp2 ) ) {
 				r = w;
 				break;
 			} else if ( w->sugar < s ) s = (r = w)->sugar;
@@ -1946,7 +2147,8 @@ int dl_redble(DL dl1,DL dl2)
 	return 1;
 }
 
-int dl_weight(DL dl)
+#if 0
+int dl_weyl_weight(DL dl)
 {
 	int n,w,i;
 
@@ -1955,6 +2157,7 @@ int dl_weight(DL dl)
 		w += (-dl->d[i]+dl->d[n+i]);
 	return w;
 }
+#endif
 
 int gbcheck(NODE f)
 {
@@ -1994,27 +2197,27 @@ int gbcheck(NODE f)
 	return 1;
 }
 
-void gbcheck_list(NODE f,LIST *gp,LIST *pp)
+void gbcheck_list(NODE f,int n,VECT *gp,LIST *pp)
 {
 	int i;
 	NODE r,g,gall,u,u0,t;
+	VECT vect;
 	LIST pair;
 	DP_pairs d,l;
 	Q q1,q2;
 
-	for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) {
+	/* we need the following settings */
+	NVars = CNVars = n;
+	setup_arrays(f,0,&r);
+	for ( gall = g = 0, d = 0; r; r = NEXT(r) ) {
 		i = (int)BDY(r);
 		d = updpairs(d,g,i);
 		g = updbase(g,i);
 		gall = append_one(gall,i);
 	}
-	for ( u0 = 0, t = gall; t; t = NEXT(t) ) {
-		NEXTNODE(u0,u);
-		BDY(u) = ps[(int)BDY(t)];
-	}
-	if ( u0 )
-		NEXT(u) = 0;
-	MKLIST(*gp,u);
+	NEWVECT(vect); vect->len = psn; vect->body = (pointer)ps;
+	*gp = vect;
+
 	for ( u0 = 0, l = d; l; l = NEXT(l) ) {
 		NEXTNODE(u0,u);
 		STOQ(l->dp1,q1);
@@ -2025,7 +2228,7 @@ void gbcheck_list(NODE f,LIST *gp,LIST *pp)
 	}
 	if ( u0 )
 		NEXT(u) = 0;
-	MKLIST(*pp,u);
+	MKLIST(*pp,u0);
 }
 
 int membercheck(NODE f,NODE x)
@@ -2059,6 +2262,7 @@ void dp_set_flag(Obj name,Obj value)
 {
 	char *n;
 	int v;
+	Q ratio;
 
 	if ( OID(name) != O_STR )
 		return;
@@ -2069,6 +2273,16 @@ void dp_set_flag(Obj name,Obj value)
 	if ( !strcmp(n,"Dist") ) {
 		Dist = (LIST)value; return;
 	}
+	if ( !strcmp(n,"Content") ) {
+		ratio = (Q)value;
+		if ( ratio ) {
+			DP_Multiple = BD(NM(ratio))[0];
+			Denominator = INT(ratio) ? 1 : BD(DN(ratio))[0];
+		} else {
+			DP_Multiple = 0;
+			Denominator = 1;
+		}
+	}
 	if ( value && OID(value) != O_N )
 		return;
 	v = QTOS((Q)value);
@@ -2108,16 +2322,28 @@ void dp_set_flag(Obj name,Obj value)
 		GenTrace = v;
 	else if ( !strcmp(n,"OXCheck") )
 		OXCheck = v;
+	else if ( !strcmp(n,"GenSyz") )
+		GenSyz = v;
+	else if ( !strcmp(n,"OneZeroHomo") )
+		OneZeroHomo = v;
 }
 
 void dp_make_flaglist(LIST *list)
 {
-	Q v;
+	Q v,nm,dn;
 	STRING name,path;
 	NODE n,n1;
 
+#if 0
 	STOQ(DP_Multiple,v); MKNODE(n,v,0); MKSTR(name,"DP_Multiple"); MKNODE(n1,name,n); n = n1;
 	STOQ(Denominator,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Denominator"); MKNODE(n1,name,n); n = n1;
+#else
+	if ( DP_Multiple ) {
+		STOQ(DP_Multiple,nm); STOQ(Denominator,dn); divq(nm,dn,&v);
+	} else
+		v = 0;
+	MKNODE(n,v,0); MKSTR(name,"Content"); MKNODE(n1,name,n); n = n1;
+#endif
 	MKNODE(n1,Dist,n); n = n1; MKSTR(name,"Dist"); MKNODE(n1,name,n); n = n1;
 	STOQ(Reverse,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Reverse"); MKNODE(n1,name,n); n = n1;
 	STOQ(Stat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Stat"); MKNODE(n1,name,n); n = n1;
@@ -2126,6 +2352,7 @@ void dp_make_flaglist(LIST *list)
 	STOQ(DP_NFStat,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"NFStat"); MKNODE(n1,name,n); n = n1;
 	STOQ(OXCheck,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"OXCheck"); MKNODE(n1,name,n); n = n1;
 	STOQ(GenTrace,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenTrace"); MKNODE(n1,name,n); n = n1;
+	STOQ(GenSyz,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"GenSyz"); MKNODE(n1,name,n); n = n1;
 	STOQ(PtozpRA,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"PtozpRA"); MKNODE(n1,name,n); n = n1;
 	STOQ(ShowMag,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"ShowMag"); MKNODE(n1,name,n); n = n1;
 	STOQ(Top,v); MKNODE(n1,v,n); n = n1; MKSTR(name,"Top"); MKNODE(n1,name,n); n = n1;
@@ -2170,7 +2397,11 @@ void dp_load(int index,DP *p)
 		sprintf(path,"%s%c%d",Demand,DELIM,index);
 		if ( !(fp = fopen(path,"rb") ) )
 			error("dp_load : cannot open a file");
-		skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp);
+		if ( PCoeffs )
+			loadvl(fp);
+		else
+			skipvl(fp);
+		loadobj(fp,(Obj *)p); fclose(fp);
 	}
 }
 
@@ -2183,7 +2414,11 @@ int dp_load_t(int index,DP *p)
 	if ( !(fp = fopen(path,"rb") ) )
 		return 0;
 	else {
-		skipvl(fp); loadobj(fp,(Obj *)p); fclose(fp); return 1;
+		if ( PCoeffs )
+			loadvl(fp);
+		else
+			skipvl(fp);
+		loadobj(fp,(Obj *)p); fclose(fp); return 1;
 	}
 }
 
@@ -2191,7 +2426,7 @@ void init_stat() {
 	init_eg(&eg_nf); init_eg(&eg_nfm); init_eg(&eg_znfm);
 	init_eg(&eg_pz); init_eg(&eg_np);
 	init_eg(&eg_ra); init_eg(&eg_mc); init_eg(&eg_gc);
-	ZR = NZR = TP = NMP = NBP = NFP = NDP = 0;
+	ZR = NZR = TP = NMP = N_BP = NFP = NDP = 0;
 }
 
 void print_stat() {
@@ -2200,7 +2435,7 @@ void print_stat() {
 	print_eg("NF",&eg_nf); print_eg("NFM",&eg_nfm); print_eg("ZNFM",&eg_znfm);
 	print_eg("PZ",&eg_pz); print_eg("NP",&eg_np);
 	print_eg("RA",&eg_ra); print_eg("MC",&eg_mc); print_eg("GC",&eg_gc);
-	fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,NBP,NMP,NFP,NDP,ZR,NZR);
+	fprintf(asir_out,"T=%d,B=%d M=%d F=%d D=%d ZR=%d NZR=%d\n",TP,N_BP,NMP,NFP,NDP,ZR,NZR);
 }
 
 /*
@@ -2335,7 +2570,7 @@ void _dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple
 				if ( GenTrace ) {
 					/* u = cr*rp + (-cred)*shift*red */ 
 					STOQ((int)BDY(l),cq);
-					node = mknode(4,cr,cq,0,0);
+					node = mknode(4,cr,cq,NULLP,NULLP);
 					mulq(cred,rc,&rcred);
 					chsgnnum((Num)rcred,(Num *)&mrcred);
 					muldc(CO,shift,(P)mrcred,(DP *)&ARG2(node));