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

Annotation of OpenXM_contrib2/asir2000/builtin/bfaux.c, Revision 1.6

1.6     ! takayama    1: /* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.5 2015/08/07 05:30:35 noro Exp $ */
1.1       noro        2: #include "ca.h"
                      3: #include "parse.h"
                      4:
1.4       noro        5: void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround();
1.6     ! takayama    6: void Pmpfr_gamma(), Pmpfr_floor();
1.1       noro        7:
                      8: struct ftab bf_tab[] = {
                      9:        {"eval",Peval,-2},
                     10:        {"setprec",Psetprec,-1},
1.3       noro       11:        {"setbprec",Psetbprec,-1},
1.4       noro       12:        {"setround",Psetround,-1},
1.1       noro       13:        {"todouble",Ptodouble,1},
1.5       noro       14:        {"mpfr_gamma",Pmpfr_gamma,-2},
1.6     ! takayama   15:        {"mpfr_floor",Pmpfr_floor,-1},
1.1       noro       16:        {0,0,0},
                     17: };
                     18:
1.4       noro       19: int mpfr_roundmode = MPFR_RNDN;
                     20:
1.1       noro       21: void Ptodouble(NODE arg,Num *rp)
                     22: {
                     23:        double r,i;
                     24:        Real real,imag;
                     25:        Num num;
                     26:
                     27:        asir_assert(ARG0(arg),O_N,"todouble");
                     28:        num = (Num)ARG0(arg);
                     29:        if ( !num ) {
                     30:                *rp = 0;
                     31:                return;
                     32:        }
                     33:        switch ( NID(num) ) {
                     34:                case N_R: case N_Q: case N_B:
                     35:                        r = ToReal(num);
                     36:                        MKReal(r,real);
                     37:                        *rp = (Num)real;
                     38:                        break;
                     39:                case N_C:
                     40:                        r = ToReal(((C)num)->r);
                     41:                        i = ToReal(((C)num)->i);
                     42:                        MKReal(r,real);
                     43:                        MKReal(i,imag);
                     44:                        reimtocplx((Num)real,(Num)imag,rp);
                     45:                        break;
                     46:                default:
                     47:                        *rp = num;
                     48:                        break;
                     49:        }
                     50: }
                     51:
1.4       noro       52: void Peval(NODE arg,Obj *rp)
1.1       noro       53: {
                     54:   int prec;
                     55:
                     56:        asir_assert(ARG0(arg),O_R,"eval");
                     57:   if ( argc(arg) == 2 ) {
1.3       noro       58:          prec = QTOS((Q)ARG1(arg))*3.32193;
1.1       noro       59:     if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                     60:     else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                     61:   } else
                     62:     prec = 0;
1.2       noro       63:        evalr(CO,(Obj)ARG0(arg),prec,rp);
1.1       noro       64: }
                     65:
1.3       noro       66: /* set/get decimal precision */
1.1       noro       67:
                     68: void Psetprec(NODE arg,Obj *rp)
                     69: {
                     70:        int p;
                     71:        Q q;
1.3       noro       72:        int prec,dprec;
                     73:
                     74:   prec = mpfr_get_default_prec();
                     75:   /* decimal precision */
                     76:   dprec = prec*0.30103;
                     77:   STOQ(dprec,q); *rp = (Obj)q;
                     78:        if ( arg ) {
                     79:                asir_assert(ARG0(arg),O_N,"setprec");
                     80:                prec = QTOS((Q)ARG0(arg))*3.32193;
                     81:                if ( p > 0 )
                     82:                        prec = p;
                     83:        }
                     84:   if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                     85:   else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                     86:        mpfr_set_default_prec(prec);
                     87: }
1.1       noro       88:
1.3       noro       89: /* set/get bit precision */
1.1       noro       90:
1.3       noro       91: void Psetbprec(NODE arg,Obj *rp)
                     92: {
                     93:        int p;
                     94:        Q q;
                     95:        int prec;
                     96:
                     97:   prec = mpfr_get_default_prec();
                     98:   STOQ(prec,q); *rp = (Obj)q;
1.1       noro       99:        if ( arg ) {
1.3       noro      100:                asir_assert(ARG0(arg),O_N,"setbprec");
                    101:                prec = QTOS((Q)ARG0(arg));
1.1       noro      102:                if ( p > 0 )
                    103:                        prec = p;
                    104:        }
                    105:   if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN;
                    106:   else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX;
                    107:        mpfr_set_default_prec(prec);
                    108: }
                    109:
