=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/bfaux.c,v retrieving revision 1.16 retrieving revision 1.17 diff -u -p -r1.16 -r1.17 --- OpenXM_contrib2/asir2000/builtin/bfaux.c 2018/03/28 05:27:22 1.16 +++ OpenXM_contrib2/asir2000/builtin/bfaux.c 2018/03/29 01:32:50 1.17 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/builtin/bfaux.c,v 1.15 2017/08/31 04:21:48 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" @@ -17,74 +17,74 @@ 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_sin",mp_sin,-2}, - {"mpfr_cos",mp_cos,-2}, - {"mpfr_tan",mp_tan,-2}, - {"mpfr_asin",mp_asin,-2}, - {"mpfr_acos",mp_acos,-2}, - {"mpfr_atan",mp_atan,-2}, - {"mpfr_sinh",mp_sinh,-2}, - {"mpfr_cosh",mp_cosh,-2}, - {"mpfr_tanh",mp_tanh,-2}, - {"mpfr_asinh",mp_asinh,-2}, - {"mpfr_acosh",mp_acosh,-2}, - {"mpfr_atanh",mp_atanh,-2}, - {"mpfr_exp",mp_exp,-2}, - {"mpfr_log",mp_log,-2}, - {"mpfr_pow",mp_pow,-3}, - {"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}, - {"rk_ratmat",Prk_ratmat,7}, - {0,0,0}, + {"eval",Peval,-2}, + {"setprec",Psetprec,-1}, + {"setbprec",Psetbprec,-1}, + {"setround",Psetround,-1}, + {"todouble",Ptodouble,1}, + {"mpfr_sin",mp_sin,-2}, + {"mpfr_cos",mp_cos,-2}, + {"mpfr_tan",mp_tan,-2}, + {"mpfr_asin",mp_asin,-2}, + {"mpfr_acos",mp_acos,-2}, + {"mpfr_atan",mp_atan,-2}, + {"mpfr_sinh",mp_sinh,-2}, + {"mpfr_cosh",mp_cosh,-2}, + {"mpfr_tanh",mp_tanh,-2}, + {"mpfr_asinh",mp_asinh,-2}, + {"mpfr_acosh",mp_acosh,-2}, + {"mpfr_atanh",mp_atanh,-2}, + {"mpfr_exp",mp_exp,-2}, + {"mpfr_log",mp_log,-2}, + {"mpfr_pow",mp_pow,-3}, + {"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}, + {"rk_ratmat",Prk_ratmat,7}, + {0,0,0}, }; int mpfr_roundmode = MPFR_RNDN; void todoublen(Num a,Num *rp) { - double r,i; - Real real,imag; + double r,i; + Real real,imag; - if ( !a ) { - *rp = 0; - return; - } - switch ( NID(a) ) { - case N_R: case N_Q: case N_B: - r = ToReal(a); - MKReal(r,real); - *rp = (Num)real; - break; - case N_C: - r = ToReal(((C)a)->r); - i = ToReal(((C)a)->i); - MKReal(r,real); - MKReal(i,imag); - reimtocplx((Num)real,(Num)imag,rp); - break; - default: - *rp = a; - break; - } + if ( !a ) { + *rp = 0; + return; + } + switch ( NID(a) ) { + case N_R: case N_Q: case N_B: + r = ToReal(a); + MKReal(r,real); + *rp = (Num)real; + break; + case N_C: + r = ToReal(((C)a)->r); + i = ToReal(((C)a)->i); + MKReal(r,real); + MKReal(i,imag); + reimtocplx((Num)real,(Num)imag,rp); + break; + default: + *rp = a; + break; + } } void todoublep(P a,P *rp) @@ -120,81 +120,81 @@ void todoubler(R a,R *rp) 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; + 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: + 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; + 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"); - } + break; + default: + error("todouble : invalid argument"); + } } void Ptodouble(NODE arg,Obj *rp) @@ -206,58 +206,58 @@ void Peval(NODE arg,Obj *rp) { int prec; - asir_assert(ARG0(arg),O_R,"eval"); + asir_assert(ARG0(arg),O_R,"eval"); if ( argc(arg) == 2 ) { - prec = QTOS((Q)ARG1(arg))*3.32193; + 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); + evalr(CO,(Obj)ARG0(arg),prec,rp); } /* set/get decimal precision */ void Psetprec(NODE arg,Obj *rp) { - int p; - Q q; - int prec,dprec; + 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 ( 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); + mpfr_set_default_prec(prec); } /* set/get bit precision */ void Psetbprec(NODE arg,Obj *rp) { - int p; - Q q; - int prec; + 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 ( 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); + mpfr_set_default_prec(prec); } void Psetround(NODE arg,Q *rp) @@ -266,7 +266,7 @@ void Psetround(NODE arg,Q *rp) STOQ(mpfr_roundmode,*rp); if ( arg ) { - asir_assert(ARG0(arg),O_N,"setround"); + asir_assert(ARG0(arg),O_N,"setround"); round = QTOS((Q)ARG0(arg)); switch ( round ) { case 0: @@ -302,12 +302,12 @@ Num tobf(Num a,int prec); void mp_pi(NODE arg,BF *rp) { int prec; - BF r; + 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); + 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; } @@ -315,29 +315,29 @@ void mp_pi(NODE arg,BF *rp) void mp_e(NODE arg,BF *rp) { int prec; - mpfr_t one; - BF r; + 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); + 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; + Num a; int prec; - BF r,re,im; + 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); + 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; @@ -354,9 +354,9 @@ 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); + 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; } @@ -455,15 +455,15 @@ void mp_factorial(NODE arg,Num *rp) void mp_pow(NODE arg,Num *rp) { - Num a,e; + Num a,e; int prec; - BF r,re,im; + 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); + 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 ) { @@ -492,9 +492,9 @@ 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); + NEWBF(r); + mpfr_init2(r->body,prec); + mpfr_pow(r->body,((BF)a)->body,((BF)e)->body,mpfr_roundmode); *rp = (Num)r; } } @@ -509,193 +509,193 @@ void mp_pow(NODE arg,Num *rp) void Pmpfr_gamma(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_gamma(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_lngamma(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_lngamma(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_digamma(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_digamma(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_zeta(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_zeta(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_eint(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_eint(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_erf(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_erf(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_erfc(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_erfc(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_j0(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_j0(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_j1(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_j1(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_y0(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_y0(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_y1(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_y1(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_li2(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_li2(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_ai(NODE arg,BF *rp) { - Num a; + Num a; int prec; - BF r; + BF r; SETPREC - mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode); + mpfr_ai(r->body,((BF)a)->body,mpfr_roundmode); *rp = r; } void Pmpfr_floor(NODE arg,Q *rp) { - Num a; + Num a; int prec; - BF r; - mpz_t t; - GZ rz; + 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); + 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; + Num a; int prec; - BF r; - mpz_t t; - GZ rz; + 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); + 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; + Num a; int prec; - BF r; - mpz_t t; - GZ rz; + 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); + 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)