=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v retrieving revision 1.3 retrieving revision 1.8 diff -u -p -r1.3 -r1.8 --- OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/05 01:23:02 1.3 +++ OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/16 03:12:09 1.8 @@ -1,17 +1,42 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.2 2015/08/05 01:10:38 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.7 2015/08/07 08:00:30 takayama Exp $ */ #include "ca.h" #include "parse.h" -void Peval(), Psetprec(), Psetbprec(), Ptodouble(); +void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround(); +void Pmpfr_ai(); +void Pmpfr_eint(), Pmpfr_erf(),Pmpfr_li2(); +void Pmpfr_zeta(); +void Pmpfr_j0(), Pmpfr_j1(); +void Pmpfr_y0(), Pmpfr_y1(); +void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma(); +void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil(); struct ftab bf_tab[] = { {"eval",Peval,-2}, {"setprec",Psetprec,-1}, {"setbprec",Psetbprec,-1}, + {"setround",Psetround,-1}, {"todouble",Ptodouble,1}, + {"mpfr_ai",Pmpfr_ai,-2}, + {"mpfr_zeta",Pmpfr_zeta,-2}, + {"mpfr_j0",Pmpfr_j0,-2}, + {"mpfr_j1",Pmpfr_j1,-2}, + {"mpfr_y0",Pmpfr_y0,-2}, + {"mpfr_y1",Pmpfr_y1,-2}, + {"mpfr_eint",Pmpfr_eint,-2}, + {"mpfr_erf",Pmpfr_erf,-2}, + {"mpfr_li2",Pmpfr_li2,-2}, + {"mpfr_gamma",Pmpfr_gamma,-2}, + {"mpfr_lngamma",Pmpfr_gamma,-2}, + {"mpfr_digamma",Pmpfr_gamma,-2}, + {"mpfr_floor",Pmpfr_floor,-2}, + {"mpfr_ceil",Pmpfr_ceil,-2}, + {"mpfr_round",Pmpfr_round,-2}, {0,0,0}, }; +int mpfr_roundmode = MPFR_RNDN; + void Ptodouble(NODE arg,Num *rp) { double r,i; @@ -43,9 +68,7 @@ void Ptodouble(NODE arg,Num *rp) } } -void Peval(arg,rp) -NODE arg; -Obj *rp; +void Peval(NODE arg,Obj *rp) { int prec; @@ -103,6 +126,43 @@ void Psetbprec(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) @@ -113,7 +173,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; } @@ -127,8 +187,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; } @@ -142,7 +202,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; } @@ -156,7 +216,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; } @@ -170,7 +230,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; } @@ -184,7 +244,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) @@ -197,7 +257,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) @@ -210,7 +270,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; } @@ -224,7 +284,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; } @@ -238,7 +298,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; } @@ -252,7 +312,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; } @@ -266,7 +326,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) @@ -279,7 +339,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) @@ -292,7 +352,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; } @@ -306,7 +366,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; } @@ -320,7 +380,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; } @@ -335,6 +395,194 @@ 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; +} + +#define SETPREC \ + (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); + + +void Pmpfr_gamma(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_lngamma(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_digamma(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_zeta(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_eint(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_erf(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_j0(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_j1(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_y0(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_y1(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_li2(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_ai(NODE arg,BF *rp) +{ + Num a; + int prec; + BF r; + + SETPREC + mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; +} + +void Pmpfr_floor(NODE arg,Q *rp) +{ + Num a; + int prec; + BF r; + mpz_t t; + GZ rz; + + SETPREC + mpfr_floor(r->body,((BF)a)->body); + mpz_init(t); + mpfr_get_z(t,r->body,mpfr_roundmode); + MPZTOGZ(t,rz); + *rp = gztoz(rz); +} + +void Pmpfr_ceil(NODE arg,Q *rp) +{ + Num a; + int prec; + BF r; + mpz_t t; + GZ rz; + + SETPREC + mpfr_ceil(r->body,((BF)a)->body); + mpz_init(t); + mpfr_get_z(t,r->body,mpfr_roundmode); + MPZTOGZ(t,rz); + *rp = gztoz(rz); +} + +void Pmpfr_round(NODE arg,Q *rp) +{ + Num a; + int prec; + BF r; + mpz_t t; + GZ rz; + + SETPREC + mpfr_round(r->body,((BF)a)->body); + mpz_init(t); + mpfr_get_z(t,r->body,mpfr_roundmode); + MPZTOGZ(t,rz); + *rp = gztoz(rz); }