=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/algnum.c,v retrieving revision 1.1.1.1 retrieving revision 1.8 diff -u -p -r1.1.1.1 -r1.8 --- OpenXM_contrib2/asir2000/builtin/algnum.c 1999/12/03 07:39:07 1.1.1.1 +++ OpenXM_contrib2/asir2000/builtin/algnum.c 2004/12/06 01:15:18 1.8 @@ -1,28 +1,74 @@ -/* $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.7 2004/12/02 13:48:43 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(); -#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 *); +void get_algtree(Obj f,VL *r); struct ftab alg_tab[] = { + {"set_field",Pset_field,1}, + {"algtodalg",Palgtodalg,1}, + {"dalgtoalg",Pdalgtoalg,1}, + {"invalg_le",Pinvalg_le,1}, {"defpoly",Pdefpoly,1}, {"newalg",Pnewalg,1}, {"mainalg",Pmainalg,1}, @@ -37,6 +83,22 @@ struct ftab alg_tab[] = { static int UCN,ACNT; +void Pset_field(NODE arg,Q *rp) +{ + setfield_dalg(BDY((LIST)ARG0(arg))); + *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 Pnewalg(arg,rp) NODE arg; Alg *rp; @@ -205,13 +267,14 @@ LIST *rp; Alg b; NODE n0,n; +#if 0 if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R ) vl = 0; else { 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); @@ -220,6 +283,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; } @@ -371,4 +437,229 @@ 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); +} + +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; + } }