=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/parse/puref.c,v retrieving revision 1.1.1.1 retrieving revision 1.14 diff -u -p -r1.1.1.1 -r1.14 --- OpenXM_contrib2/asir2000/parse/puref.c 1999/12/03 07:39:12 1.1.1.1 +++ OpenXM_contrib2/asir2000/parse/puref.c 2018/03/28 05:27:22 1.14 @@ -1,12 +1,60 @@ -/* $OpenXM: OpenXM/src/asir99/parse/puref.c,v 1.1.1.1 1999/11/10 08:12:34 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/parse/puref.c,v 1.13 2018/03/27 06:29:19 noro Exp $ +*/ #include "ca.h" #include "parse.h" +void instoobj(PFINS ins,Obj *rp); + NODE pflist; -void searchpf(name,fp) -char *name; -FUNC *fp; +void searchpf(char *name,FUNC *fp) { NODE node; PF pf; @@ -17,15 +65,13 @@ FUNC *fp; pf = (PF)node->body; *fp = t = (FUNC)MALLOC(sizeof(struct oFUNC)); t->name = name; t->id = A_PURE; t->argc = pf->argc; - t->f.puref = pf; + t->f.puref = pf; t->fullname = name; return; } *fp = 0; } -void searchc(name,fp) -char *name; -FUNC *fp; +void searchc(char *name,FUNC *fp) { NODE node; PF pf; @@ -37,20 +83,14 @@ FUNC *fp; pf = (PF)node->body; *fp = t = (FUNC)MALLOC(sizeof(struct oFUNC)); t->name = name; t->id = A_PURE; t->argc = pf->argc; - t->f.puref = pf; + t->f.puref = pf; t->fullname = name; return; } *fp = 0; } -void mkpf(name,body,argc,args,parif,libmf,simp,pfp) -char *name; -Obj body; -int argc; -V *args; -int (*parif)(),(*simp)(); -double (*libmf)(); -PF *pfp; +void mkpf(char *name,Obj body,int argc,V *args, + int (*parif)(),double (*libmf)(), int (*simp)(),PF *pfp) { PF pf; NODE node; @@ -75,10 +115,7 @@ PF *pfp; each arg. */ -void mkpfins(pf,args,vp) -PF pf; -V *args; -V *vp; +void mkpfins(PF pf,V *args,V *vp) { V v; PFINS ins; @@ -99,10 +136,7 @@ V *vp; /* the same as above. Argements are given as an array of Obj */ -void _mkpfins(pf,args,vp) -PF pf; -Obj *args; -V *vp; +void _mkpfins(PF pf,Obj *args,V *vp) { V v; PFINS ins; @@ -122,11 +156,7 @@ V *vp; /* the same as above. darray is also given */ -void _mkpfins_with_darray(pf,args,darray,vp) -PF pf; -Obj *args; -int *darray; -V *vp; +void _mkpfins_with_darray(PF pf,Obj *args,int *darray,V *vp) { V v; PFINS ins; @@ -144,9 +174,7 @@ V *vp; appendpfins(v,vp); } -void appendpfins(v,vp) -V v; -V *vp; +void appendpfins(V v,V *vp) { PF fdef; PFAD ad,tad; @@ -157,7 +185,7 @@ V *vp; for ( node = fdef->ins; node; node = NEXT(node) ) { for ( i = 0, tad = ((PFINS)((V)node->body)->priv)->ad; i < fdef->argc; i++ ) - if ( (ad[i].d != tad[i].d) || compr(CO,ad[i].arg,tad[i].arg) ) + if ( (ad[i].d != tad[i].d) || !equalr(CO,ad[i].arg,tad[i].arg) ) break; if ( i == fdef->argc ) { *vp = (V)node->body; @@ -168,9 +196,7 @@ V *vp; fdef->ins = node; appendvar(CO,v); *vp = v; } -void duppfins(v,vp) -V v; -V *vp; +void duppfins(V v,V *vp) { V tv; PFINS tins; @@ -183,10 +209,7 @@ V *vp; *vp = tv; } -void derivvar(vl,pf,v,a) -VL vl; -V pf,v; -Obj *a; +void derivvar(VL vl,V pf,V v,Obj *a) { Obj t,s,u,w,u1; P p; @@ -216,10 +239,7 @@ Obj *a; } } -void derivr(vl,a,v,b) -VL vl; -V v; -Obj a,*b; +void derivr(VL vl,Obj a,V v,Obj *b) { VL tvl,svl; Obj r,s,t,u,nm,dn,dnm,ddn,m; @@ -266,14 +286,61 @@ Obj a,*b; } } -void substr(vl,partial,a,v,b,c) -VL vl; -int partial; -Obj a; -V v; -Obj b; -Obj *c; +void simple_derivr(VL vl,Obj a,V v,Obj *b) { + Obj r,s,t,u,nm,dn; + + if ( !a || NUM(a) ) + *b = 0; + else + switch ( OID(a) ) { + case O_P: + pderivr(vl,a,v,b); break; + case O_R: + nm = (Obj)NM((R)a); dn = (Obj)DN((R)a); + /* (nm/dn)' = nm'/dn - dn'/dn*nm/dn */ + pderivr(vl,nm,v,&s); divr(vl,s,dn,&u); reductr(vl,u,&t); + pderivr(vl,dn,v,&s); divr(vl,s,dn,&u); reductr(vl,u,&s); mulr(vl,s,a,&u); + subr(vl,t,u,&s); reductr(vl,s,b); + break; + default: + error("simple_derivr : invalid argument"); + } +} + +int obj_is_dependent(Obj a,V v) +{ + if ( !a || OID(a) <= O_N ) return 0; + else if ( OID(a) == O_P ) return poly_is_dependent((P)a,v); + else if ( OID(a) == O_R ) return poly_is_dependent(NM((R)a),v) + || poly_is_dependent(DN((R)a),v); + else + error("obj_is_dependent : not implemented"); +} + +int poly_is_dependent(P p,V v) +{ + DCP dc; + + if ( !p || OID(p) <= O_N ) return 0; + else if ( v == VR(p) ) return 1; + else { + for ( dc = DC(p); dc; dc = NEXT(dc) ) + if ( poly_is_dependent(COEF(dc),v) ) return 1; + return 0; + } +} + +void gen_pwrr(VL vl,Obj a,Obj d,Obj *r) +{ + if ( INT(d) ) + pwrr(vl,a,d,r); + else + mkpow(vl,a,d,r); +} + +void substr(VL vl,int partial,Obj a,V v,Obj b,Obj *c) +{ Obj nm,dn,t; if ( !a ) @@ -295,13 +362,7 @@ Obj *c; } } -void substpr(vl,partial,p,v0,p0,pr) -VL vl; -int partial; -V v0; -Obj p; -Obj p0; -Obj *pr; +void substpr(VL vl,int partial,Obj p,V v0,Obj p0,Obj *pr) { P x; Obj t,m,c,s,a; @@ -332,7 +393,8 @@ Obj *pr; for ( c = 0, dc = DC((P)p); dc; dc = NEXT(dc) ) { substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t); if ( DEG(dc) ) { - pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m); + gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s); + mulr(vl,s,t,&m); addr(vl,m,c,&a); c = a; } else { addr(vl,t,c,&a); c = a; @@ -346,7 +408,8 @@ Obj *pr; else c = (Obj)COEF(dc); for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) { - subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); + subq(d,DEG(dc),(Q *)&t); + gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); if ( !partial ) substpr(vl,partial,(Obj)COEF(dc),v0,p0,&t); else @@ -354,18 +417,15 @@ Obj *pr; addr(vl,m,t,&c); } if ( d ) { - pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); + gen_pwrr(vl,p0,(Obj)d,&t); + mulr(vl,t,c,&m); c = m; } *pr = c; } } -void evalr(vl,a,prec,c) -VL vl; -Obj a; -int prec; -Obj *c; +void evalr(VL vl,Obj a,int prec,Obj *c) { Obj nm,dn; @@ -387,11 +447,7 @@ Obj *c; } } -void evalp(vl,p,prec,pr) -VL vl; -P p; -int prec; -P *pr; +void evalp(VL vl,P p,int prec,P *pr) { P t; DCP dc,dcr0,dcr; @@ -419,11 +475,7 @@ P *pr; } } -void evalv(vl,v,prec,rp) -VL vl; -V v; -int prec; -Obj *rp; +void evalv(VL vl,V v,int prec,Obj *rp) { PFINS ins,tins; PFAD ad,tad; @@ -444,13 +496,11 @@ Obj *rp; } } -void evalins(ins,prec,rp) -PFINS ins; -int prec; -Obj *rp; +void evalins(PFINS ins,int prec,Obj *rp) { PF pf; - PFAD ad; + PFINS tins; + PFAD ad,tad; int i; Q q; V v; @@ -458,14 +508,18 @@ Obj *rp; NODE n0,n; pf = ins->pf; ad = ins->ad; + tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD)); + tins->pf = pf; tad = tins->ad; + for ( i = 0; i < pf->argc; i++ ) { + tad[i].d = ad[i].d; evalr(CO,ad[i].arg,prec,&tad[i].arg); + } for ( i = 0; i < pf->argc; i++ ) - if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) ) - break; + if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break; if ( (i != pf->argc) || !pf->pari ) { - instov(ins,&v); MKV(v,x); *rp = (Obj)x; + instoobj(tins,rp); } else { for ( n0 = 0, i = 0; i < pf->argc; i++ ) { - NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg; + NEXTNODE(n0,n); BDY(n) = (pointer)tad[i].arg; } if ( prec ) { NEXTNODE(n0,n); STOQ(prec,q); BDY(n) = (pointer)q; @@ -476,28 +530,31 @@ Obj *rp; } } -void devalins(PFINS,Obj *); -void devalv(VL,V,Obj *); -void devalp(VL,P,P *); -void devalr(VL,Obj,Obj *); - -void devalr(vl,a,c) -VL vl; -Obj a; -Obj *c; +void devalr(VL vl,Obj a,Obj *c) { Obj nm,dn; double d; - Real r; + Real r,re,im; + C z; + int nid; if ( !a ) *c = 0; else { switch ( OID(a) ) { case O_N: - d = ToReal(a); - MKReal(d,r); - *c = (Obj)r; + nid = NID((Num)a); + if ( nid == N_C ) { + d = ToReal(((C)a)->r); MKReal(d,re); + d = ToReal(((C)a)->i); MKReal(d,im); + reimtocplx(re,im,&z); + *c = (Obj)z; + } else if ( nid == N_Q || nid == N_R || nid == N_B ) { + d = ToReal(a); + MKReal(d,r); + *c = (Obj)r; + } else + error("devalr : unsupported"); break; case O_P: devalp(vl,(P)a,(P *)c); break; @@ -512,10 +569,7 @@ Obj *c; } } -void devalp(vl,p,pr) -VL vl; -P p; -P *pr; +void devalp(VL vl,P p,P *pr) { P t; DCP dc,dcr0,dcr; @@ -558,16 +612,12 @@ P *pr; } } -void devalv(vl,v,rp) -VL vl; -V v; -Obj *rp; +void devalv(VL vl,V v,Obj *rp) { PFINS ins,tins; PFAD ad,tad; PF pf; P t; - Obj s; int i; if ( (vid)v->attr != V_PF ) { @@ -583,76 +633,125 @@ Obj *rp; } } -void devalins(ins,rp) -PFINS ins; -Obj *rp; +void devalins(PFINS ins,Obj *rp) { + PFINS tins; PF pf; - PFAD ad; + PFAD ad,tad; int i; Real r; double d; - Q q; V v; P x; pf = ins->pf; ad = ins->ad; + tins = (PFINS)CALLOC(1,sizeof(PF)+pf->argc*sizeof(struct oPFAD)); + tins->pf = pf; tad = tins->ad; + for ( i = 0; i < pf->argc; i++ ) { + tad[i].d = ad[i].d; devalr(CO,ad[i].arg,&tad[i].arg); + } for ( i = 0; i < pf->argc; i++ ) - if ( ad[i].d || (ad[i].arg && !NUM(ad[i].arg)) ) - break; + if ( tad[i].d || (tad[i].arg && !NUM(tad[i].arg)) ) break; if ( (i != pf->argc) || !pf->libm ) { - instov(ins,&v); MKV(v,x); *rp = (Obj)x; + instoobj(tins,rp); } else { + for ( i = 0; i < pf->argc; i++ ) + if ( tad[i].arg && NID((Num)tad[i].arg) == N_C ) + error("devalins : not supported"); switch ( pf->argc ) { case 0: d = (*pf->libm)(); break; case 1: - d = (*pf->libm)(ToReal(ad[0].arg)); break; + d = (*pf->libm)(ToReal(tad[0].arg)); break; case 2: - d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg)); break; + d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg)); break; case 3: - d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg), - ToReal(ad[2].arg)); break; + d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg), + ToReal(tad[2].arg)); break; case 4: - d = (*pf->libm)(ToReal(ad[0].arg),ToReal(ad[1].arg), - ToReal(ad[2].arg),ToReal(ad[3].arg)); break; + d = (*pf->libm)(ToReal(tad[0].arg),ToReal(tad[1].arg), + ToReal(tad[2].arg),ToReal(tad[3].arg)); break; default: - error("devalv : not supported"); + error("devalins : not supported"); } MKReal(d,r); *rp = (Obj)r; } } -void simplify_ins(ins,rp) -PFINS ins; -Obj *rp; +extern int evalef,bigfloat; + +void simplify_elemfunc_ins(PFINS ins,Obj *rp) { + if ( evalef ) { + if ( bigfloat ) evalins(ins,0,rp); + else devalins(ins,rp); + } else instoobj(ins,rp); +} + +void simplify_factorial_ins(PFINS ins,Obj *rp) +{ + PFAD ad; + Obj a; + Q q; + + ad = ins->ad; + a = ad[0].arg; + if ( !ad[0].d && INT(a) && ( !a || (PL(NM((Q)a)) == 1 && SGN((Q)a) > 0) ) ) { + factorial(QTOS((Q)a),&q); + *rp = (Obj)q; + } else simplify_elemfunc_ins(ins,rp); +} + +void simplify_abs_ins(PFINS ins,Obj *rp) +{ + PFAD ad; + Obj a; + Q q; + double t; + Real r; + struct oNODE arg0; + + ad = ins->ad; + a = ad[0].arg; + if ( !ad[0].d && NUM(a) && (!a || RATN(a)) ) { + if ( !a || SGN((Q)a) > 0 ) *rp = (Obj)a; + else { + chsgnq((Q)a,&q); *rp = (Obj)q; + } + } else if ( !ad[0].d && REAL(a) ) { + t = fabs(((Real)a)->body); + MKReal(t,r); *rp = (Obj)r; + } else if ( !ad[0].d && BIGFLOAT(a) ) { + arg0.body = (pointer)a; arg0.next = 0; + mp_abs(&arg0,rp); + } else simplify_elemfunc_ins(ins,rp); +} + +void simplify_ins(PFINS ins,Obj *rp) +{ V v; P t; if ( ins->pf->simplify ) (*ins->pf->simplify)(ins,rp); else { - instov(ins,&v); MKV(v,t); *rp = (Obj)t; + instoobj(ins,rp); } } -void instov(ins,vp) -PFINS ins; -V *vp; +void instoobj(PFINS ins,Obj *rp) { - V v; + V v,newv; + P t; NEWV(v); NAME(v) = 0; v->attr = (pointer)V_PF; v->priv = (pointer)ins; - appendpfins(v,vp); + appendpfins(v,&newv); + MKV(newv,t); + *rp = (Obj)t; } -void substfr(vl,a,u,f,c) -VL vl; -Obj a; -PF u,f; -Obj *c; +void substfr(VL vl,Obj a,PF u,PF f,Obj *c) { Obj nm,dn; @@ -674,11 +773,7 @@ Obj *c; } } -void substfp(vl,p,u,f,pr) -VL vl; -Obj p; -PF u,f; -Obj *pr; +void substfp(VL vl,Obj p,PF u,PF f,Obj *pr) { V v; DCP dc; @@ -697,7 +792,8 @@ Obj *pr; for ( c = 0; dc; dc = NEXT(dc) ) { substfp(vl,(Obj)COEF(dc),u,f,&t); if ( DEG(dc) ) { - pwrp(vl,x,DEG(dc),(P *)&s); mulr(vl,s,t,&m); + gen_pwrr(vl,(Obj)x,(Obj)DEG(dc),&s); + mulr(vl,s,t,&m); addr(vl,m,c,&a); c = a; } else { addr(vl,t,c,&a); c = a; @@ -707,11 +803,12 @@ Obj *pr; substfv(vl,v,u,f,&p0); substfp(vl,(Obj)COEF(dc),u,f,&c); for ( d = DEG(dc), dc = NEXT(dc); dc; d = DEG(dc), dc = NEXT(dc) ) { - subq(d,DEG(dc),(Q *)&t); pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); + subq(d,DEG(dc),(Q *)&t); + gen_pwrr(vl,p0,t,&s); mulr(vl,s,c,&m); substfp(vl,(Obj)COEF(dc),u,f,&t); addr(vl,m,t,&c); } if ( d ) { - pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); + gen_pwrr(vl,p0,(Obj)d,&t); mulr(vl,t,c,&m); c = m; } } @@ -719,11 +816,7 @@ Obj *pr; } } -void substfv(vl,v,u,f,c) -VL vl; -V v; -PF u,f; -Obj *c; +void substfv(VL vl,V v,PF u,PF f,Obj *c) { P t; Obj r,s,w; @@ -754,6 +847,6 @@ Obj *c; for ( i = 0; i < f->argc; i++ ) { tad[i].d = ad[i].d; substfr(vl,ad[i].arg,u,f,&tad[i].arg); } - instov(tins,&v); MKV(v,t); *c = (Obj)t; + instoobj(tins,c); } }