/* $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(), 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;
Real real,imag;
Num num;
asir_assert(ARG0(arg),O_N,"todouble");
num = (Num)ARG0(arg);
if ( !num ) {
*rp = 0;
return;
}
switch ( NID(num) ) {
case N_R: case N_Q: case N_B:
r = ToReal(num);
MKReal(r,real);
*rp = (Num)real;
break;
case N_C:
r = ToReal(((C)num)->r);
i = ToReal(((C)num)->i);
MKReal(r,real);
MKReal(i,imag);
reimtocplx((Num)real,(Num)imag,rp);
break;
default:
*rp = num;
break;
}
}
void Peval(NODE arg,Obj *rp)
{
int prec;
asir_assert(ARG0(arg),O_R,"eval");
if ( argc(arg) == 2 ) {
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
prec = 0;
evalr(CO,(Obj)ARG0(arg),prec,rp);
}
/* 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);
}
/* set/get bit precision */
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,"setbprec");
prec = QTOS((Q)ARG0(arg));
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);
}
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;
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_roundmode);
if ( !cmpbf((Num)r,0) ) r = 0;
*rp = r;
}
void mp_e(NODE arg,BF *rp)
{
int prec;
mpfr_t one;
BF r;
prec = arg ? QTOS((Q)ARG0(arg)) : 0;
NEWBF(r);
prec ? mpfr_init2(r->body,prec) : mpfr_init(r->body);
mpfr_init(one);
mpfr_set_ui(one,1,mpfr_roundmode);
mpfr_exp(r->body,one,mpfr_roundmode);
if ( !cmpbf((Num)r,0) ) r = 0;
*rp = r;
}
void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)(),Num *rp)
{
Num a;
int prec;
BF r,re,im;
C c;
mpc_t mpc,a1;
prec = NEXT(arg) ? QTOS((Q)ARG1(arg)) : mpfr_get_default_prec();
a = tobf(ARG0(arg),prec);
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_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;
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_erfc(NODE arg,BF *rp)
{
Num a;
int prec;
BF r;
SETPREC
mpfr_erfc(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);
}