=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v retrieving revision 1.2 retrieving revision 1.6 diff -u -p -r1.2 -r1.6 --- OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/05 01:10:38 1.2 +++ OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/07 06:15:00 1.6 @@ -1,16 +1,23 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.1 2015/08/04 06:55:02 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.5 2015/08/07 05:30:35 noro Exp $ */ #include "ca.h" #include "parse.h" -void Peval(), Psetprec(), Ptodouble(); +void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround(); +void Pmpfr_gamma(), Pmpfr_floor(); struct ftab bf_tab[] = { {"eval",Peval,-2}, {"setprec",Psetprec,-1}, + {"setbprec",Psetbprec,-1}, + {"setround",Psetround,-1}, {"todouble",Ptodouble,1}, + {"mpfr_gamma",Pmpfr_gamma,-2}, + {"mpfr_floor",Pmpfr_floor,-1}, {0,0,0}, }; +int mpfr_roundmode = MPFR_RNDN; + void Ptodouble(NODE arg,Num *rp) { double r,i; @@ -42,15 +49,13 @@ void Ptodouble(NODE arg,Num *rp) } } -void Peval(arg,rp) -NODE arg; -Obj *rp; +void Peval(NODE arg,Obj *rp) { int prec; asir_assert(ARG0(arg),O_R,"eval"); if ( argc(arg) == 2 ) { - prec = QTOS((Q)ARG1(arg)); + prec = QTOS((Q)ARG1(arg))*3.32193; if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN; else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX; } else @@ -58,19 +63,42 @@ Obj *rp; evalr(CO,(Obj)ARG0(arg),prec,rp); } -/* bit precision */ +/* set/get decimal precision */ void Psetprec(NODE arg,Obj *rp) { int p; Q q; + int prec,dprec; + + prec = mpfr_get_default_prec(); + /* decimal precision */ + dprec = prec*0.30103; + STOQ(dprec,q); *rp = (Obj)q; + if ( arg ) { + asir_assert(ARG0(arg),O_N,"setprec"); + prec = QTOS((Q)ARG0(arg))*3.32193; + if ( p > 0 ) + prec = p; + } + if ( prec < MPFR_PREC_MIN ) prec = MPFR_PREC_MIN; + else if ( prec > MPFR_PREC_MAX ) prec = MPFR_PREC_MAX; + mpfr_set_default_prec(prec); +} - long prec = mpfr_get_default_prec(); +/* set/get bit precision */ - STOQ(prec,q); *rp = (Obj)q; +void Psetbprec(NODE arg,Obj *rp) +{ + int p; + Q q; + int prec; + + prec = mpfr_get_default_prec(); + STOQ(prec,q); *rp = (Obj)q; if ( arg ) { - asir_assert(ARG0(arg),O_N,"setprec"); - p = QTOS((Q)ARG0(arg)); + asir_assert(ARG0(arg),O_N,"setbprec"); + prec = QTOS((Q)ARG0(arg)); if ( p > 0 ) prec = p; } @@ -79,6 +107,43 @@ void Psetprec(NODE arg,Obj *rp) mpfr_set_default_prec(prec); } +void Psetround(NODE arg,Q *rp) +{ + int round; + + STOQ(mpfr_roundmode,*rp); + if ( arg ) { + asir_assert(ARG0(arg),O_N,"setround"); + round = QTOS((Q)ARG0(arg)); + switch ( round ) { + case 0: + mpfr_roundmode = MPFR_RNDN; + break; + case 1: + mpfr_roundmode = MPFR_RNDZ; + break; + case 2: + mpfr_roundmode = MPFR_RNDU; + break; + case 3: + mpfr_roundmode = MPFR_RNDD; + break; + case 4: + mpfr_roundmode = MPFR_RNDA; + break; + case 5: + mpfr_roundmode = MPFR_RNDF; + break; + case 6: + mpfr_roundmode = MPFR_RNDNA; + break; + default: + error("setround : invalid rounding mode"); + break; + } + } +} + Num tobf(Num a,int prec); void mp_pi(NODE arg,BF *rp) @@ -89,7 +154,7 @@ void mp_pi(NODE arg,BF *rp) prec = arg ? QTOS((Q)ARG0(arg)) : 0; NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_const_pi(r->body,MPFR_RNDN); + mpfr_const_pi(r->body,mpfr_roundmode); *rp = r; } @@ -103,8 +168,8 @@ void mp_e(NODE arg,BF *rp) NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); mpfr_init(one); - mpfr_set_ui(one,1,MPFR_RNDN); - mpfr_exp(r->body,one,MPFR_RNDN); + mpfr_set_ui(one,1,mpfr_roundmode); + mpfr_exp(r->body,one,mpfr_roundmode); *rp = r; } @@ -118,7 +183,7 @@ void mp_sin(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_sin(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_sin(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -132,7 +197,7 @@ void mp_cos(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_cos(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_cos(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -146,7 +211,7 @@ void mp_tan(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_tan(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_tan(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -160,7 +225,7 @@ void mp_asin(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_asin(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_asin(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void mp_acos(NODE arg,BF *rp) @@ -173,7 +238,7 @@ void mp_acos(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_acos(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_acos(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void mp_atan(NODE arg,BF *rp) @@ -186,7 +251,7 @@ void mp_atan(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_atan(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_atan(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -200,7 +265,7 @@ void mp_sinh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_sinh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_sinh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -214,7 +279,7 @@ void mp_cosh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_cosh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_cosh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -228,7 +293,7 @@ void mp_tanh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_tanh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_tanh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -242,7 +307,7 @@ void mp_asinh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_asinh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_asinh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void mp_acosh(NODE arg,BF *rp) @@ -255,7 +320,7 @@ void mp_acosh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_acosh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_acosh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void mp_atanh(NODE arg,BF *rp) @@ -268,7 +333,7 @@ void mp_atanh(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_atanh(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_atanh(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -282,7 +347,7 @@ void mp_exp(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_exp(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_exp(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -296,7 +361,7 @@ void mp_log(NODE arg,BF *rp) a = tobf(ARG0(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_log(r->body,((BF)a)->body,MPFR_RNDN); + mpfr_log(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } @@ -311,6 +376,36 @@ void mp_pow(NODE arg,BF *rp) e = tobf(ARG1(arg),prec); NEWBF(r); prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); - mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,MPFR_RNDN); + mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode); *rp = r; +} + +void Pmpfr_gamma(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; + prec *= 3.32193; + a = tobf(ARG0(arg),prec); + NEWBF(r); + prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); + mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_floor(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; + prec *= 3.32193; + a = tobf(ARG0(arg),prec); + NEWBF(r); + prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body); + mpfr_floor(r->body,((BF)a)->body); + *rp = r; }