1.4       noro      110: void Psetround(NODE arg,Q *rp)
                    111: {
                    112:   int round;
                    113:
                    114:   STOQ(mpfr_roundmode,*rp);
                    115:   if ( arg ) {
                    116:                asir_assert(ARG0(arg),O_N,"setround");
                    117:     round = QTOS((Q)ARG0(arg));
                    118:     switch ( round ) {
                    119:     case 0:
                    120:       mpfr_roundmode = MPFR_RNDN;
                    121:       break;
                    122:     case 1:
                    123:       mpfr_roundmode = MPFR_RNDZ;
                    124:       break;
                    125:     case 2:
                    126:       mpfr_roundmode = MPFR_RNDU;
                    127:       break;
                    128:     case 3:
                    129:       mpfr_roundmode = MPFR_RNDD;
                    130:       break;
                    131:     case 4:
                    132:       mpfr_roundmode = MPFR_RNDA;
                    133:       break;
                    134:     case 5:
                    135:       mpfr_roundmode = MPFR_RNDF;
                    136:       break;
                    137:     case 6:
                    138:       mpfr_roundmode = MPFR_RNDNA;
                    139:       break;
                    140:     default:
                    141:       error("setround : invalid rounding mode");
                    142:       break;
                    143:     }
                    144:   }
                    145: }
                    146:
