diff options
Diffstat (limited to 'mpfr/asinh.c')
| -rw-r--r-- | mpfr/asinh.c | 137 |
1 files changed, 63 insertions, 74 deletions
diff --git a/mpfr/asinh.c b/mpfr/asinh.c index 048a01db0..8ca8f41d7 100644 --- a/mpfr/asinh.c +++ b/mpfr/asinh.c @@ -1,4 +1,4 @@ -/* mpfr_asinh -- Inverse Hyperbolic Sinus of Unsigned Integer Number +/* mpfr_asinh -- inverse hyperbolic sine Copyright 2001, 2002 Free Software Foundation. @@ -26,17 +26,18 @@ MA 02111-1307, USA. */ /* The computation of asinh is done by - asinh= ln(x+sqrt(x^2+1)) + asinh = ln(x + sqrt(x^2 + 1)) */ int -mpfr_asinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) +mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int inexact; - mpfr_t x; - int flag_neg=0; - mp_prec_t Nx; + int neg = 0; + mp_prec_t Nx, Ny, Nt; + mpfr_t t, te, ti; /* auxiliary variables */ + long int err; - if (MPFR_IS_NAN(xt)) + if (MPFR_IS_NAN(x)) { MPFR_SET_NAN(y); MPFR_RET_NAN; @@ -44,87 +45,75 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) MPFR_CLEAR_NAN(y); - if (MPFR_IS_INF(xt)) + if (MPFR_IS_INF(x)) { MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y, xt); + MPFR_SET_SAME_SIGN(y, x); MPFR_RET(0); } MPFR_CLEAR_INF(y); - if (MPFR_IS_ZERO(xt)) + if (MPFR_IS_ZERO(x)) { MPFR_SET_ZERO(y); /* asinh(0) = 0 */ - MPFR_SET_SAME_SIGN(y, xt); + MPFR_SET_SAME_SIGN(y, x); MPFR_RET(0); } - Nx = MPFR_PREC(xt); /* Precision of input variable */ - mpfr_init2(x, Nx); - mpfr_set(x, xt, GMP_RNDN); + Nx = MPFR_PREC(x); /* Precision of input variable */ + Ny = MPFR_PREC(y); /* Precision of output variable */ - if (MPFR_SIGN(x) < 0) - { - MPFR_CHANGE_SIGN(x); - flag_neg=1; - } + neg = MPFR_SIGN(x) < 0; /* General case */ - { - /* Declaration of the intermediary variable */ - mpfr_t t, te,ti; - /* Declaration of the size variable */ - mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */ - mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */ - - mp_prec_t Nt; /* Precision of the intermediary variable */ - long int err; /* Precision of error */ - - /* compute the precision of intermediary variable */ - Nt=MAX(Nx,Ny); - /* the optimal number of bits : see algorithms.ps */ - Nt=Nt+4+_mpfr_ceil_log2(Nt); - - /* initialise of intermediary variable */ - mpfr_init(t); - mpfr_init(te); - mpfr_init(ti); - - /* First computation of cosh */ - do { - - /* reactualisation of the precision */ - mpfr_set_prec(t,Nt); - mpfr_set_prec(te,Nt); - mpfr_set_prec(ti,Nt); - - /* compute asinh */ - mpfr_mul(te,x,x,GMP_RNDD); /* (x^2) */ - mpfr_add_ui(ti,te,1,GMP_RNDD); /* (x^2+1) */ - mpfr_sqrt(t,ti,GMP_RNDN); /* sqrt(x^2+1) */ - mpfr_add(t,t,x,GMP_RNDN); /* sqrt(x^2+1)+x */ - mpfr_log(t,t,GMP_RNDN); /* ln(sqrt(x^2+1)+x)*/ - - /* estimation of the error see- algorithms.ps*/ - /*err=Nt-_mpfr_ceil_log2(1+pow(2,3-MPFR_EXP(t)));*/ - err=Nt-(MAX(3-MPFR_EXP(t),0)+1); - - /* actualisation of the precision */ - Nt += 10; - - } while ((err < 0) || (!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny) || (MPFR_IS_ZERO(t)))); - - if(flag_neg) - MPFR_CHANGE_SIGN(t); - - inexact = mpfr_set(y,t,rnd_mode); - - mpfr_clear(t); - mpfr_clear(ti); - mpfr_clear(te); - } - mpfr_clear(x); - MPFR_RET(inexact); + /* compute the precision of intermediary variable */ + Nt = MAX(Nx, Ny); + + /* the optimal number of bits : see algorithms.ps */ + Nt = Nt + 4 + __gmpfr_ceil_log2 (Nt); + + /* initialize intermediary variables */ + mpfr_init2 (t, 2); + mpfr_init2 (te, 2); + mpfr_init2 (ti, 2); + + mpfr_save_emin_emax (); + + /* First computation of asinh */ + do { + + /* reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + mpfr_set_prec (ti, Nt); + + /* compute asinh */ + mpfr_mul (te, x, x, GMP_RNDD); /* x^2 */ + mpfr_add_ui (ti, te, 1, GMP_RNDD); /* x^2+1 */ + mpfr_sqrt (t, ti, GMP_RNDN); /* sqrt(x^2+1) */ + ((neg) ? mpfr_sub : mpfr_add) (t, t, x, GMP_RNDN); /* sqrt(x^2+1)+x */ + mpfr_log (t, t, GMP_RNDN); /* ln(sqrt(x^2+1)+x)*/ + + /* estimation of the error -- see algorithms.ps */ + err = Nt - (MAX(3 - MPFR_EXP(t), 0) + 1); + + /* actualisation of the precision */ + Nt += 10; + + } while ((err < 0) || (!mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny) || (MPFR_IS_ZERO(t)))); + + mpfr_restore_emin_emax (); + + if (neg) + MPFR_CHANGE_SIGN(t); + + inexact = mpfr_set (y, t, rnd_mode); + + mpfr_clear (t); + mpfr_clear (ti); + mpfr_clear (te); + + return mpfr_check_range (y, inexact, rnd_mode); } |