===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.111
retrieving revision 1.118
diff -u -p -r1.111 -r1.118
--- OpenXM_contrib2/asir2000/engine/nd.c	2004/09/21 07:19:01	1.111
+++ OpenXM_contrib2/asir2000/engine/nd.c	2004/12/04 09:39:27	1.118
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.110 2004/09/21 05:23:14 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.117 2004/12/03 08:57:30 noro Exp $ */
 
 #include "nd.h"
 
@@ -8,6 +8,7 @@ NM _nm_free_list;
 ND _nd_free_list;
 ND_pairs _ndp_free_list;
 
+static int nd_nalg;
 #if 0
 static int ndv_alloc;
 #endif
@@ -41,6 +42,10 @@ static int nd_found,nd_create,nd_notfirst;
 static int nmv_adv;
 static int nd_demand;
 
+UINT *nd_det_compute_bound(NDV **dm,int n,int j);
+void nd_det_reconstruct(NDV **dm,int n,int j,NDV d);
+int nd_monic(int m,ND *p);
+
 void nd_free_private_storage()
 {
 	_nm_free_list = 0;
@@ -724,13 +729,11 @@ int ndl_disjoint(UINT *d1,UINT *d2)
 #endif
 }
 
-int ndl_check_bound2(int index,UINT *d2)
+int ndl_check_bound(UINT *d1,UINT *d2)
 {
 	UINT u2;
-	UINT *d1;
 	int i,j,ind,k;
 
-	d1 = nd_bound[index];
 	ind = 0;
 #if USE_UNROLL
 	switch ( nd_bpe ) {
@@ -819,6 +822,11 @@ int ndl_check_bound2(int index,UINT *d2)
 #endif
 }
 
+int ndl_check_bound2(int index,UINT *d2)
+{
+	return ndl_check_bound(nd_bound[index],d2);
+}
+
 INLINE int ndl_hash_value(UINT *d)
 {
 	int i;
@@ -1574,7 +1582,7 @@ NODE nd_gb(int m,int ishomo,int checkonly)
 	NODE r,g,t;
 	ND_pairs d;
 	ND_pairs l;
-	ND h,nf,s,head;
+	ND h,nf,s,head,nf1;
 	NDV nfv;
 	Q q,num,den;
 	union oNDC dn;
@@ -1613,6 +1621,10 @@ again:
 			if ( checkonly ) return 0;
 			if ( DP_Print ) { printf("+"); fflush(stdout); }
 			nd_removecont(m,nf);
+			if ( !m && nd_nalg ) {
+				nd_monic(0,&nf);
+				nd_removecont(m,nf);
+			}
 			nfv = ndtondv(m,nf); nd_free(nf);
 			nh = ndv_newps(m,nfv,0);
 			d = update_pairs(d,g,nh);
@@ -1679,6 +1691,9 @@ void do_diagonalize_trace(int sugar,int m)
 	}
 }
 
+static struct oEGT eg_invdalg;
+struct oEGT eg_le;
+
 NODE nd_gb_trace(int m,int ishomo)
 {
 	int i,nh,sugar,stat;
@@ -1689,7 +1704,11 @@ NODE nd_gb_trace(int m,int ishomo)
 	NDV nfv,nfqv;
 	Q q,den,num;
 	union oNDC dn;
+	struct oEGT eg_monic,egm0,egm1;
 
+	init_eg(&eg_monic);
+	init_eg(&eg_invdalg);
+	init_eg(&eg_le);
 	g = 0; d = 0;
 	for ( i = 0; i < nd_psn; i++ ) {
 		d = update_pairs(d,g,i);
@@ -1737,8 +1756,17 @@ again:
 				if ( !rem(NM(HCQ(nfq)),m) ) return 0;
 
 				if ( DP_Print ) { printf("+"); fflush(stdout); }
-				nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf);
-				nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq);
+				if ( nd_nalg ) {
+					/* m|DN(HC(nf)^(-1)) => failure */
+					get_eg(&egm0);
+					if ( !nd_monic(m,&nfq) ) return 0;
+					get_eg(&egm1); add_eg(&eg_monic,&egm0,&egm1);
+					nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq);
+					nfv = ndv_dup(0,nfqv); ndv_mod(m,nfv); nd_free(nf);
+				} else {
+					nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq);
+					nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf);
+				}
 				nh = ndv_newps(0,nfv,nfqv);
 				d = update_pairs(d,g,nh);
 				g = update_base(g,nh);
@@ -1756,6 +1784,11 @@ again:
 	else
 		for ( t = g; t; t = NEXT(t) )
 			BDY(t) = (pointer)nd_ps_trace[(int)BDY(t)];
+	if ( nd_nalg ) {
+		print_eg("monic",&eg_monic);
+		print_eg("invdalg",&eg_invdalg);
+		print_eg("le",&eg_le);
+	}
 	return g;
 }
 
@@ -2616,6 +2649,7 @@ void nd_mul_c(int mod,ND p,int mul)
 	int c,c1;
 
 	if ( !p ) return;
