===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/dalg.c,v
retrieving revision 1.13
retrieving revision 1.17
diff -u -p -r1.13 -r1.17
--- OpenXM_contrib2/asir2000/engine/dalg.c	2006/01/05 00:21:20	1.13
+++ OpenXM_contrib2/asir2000/engine/dalg.c	2017/08/31 02:36:21	1.17
@@ -1,5 +1,5 @@
 /*
- * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.12 2005/10/12 14:43:36 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.16 2013/11/05 02:55:03 noro Exp $
 */
 
 #include "ca.h"
@@ -36,7 +36,7 @@ void setfield_dalg(NODE alist)
 	current_numberfield = nf;
 	vl = 0;
 	for ( t = alist; t; t = NEXT(t) ) {
-		clctalg(BDY((Alg)BDY(t)),&vl1);
+		clctalg((P)BDY((Alg)BDY(t)),&vl1);
 		mergev(ALG,vl,vl1,&vl2); vl = vl2;
 	}
 	for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ ); 
@@ -121,11 +121,11 @@ void qtodalg(Q q,DAlg *r)
 		*r = (DAlg)q;
 	else if ( NID(q) == N_Q ) {
 		if ( INT(q) ) {
-			muldc(CO,nf->one->nm,(P)q,&nm);
+			muldc(CO,nf->one->nm,(Obj)q,&nm);
 			MKDAlg(nm,ONE,*r);
 		} else {
 			NTOQ(NM(q),SGN(q),t);
-			muldc(CO,nf->one->nm,(P)t,&nm);
+			muldc(CO,nf->one->nm,(Obj)t,&nm);
 			NTOQ(DN(q),1,t);
 			MKDAlg(nm,t,*r);
 		}
@@ -330,12 +330,12 @@ void algtodalg(Alg a,DAlg *r)
 		case N_Q:
 			c = (Q)a;
 			if ( INT(c) ) {
-				muldc(CO,nf->one->nm,(P)c,&dp);
+				muldc(CO,nf->one->nm,(Obj)c,&dp);
 				MKDAlg(dp,ONE,*r);
 			} else {
 				NTOQ(NM(c),SGN(c),c1);
 				NTOQ(DN(c),1,d1);
-				muldc(CO,nf->one->nm,(P)c1,&dp);
+				muldc(CO,nf->one->nm,(Obj)c1,&dp);
 				MKDAlg(dp,d1,*r);
 			}
 			break;
@@ -351,7 +351,7 @@ void algtodalg(Alg a,DAlg *r)
 				mulpq(p,(P)nm,&p1); p = p1;
 			}
 			current_spec = dp_current_spec; initd(nf->spec);
-			get_vars(p,&vl);
+			get_vars((Obj)p,&vl);
 			for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
 				v = tvl->v;
 				for ( svl = nf->vl; svl; svl = NEXT(svl) )
@@ -378,10 +378,13 @@ void dalgtoalg(DAlg da,Alg *r)
 
 	if ( !(nf=current_numberfield) )
 		error("dalgtoalg : current_numberfield is not set");
-	dtop(ALG,nf->vl,da->nm,&p);
-	invq(da->dn,&inv);
-	mulpq(p,(P)inv,&p1);
-	MKAlg(p1,*r);
+	if ( !da ) *r = 0;
+	else {
+		dtop(ALG,nf->vl,da->nm,(Obj *)&p);
+		invq(da->dn,&inv);
+		mulpq(p,(P)inv,&p1);
+		MKAlg(p1,*r);
+	}
 }
 
 void simpdalg(DAlg da,DAlg *r)
@@ -399,7 +402,7 @@ void simpdalg(DAlg da,DAlg *r)
 		return;
 	}
 	current_spec = dp_current_spec; initd(nf->spec);
-	dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
+	dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,(P *)&dn);
 	if ( !nm ) *r = 0;
 	else {
 		initd(current_spec);
@@ -432,7 +435,7 @@ void adddalg(DAlg a,DAlg b,DAlg *c)
 		divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
 		NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
 		/* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
-		muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
+		muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
 		current_spec = dp_current_spec; initd(nf->spec);
 		addd(CO,ta,tb,&nm);
 		initd(current_spec);
@@ -468,7 +471,7 @@ void subdalg(DAlg a,DAlg b,DAlg *c)
 		divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
 		NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
 		/* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
-		muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);
+		muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
 		current_spec = dp_current_spec; initd(nf->spec);
 		subd(CO,ta,tb,&nm);
 		initd(current_spec);
@@ -540,7 +543,7 @@ void rmcontdalg(DAlg a, DAlg *r)
 		gcdn(NM(cont),NM(a->dn),&gn);
 		divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
 		divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
-		muldc(CO,u,(P)c,&u1);
+		muldc(CO,u,(Obj)c,&u1);
 		MKDAlg(u1,d,*r);
 	}
 }
@@ -571,7 +574,7 @@ int invdalg(DAlg a,DAlg *c)
 		error("invdalg : division by 0");
 	else if ( NID(a) == N_Q ) {
 		invq((Q)a,&dn); *c = (DAlg)dn;
-		return;
+		return 1;
 	}
 	dim = nf->dim;
 	mb = nf->mb;
@@ -620,7 +623,7 @@ int invdalg(DAlg a,DAlg *c)
 		for ( i = dim-1, mp0 = 0; i >= 0; i-- )
 			if ( solmat[i][0] ) {
 				NEXTMP(mp0,mp);
-				mp->c = (P)solmat[i][0];
+				mp->c = (Obj)solmat[i][0];
 				mp->dl = BDY(mb[i])->dl;
 			}
 		NEXT(mp) = 0; MKDP(n,mp0,u);
@@ -659,7 +662,7 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 		error("invdalg : division by 0");
 	else if ( NID(a) == N_Q ) {
 		invq((Q)a,&dn); *c = (DAlg)dn;
-		return;
+		return 0;
 	}
 	dim = nf->dim;
 	mb = nf->mb;
@@ -667,7 +670,7 @@ 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); }
@@ -708,7 +711,7 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 		}
 	}
 	get_eg(&eg0);
-	rank = generic_gauss_elim_hensel_dalg(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 */
@@ -716,7 +719,7 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 		for ( i = dim-1, mp0 = 0; i >= 0; i-- )
 			if ( solmat[i][0] ) {
 				NEXTMP(mp0,mp);
-				mp->c = (P)solmat[i][0];
+				mp->c = (Obj)solmat[i][0];
 				mp->dl = BDY(mb[i])->dl;
 			}
 		NEXT(mp) = 0; MKDP(n,mp0,u);
@@ -734,17 +737,140 @@ NODE inv_or_split_dalg(DAlg a,DAlg *c)
 			m = mb[cinfo[k]];
 			mp0 = 0;
 			NEXTMP(mp0,mp);
-			chsgnq(dnsol,&dn1); mp->c = (P)dn1;
+			chsgnq(dnsol,&dn1); mp->c = (Obj)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->c = (Obj)solmat[l][k];
 					mp->dl = BDY(mb[rinfo[l]])->dl;
 				}
 			}
 			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,(P *)&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,(P *)&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,(P *)&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 = (Obj)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 = (Obj)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 = (Obj)solmat[l][k];
+					mp->dl = BDY(mb[rinfo[l]])->dl;
+				}
+			}
+			NEXT(mp) = 0; MKDP(nv,mp0,u);
 			NEXTNODE(nd0,nd);
 			BDY(nd) = (pointer)u;
 			NEXT(nd) = 0;