===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dalg.c,v
retrieving revision 1.8
retrieving revision 1.14
diff -u -p -r1.8 -r1.14
--- OpenXM_contrib2/asir2000/engine/dalg.c	2005/01/23 14:03:47	1.8
+++ OpenXM_contrib2/asir2000/engine/dalg.c	2006/10/26 10:49:17	1.14
@@ -1,5 +1,5 @@
 /*
- * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.7 2004/12/10 02:45:02 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.13 2006/01/05 00:21:20 noro Exp $
 */
 
 #include "ca.h"
@@ -65,6 +65,48 @@ void setfield_dalg(NODE alist)
 		mb[dim-i-1] = (DP)BDY(t);	
 }
 
+void setfield_gb(NODE gb,VL vl,struct order_spec *spec)
+{
+	NumberField nf;
+	VL vl1,vl2;
+	int n,i,dim;
+	Alg *gen;
+	P *defpoly;
+	P p;
+	Q c,iq,two;
+	DP *ps,*mb;
+	DP one;
+	NODE t,b,b1,b2,hlist,mblist;
+	struct order_spec *current_spec;
+
+	nf = (NumberField)MALLOC(sizeof(struct oNumberField));
+	current_numberfield = nf;
+	for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );
+	nf->n = n;
+	nf->psn = length(gb);
+	nf->vl = vl;
+	nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));
+	nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));
+	current_spec = dp_current_spec;
+	nf->spec = spec;
+	initd(nf->spec);	
+	for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {
+		ptozp((P)BDY(t),1,&c,&defpoly[i]);
+		ptod(CO,vl,defpoly[i],&ps[i]);
+		STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
+		MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
+	}
+	ptod(ALG,vl,(P)ONE,&one);
+	MKDAlg(one,ONE,nf->one);
+	nf->ind = b;
+	dp_mbase(hlist,&mblist);
+	initd(current_spec);
+	nf->dim = dim = length(mblist);
+	nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
+	for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
+		mb[dim-i-1] = (DP)BDY(t);	
+}
+
 void qtodalg(Q q,DAlg *r)
 {
 	NumberField nf;
@@ -294,7 +336,7 @@ void algtodalg(Alg a,DAlg *r)
 				NTOQ(NM(c),SGN(c),c1);
 				NTOQ(DN(c),1,d1);
 				muldc(CO,nf->one->nm,(P)c1,&dp);
-				MKDAlg(dp,c1,*r);
+				MKDAlg(dp,d1,*r);
 			}
 			break;
 		case N_A:
@@ -358,10 +400,13 @@ void simpdalg(DAlg da,DAlg *r)
 	}
 	current_spec = dp_current_spec; initd(nf->spec);
 	dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
-	initd(current_spec);
-	mulq(da->dn,dn,&dn1);
-	MKDAlg(nm,dn1,d);
-	rmcontdalg(d,r);
+	if ( !nm ) *r = 0;
+	else {
+		initd(current_spec);
+		mulq(da->dn,dn,&dn1);
+		MKDAlg(nm,dn1,d);
+		rmcontdalg(d,r);
+	}
 }
 
 void adddalg(DAlg a,DAlg b,DAlg *c)
@@ -606,6 +651,7 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 	struct order_spec *current_spec;
 	struct oEGT eg0,eg1;
 	extern struct oEGT eg_le;
+	extern int DP_Print;
 
 	if ( !(nf=current_numberfield) )
 		error("invdalg : current_numberfield is not set");
@@ -621,24 +667,30 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 	ln = ONEN;
 	dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
 	MKDAlg(u,ONE,a0);
-	simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
+	simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
 	current_spec = dp_current_spec; initd(nf->spec);
 	for ( i = 0; i < dim; i++ ) {
+		if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
 		m = mb[i];
 		for ( j = i-1; j >= 0; j-- )
 			if ( dp_redble(m,mb[j]) )
 				break;
 		if ( j >= 0 ) {
 			dp_subd(m,mb[j],&d);
-			muld(CO,d,simp[j]->nm,&u);
-			MKDAlg(u,simp[j]->dn,t);
-			simpdalg(t,&simp[i]);
+			if ( simp[j] ) {
+				muld(CO,d,simp[j]->nm,&u);
+				MKDAlg(u,simp[j]->dn,t);
+				simpdalg(t,&simp[i]);
+			} else
+				simp[i] = 0;
 		} else {
 			MKDAlg(m,ONE,t);
 			muldalg(t,a0,&simp[i]);
 		}
-		gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
-		muln(NM(simp[i]->dn),qn,&ln);
+		if ( simp[i] ) {
+			gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
+			muln(NM(simp[i]->dn),qn,&ln);
+		}
 	}
 	initd(current_spec);
 	NTOQ(ln,1,dn);
