[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.5

1.2       noro        1: /*
                      2:  * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
                      3:  * All rights reserved.
                      4:  *
                      5:  * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
                      6:  * non-exclusive and royalty-free license to use, copy, modify and
                      7:  * redistribute, solely for non-commercial and non-profit purposes, the
                      8:  * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
                      9:  * conditions of this Agreement. For the avoidance of doubt, you acquire
                     10:  * only a limited right to use the SOFTWARE hereunder, and FLL or any
                     11:  * third party developer retains all rights, including but not limited to
                     12:  * copyrights, in and to the SOFTWARE.
                     13:  *
                     14:  * (1) FLL does not grant you a license in any way for commercial
                     15:  * purposes. You may use the SOFTWARE only for non-commercial and
                     16:  * non-profit purposes only, such as academic, research and internal
                     17:  * business use.
                     18:  * (2) The SOFTWARE is protected by the Copyright Law of Japan and
                     19:  * international copyright treaties. If you make copies of the SOFTWARE,
                     20:  * with or without modification, as permitted hereunder, you shall affix
                     21:  * to all such copies of the SOFTWARE the above copyright notice.
                     22:  * (3) An explicit reference to this SOFTWARE and its copyright owner
                     23:  * shall be made on your publication or presentation in any form of the
                     24:  * results obtained by use of the SOFTWARE.
                     25:  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3       noro       26:  * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2       noro       27:  * for such modification or the source code of the modified part of the
                     28:  * SOFTWARE.
                     29:  *
                     30:  * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
                     31:  * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
                     32:  * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
                     33:  * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
                     34:  * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
                     35:  * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
                     36:  * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
                     37:  * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
                     38:  * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
                     39:  * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
                     40:  * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
                     41:  * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
                     42:  * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
                     43:  * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
                     44:  * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
                     45:  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
                     46:  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
                     47:  *
1.5     ! noro       48:  * $OpenXM: OpenXM_contrib2/asir2000/builtin/pf.c,v 1.4 2001/04/20 02:34:21 noro Exp $
1.2       noro       49: */
1.1       noro       50: #include "ca.h"
                     51: #include "math.h"
                     52: #include "parse.h"
                     53: #if 0
                     54: #include <alloca.h>
                     55: #endif
                     56:
                     57: double const_pi(),const_e();
                     58:
                     59: void make_ihyp(void);
                     60: void make_hyp(void);
                     61: void make_itri(void);
                     62: void make_tri(void);
                     63: void make_exp(void);
                     64: void simplify_pow(PFINS,Obj *);
                     65:
                     66: void Pfunctor(),Pargs(),Pfunargs(),Pvtype(),Pcall(),Pdeval();
                     67: void Pregister_handler();
1.4       noro       68: void Peval_quote();
1.1       noro       69:
                     70: struct ftab puref_tab[] = {
                     71:        {"functor",Pfunctor,1},
                     72:        {"args",Pargs,1},
                     73:        {"funargs",Pfunargs,1},
                     74:        {"register_handler",Pregister_handler,1},
                     75:        {"call",Pcall,2},
                     76:        {"vtype",Pvtype,1},
                     77:        {"deval",Pdeval,1},
1.4       noro       78:        {"eval_quote",Peval_quote,1},
1.1       noro       79:        {0,0,0},
                     80: };
                     81:
                     82: #if PARI
                     83: int p_pi(),p_e();
                     84: int p_log(),p_exp(),p_pow();
                     85: int p_sin(),p_cos(),p_tan(),p_asin(),p_acos(),p_atan();
                     86: int p_sinh(),p_cosh(),p_tanh(),p_asinh(),p_acosh(),p_atanh();
                     87: #else
                     88: int p_pi,p_e;
                     89: int p_log,p_exp,p_pow;
                     90: int p_sin,p_cos,p_tan,p_asin,p_acos,p_atan;
                     91: int p_sinh,p_cosh,p_tanh,p_asinh,p_acosh,p_atanh;
                     92: #endif
                     93:
                     94: static V *uarg,*darg;
                     95: static P x,y;
                     96: static PF pidef,edef;
                     97: static PF logdef,expdef,powdef;
                     98: static PF sindef,cosdef,tandef;
                     99: static PF asindef,acosdef,atandef;
                    100: static PF sinhdef,coshdef,tanhdef;
                    101: static PF asinhdef,acoshdef,atanhdef;
                    102:
                    103: #define OALLOC(p,n) ((p)=(Obj *)CALLOC((n),sizeof(Obj)))
                    104:
                    105: double const_pi() { return 3.14159265358979323846264338327950288; }
                    106: double const_e() { return 2.718281828459045235360287471352662497; }
                    107:
                    108: void pf_init() {
                    109:        uarg = (V *)CALLOC(1,sizeof(V));
                    110:        uarg[0] = &oVAR[26]; MKV(uarg[0],x);
                    111:
                    112:        darg = (V *)CALLOC(2,sizeof(V));
                    113:        darg[0] = &oVAR[26];
                    114:        darg[1] = &oVAR[27]; MKV(darg[1],y);
                    115:
                    116:        mkpf("@pi",0,0,0,(int (*)())p_pi,const_pi,0,&pidef);
                    117:        mkpf("@e",0,0,0,(int (*)())p_e,const_e,0,&edef);
                    118:
                    119:        mkpf("log",0,1,uarg,(int (*)())p_log,log,0,&logdef);
                    120:        mkpf("exp",0,1,uarg,(int (*)())p_exp,exp,0,&expdef);
                    121:        mkpf("pow",0,2,darg,(int (*)())p_pow,pow,(int (*)())simplify_pow,&powdef);
                    122:
                    123:        mkpf("sin",0,1,uarg,(int (*)())p_sin,sin,0,&sindef);
                    124:        mkpf("cos",0,1,uarg,(int (*)())p_cos,cos,0,&cosdef);
                    125:        mkpf("tan",0,1,uarg,(int (*)())p_tan,tan,0,&tandef);
                    126:        mkpf("asin",0,1,uarg,(int (*)())p_asin,asin,0,&asindef);
                    127:        mkpf("acos",0,1,uarg,(int (*)())p_acos,acos,0,&acosdef);
                    128:        mkpf("atan",0,1,uarg,(int (*)())p_atan,atan,0,&atandef);
                    129:
                    130:        mkpf("sinh",0,1,uarg,(int (*)())p_sinh,sinh,0,&sinhdef);
                    131:        mkpf("cosh",0,1,uarg,(int (*)())p_cosh,cosh,0,&coshdef);
                    132:        mkpf("tanh",0,1,uarg,(int (*)())p_tanh,tanh,0,&tanhdef);
                    133: #if !defined(VISUAL)
                    134:        mkpf("asinh",0,1,uarg,(int (*)())p_asinh,asinh,0,&asinhdef);
                    135:        mkpf("acosh",0,1,uarg,(int (*)())p_acosh,acosh,0,&acoshdef);
                    136:        mkpf("atanh",0,1,uarg,(int (*)())p_atanh,atanh,0,&atanhdef);
                    137: #endif
                    138:        make_exp();
                    139:        make_tri();
                    140:        make_itri();
                    141:        make_hyp();
                    142: #if !defined(VISUAL)
                    143:        make_ihyp();
                    144: #endif
                    145: }
                    146:
                    147: void make_exp() {
                    148:        V v;
                    149:        P u,vexp,vlog,vpow;
                    150:        Obj *args;
                    151:
                    152:        mkpfins(expdef,uarg,&v); MKV(v,vexp);
                    153:        mkpfins(powdef,darg,&v); MKV(v,vpow);
                    154:        mkpfins(logdef,uarg,&v); MKV(v,vlog);
                    155:
                    156:        /* d/dx(log(x)) = 1/x */
                    157:        OALLOC(logdef->deriv,1); divr(CO,(Obj)ONE,(Obj)x,&logdef->deriv[0]);
                    158:
                    159:        /* d/dx(exp(x)) = exp(x) */
                    160:        OALLOC(expdef->deriv,1); expdef->deriv[0] = (Obj)vexp;
                    161:
                    162:        /* d/dy(x^y) = log(x)*x^y */
                    163:        OALLOC(powdef->deriv,2); mulp(CO,vpow,vlog,(P *)&powdef->deriv[1]);
                    164:
                    165:        /* d/dx(x^y) = y*x^(y-1) */
                    166:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    167:        args[0] = (Obj)x; subp(CO,y,(P)ONE,(P *)&args[1]);
                    168:        _mkpfins(powdef,args,&v); MKV(v,u);
                    169:        mulr(CO,(Obj)u,(Obj)y,&powdef->deriv[0]);
                    170: }
                    171:
                    172: void make_tri() {
                    173:        V v;
                    174:        P vcos,vsin,vtan,t;
                    175:
                    176:        mkpfins(cosdef,uarg,&v); MKV(v,vcos);
                    177:        mkpfins(sindef,uarg,&v); MKV(v,vsin);
                    178:        mkpfins(tandef,uarg,&v); MKV(v,vtan);
                    179:
                    180:        /* d/dx(sin(x)) = cos(x) */
                    181:        OALLOC(sindef->deriv,1); sindef->deriv[0] = (Obj)vcos;
                    182:
                    183:        /* d/dx(cos(x)) = -sin(x) */
                    184:        OALLOC(cosdef->deriv,1); chsgnp(vsin,(P *)&cosdef->deriv[0]);
                    185:
                    186:        /* d/dx(tan(x)) = 1+tan(x)^2 */
                    187:        OALLOC(tandef->deriv,1);
                    188:        mulr(CO,(Obj)vtan,(Obj)vtan,(Obj *)&t); addp(CO,(P)ONE,t,(P *)&tandef->deriv[0]);
                    189: }
                    190:
                    191: void make_itri() {
                    192:        P t,xx;
                    193:        Q mtwo;
                    194:        V v;
                    195:        Obj *args;
                    196:
                    197:        /* d/dx(asin(x)) = (1-x^2)^(-1/2) */
                    198:        OALLOC(asindef->deriv,1);
                    199:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    200:        mulp(CO,x,x,&xx); subp(CO,(P)ONE,xx,(P *)&args[0]);
                    201:        STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
                    202:        _mkpfins(powdef,args,&v); MKV(v,t);
                    203:        asindef->deriv[0] = (Obj)t;
                    204:
                    205:        /* d/dx(acos(x)) = -(1-x^2)^(-1/2) */
                    206:        OALLOC(acosdef->deriv,1); chsgnp((P)asindef->deriv[0],(P *)&acosdef->deriv[0]);
                    207:
                    208:        /* d/dx(atan(x)) = 1/(x^2+1) */
                    209:        OALLOC(atandef->deriv,1);
                    210:        addp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atandef->deriv[0]);
                    211: }
                    212:
                    213: void make_hyp() {
                    214:        V v;
                    215:        P vcosh,vsinh,vtanh,t;
                    216:
                    217:        mkpfins(coshdef,uarg,&v); MKV(v,vcosh);
                    218:        mkpfins(sinhdef,uarg,&v); MKV(v,vsinh);
                    219:        mkpfins(tanhdef,uarg,&v); MKV(v,vtanh);
                    220:
                    221:        /* d/dx(sinh(x)) = cosh(x) */
                    222:        OALLOC(sinhdef->deriv,1); sinhdef->deriv[0] = (Obj)vcosh;
                    223:
                    224:        /* d/dx(cosh(x)) = sinh(x) */
                    225:        OALLOC(coshdef->deriv,1); coshdef->deriv[0] = (Obj)vsinh;
                    226:
                    227:        /* d/dx(tanh(x)) = 1-tanh(x)^2 */
                    228:        OALLOC(tanhdef->deriv,1);
                    229:        mulr(CO,(Obj)vtanh,(Obj)vtanh,(Obj *)&t); subp(CO,(P)ONE,t,(P *)&tanhdef->deriv[0]);
                    230: }
                    231:
                    232: void make_ihyp() {
                    233:        P t,xx;
                    234:        Q mtwo;
                    235:        V v;
                    236:        Obj *args;
                    237:
                    238:        /* d/dx(asinh(x)) = (1+x^2)^(-1/2) */
                    239:        OALLOC(asinhdef->deriv,1);
                    240:        args = (Obj *)ALLOCA(2*sizeof(Obj));
                    241:        mulp(CO,x,x,&xx); addp(CO,(P)ONE,xx,(P *)&args[0]);
                    242:        STOQ(-2,mtwo); divq(ONE,mtwo,(Q *)&args[1]);
                    243:        _mkpfins(powdef,args,&v); MKV(v,t);
                    244:        asinhdef->deriv[0] = (Obj)t;
                    245:
                    246:        /* d/dx(acosh(x)) = (x^2-1)^(-1/2) */
                    247:        OALLOC(acoshdef->deriv,1);
                    248:        subp(CO,xx,(P)ONE,(P *)&args[0]);
                    249:        _mkpfins(powdef,args,&v); MKV(v,t);
                    250:        acoshdef->deriv[0] = (Obj)t;
                    251:
                    252:        /* d/dx(atanh(x)) = 1/(1-x^2) */
                    253:        OALLOC(atanhdef->deriv,1);
                    254:        subp(CO,(P)ONE,xx,&t); divr(CO,(Obj)ONE,(Obj)t,&atanhdef->deriv[0]);
                    255: }
                    256:
                    257: void mkpow(vl,a,e,r)
                    258: VL vl;
                    259: Obj a;
                    260: Obj e;
                    261: Obj *r;
                    262: {
                    263:        PFINS ins;
                    264:        PFAD ad;
                    265:
                    266:        ins = (PFINS)CALLOC(1,sizeof(PF)+2*sizeof(struct oPFAD));
                    267:        ins->pf = powdef; ad = ins->ad;
                    268:        ad[0].d = 0; ad[0].arg = a; ad[1].d = 0; ad[1].arg = e;
                    269:        simplify_ins(ins,r);
                    270: }
                    271:
                    272: void simplify_pow(ins,rp)
                    273: PFINS ins;
                    274: Obj *rp;
                    275: {
                    276:        PF pf;
                    277:        PFAD ad;
                    278:        Obj a0,a1;
                    279:        V v;
                    280:        P t;
                    281:
                    282:        pf = ins->pf; ad = ins->ad; a0 = ad[0].arg; a1 = ad[1].arg;
                    283:        if ( !a1 )
                    284:                *rp = (Obj)ONE;
                    285:        else if ( !a0 )
                    286:                *rp = 0;
                    287:        else if ( NUM(a1) && INT(a1) )
                    288:                arf_pwr(CO,a0,a1,rp);
                    289:        else {
                    290:                instov(ins,&v); MKV(v,t); *rp = (Obj)t;
                    291:        }
                    292: }
                    293:
                    294: #define ISPFINS(p)\
                    295: (p)&&(ID(p) == O_P)&&((int)VR((P)p)->attr!=V_PF)&&\
                    296: UNIQ(DEG(DC((P)p)))&&UNIQ(COEF(DC((P)p)))
                    297:
                    298: void Pfunctor(arg,rp)
                    299: NODE arg;
                    300: P *rp;
                    301: {
                    302:        P p;
                    303:        FUNC t;
                    304:        PF pf;
                    305:        PFINS ins;
                    306:
                    307:        p = (P)ARG0(arg);
                    308:        if ( !ISPFINS(p) )
                    309:                *rp = 0;
                    310:        else {
                    311:                ins = (PFINS)VR(p)->priv; 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,rp);
                    316:        }
                    317: }
                    318:
                    319: void Pargs(arg,rp)
                    320: NODE arg;
                    321: LIST *rp;
                    322: {
                    323:        P p;
                    324:        PF pf;
                    325:        PFAD ad;
                    326:        PFINS ins;
                    327:        NODE n,n0;
                    328:        int i;
                    329:
                    330:        p = (P)ARG0(arg);
                    331:        if ( !ISPFINS(p) )
                    332:                *rp = 0;
                    333:        else {
                    334:                ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
                    335:                for ( i = 0, n0 = 0; i < pf->argc; i++ ) {
                    336:                        NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
                    337:                }
                    338:                if ( n0 )
                    339:                        NEXT(n) = 0;
                    340:                MKLIST(*rp,n0);
                    341:        }
                    342: }
                    343:
                    344: void Pfunargs(arg,rp)
                    345: NODE arg;
                    346: LIST *rp;
                    347: {
                    348:        P p;
                    349:        P f;
                    350:        FUNC t;
                    351:        PF pf;
                    352:        PFINS ins;
                    353:        PFAD ad;
                    354:        NODE n,n0;
                    355:        int i;
                    356:
                    357:        p = (P)ARG0(arg);
                    358:        if ( !ISPFINS(p) )
                    359:                *rp = 0;
                    360:        else {
                    361:                ins = (PFINS)VR(p)->priv; ad = ins->ad; pf = ins->pf;
                    362:                t = (FUNC)MALLOC(sizeof(struct oFUNC));
                    363:                t->name = pf->name; t->id = A_PURE; t->argc = pf->argc;
                    364:                t->f.puref = pf;
                    365:                makesrvar(t,&f);
1.5     ! noro      366:                n = n0 = 0; NEXTNODE(n0,n); BDY(n) = (pointer)f;
1.1       noro      367:                for ( i = 0; i < pf->argc; i++ ) {
                    368:                        NEXTNODE(n0,n); BDY(n) = (pointer)ad[i].arg;
                    369:                }
                    370:                NEXT(n) = 0;
                    371:                MKLIST(*rp,n0);
                    372:        }
                    373: }
                    374:
                    375: void Pvtype(arg,rp)
                    376: NODE arg;
                    377: Q *rp;
                    378: {
                    379:        P p;
                    380:
                    381:        p = (P)ARG0(arg);
                    382:        if ( !p || ID(p) != O_P )
                    383:                *rp = 0;
                    384:        else
                    385:                STOQ((int)VR(p)->attr,*rp);
                    386: }
                    387:
                    388: extern FUNC registered_handler;
                    389:
                    390: void Pregister_handler(arg,rp)
                    391: NODE arg;
                    392: Q *rp;
                    393: {
                    394:        P p;
                    395:        V v;
                    396:        FUNC func;
                    397:
                    398:        p = (P)ARG0(arg);
                    399:        if ( !p )
                    400:                registered_handler = 0;
                    401:        else if ( OID(p) != 2 )
                    402:                error("register_hanlder : invalid argument");
                    403:        v = VR(p);
                    404:        if ( (int)v->attr != V_SR )
                    405:                error("register_hanlder : no such function");
                    406:        else {
                    407:                func = (FUNC)v->priv;
                    408:                if ( func->argc )
                    409:                        error("register_hanlder : the function must be with no argument");
                    410:                else {
                    411:                        registered_handler = func;
                    412:                        *rp = ONE;
                    413:                }
                    414:        }
                    415: }
                    416:
                    417: void Pcall(arg,rp)
                    418: NODE arg;
                    419: Obj *rp;
                    420: {
                    421:        P p;
                    422:        V v;
                    423:
                    424:        p = (P)ARG0(arg);
                    425:        if ( !p || OID(p) != 2 )
                    426:                error("call : invalid argument");
                    427:        v = VR(p);
                    428:        if ( (int)v->attr != V_SR )
                    429:                error("call : no such function");
                    430:
                    431:        else
                    432:                *rp = (Obj)bevalf((FUNC)v->priv,BDY((LIST)ARG1(arg)));
                    433: }
                    434:
                    435: void Pdeval(arg,rp)
                    436: NODE arg;
                    437: Obj *rp;
                    438: {
                    439:        asir_assert(ARG0(arg),O_R,"deval");
                    440:        devalr(CO,(Obj)ARG0(arg),rp);
                    441: }
                    442:
1.4       noro      443: void Peval_quote(arg,rp)
                    444: NODE arg;
                    445: Obj *rp;
                    446: {
                    447:        asir_assert(ARG0(arg),O_QUOTE,"eval_quote");
                    448:        *rp = eval((FNODE)BDY((QUOTE)ARG0(arg)));
                    449: }

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