+	if ( mul == 1 ) return;
 	if ( mod == -1 )
 		for ( m = BDY(p); m; m = NEXT(m) )
 			CM(m) = _mulsf(CM(m),mul);
@@ -2631,6 +2665,7 @@ void nd_mul_c_q(ND p,Q mul)
 	Q c;
 
 	if ( !p ) return;
+	if ( UNIQ(mul) ) return;
 	for ( m = BDY(p); m; m = NEXT(m) ) {
 		mulq(CQ(m),mul,&c); CQ(m) = c;
 	}
@@ -3709,6 +3744,10 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec
 	int stat,nvar,max,e;
 	union oNDC dn;
 
+	if ( !f ) {
+		*rp = 0;
+		return;
+	}
 	pltovl(v,&vv);
 	for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ );
 
@@ -4887,6 +4926,7 @@ void nd_det(int mod,MAT f,P *rp)
 	NDV d,s,mij,mjj;
 	ND u;
 	NMV nmv;
+	UINT *bound;
 	PGeoBucket bucket;
 	struct order_spec *ord;
 
@@ -4904,7 +4944,7 @@ void nd_det(int mod,MAT f,P *rp)
 				e = getdeg(tv->v,(P)m[i][j]);
 				max = MAX(e,max);
 			}
-	nd_setup_parameters(nvar,1024);
+	nd_setup_parameters(nvar,max);
 	dm = (NDV **)almat_pointer(n,n);
 	for ( i = 0, max = 0; i < n; i++ )
 		for ( j = 0; j < n; j++ ) {
@@ -4938,15 +4978,19 @@ void nd_det(int mod,MAT f,P *rp)
 			}
 			sgn = -sgn;
 		}
+		bound = nd_det_compute_bound(dm,n,j);
+		if ( ndl_check_bound(bound,bound) )
+			nd_det_reconstruct(dm,n,j,d);
+
 		for ( i = j+1, mj = dm[j], mjj = mj[j]; i < n; i++ ) {
-			if ( DP_Print ) fprintf(stderr,"	i=%d\n		",i);
+/*			if ( DP_Print ) fprintf(stderr,"	i=%d\n		",i); */
 			mi = dm[i]; mij = mi[j];
 			if ( mod )
 				ndv_mul_c(mod,mij,mod-1);
 			else
 				ndv_mul_c_q(mij,mone);
 			for ( k = j+1; k < n; k++ ) {
-				if ( DP_Print ) fprintf(stderr,"k=%d ",k);
+/*				if ( DP_Print ) fprintf(stderr,"k=%d ",k); */
 				bucket = create_pbucket();
 				if ( mi[k] ) {
 					nmv = BDY(mjj); len = LEN(mjj);
@@ -4965,7 +5009,7 @@ void nd_det(int mod,MAT f,P *rp)
 				u = nd_quo(mod,bucket,d);
 				mi[k] = ndtondv(mod,u);
 			}
-			if ( DP_Print ) fprintf(stderr,"\n",k);
+/*			if ( DP_Print ) fprintf(stderr,"\n",k); */
 		}
 		d = mjj;
 	}
@@ -5025,4 +5069,184 @@ ND ndv_mul_nmv_trunc(int mod,NMV m0,NDV p,UINT *d)
 			return r;
 		}
 	}
