===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/algnum.c,v
retrieving revision 1.5
retrieving revision 1.14
diff -u -p -r1.5 -r1.14
--- OpenXM_contrib2/asir2000/builtin/algnum.c	2001/10/09 01:36:05	1.5
+++ OpenXM_contrib2/asir2000/builtin/algnum.c	2013/11/17 17:34:59	1.14
@@ -45,13 +45,19 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.4 2000/12/05 01:24:49 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.13 2007/02/15 02:06:20 noro Exp $ 
 */
 #include "ca.h"
 #include "parse.h"
 
 void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();
 void Palg(), Palgv(), Pgetalgtree();
+void Pinvalg_le();
+void Pset_field(),Palgtodalg(),Pdalgtoalg();
+void Pinv_or_split_dalg();
+void Pdalgtoup();
+void Pget_field_defpoly();
+void Pget_field_generator();
 
 void mkalg(P,Alg *);
 int cmpalgp(P,P);
@@ -60,8 +66,23 @@ void algtorat(Num,Obj *);
 void rattoalg(Obj,Alg *);
 void ptoalgp(P,P *);
 void clctalg(P,VL *);
+void get_algtree(Obj f,VL *r);
+void Pinvalg_chrem();
+void Pdalgtodp();
+void Pdptodalg();
 
 struct ftab alg_tab[] = {
+	{"set_field",Pset_field,-3},
+	{"get_field_defpoly",Pget_field_defpoly,1},
+	{"get_field_generator",Pget_field_generator,1},
+	{"algtodalg",Palgtodalg,1},
+	{"dalgtoalg",Pdalgtoalg,1},
+	{"dalgtodp",Pdalgtodp,1},
+	{"dalgtoup",Pdalgtoup,1},
+	{"dptodalg",Pdptodalg,1},
+	{"inv_or_split_dalg",Pinv_or_split_dalg,1},
+	{"invalg_chrem",Pinvalg_chrem,2},
+	{"invalg_le",Pinvalg_le,1},
 	{"defpoly",Pdefpoly,1},
 	{"newalg",Pnewalg,1},
 	{"mainalg",Pmainalg,1},
@@ -76,6 +97,174 @@ struct ftab alg_tab[] = {
 
 static int UCN,ACNT;
 
+void Pset_field(NODE arg,Q *rp)
+{
+	int ac;
+	NODE a0,a1;
+	VL vl0,vl;
+	struct order_spec *spec;
+
+	if ( (ac = argc(arg)) == 1 )
+		setfield_dalg(BDY((LIST)ARG0(arg)));
+	else if ( ac == 3 ) {
+		a0 = BDY((LIST)ARG0(arg));
+		a1 = BDY((LIST)ARG1(arg));
+		for ( vl0 = 0; a1; a1 = NEXT(a1) ) {
+			NEXTVL(vl0,vl);
+			vl->v = VR((P)BDY(a1));
+		}
+		if ( vl0 ) NEXT(vl) = 0;
+		create_order_spec(0,ARG2(arg),&spec);
+		setfield_gb(a0,vl0,spec);
+	}
+	*rp = 0;
+}
+
+void Palgtodalg(NODE arg,DAlg *rp)
+{
+	algtodalg((Alg)ARG0(arg),rp);
+}
+
+void Pdalgtoalg(NODE arg,Alg *rp)
+{
+	dalgtoalg((DAlg)ARG0(arg),rp);
+}
+
+void Pdalgtodp(NODE arg,LIST *r)
+{
+	NODE b;
+	DP nm;
+	Q dn;
+	DAlg da;
+
+	da = (DAlg)ARG0(arg);
+	nm = da->nm;
+	dn = da->dn;
+	b = mknode(2,nm,dn);
+	MKLIST(*r,b);
+}
+
+void Pdptodalg(NODE arg,DAlg *r)
+{
+	DP d,nm,nm1;
+	MP m;
+	Q c,a;
+	DAlg t;
+	
+	d = (DP)ARG0(arg);
+	if ( !d ) *r = 0;
+	else {
+		for ( m = BDY(d); m; m = NEXT(m) )
+			if ( !INT((Q)m->c) ) break;
+		if ( !m ) {
+			MKDAlg(d,(Q)ONE,t);
+		} else {
+			dp_ptozp(d,&nm);
+			divq((Q)BDY(d)->c,(Q)BDY(nm)->c,&c);
+			NTOQ(NM(c),SGN(c),a);
+			muldc(CO,nm,(P)a,&nm1);
+			NTOQ(DN(c),1,a);
+			MKDAlg(nm1,a,t);
+		}
+		simpdalg(t,r);
+	}
+}
+
+void Pdalgtoup(NODE arg,LIST *r)
+{
+	NODE b;
+	int pos;
+	P up;
+	DP nm;
+	Q dn,q;
+
+	pos = dalgtoup((DAlg)ARG0(arg),&up,&dn);
+	STOQ(pos,q);
+	b = mknode(3,up,dn,q);
+	MKLIST(*r,b);
+}
+
+NODE inv_or_split_dalg(DAlg,DAlg *);
+NumberField	get_numberfield();
+
+void Pget_field_defpoly(NODE arg,DAlg *r)
+{
+	NumberField nf;
+	DP d;
+
+	nf = get_numberfield();
+	d = nf->ps[QTOS((Q)ARG0(arg))];
+	MKDAlg(d,ONE,*r);
+}
+
+void Pget_field_generator(NODE arg,DAlg *r)
+{
+	int index,n,i;
+	DL dl;
+	MP m;
+	DP d;
+
+	index = QTOS((Q)ARG0(arg));
+	n = get_numberfield()->n;
+	NEWDL(dl,n);
+	for ( i = 0; i < n; i++ ) dl->d[i] = 0;
+	dl->d[index] = 1; dl->td = 1;
+	NEWMP(m); m->dl = dl; m->c = (P)ONE; NEXT(m) = 0;
+	MKDP(n,m,d);
+	MKDAlg(d,ONE,*r);
+}
+
+
+void Pinv_or_split_dalg(NODE arg,Obj *rp)
+{
+	NODE gen,t,nd0,nd;
+	LIST list;
+	int l,i,j,n;
+	DP *ps,*ps1,*psw;
+	NumberField nf;
+	DAlg inv;
+	extern struct order_spec *dp_current_spec;
+	struct order_spec *current_spec;
+
+	gen = inv_or_split_dalg((DAlg)ARG0(arg),&inv);
+	if ( !gen )
+		*rp = (Obj)inv;
+	else {
+		nf = get_numberfield();
+		current_spec = dp_current_spec; initd(nf->spec);
+		l = length(gen);
+		n = nf->n;
+		ps = nf->ps;
+		psw = (DP *)ALLOCA((n+l)*sizeof(DP));
+		for ( i = j = 0; i < n; i++ ) {
+			for ( t = gen; t; t = NEXT(t) )
+				if ( dp_redble(ps[i],(DP)BDY(t)) ) break;
+			if ( !t )
+				psw[j++] = ps[i];
+		}
+		nd0  = 0;
+		/* gen[0] < gen[1] < ... */
+		/* psw[0] > psw[1] > ... */
+		for ( i = j-1, t = gen; i >= 0 && t; ) {
+			NEXTNODE(nd0,nd);
+			if ( compd(CO,psw[i],(DP)BDY(t)) > 0 ) {
+				BDY(nd) = BDY(t); t = NEXT(t);
+			} else
+				BDY(nd) = (pointer)psw[i--];
+		}
+		for ( ; i >= 0; i-- ) {
+			NEXTNODE(nd0,nd); BDY(nd) = (pointer)psw[i];
+		}
+		for ( ; t; t = NEXT(t) ) {
+			NEXTNODE(nd0,nd); BDY(nd) = BDY(t);
+		}
+		NEXT(nd) = 0;
+		MKLIST(list,nd0);
+		initd(current_spec);
+		*rp = (Obj)list;
+	}
+}
+
 void Pnewalg(arg,rp)
 NODE arg;
 Alg *rp;
@@ -244,6 +433,7 @@ LIST *rp;
 	Alg b;
 	NODE n0,n;
 
+#if 0
 	if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
 		vl = 0;
 	else {
@@ -259,6 +449,9 @@ LIST *rp;
 				vl = 0; break;
 		}
 	}
+#else
+	get_algtree((Obj)ARG0(arg),&vl);
+#endif
 	for ( n0 = 0; vl; vl = NEXT(vl) ) {
 		NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
 	}
@@ -410,4 +603,310 @@ P p,*r;
 		}
 		NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
 	}
