[BACK]Return to pf.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/pf.c, Revision 1.1.1.1

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/builtin/pf.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
                      2: #include "ca.h"
                      3: #include "math.h"
                      4: #include "parse.h"
                      5: #if 0
                      6: #include <alloca.h>
                      7: #endif
                      8:
                      9: double const_pi(),const_e();
                     10:
                     11: void make_ihyp(void);
                     12: void make_hyp(void);
                     13: void make_itri(void);
                     14: void make_tri(void);
                     15: void make_exp(void);
                     16: void simplify_pow(PFINS,Obj *);
                     17:
                     18: void Pfunctor(),Pargs(),Pfunargs(),Pvtype(),Pcall(),Pdeval();
                     19: void Pregister_handler();
                     20:
                     21: struct ftab puref_tab[] = {
                     22:        {"functor",Pfunctor,1},
                     23:        {"args",Pargs,1},
                     24:        {"funargs",Pfunargs,1},
                     25:        {"register_handler",Pregister_handler,1},
                     26:        {"call",Pcall,2},
                     27:        {"vtype",Pvtype,1},
                     28:        {"deval",Pdeval,1},
                     29:        {0,0,0},
                     30: };
                     31:
                     32: #if PARI
                     33: int p_pi(),p_e();
                     34: int p_log(),p_exp(),p_pow();
                     35: int p_sin(),p_cos(),p_tan(),p_asin(),p_acos(),p_atan();
                     36: int p_sinh(),p_cosh(),p_tanh(),p_asinh(),p_acosh(),p_atanh();
                     37: #else
                     38: int p_pi,p_e;
                     39: int p_log,p_exp,p_pow;
                     40: int p_sin,p_cos,p_tan,p_asin,p_acos,p_atan;
                     41: int p_sinh,p_cosh,p_tanh,p_asinh,p_acosh,p_atanh;
                     42: #endif
                     43:
                     44: static V *uarg,*darg;
                     45: static P x,y;
                     46: static PF pidef,edef;
                     47: static PF logdef,expdef,powdef;
                     48: static PF sindef,cosdef,tandef;
                     49: static PF asindef,acosdef,atandef;
                     50: static PF sinhdef,coshdef,tanhdef;
                     51: static PF asinhdef,acoshdef,atanhdef;
                     52:
                     53: #define OALLOC(p,n) ((p)=(Obj *)CALLOC((n),sizeof(Obj)))
                     54:
                     55: double const_pi() { return 3.14159265358979323846264338327950288; }
                     56: double const_e() { return 2.718281828459045235360287471352662497; }
                     57:
                     58: void pf_init() {
                     59:        uarg = (V *)CALLOC(1,sizeof(V));
                     60:        uarg[0] = &oVAR[26]; MKV(uarg[0],x);
                     61:
                     62:        darg = (V *)CALLOC(2,sizeof(V));
                     63:        darg[0] = &oVAR[26];
                     64:        darg[1] = &oVAR[27]; MKV(darg[1],y);
                     65:
                     66:        mkpf("@pi",0,0,0,(int (*)())p_pi,const_pi,0,&pidef);
                     67:        mkpf("@e",0,0,0,(int (*)())p_e,const_e,0,&edef);
                     68:
                     69:        mkpf("log",0,1,uarg,(int (*)())p_log,log,0,&logdef);
                     70:        mkpf("exp",0,1,uarg,(int (*)())p_exp,exp,0,&expdef);
                     71:        mkpf("pow",0,2,darg,(int (*)())p_pow,pow,(int (*)())simplify_pow,&powdef);
                     72:
                     73:        mkpf("sin",0,1,uarg,(int (*)())p_sin,sin,0,&sindef);
                     74:        mkpf("cos",0,1,uarg,(int (*)())p_cos,cos,0,&cosdef);
                     75:        mkpf("tan",0,1,uarg,(int (*)())p_tan,tan,0,&tandef);
                     76:        mkpf("asin",0,1,uarg,(int (*)())p_asin,asin,0,&asindef);
                     77:        mkpf("acos",0,1,uarg,(int (*)())p_acos,acos,0,&acosdef);
                     78:        mkpf("atan",0,1,uarg,(int (*)())p_atan,atan,0,&atandef);
                     79:
                     80:        mkpf("sinh",0,1,uarg,(int (*)())p_sinh,sinh,0,&sinhdef);
                     81:        mkpf("cosh",0,1,uarg,(int (*)())p_cosh,cosh,0,&coshdef);
                     82:        mkpf("tanh",0,1,uarg,(int (*)())p_tanh,tanh,0,&tanhdef);
                     83: #if !defined(VISUAL)
                     84:        mkpf("asinh",0,1,uarg,(int (*)())p_asinh,asinh,0,&asinhdef);
                     85:        mkpf("acosh",0,1,uarg,(int (*)())p_acosh,acosh,0,&acoshdef);
                     86:        mkpf("atanh",0,1,uarg,(int (*)())p_atanh,atanh,0,&atanhdef);
                     87: #endif
                     88:        make_exp();
                     89:        make_tri();
                     90:        make_itri();
                     91:        make_hyp();
                     92: #if !defined(VISUAL)
                     93:        make_ihyp();
                     94: #endif
                     95: }
                     96:
                     97: void make_exp() {
                     98:        V v;
                     99:        P u,vexp,vlog,vpow;
                    100:        Obj *args;
                    101:
                    102:        mkpfins(expdef,uarg,&v); MKV(v,vexp);
                    103:        mkpfins(powdef,darg,&v); MKV(v,vpow);
                    104:        mkpfins(logdef,uarg,&v); MKV(v,vlog);
                    105:
                    106:        /* d/dx(log(x)) = 1/x */
                    107:        OALLOC(logdef->deriv,1); divr(CO,(Obj)ONE,(Obj)x,&logdef->deriv[0]);
                    108:
                    109:        /* d/dx(exp(x)) = exp(x) */
                    110:        OALLOC(expdef->deriv,1); expdef->deriv[0] = (Obj)vexp;
                    111:
                    112:        /* d/dy(x^y) = log(x)*x^y */
                    113:        OALLOC(powdef->deriv,2); mulp(CO,vpow,vlog,(P *)&powdef->deriv[1]);
                    114:
                    115:        /* d/dx(x^y) = y*x^(y-1) */
                    116:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    117:        args[0] = (Obj)x; subp(CO,y,(P)ONE,(P *)&args[1]);
                    118:        _mkpfins(powdef,args,&v); MKV(v,u);
                    119:        mulr(CO,(Obj)u,(Obj)y,&powdef->deriv[0]);
                    120: }
                    121:
                    122: void make_tri() {
                    123:        V v;
                    124:        P vcos,vsin,vtan,t;
                    125:
                    126:        mkpfins(cosdef,uarg,&v); MKV(v,vcos);
                    127:        mkpfins(sindef,uarg,&v); MKV(v,vsin);
                    128:        mkpfins(tandef,uarg,&v); MKV(v,vtan);
                    129:
                    130:        /* d/dx(sin(x)) = cos(x) */
                    131:        OALLOC(sindef->deriv,1); sindef->deriv[0] = (Obj)vcos;
                    132:
                    133:        /* d/dx(cos(x)) = -sin(x) */
                    134:        OALLOC(cosdef->deriv,1); chsgnp(vsin,(P *)&cosdef->deriv[0]);
                    135:
                    136:        /* d/dx(tan(x)) = 1+tan(x)^2 */
                    137:        OALLOC(tandef->deriv,1);
                    138:        mulr(CO,(Obj)vtan,(Obj)vtan,(Obj *)&t); addp(CO,(P)ONE,t,(P *)&tandef->deriv[0]);
                    139: }
                    140:
                    141: void make_itri() {
                    142:        P t,xx;
                    143:        Q mtwo;
                    144:        V v;
                    145:        Obj *args;
                    146:
                    147:        /* d/dx(asin(x)) = (1-x^2)^(-1/2) */
                    148:        OALLOC(asindef->deriv,1);
                    149:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    150:        mulp(CO,x,x,&xx); subp(CO,(P)ONE,xx,(P *)&args[0]);
                    151:        STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
                    152:        _mkpfins(powdef,args,&v); MKV(v,t);
                    153:        asindef->deriv[0] = (Obj)t;
                    154:
                    155:        /* d/dx(acos(x)) = -(1-x^2)^(-1/2) */
                    156:        OALLOC(acosdef->deriv,1); chsgnp((P)asindef->deriv[0],(P *)&acosdef->deriv[0]);
                    157:
                    158:        /* d/dx(atan(x)) = 1/(x^2+1) */
                    159:        OALLOC(atandef->deriv,1);
                    160:        addp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atandef->deriv[0]);
                    161: }
                    162:
                    163: void make_hyp() {
                    164:        V v;
                    165:        P vcosh,vsinh,vtanh,t;
                    166:
                    167:        mkpfins(coshdef,uarg,&v); MKV(v,vcosh);
                    168:        mkpfins(sinhdef,uarg,&v); MKV(v,vsinh);
                    169:        mkpfins(tanhdef,uarg,&v); MKV(v,vtanh);
                    170:
                    171:        /* d/dx(sinh(x)) = cosh(x) */
                    172:        OALLOC(sinhdef->deriv,1); sinhdef->deriv[0] = (Obj)vcosh;
                    173:
                    174:        /* d/dx(cosh(x)) = sinh(x) */
                    175:        OALLOC(coshdef->deriv,1); coshdef->deriv[0] = (Obj)vsinh;
                    176:
                    177:        /* d/dx(tanh(x)) = 1-tanh(x)^2 */
                    178:        OALLOC(tanhdef->deriv,1);
                    179:        mulr(CO,(Obj)vtanh,(Obj)vtanh,(Obj *)&t); subp(CO,(P)ONE,t,(P *)&tanhdef->deriv[0]);
                    180: }
                    181:
                    182: void make_ihyp() {
                    183:        P t,xx;
                    184:        Q mtwo;
                    185:        V v;
                    186:        Obj *args;
                    187:
                    188:        /* d/dx(asinh(x)) = (1+x^2)^(-1/2) */
                    189:        OALLOC(asinhdef->deriv,1);
                    190:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    191:        mulp(CO,x,x,&xx); addp(CO,(P)ONE,xx,(P *)&args[0]);
                    192:        STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
                    193:        _mkpfins(powdef,args,&v); MKV(v,t);
                    194:        asinhdef->deriv[0] = (Obj)t;
                    195:
                    196:        /* d/dx(acosh(x)) = (x^2-1)^(-1/2) */
                    197:        OALLOC(acoshdef->deriv,1);
                    198:        subp(CO,xx,(P)ONE,(P *)&args[0]);
                    199:        _mkpfins(powdef,args,&v); MKV(v,t);
                    200:        acoshdef->deriv[0] = (Obj)t;
                    201:
                    202:        /* d/dx(atanh(x)) = 1/(1-x^2) */
                    203:        OALLOC(atanhdef->deriv,1);
                    204:        subp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atanhdef->deriv[0]);
                    205: }
                    206:
                    207: void mkpow(vl,a,e,r)
                    208: VL vl;
                    209: Obj a;
                    210: Obj e;
                    211: Obj *r;
                    212: {
                    213:        PFINS ins;
                    214:        PFAD ad;
                    215:
                    216:        ins = (PFINS)CALLOC(1,sizeof(PF)+2*sizeof(struct oPFAD));
                    217:        ins->pf = powdef; ad = ins->ad;
                    218:        ad[0].d = 0; ad[0].arg = a; ad[1].d = 0; ad[1].arg = e;
                    219:        simplify_ins(ins,r);
                    220: }
                    221:
                    222: void simplify_pow(ins,rp)
                    223: PFINS ins;
                    224: Obj *rp;
                    225: {
                    226:        PF pf;
                    227:        PFAD ad;
                    228:        Obj a0,a1;
                    229:        V v;
                    230:        P t;
                    231:
                    232:        pf = ins->pf; ad = ins->ad; a0 = ad[0].arg; a1 = ad[1].arg;
                    233:        if ( !a1 )
                    234:                *rp = (Obj)ONE;
                    235:        else if ( !a0 )
                    236:                *rp = 0;
                    237:        else if ( NUM(a1) && INT(a1) )
                    238:                arf_pwr(CO,a0,a1,rp);
                    239:        else {
                    240:                instov(ins,&v); MKV(v,t); *rp = (Obj)t;
                    241:        }
                    242: }
                    243:
                    244: #define ISPFINS(p)\
                    245: (p)&&(ID(p) == O_P)&&((int)VR((P)p)->attr!=V_PF)&&\
                    246: UNIQ(DEG(DC((P)p)))&&UNIQ(COEF(DC((P)p)))
                    247:
                    248: void Pfunctor(arg,rp)
                    249: NODE arg;
                    250: P *rp;
                    251: {
                    252:        P p;
                    253:        FUNC t;
                    254:        PF pf;
                    255:        PFINS ins;
                    256:
                    257:        p = (P)ARG0(arg);
                    258:        if ( !ISPFINS(p) )
                    259:                *rp = 0;
                    260:        else {
                    261:                ins = (PFINS)VR(p)->priv; pf = ins->pf;
                    262:                t = (FUNC)MALLOC(sizeof(struct oFUNC));
                    263:                t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
                    264:                t->f.puref = pf;
                    265:                makesrvar(t,rp);
                    266:        }
                    267: }
                    268:
                    269: void Pargs(arg,rp)
                    270: NODE arg;
                    271: LIST *rp;
                    272: {
                    273:        P p;
                    274:        PF pf;
                    275:        PFAD ad;
                    276:        PFINS ins;
                    277:        NODE n,n0;
                    278:        int i;
                    279:
                    280:        p = (P)ARG0(arg);
                    281:        if ( !ISPFINS(p) )
                    282:                *rp = 0;
                    283:        else {
                    284:                ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
                    285:                for ( i = 0, n0 = 0; i < pf->argc; i++ ) {
                    286:                        NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
                    287:                }
                    288:                if ( n0 )
                    289:                        NEXT(n) = 0;
                    290:                MKLIST(*rp,n0);
                    291:        }
                    292: }
                    293:
                    294: void Pfunargs(arg,rp)
                    295: NODE arg;
                    296: LIST *rp;
                    297: {
                    298:        P p;
                    299:        P f;
                    300:        FUNC t;
                    301:        PF pf;
                    302:        PFINS ins;
                    303:        PFAD ad;
                    304:        NODE n,n0;
                    305:        int i;
                    306:
                    307:        p = (P)ARG0(arg);
                    308:        if ( !ISPFINS(p) )
                    309:                *rp = 0;
                    310:        else {
                    311:                ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
                    312:                t = (FUNC)MALLOC(sizeof(struct oFUNC));
                    313:                t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
                    314:                t->f.puref = pf;
                    315:                makesrvar(t,&f);
                    316:                n0 = 0; NEXTNODE(n0,n); BDY(n) = (pointer)f;
                    317:                for ( i = 0; i < pf->argc; i++ ) {
                    318:                        NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
                    319:                }
                    320:                NEXT(n) = 0;
                    321:                MKLIST(*rp,n0);
                    322:        }
                    323: }
                    324:
                    325: void Pvtype(arg,rp)
                    326: NODE arg;
                    327: Q *rp;
                    328: {
                    329:        P p;
                    330:
                    331:        p = (P)ARG0(arg);
                    332:        if ( !p || ID(p) != O_P )
                    333:                *rp = 0;
                    334:        else
                    335:                STOQ((int)VR(p)->attr,*rp);
                    336: }
                    337:
                    338: extern FUNC registered_handler;
                    339:
                    340: void Pregister_handler(arg,rp)
                    341: NODE arg;
                    342: Q *rp;
                    343: {
                    344:        P p;
                    345:        V v;
                    346:        FUNC func;
                    347:
                    348:        p = (P)ARG0(arg);
                    349:        if ( !p )
                    350:                registered_handler = 0;
                    351:        else if ( OID(p) != 2 )
                    352:                error("register_hanlder : invalid argument");
                    353:        v = VR(p);
                    354:        if ( (int)v->attr != V_SR )
                    355:                error("register_hanlder : no such function");
                    356:        else {
                    357:                func = (FUNC)v->priv;
                    358:                if ( func->argc )
                    359:                        error("register_hanlder : the function must be with no argument");
                    360:                else {
                    361:                        registered_handler = func;
                    362:                        *rp = ONE;
                    363:                }
                    364:        }
                    365: }
                    366:
                    367: void Pcall(arg,rp)
                    368: NODE arg;
                    369: Obj *rp;
                    370: {
                    371:        P p;
                    372:        V v;
                    373:
                    374:        p = (P)ARG0(arg);
                    375:        if ( !p || OID(p) != 2 )
                    376:                error("call : invalid argument");
                    377:        v = VR(p);
                    378:        if ( (int)v->attr != V_SR )
                    379:                error("call : no such function");
                    380:
                    381:        else
                    382:                *rp = (Obj)bevalf((FUNC)v->priv,BDY((LIST)ARG1(arg)));
                    383: }
                    384:
                    385: void Pdeval(arg,rp)
                    386: NODE arg;
                    387: Obj *rp;
                    388: {
                    389:        asir_assert(ARG0(arg),O_R,"deval");
                    390:        devalr(CO,(Obj)ARG0(arg),rp);
                    391: }
                    392:

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>