[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     ! 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>