@@ -646,15 +698,17 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 	mat = (Q **)BDY(mobj);
 	mulq(dn,a->dn,&mat[0][dim]);
 	for ( j = 0; j < dim; j++ ) {
-		divq(dn,simp[j]->dn,&mul);
-		for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
-			if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
-				mulq(mul,(Q)mp->c,&mat[i][j]);
-				mp = NEXT(mp);
-			}
+		if ( simp[j] ) {
+			divq(dn,simp[j]->dn,&mul);
+			for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
+				if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
+					mulq(mul,(Q)mp->c,&mat[i][j]);
+					mp = NEXT(mp);
+				}
+		}
 	}
 	get_eg(&eg0);
-	rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
+	rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
 	get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
 	if ( cinfo[0] == dim ) {
 		/* the input is invertible */
@@ -672,21 +726,135 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 		return 0;
 	} else {
 		/* the input is not invertible */
-		nparam = (dim+1)-rank;
-		/* the index 'dim' should not be in cinfo[] */
+		nparam = sol->col;
 		solmat = (Q **)BDY(sol);
-		for ( k = 0; k < nparam; k++ )
-			if ( cinfo[k] == dim )
-				error("invdalg : cannot happen");
 		nd0 = 0;
 		for ( k = 0; k < nparam; k++ ) {
+			/* construct a new basis element */
 			m = mb[cinfo[k]];
-			for ( ndt = nd0; ndt; ndt = NEXT(ndt) ) {
-				if ( dp_redble(m,(DP)BDY(ndt)) ) break;
+			mp0 = 0;
+			NEXTMP(mp0,mp);
+			chsgnq(dnsol,&dn1); mp->c = (P)dn1;
+			mp->dl = BDY(m)->dl;
+			/* skip the last parameter */
+			for ( l = rank-2; l >= 0; l-- ) {
+				if ( solmat[l][k] ) {
+					NEXTMP(mp0,mp);
+					mp->c = (P)solmat[l][k];
+					mp->dl = BDY(mb[rinfo[l]])->dl;
+				}
 			}
-			/* skip a redundunt basis element */
-			if ( ndt ) continue;
+			NEXT(mp) = 0; MKDP(n,mp0,u);
+			NEXTNODE(nd0,nd);
+			BDY(nd) = (pointer)u;
+			NEXT(nd) = 0;
+		}
+		NEXT(nd) = 0;
+		return nd0;
+	}
+}
+
+NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv)
+{
+	int dim,n,i,j,k,l,nv;
+	DP *mb,*ps;
+	DP m,d,u,nm;
+	N ln,gn,qn;
+	DAlg *simp;
+	DAlg a0,r;
+	Q dn,dnsol,mul,nmc,dn1,iq;
+	MAT mobj,sol;
+	Q **mat,**solmat;
+	MP mp0,mp;
+	int *rinfo,*cinfo;
+	int rank,nparam;
+	NODE nd0,nd,ndt,ind,indt,t,mblist;
+	struct oEGT eg0,eg1;
+	extern struct oEGT eg_le;
+	extern int DP_Print;
+	initd(spec);
+	dp_ptozp(f,&u); f = u;
+
+	n = length(gb);
+	ps = (DP *)MALLOC(n*sizeof(DP));
+	for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) {
+		ps[i] = (DP)BDY(t);
+		NEXTNODE(ind,indt);
+		STOQ(i,iq); BDY(indt) = iq;
+	}
+	if ( ind ) NEXT(indt) = 0;
+	dp_true_nf(ind,f,ps,1,&nm,&dn);
+	if ( !nm ) error("dp_inv_or_split : input is 0");
+	f = nm;
+
+	dp_mbase(gb,&mblist);
+	dim = length(mblist);
+	mb = (DP *)MALLOC(dim*sizeof(DP));
+	for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) )
+		mb[dim-i-1] = (DP)BDY(t);
+	nv = mb[0]->nv;
+	ln = ONEN;
+	simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
+	for ( i = 0; i < dim; i++ ) {
+		if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
+		m = mb[i];
+		for ( j = i-1; j >= 0; j-- )
+			if ( dp_redble(m,mb[j]) )
+				break;
+		if ( j >= 0 ) {
+			dp_subd(m,mb[j],&d);
+			if ( simp[j] ) {
+				muld(CO,d,simp[j]->nm,&u);
+				dp_true_nf(ind,u,ps,1,&nm,&dn);
+				mulq(simp[j]->dn,dn,&dn1);
+				MKDAlg(nm,dn1,simp[i]);
+			} else
+				simp[i] = 0;
+		} else {
+			dp_true_nf(ind,f,ps,1,&nm,&dn);
+			MKDAlg(nm,dn,simp[i]);
+		}
+		if ( simp[i] ) {
+			gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
+			muln(NM(simp[i]->dn),qn,&ln);
+		}
+	}
+	NTOQ(ln,1,dn);
+	MKMAT(mobj,dim,dim+1);
+	mat = (Q **)BDY(mobj);
+	mat[0][dim] = dn;
+	for ( j = 0; j < dim; j++ ) {
+		if ( simp[j] ) {
+			divq(dn,simp[j]->dn,&mul);
+			for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
+				if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) {
+					mulq(mul,(Q)mp->c,&mat[i][j]);
+					mp = NEXT(mp);
+				}
+		}
+	}
+	get_eg(&eg0);
+	rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
+	get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
+	if ( cinfo[0] == dim ) {
+		/* the input is invertible */
+		solmat = (Q **)BDY(sol);
+		for ( i = dim-1, mp0 = 0; i >= 0; i-- )
+			if ( solmat[i][0] ) {
+				NEXTMP(mp0,mp);
+				mp->c = (P)solmat[i][0];
+				mp->dl = BDY(mb[i])->dl;
+			}
+		NEXT(mp) = 0; MKDP(nv,mp0,*inv);
+		return 0;
+	} else {
+		/* the input is not invertible */
+		nparam = sol->col;
+		solmat = (Q **)BDY(sol);
+		nd0 = 0;
+		for ( k = 0; k < nparam; k++ ) {
 			/* construct a new basis element */
+			m = mb[cinfo[k]];
 			mp0 = 0;
 			NEXTMP(mp0,mp);
 			chsgnq(dnsol,&dn1); mp->c = (P)dn1;
@@ -699,7 +867,7 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 					mp->dl = BDY(mb[rinfo[l]])->dl;
 				}
 			}
