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

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

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