+}
+
+void nd_det_reconstruct(NDV **dm,int n,int j,NDV d)
+{
+	int i,obpe,oadv,h,k,l;
+	static NM prev_nm_free_list;
+	EPOS oepos;
+
+	obpe = nd_bpe;
+	oadv = nmv_adv;
+	oepos = nd_epos;
+	if ( obpe < 2 ) nd_bpe = 2;
+	else if ( obpe < 3 ) nd_bpe = 3;
+	else if ( obpe < 4 ) nd_bpe = 4;
+	else if ( obpe < 5 ) nd_bpe = 5;
+	else if ( obpe < 6 ) nd_bpe = 6;
+	else if ( obpe < 8 ) nd_bpe = 8;
+	else if ( obpe < 10 ) nd_bpe = 10;
+	else if ( obpe < 16 ) nd_bpe = 16;
+	else if ( obpe < 32 ) nd_bpe = 32;
+	else error("nd_det_reconstruct : exponent too large");
+
+	nd_setup_parameters(nd_nvar,0);
+	prev_nm_free_list = _nm_free_list;
+	_nm_free_list = 0;
+	for ( k = j; k < n; k++ )
+		for (l = j; l < n; l++ )
+			ndv_realloc(dm[k][l],obpe,oadv,oepos);
+	ndv_realloc(d,obpe,oadv,oepos);
+	prev_nm_free_list = 0;
+#if 0
+	GC_gcollect();
+#endif
+}
+
+UINT *nd_det_compute_bound(NDV **dm,int n,int j)
+{
+	UINT *d0,*d1,*d,*t,*r;
+	int k,l;
+
+	d0 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	d1 = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	for ( k = 0; k < nd_wpd; k++ ) d0[k] = 0;
+	for ( k = j; k < n; k++ )
+		for ( l = j; l < n; l++ )
+			if ( dm[k][l] ) {
+				d = ndv_compute_bound(dm[k][l]);
+				ndl_lcm(d,d0,d1);
+				t = d1; d1 = d0; d0 = t;
+			}
+	r = (UINT *)ALLOCA(nd_wpd*sizeof(UINT));
+	for ( k = 0; k < nd_wpd; k++ ) r[k] = d0[k];
+	return r;
+}
+
+DL nd_separate_d(UINT *d,UINT *trans)
+{
+	int n,ntrans,td,i,e;
+	DL a;
+
+	n = nd_nvar; ntrans = n-nd_nalg;
+	ndl_zero(trans);
+	td = 0;
+	for ( i = 0; i < ntrans; i++ ) {
+		e = GET_EXP(d,i);
+		PUT_EXP(trans,i,e);
+		td += MUL_WEIGHT(e,i);
+	}
+	TD(trans) = td;
+	if ( nd_blockmask) ndl_weight_mask(trans);
+	NEWDL(a,nd_nalg);
+	td = 0;
+	for ( ; i < n; i++ ) {
+		e = GET_EXP(d,i);
+		a->d[i-ntrans] = e;
+		td += e;
+	}
+	a->td = td;
+	return a;
+}
+
+NumberField get_numberfield();
+
+int nd_monic(int mod,ND *p)
+{
+	UINT *trans,*t;
+	DL alg;
+	MP mp0,mp;
+	NM m,m0,m1,ma0,ma,mb,mr0,mr;
+	ND r;
+	DL dl;
+	DP nm;
+	NDV ndv;
+	DAlg inv,cd;
+	ND s,c;
+	Q l,mul;
+	N ln;
+	int n,ntrans,i,e,td,is_lc,len;
+	NumberField nf;
+	struct oEGT eg0,eg1;
+
+	if ( !(nf = get_numberfield()) )
+		error("nd_monic : current_numberfield is not set");
+
+	n = nd_nvar; ntrans = n-nd_nalg;
+	/* Q coef -> DAlg coef */
+	NEWNM(ma0); ma = ma0;
+	m = BDY(*p); 
+	is_lc = 1;
+	while ( 1 ) {
+		NEWMP(mp0); mp = mp0;
+		mp->c = (P)CQ(m);
+		mp->dl = nd_separate_d(DL(m),DL(ma));
+		NEWNM(mb);
+		for ( m = NEXT(m); m; m = NEXT(m) ) {
+			alg = nd_separate_d(DL(m),DL(mb));
+			if ( !ndl_equal(DL(ma),DL(mb)) )
+				break;
+			NEXTMP(mp0,mp); mp->c = (P)CQ(m); mp->dl = alg;
+		}
+		NEXT(mp) = 0;
+		MKDP(nd_nalg,mp0,nm);
+		MKDAlg(nm,ONE,cd);
+		if ( is_lc == 1 ) {
+			/* if the lc is a rational number, we have nothing to do */
+			if ( !mp0->dl->td )
+				return 1;
+
+			get_eg(&eg0);
+			invdalg(cd,&inv);
+			get_eg(&eg1); add_eg(&eg_invdalg,&eg0,&eg1);
+			/* check the validity of inv */
+			if ( mod && !rem(NM(inv->dn),mod) )
+				return 0;
+			CA(ma) = nf->one;
+			is_lc = 0;
+			ln = ONEN;
+		} else {
+			muldalg(cd,inv,&CA(ma));	
+			lcmn(ln,NM(CA(ma)->dn),&ln);
+		}
+		if ( m ) {
+			NEXT(ma) = mb; ma = mb;
+		} else {
+			NEXT(ma) = 0;
+			break;
+		}
+	}
+	/* l = lcm(denoms) */
+	NTOQ(ln,1,l);
+	for ( mr0 = 0, m = ma0; m; m = NEXT(m) ) {
+		divq(l,CA(m)->dn,&mul);
+		for ( mp = BDY(CA(m)->nm); mp; mp = NEXT(mp) ) {
+			NEXTNM(mr0,mr);
+			mulq((Q)mp->c,mul,&CQ(mr));
+			dl = mp->dl;
+			td = TD(DL(m));
+			ndl_copy(DL(m),DL(mr));
+			for ( i = ntrans; i < n; i++ ) {
+				e = dl->d[i-ntrans];
+				PUT_EXP(DL(mr),i,e);
+				td += MUL_WEIGHT(e,i);
+			}
+			TD(DL(mr)) = td;
+			if ( nd_blockmask) ndl_weight_mask(DL(mr));
+		}
+	}
+	NEXT(mr) = 0;
+	for ( len = 0, mr = mr0; mr; mr = NEXT(mr), len++ );
+	MKND(NV(*p),mr0,len,r);
+	/* XXX */
+	SG(r) = SG(*p);
+	nd_free(*p);
+	*p = r;
+	return 1;
+}
+
+void nd_set_nalg(int nalg)
+{
+	nd_nalg = nalg;
 }