===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/algnum.c,v
retrieving revision 1.1.1.1
retrieving revision 1.6
diff -u -p -r1.1.1.1 -r1.6
--- OpenXM_contrib2/asir2000/builtin/algnum.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/algnum.c	2004/12/01 08:49:42	1.6
@@ -1,28 +1,69 @@
-/* $OpenXM: OpenXM/src/asir99/builtin/algnum.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
+/*
+ * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
+ * All rights reserved.
+ * 
+ * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
+ * non-exclusive and royalty-free license to use, copy, modify and
+ * redistribute, solely for non-commercial and non-profit purposes, the
+ * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
+ * conditions of this Agreement. For the avoidance of doubt, you acquire
+ * only a limited right to use the SOFTWARE hereunder, and FLL or any
+ * third party developer retains all rights, including but not limited to
+ * copyrights, in and to the SOFTWARE.
+ * 
+ * (1) FLL does not grant you a license in any way for commercial
+ * purposes. You may use the SOFTWARE only for non-commercial and
+ * non-profit purposes only, such as academic, research and internal
+ * business use.
+ * (2) The SOFTWARE is protected by the Copyright Law of Japan and
+ * international copyright treaties. If you make copies of the SOFTWARE,
+ * with or without modification, as permitted hereunder, you shall affix
+ * to all such copies of the SOFTWARE the above copyright notice.
+ * (3) An explicit reference to this SOFTWARE and its copyright owner
+ * shall be made on your publication or presentation in any form of the
+ * results obtained by use of the SOFTWARE.
+ * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
+ * for such modification or the source code of the modified part of the
+ * SOFTWARE.
+ * 
+ * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
+ * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
+ * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
+ * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
+ * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
+ * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
+ * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
+ * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
+ * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
+ * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
+ * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
+ * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
+ * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
+ * 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.5 2001/10/09 01:36:05 noro Exp $ 
+*/
 #include "ca.h"
 #include "parse.h"
 
 void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();
 void Palg(), Palgv(), Pgetalgtree();
+void Pinvalg_le();
 
-#if defined(THINK_C)
 void mkalg(P,Alg *);
 int cmpalgp(P,P);
 void algptop(P,P *);
 void algtorat(Num,Obj *);
 void rattoalg(Obj,Alg *);
 void ptoalgp(P,P *);
-#else
-void mkalg();
-int cmpalgp();
-void algptop();
-void algtorat();
-void rattoalg();
-void ptoalgp();
-void clctalg();
-#endif
+void clctalg(P,VL *);
 
 struct ftab alg_tab[] = {
+	{"invalg_le",Pinvalg_le,1},
 	{"defpoly",Pdefpoly,1},
 	{"newalg",Pnewalg,1},
 	{"mainalg",Pmainalg,1},
@@ -211,7 +252,7 @@ LIST *rp;
 		t = BDY((Alg)a);
 		switch ( OID(t) ) {
 			case O_P:
-				clctalg(t,&vl); break;
+				clctalg((P)t,&vl); break;
 			case O_R:
 				clctalg(NM((R)t),&vl1);
 				clctalg(DN((R)t),&vl2);
@@ -371,4 +412,148 @@ P p,*r;
 		}
 		NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
 	}
+}
+
+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);
 }