diff options
Diffstat (limited to 'cosh.c')
-rw-r--r-- | cosh.c | 158 |
1 files changed, 74 insertions, 84 deletions
@@ -23,99 +23,89 @@ MA 02111-1307, USA. */ #include "mpfr-impl.h" /* The computation of cosh is done by - - cosh= 1/2[e^(x)+e^(-x)] - */ + cosh= 1/2[e^(x)+e^(-x)] */ int mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) { - /****** Declaration ******/ - mpfr_t x; - mp_prec_t Nxt = MPFR_PREC(xt); - int inexact =0; + mpfr_t x; + mp_prec_t Nxt = MPFR_PREC(xt); + int inexact =0; + + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) + { + if (MPFR_IS_NAN(xt)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(xt)) + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + MPFR_RET(0); + } + else + { + MPFR_ASSERTD(MPFR_IS_ZERO(xt)); + return mpfr_set_ui (y, 1, rnd_mode); /* cosh(0) = 1 */ + } + } + + mpfr_init2 (x, Nxt); + mpfr_abs (x, xt, GMP_RNDN); - if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) + /* General case */ + { + /* Declaration of the intermediary variable */ + mpfr_t t, te; + + /* Declaration of the size variable */ + mp_prec_t Nx = Nxt; /* Precision of input variable */ + mp_prec_t Ny = MPFR_PREC(y); /* Precision of output 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 + 3 + __gmpfr_ceil_log2 (Nt); + + /* initialise of intermediary variables */ + mpfr_init2 (t, Nt); + mpfr_init2 (te, Nt); + + /* First computation of cosh */ + do { - if (MPFR_IS_NAN(xt)) - { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - else if (MPFR_IS_INF(xt)) - { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - MPFR_RET(0); - } - else - { - MPFR_ASSERTD(MPFR_IS_ZERO(xt)); - return mpfr_set_ui (y, 1, rnd_mode); /* cosh(0) = 1 */ - } + /* Reactualisation of the precision */ + mpfr_set_prec (t, Nt); + mpfr_set_prec (te, Nt); + + /* Compute cosh */ + mpfr_exp (te, x, GMP_RNDD); /* exp(x) */ + mpfr_ui_div (t, 1, te, GMP_RNDU); /* 1/exp(x) */ + mpfr_add (t, te, t, GMP_RNDN); /* exp(x) + 1/exp(x)*/ + mpfr_div_2ui (t, t, 1, GMP_RNDN); /* 1/2(exp(x) + 1/exp(x))*/ + + /* Estimation of the error */ + err = Nt - 3; + + /* Actualisation of the precision */ + Nt += 10; } - - mpfr_init2 (x, Nxt); - mpfr_abs (x, xt, GMP_RNDN); - - /* General case */ - { - /* Declaration of the intermediary variable */ - mpfr_t t, te,ti; - - /* Declaration of the size variable */ - mp_prec_t Nx = Nxt; /* 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 + 3 + __gmpfr_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 cosh */ - mpfr_exp (te, x, GMP_RNDD); /* exp(x) */ - mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */ - mpfr_add (t, te, ti, GMP_RNDN); /* exp(x) + 1/exp(x)*/ - mpfr_div_2ui (t, t, 1, GMP_RNDN); /* 1/2(exp(x) + 1/exp(x))*/ - - /* estimation of the error */ - err = Nt - 3; - - /* actualisation of the precision */ - Nt += 10; - - } - while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, + while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, Ny + (rnd_mode == GMP_RNDN))); - inexact = mpfr_set (y, t, rnd_mode); - - mpfr_clear (t); - mpfr_clear (ti); - mpfr_clear (te); - } - - mpfr_clear (x); + inexact = mpfr_set (y, t, rnd_mode); - return inexact; + mpfr_clear (te); + mpfr_clear (t); + } + + mpfr_clear (x); + + return inexact; } |