diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2006-08-15 11:37:20 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2006-08-15 11:37:20 +0000 |
commit | 7c74c3ec0bc3f1d3a1ee664101498305b05d1236 (patch) | |
tree | 5f39e8bd23cc15c7cb8afe063168ea2893ad6450 | |
parent | 9dbf737924b7d31ba2254c1eb0af9dcaa8306bf3 (diff) | |
download | mpfr-7c74c3ec0bc3f1d3a1ee664101498305b05d1236.tar.gz |
Better fix for exp2 exponent range bug.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4119 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | exp2.c | 84 |
1 files changed, 40 insertions, 44 deletions
@@ -86,53 +86,49 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_set_ui (y, 1, GMP_RNDZ); inexact = mpfr_mul_2si (y, y, xd, rnd_mode); + } + else /* General case */ + { + /* Declaration of the intermediary variable */ + mpfr_t t; + + /* Declaration of the size variable */ + mp_prec_t Ny = MPFR_PREC(y); /* target precision */ + mp_prec_t Nt; /* working precision */ + mp_exp_t err; /* error */ + MPFR_ZIV_DECL (loop); + + /* compute the precision of intermediary variable */ + /* the optimal number of bits : see algorithms.tex */ + Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); + + /* initialise of intermediary variable */ + mpfr_init2 (t, Nt); - MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (y, inexact, rnd_mode); + /* First computation */ + MPFR_ZIV_INIT (loop, Nt); + for (;;) + { + /* compute exp(x*ln(2))*/ + mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ + mpfr_mul (t, x, t, GMP_RNDU); /* x*ln(2) */ + err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ + mpfr_exp (t, t, GMP_RNDN); /* exp(x*ln(2))*/ + + if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) + break; + + /* Actualisation of the precision */ + MPFR_ZIV_NEXT (loop, Nt); + mpfr_set_prec (t, Nt); + } + MPFR_ZIV_FREE (loop); + + inexact = mpfr_set (y, t, rnd_mode); + + mpfr_clear (t); } - /* General case */ - { - /* Declaration of the intermediary variable */ - mpfr_t t; - - /* Declaration of the size variable */ - mp_prec_t Ny = MPFR_PREC(y); /* target precision */ - mp_prec_t Nt; /* working precision */ - mp_exp_t err; /* error */ - MPFR_ZIV_DECL (loop); - - /* compute the precision of intermediary variable */ - /* the optimal number of bits : see algorithms.tex */ - Nt = Ny + 5 + MPFR_INT_CEIL_LOG2 (Ny); - - /* initialise of intermediary variable */ - mpfr_init2 (t, Nt); - - /* First computation */ - MPFR_ZIV_INIT (loop, Nt); - for (;;) - { - /* compute exp(x*ln(2))*/ - mpfr_const_log2 (t, GMP_RNDU); /* ln(2) */ - mpfr_mul (t, x, t, GMP_RNDU); /* x*ln(2) */ - err = Nt - (MPFR_GET_EXP (t) + 2); /* Estimate of the error */ - mpfr_exp (t, t, GMP_RNDN); /* exp(x*ln(2))*/ - - if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) - break; - - /* Actualisation of the precision */ - MPFR_ZIV_NEXT (loop, Nt); - mpfr_set_prec (t, Nt); - } - MPFR_ZIV_FREE (loop); - - inexact = mpfr_set (y, t, rnd_mode); - - mpfr_clear (t); - } MPFR_SAVE_EXPO_FREE (expo); - return mpfr_check_range (y, inexact, rnd_mode); } |