+}
+
+void Pinvalg_chrem(NODE arg,LIST *r)
+{
+	NODE n;
+
+	inva_chrem((P)ARG0(arg),(P)ARG1(arg),&n);
+	MKLIST(*r,n);
+}
+
+void invalg_le(Alg a,LIST *r);
+
+void Pinvalg_le(NODE arg,LIST *r)
+{
+	invalg_le((Alg)ARG0(arg),r);
+}
+
+typedef struct oMono_nf {
+	DP mono;
+	DP nf;
+	Q dn;
+} *Mono_nf;
+
+void invalg_le(Alg a,LIST *r)
+{
+	Alg inv;
+	MAT mobj,sol;
+	int *rinfo,*cinfo;
+	P p,dn,dn1,ap;
+	VL vl,tvl;
+	Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol;
+	int i,j,n,len,k;
+	MP mp,mp0;
+	DP dp,nm,nm1,m,d,u,u1;
+	NODE b,b1,hlist,mblist,t,s,rev0,rev,hist;
+	DP *ps;
+	struct order_spec *spec;
+	Mono_nf h,h1;
+	N nq,nr,nl,ng;
+	Q **mat,**solmat;
+	Q *w;
+	int *wi;
+
+	ap = (P)BDY(a);
+	asir_assert(ap,O_P,"invalg_le");
+
+	/* collecting algebraic numbers */
+	clctalg(ap,&vl);
+
+	/* setup */
+	ptozp(ap,1,&c,&p);
+	STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec);
+	for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
+	ps = (DP *)ALLOCA(n*sizeof(DP));
+
+	/* conversion to DP */
+	for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) {
+		ptod(ALG,vl,tvl->v->attr,&ps[i]);
+	}
+	ptod(ALG,vl,p,&dp);
+	/* index list */
+	for ( b = 0, i = 0; i < n; i++ ) {
+		STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
+	}
+	/* simplification */
+	dp_true_nf(b,dp,ps,1,&nm,&dn);
+
+	/* construction of NF table */
+
+	/* stdmono: <<0,...,0>> < ... < max */
+	for ( hlist = 0, i = 0; i < n; i++ ) {
+		MKNODE(b1,(pointer)ps[i],hlist); hlist = b1;
+	}
+	dp_mbase(hlist,&rev0);
+	for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) {
+		MKNODE(b1,BDY(rev),mblist); mblist = b1;
+	}
+	dn0 = ONE;
+	for ( hist = 0, t = mblist; t; t = NEXT(t) ) {
+		/* searching a predecessor */
+		for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) {
+			h = (Mono_nf)BDY(s);
+			if ( dp_redble(m,h->mono) )
+				break;
+		}
+		h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf));
+		if ( s ) {
+			dp_subd(m,h->mono,&d);
+			muld(CO,d,h->nf,&u);
+			dp_true_nf(b,u,ps,1,&nm1,&dn1);
+			mulq(h->dn,(Q)dn1,&h1->dn);
+		} else {
+			muld(CO,m,nm,&u);
+			dp_true_nf(b,u,ps,1,&nm1,&dn1);
+			h1->dn = (Q)dn1;
+		}
+		h1->mono = m;
+		h1->nf = nm1;
+		MKNODE(b1,(pointer)h1,hist); hist = b1;
+
+		/* dn0 = LCM(dn0,h1->dn) */
+		gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr);
+		muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0);
+	}
+	/* create a matrix */
+	len = length(mblist);
+	MKMAT(mobj,len,len+1);
+	mat = (Q **)BDY(mobj);
+	mat[len-1][len] = dn0;
+	for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) {
+		h = (Mono_nf)BDY(t);
+		nm1 = h->nf;
+		divq((Q)dn0,h->dn,&mul);
+		for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) )
+			if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) {
+				mulq(mul,(Q)mp->c,&mat[i][j]);
+				mp = NEXT(mp);
+			}
+	}
+#if 0
+	w = (Q *)ALLOCA((len+1)*sizeof(Q));
+	wi = (int *)ALLOCA((len+1)*sizeof(int));
+	for ( i = 0; i < len; i++ ) {
+		for ( j = 0, k = 0; j <= len; j++ )
+			if ( mat[i][j] ) {
+				w[k] = mat[i][j];
+				wi[k] = j;
+				k++;
+			}
+		removecont_array(w,k);
+		for ( j = 0; j < k; j++ )
+			mat[i][wi[j]] = w[j];
+	}
+#endif
+	generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
+	solmat = (Q **)BDY(sol);
+	for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) )
+		if ( solmat[i][0] ) {
+			NEXTMP(mp0,mp);
+			mp->c = (P)solmat[i][0];
+			mp->dl = BDY((DP)BDY(t))->dl;
+		}
+	NEXT(mp) = 0; MKDP(n,mp0,u);
+	dp_ptozp(u,&u1);
+	divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont);
+	dtop(ALG,vl,u1,&ap);
+	MKAlg(ap,inv);
+	mulq(dnsol,(Q)dn,&c1);
+	mulq(c1,c,&c2);
+	divq(c2,cont,&c3);
+	b = mknode(2,inv,c3);
+	MKLIST(*r,b);
+}
+
+void get_algtree(Obj f,VL *r)
+{
+	VL vl1,vl2,vl3;
+	Obj t;
+	DCP dc;
+	NODE b;
+	pointer *a;
+	pointer **m;
+	int len,row,col,i,j,l;
+
+	if ( !f ) *r = 0;
+	else
+		switch ( OID(f) ) {
+			case O_N:
+				if ( NID((Num)f) != N_A ) *r = 0;
+				else  {
+					t = BDY((Alg)f);
+					switch ( OID(t) ) {
+						case O_P:
+							clctalg((P)t,r); break;
+						case O_R:
+							clctalg(NM((R)t),&vl1);
+							clctalg(DN((R)t),&vl2);
+							mergev(ALG,vl1,vl2,r); break;
+						default:
+							*r = 0; break;
+					}
+				}
+				break;
+			case O_P:
+				vl1 = 0;
+				for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
+					get_algtree((Obj)COEF(dc),&vl2);
+					mergev(ALG,vl1,vl2,&vl3);
+					vl1 = vl3;
+				}
+				*r = vl1;
+				break;
+			case O_R:
+				get_algtree((Obj)NM((R)f),&vl1);
+				get_algtree((Obj)DN((R)f),&vl2);
+				mergev(ALG,vl1,vl2,r);
+				break;
+			case O_LIST:
+				vl1 = 0;
+				for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
+					get_algtree((Obj)BDY(b),&vl2);
+					mergev(ALG,vl1,vl2,&vl3);
+					vl1 = vl3;
+				}
+				*r = vl1;
+				break;
+			case O_VECT:
+				vl1 = 0;
+				l = ((VECT)f)->len;
+				a = BDY((VECT)f);
+				for ( i = 0; i < l; i++ ) {
+					get_algtree((Obj)a[i],&vl2);
+					mergev(ALG,vl1,vl2,&vl3);
+					vl1 = vl3;
+				}
+				*r = vl1;
+				break;
+			case O_MAT:
+				vl1 = 0;
+				row = ((MAT)f)->row; col = ((MAT)f)->col;
+				m = BDY((MAT)f);
+				for ( i = 0; i < row; i++ )
+					for ( j = 0; j < col; j++ ) {
+						get_algtree((Obj)m[i][j],&vl2);
+						mergev(ALG,vl1,vl2,&vl3);
+						vl1 = vl3;
+					}
+				*r = vl1;
+				break;
+			default:
+				*r = 0;
+				break;
+		}
+}
+
+void algobjtorat(Obj f,Obj *r)
+{
+	Obj t;
+	DCP dc,dcr,dcr0;
+	P p,nm,dn;
+	R rat;
+	NODE b,s,s0;
+	VECT v;
+	MAT mat;
+	LIST list;
+	pointer *a;
+	pointer **m;
+	int len,row,col,i,j,l;
+
+	if ( !f ) *r = 0;
+	else
+		switch ( OID(f) ) {
+			case O_N:
+				algtorat((Num)f,r);
+				break;
+			case O_P:
+				dcr0 = 0;
+				for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
+					NEXTDC(dcr0,dcr);
+					algobjtorat((Obj)COEF(dc),&t);
+					COEF(dcr) = (P)t;
+					DEG(dcr) = DEG(dc);
+				}
+				NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p;
+				break;
+			case O_R:
+				algobjtorat((Obj)NM((R)f),&t); nm = (P)t;
+				algobjtorat((Obj)DN((R)f),&t); dn = (P)t;
+				MKRAT(nm,dn,0,rat); *r = (Obj)rat;
+				break;
+			case O_LIST:
+				s0 = 0;
+				for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
+					NEXTNODE(s0,s);
+					algobjtorat((Obj)BDY(b),&t);
+					BDY(s) = (pointer)t;
+				}
+				NEXT(s) = 0;
+				MKLIST(list,s0);
+				*r = (Obj)list;
+				break;
+			case O_VECT:
+				l = ((VECT)f)->len;
+				a = BDY((VECT)f);
+				MKVECT(v,l);
+				for ( i = 0; i < l; i++ ) {
+					algobjtorat((Obj)a[i],&t);
+					BDY(v)[i] = (pointer)t;
+				}
+				*r = (Obj)v;
+				break;
+			case O_MAT:
+				row = ((MAT)f)->row; col = ((MAT)f)->col;
+				m = BDY((MAT)f);
+				MKMAT(mat,row,col);
+				for ( i = 0; i < row; i++ )
+					for ( j = 0; j < col; j++ ) {
+						algobjtorat((Obj)m[i][j],&t);
+						BDY(mat)[i][j] = (pointer)t;
+					}
+				*r = (Obj)mat;
+				break;
+			default:
+				*r = f;
+				break;
+		}
 }