1.1       noro      147: Num tobf(Num a,int prec);
                    148:
                    149: void mp_pi(NODE arg,BF *rp)
                    150: {
                    151:     int prec;
                    152:        BF r;
                    153:
                    154:        prec = arg ? QTOS((Q)ARG0(arg)) : 0;
                    155:        NEWBF(r);
                    156:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      157:        mpfr_const_pi(r->body,mpfr_roundmode);
1.1       noro      158:     *rp = r;
                    159: }
                    160:
                    161: void mp_e(NODE arg,BF *rp)
                    162: {
                    163:     int prec;
                    164:        mpfr_t one;
                    165:        BF r;
                    166:
                    167:        prec = arg ? QTOS((Q)ARG0(arg)) : 0;
                    168:        NEWBF(r);
                    169:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
                    170:        mpfr_init(one);
1.4       noro      171:        mpfr_set_ui(one,1,mpfr_roundmode);
                    172:        mpfr_exp(r->body,one,mpfr_roundmode);
1.1       noro      173:     *rp = r;
                    174: }
                    175:
                    176: void mp_sin(NODE arg,BF *rp)
                    177: {
                    178:        Num a;
                    179:     int prec;
                    180:        BF r;
                    181:
                    182:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    183:        a = tobf(ARG0(arg),prec);
                    184:        NEWBF(r);
                    185:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      186:        mpfr_sin(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      187:     *rp = r;
                    188: }
                    189:
                    190: void mp_cos(NODE arg,BF *rp)
                    191: {
                    192:        Num a;
                    193:     int prec;
                    194:        BF r;
                    195:
                    196:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    197:        a = tobf(ARG0(arg),prec);
                    198:        NEWBF(r);
                    199:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      200:        mpfr_cos(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      201:     *rp = r;
                    202: }
                    203:
                    204: void mp_tan(NODE arg,BF *rp)
                    205: {
                    206:        Num a;
                    207:     int prec;
                    208:        BF r;
                    209:
                    210:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    211:        a = tobf(ARG0(arg),prec);
                    212:        NEWBF(r);
                    213:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      214:        mpfr_tan(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      215:     *rp = r;
                    216: }
                    217:
                    218: void mp_asin(NODE arg,BF *rp)
                    219: {
                    220:        Num a;
                    221:     int prec;
                    222:        BF r;
                    223:
                    224:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    225:        a = tobf(ARG0(arg),prec);
                    226:        NEWBF(r);
                    227:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      228:        mpfr_asin(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      229:     *rp = r;
                    230: }
                    231: void mp_acos(NODE arg,BF *rp)
                    232: {
                    233:        Num a;
                    234:     int prec;
                    235:        BF r;
                    236:
                    237:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    238:        a = tobf(ARG0(arg),prec);
                    239:        NEWBF(r);
                    240:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      241:        mpfr_acos(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      242:     *rp = r;
                    243: }
                    244: void mp_atan(NODE arg,BF *rp)
                    245: {
                    246:        Num a;
                    247:     int prec;
                    248:        BF r;
                    249:
                    250:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    251:        a = tobf(ARG0(arg),prec);
                    252:        NEWBF(r);
                    253:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      254:        mpfr_atan(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      255:     *rp = r;
                    256: }
                    257:
                    258: void mp_sinh(NODE arg,BF *rp)
                    259: {
                    260:        Num a;
                    261:     int prec;
                    262:        BF r;
                    263:
                    264:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    265:        a = tobf(ARG0(arg),prec);
                    266:        NEWBF(r);
                    267:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      268:        mpfr_sinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      269:     *rp = r;
                    270: }
                    271:
                    272: void mp_cosh(NODE arg,BF *rp)
                    273: {
                    274:        Num a;
                    275:     int prec;
                    276:        BF r;
                    277:
                    278:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    279:        a = tobf(ARG0(arg),prec);
                    280:        NEWBF(r);
                    281:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      282:        mpfr_cosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      283:     *rp = r;
                    284: }
                    285:
                    286: void mp_tanh(NODE arg,BF *rp)
                    287: {
                    288:        Num a;
                    289:     int prec;
                    290:        BF r;
                    291:
                    292:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    293:        a = tobf(ARG0(arg),prec);
                    294:        NEWBF(r);
                    295:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      296:        mpfr_tanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      297:     *rp = r;
                    298: }
                    299:
                    300: void mp_asinh(NODE arg,BF *rp)
                    301: {
                    302:        Num a;
                    303:     int prec;
                    304:        BF r;
                    305:
                    306:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    307:        a = tobf(ARG0(arg),prec);
                    308:        NEWBF(r);
                    309:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      310:        mpfr_asinh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      311:     *rp = r;
                    312: }
                    313: void mp_acosh(NODE arg,BF *rp)
                    314: {
                    315:        Num a;
                    316:     int prec;
                    317:        BF r;
                    318:
                    319:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    320:        a = tobf(ARG0(arg),prec);
                    321:        NEWBF(r);
                    322:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      323:        mpfr_acosh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      324:     *rp = r;
                    325: }
                    326: void mp_atanh(NODE arg,BF *rp)
                    327: {
                    328:        Num a;
                    329:     int prec;
                    330:        BF r;
                    331:
                    332:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    333:        a = tobf(ARG0(arg),prec);
                    334:        NEWBF(r);
                    335:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      336:        mpfr_atanh(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      337:     *rp = r;
                    338: }
                    339:
                    340: void mp_exp(NODE arg,BF *rp)
                    341: {
                    342:        Num a;
                    343:     int prec;
                    344:        BF r;
                    345:
                    346:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    347:        a = tobf(ARG0(arg),prec);
                    348:        NEWBF(r);
                    349:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      350:        mpfr_exp(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      351:     *rp = r;
                    352: }
                    353:
                    354: void mp_log(NODE arg,BF *rp)
                    355: {
                    356:        Num a;
                    357:     int prec;
                    358:        BF r;
                    359:
                    360:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    361:        a = tobf(ARG0(arg),prec);
                    362:        NEWBF(r);
                    363:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      364:        mpfr_log(r->body,((BF)a)->body,mpfr_roundmode);
1.1       noro      365:     *rp = r;
                    366: }
                    367:
                    368: void mp_pow(NODE arg,BF *rp)
                    369: {
                    370:        Num a,e;
                    371:     int prec;
                    372:        BF r;
                    373:
                    374:        prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : 0;
                    375:        a = tobf(ARG0(arg),prec);
                    376:        e = tobf(ARG1(arg),prec);
                    377:        NEWBF(r);
                    378:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
1.4       noro      379:        mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode);
1.1       noro      380:     *rp = r;
                    381: }
1.5       noro      382:
                    383: void Pmpfr_gamma(NODE arg,BF *rp)
                    384: {
                    385:        Num a;
                    386:   int prec;
                    387:        BF r;
                    388:
                    389:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
                    390:   prec *= 3.32193;
                    391:        a = tobf(ARG0(arg),prec);
                    392:        NEWBF(r);
                    393:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
                    394:        mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode);
                    395:   *rp = r;
                    396: }
1.6     ! takayama  397:
        !           398: void Pmpfr_floor(NODE arg,BF *rp)
        !           399: {
        !           400:        Num a;
        !           401:   int prec;
        !           402:        BF r;
        !           403:
        !           404:        prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0;
        !           405:   prec *= 3.32193;
        !           406:        a = tobf(ARG0(arg),prec);
        !           407:        NEWBF(r);
        !           408:        prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
        !           409:        mpfr_floor(r->body,((BF)a)->body);
        !           410:   *rp = r;
        !           411: }

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