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

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

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