OpenXM_contrib2/asir2000/builtin/bfaux.c - diff

Return to bfaux.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/bfaux.c between version 1.13 and 1.17

version 1.13, 2017/03/09 00:46:44

version 1.17, 2018/03/29 01:32:50

Line 1

/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.12 2016/03/14 04:15:05 noro Exp $ */

/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.16 2018/03/28 05:27:22 noro Exp $ */

#include "ca.h"

#include "parse.h"

Line 11 void Pmpfr_y0(), Pmpfr_y1();

void Pmpfr_gamma(), Pmpfr_lngamma(), Pmpfr_digamma();

void Pmpfr_floor(), Pmpfr_round(), Pmpfr_ceil();

void Prk_ratmat();

void mp_sin(),mp_cos(),mp_tan(),mp_asin(),mp_acos(),mp_atan();

void mp_sinh(),mp_cosh(),mp_tanh(),mp_asinh(),mp_acosh(),mp_atanh();

void mp_exp(),mp_log(),mp_pow();

void mp_factorial(),mp_abs();

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_sin",mp_sin,-2},

{"mpfr_zeta",Pmpfr_zeta,-2},

{"mpfr_cos",mp_cos,-2},

{"mpfr_j0",Pmpfr_j0,-2},

{"mpfr_tan",mp_tan,-2},

{"mpfr_j1",Pmpfr_j1,-2},

{"mpfr_asin",mp_asin,-2},

{"mpfr_y0",Pmpfr_y0,-2},

{"mpfr_acos",mp_acos,-2},

{"mpfr_y1",Pmpfr_y1,-2},

{"mpfr_atan",mp_atan,-2},

{"mpfr_eint",Pmpfr_eint,-2},

{"mpfr_sinh",mp_sinh,-2},

{"mpfr_erf",Pmpfr_erf,-2},

{"mpfr_cosh",mp_cosh,-2},

{"mpfr_erfc",Pmpfr_erfc,-2},

{"mpfr_tanh",mp_tanh,-2},

{"mpfr_li2",Pmpfr_li2,-2},

{"mpfr_asinh",mp_asinh,-2},

{"mpfr_gamma",Pmpfr_gamma,-2},

{"mpfr_acosh",mp_acosh,-2},

{"mpfr_lngamma",Pmpfr_gamma,-2},

{"mpfr_atanh",mp_atanh,-2},

{"mpfr_digamma",Pmpfr_gamma,-2},

{"mpfr_exp",mp_exp,-2},

{"mpfr_floor",Pmpfr_floor,-2},

{"mpfr_log",mp_log,-2},

{"mpfr_ceil",Pmpfr_ceil,-2},

{"mpfr_pow",mp_pow,-3},

{"mpfr_round",Pmpfr_round,-2},

{"mpfr_ai",Pmpfr_ai,-2},

{"rk_ratmat",Prk_ratmat,7},

{"mpfr_zeta",Pmpfr_zeta,-2},

{0,0,0},

{"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},

{"rk_ratmat",Prk_ratmat,7},

{0,0,0},

};

int mpfr_roundmode = MPFR_RNDN;

void Ptodouble(NODE arg,Num *rp)

void todoublen(Num a,Num *rp)