-			NEXT(mp) = 0; MKDP(n,mp0,u);
+			NEXT(mp) = 0; MKDP(nv,mp0,u);
 			NEXTNODE(nd0,nd);
 			BDY(nd) = (pointer)u;
 			NEXT(nd) = 0;
@@ -775,3 +943,48 @@ int cmpdalg(DAlg a,DAlg b)
 	else
 		return SGN((Q)BDY(c->nm)->c);
 }
+
+/* convert da to a univariate poly; return the position of variable */
+
+int dalgtoup(DAlg da,P *up,Q *dn)
+{
+	int nv,i,hi,current_d;
+	DCP dc0,dc;
+	MP h,mp0,mp,t;
+	DL hd,d;
+	DP c;
+	DAlg cc;
+	P v;
+
+	nv = da->nm->nv;
+	h = BDY(da->nm);
+	*dn = da->dn;
+	hd = h->dl;
+	for ( i = 0; i < nv; i++ )
+		if ( hd->d[i] ) break;
+	hi = i;
+	current_d = hd->d[i];
+	dc0 = 0;
+	mp0 = 0;
+	for ( t = h; t; t = NEXT(t) ) {
+		NEWDL(d,nv);
+		for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
+		for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
+		d->td = t->dl->td - t->dl->d[hi];
+		if ( t->dl->d[hi] != current_d ) {
+			NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
+			NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
+			current_d = t->dl->d[hi];
+			mp0 = 0;
+		}
+		NEXTMP(mp0,mp);
+		mp->c = t->c; mp->dl = d;
+	}
+	NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
+	NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
+	NEXT(dc) = 0;
+	makevar("x",&v);
+	MKP(VR(v),dc0,*up);
+	return hi;
+}
+