=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v retrieving revision 1.2 retrieving revision 1.11 diff -u -p -r1.2 -r1.11 --- OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/05 01:10:38 1.2 +++ OpenXM_contrib2/asir2000/builtin/bfaux.c 2015/08/25 18:41:03 1.11 @@ -1,16 +1,43 @@ -/* $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.10 2015/08/20 08:42:07 noro Exp $ */ #include "ca.h" #include "parse.h" -void Peval(), Psetprec(), Ptodouble(); +void Peval(), Psetprec(), Psetbprec(), Ptodouble(), Psetround(); +void Pmpfr_ai(); +void Pmpfr_eint(), Pmpfr_erf(), Pmpfr_erfc(), 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_erfc",Pmpfr_erfc,-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; @@ -42,15 +69,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,18 +83,41 @@ 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"); + p = 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"); + asir_assert(ARG0(arg),O_N,"setbprec"); p = QTOS((Q)ARG0(arg)); if ( p > 0 ) prec = p; @@ -79,23 +127,61 @@ 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) { - int prec; + int prec; BF r; 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); - *rp = r; + mpfr_const_pi(r->body,mpfr_roundmode); + if ( !cmpbf((Num)r,0) ) r = 0; + *rp = r; } void mp_e(NODE arg,BF *rp) { - int prec; + int prec; mpfr_t one; BF r; @@ -103,214 +189,357 @@ 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); - *rp = r; + mpfr_set_ui(one,1,mpfr_roundmode); + mpfr_exp(r->body,one,mpfr_roundmode); + if ( !cmpbf((Num)r,0) ) r = 0; + *rp = r; } -void mp_sin(NODE arg,BF *rp) +void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp) { Num a; - int prec; - BF r; + int prec; + BF r,re,im; + C c; + mpc_t mpc,a1; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; + prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec(); 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); - *rp = r; + if ( a && NID(a)==N_C ) { + mpc_init2(mpc,prec); mpc_init2(a1,prec); + re = (BF)((C)a)->r; im = (BF)((C)a)->i; + mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode); + (*mpc_f)(mpc,a1,mpfr_roundmode); + MPFRTOBF(mpc_realref(mpc),re); + MPFRTOBF(mpc_imagref(mpc),im); + if ( !cmpbf((Num)re,0) ) re = 0; + if ( !cmpbf((Num)im,0) ) im = 0; + if ( !im ) + *rp = (Num)re; + else { + NEWC(c); c->r = (Num)re; c->i = (Num)im; + *rp = (Num)c; + } + } else { + NEWBF(r); + mpfr_init2(r->body,prec); + (*mpfr_f)(r->body,((BF)a)->body,mpfr_roundmode); + if ( !cmpbf((Num)r,0) ) r = 0; + *rp = (Num)r; + } } -void mp_cos(NODE arg,BF *rp) +void mp_sin(NODE arg,Num *rp) { + mpfr_or_mpc(arg,mpfr_sin,mpc_sin,rp); +} + +void mp_cos(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_cos,mpc_cos,rp); +} + +void mp_tan(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_tan,mpc_tan,rp); +} + +void mp_asin(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_asin,mpc_asin,rp); +} + +void mp_acos(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_acos,mpc_acos,rp); +} + +void mp_atan(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_atan,mpc_atan,rp); +} + +void mp_sinh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_sinh,mpc_sinh,rp); +} + +void mp_cosh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_cosh,mpc_cosh,rp); +} + +void mp_tanh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_tanh,mpc_tanh,rp); +} + +void mp_asinh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_asinh,mpc_asinh,rp); +} + +void mp_acosh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_acosh,mpc_acosh,rp); +} + +void mp_atanh(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_atanh,mpc_atanh,rp); +} + +void mp_exp(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_exp,mpc_exp,rp); +} + +void mp_log(NODE arg,Num *rp) +{ + mpfr_or_mpc(arg,mpfr_log,mpc_log,rp); +} + +void mp_pow(NODE arg,Num *rp) +{ + Num a,e; + int prec; + BF r,re,im; + C c; + mpc_t mpc,a1,e1; + + prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : mpfr_get_default_prec(); + a = tobf(ARG0(arg),prec); + e = tobf(ARG1(arg),prec); + if ( NID(a) == N_C || NID(e) == N_C || MPFR_SIGN(((BF)a)->body) < 0 ) { + mpc_init2(mpc,prec); mpc_init2(a1,prec); mpc_init2(e1,prec); + if ( NID(a) == N_C ) { + re = (BF)((C)a)->r; im = (BF)((C)a)->i; + mpc_set_fr_fr(a1,re->body,im->body,mpfr_roundmode); + } else { + re = (BF)a; + mpc_set_fr(a1,re->body,mpfr_roundmode); + } + if ( NID(e) == N_C ) { + re = (BF)((C)e)->r; im = (BF)((C)e)->i; + mpc_set_fr_fr(e1,re->body,im->body,mpfr_roundmode); + } else { + re = (BF)e; + mpc_set_fr(e1,re->body,mpfr_roundmode); + } + mpc_pow(mpc,a1,e1,mpfr_roundmode); + MPFRTOBF(mpc_realref(mpc),re); + MPFRTOBF(mpc_imagref(mpc),im); + if ( !cmpbf((Num)re,0) ) re = 0; + if ( !cmpbf((Num)im,0) ) im = 0; + if ( !im ) + *rp = (Num)re; + else { + NEWC(c); c->r = (Num)re; c->i = (Num)im; + *rp = (Num)c; + } + } else { + NEWBF(r); + mpfr_init2(r->body,prec); + mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode); + *rp = (Num)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; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_tan(NODE arg,BF *rp) +void Pmpfr_lngamma(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_asin(NODE arg,BF *rp) +void Pmpfr_digamma(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_acos(NODE arg,BF *rp) + +void Pmpfr_zeta(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_atan(NODE arg,BF *rp) + +void Pmpfr_eint(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_sinh(NODE arg,BF *rp) +void Pmpfr_erf(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_cosh(NODE arg,BF *rp) +void Pmpfr_erfc(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_tanh(NODE arg,BF *rp) +void Pmpfr_j0(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_asinh(NODE arg,BF *rp) +void Pmpfr_j1(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_acosh(NODE arg,BF *rp) + +void Pmpfr_y0(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_atanh(NODE arg,BF *rp) + +void Pmpfr_y1(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_exp(NODE arg,BF *rp) +void Pmpfr_li2(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_log(NODE arg,BF *rp) +void Pmpfr_ai(NODE arg,BF *rp) { Num a; - int prec; + int prec; BF r; - prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : 0; - 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); - *rp = r; + SETPREC + mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode); + *rp = r; } -void mp_pow(NODE arg,BF *rp) +void Pmpfr_floor(NODE arg,Q *rp) { - Num a,e; - int prec; + Num a; + int prec; BF r; + mpz_t t; + GZ rz; - prec = NEXT(NEXT(arg)) ? QTOS((Q)ARG2(arg)) : 0; - a = tobf(ARG0(arg),prec); - 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); - *rp = r; + 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); }