{

double r,i;

Real real,imag;

Num num;

asir_assert(ARG0(arg),O_N,"todouble");

if ( !a ) {

num = (Num)ARG0(arg);

*rp = 0;

if ( !num ) {

return;

*rp = 0;

}

return;

switch ( NID(a) ) {

}

case N_R: case N_Q: case N_B:

switch ( NID(num) ) {

r = ToReal(a);

case N_R: case N_Q: case N_B:

MKReal(r,real);

r = ToReal(num);

*rp = (Num)real;

MKReal(r,real);

break;

*rp = (Num)real;

case N_C:

break;

r = ToReal(((C)a)->r);

case N_C:

i = ToReal(((C)a)->i);

r = ToReal(((C)num)->r);

MKReal(r,real);

i = ToReal(((C)num)->i);

MKReal(i,imag);

MKReal(r,real);

reimtocplx((Num)real,(Num)imag,rp);

MKReal(i,imag);

break;

reimtocplx((Num)real,(Num)imag,rp);

default:

break;

*rp = a;

default:

break;

*rp = num;

}

break;

}

void todoublep(P a,P *rp)

{

DCP dc,dcr,dcr0;

if ( !a ) *rp = 0;

else if ( OID(a) == O_N ) todoublen((Num)a,(Num *)rp);

else {

for ( dcr0 = 0, dc = DC(a); dc; dc = NEXT(dc) ) {

NEXTDC(dcr0,dcr);

DEG(dcr) = DEG(dc);

todoublep(COEF(dc),&COEF(dcr));

}

NEXT(dcr) = 0;

MKP(VR(a),dcr0,*rp);

}

void todoubler(R a,R *rp)

{

R b;

if ( !a ) *rp = 0;

else if ( OID(a) <= O_P ) todoublep((P)a,(P *)rp);

else {

NEWR(b);

todoublep(a->nm,&b->nm);

todoublep(a->dn,&b->dn);

*rp = b;

}

void todouble(Obj a,Obj *b)

{

Obj t;

LIST l;

V v;

int row,col,len;

VECT vect;

MAT mat;

int i,j;

NODE n0,n,nd;

MP m,mp,mp0;

DP d;

if ( !a ) {

*b = 0;

return;

}

switch ( OID(a) ) {

case O_N:

todoublen((Num)a,(Num *)b);

break;

case O_P:

todoublep((P)a,(P *)b);

break;

case O_R:

todoubler((R)a,(R *)b);

break;

case O_LIST:

n0 = 0;

for ( nd = BDY((LIST)a); nd; nd = NEXT(nd) ) {

NEXTNODE(n0,n);

todouble((Obj)BDY(nd),(Obj *)&BDY(n));

}

if ( n0 )

NEXT(n) = 0;

MKLIST(l,n0);

*b = (Obj)l;

break;

case O_VECT:

len = ((VECT)a)->len;

MKVECT(vect,len);

for ( i = 0; i < len; i++ ) {

todouble((Obj)BDY((VECT)a)[i],(Obj *)&BDY(vect)[i]);

}

*b = (Obj)vect;

break;

case O_MAT:

row = ((MAT)a)->row;

col = ((MAT)a)->col;

MKMAT(mat,row,col);

for ( i = 0; i < row; i++ )

for ( j = 0; j < col; j++ ) {

todouble((Obj)BDY((MAT)a)[i][j],(Obj *)&BDY(mat)[i][j]);

}

*b = (Obj)mat;

break;

case O_DP:

mp0 = 0;

for ( m = BDY((DP)a); m; m = NEXT(m) ) {

todouble(C(m),&t);

if ( t ) {

NEXTMP(mp0,mp);

C(mp) = t;

mp->dl = m->dl;

}

if ( mp0 ) {

MKDP(NV((DP)a),mp0,d);

d->sugar = ((DP)a)->sugar;

*b = (Obj)d;

} else

*b = 0;

break;

default:

error("todouble : invalid argument");

}

void Ptodouble(NODE arg,Obj *rp)

{

todouble((Obj)ARG0(arg),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))*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");

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);

}

/* 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");

p = 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)

Line 135 void Psetround(NODE arg,Q *rp)

Line 266 void Psetround(NODE arg,Q *rp)

STOQ(mpfr_roundmode,*rp);

if ( arg ) {

asir_assert(ARG0(arg),O_N,"setround");

round = QTOS((Q)ARG0(arg));

switch ( round ) {

case 0:

Line 171 Num tobf(Num a,int prec);

Line 302 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;

}

Line 184 void mp_pi(NODE arg,BF *rp)

Line 315 void mp_pi(NODE arg,BF *rp)

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;

Line 223 void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)

Line 354 void mpfr_or_mpc(NODE arg,int (*mpfr_f)(),int (*mpc_f)

*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;

}

Line 301 void mp_log(NODE arg,Num *rp)

Line 432 void mp_log(NODE arg,Num *rp)

mpfr_or_mpc(arg,mpfr_log,mpc_log,rp);

}

void mp_abs(NODE arg,Num *rp)

{

mpfr_or_mpc(arg,mpfr_abs,mpc_abs,rp);

}

void mp_factorial(NODE arg,Num *rp)

{

struct oNODE arg0;

Num a,a1;

a = (Num)ARG0(arg);

if ( !a ) *rp = (Num)ONE;

else if ( INT(a) ) Pfac(arg,rp);

else {

addnum(0,a,(Num)ONE,&a1);

arg0.body = (pointer)a1;

arg0.next = arg->next;

Pmpfr_gamma(&arg0,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 ) {

Line 340 void mp_pow(NODE arg,Num *rp)

Line 492 void mp_pow(NODE arg,Num *rp)

*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;

}

Line 357 void mp_pow(NODE arg,Num *rp)

Line 509 void mp_pow(NODE arg,Num *rp)

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);

}

double **almat_double(int n)

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