From 4e895d89500fe926a7e0e992670b1e2fabb25e10 Mon Sep 17 00:00:00 2001 From: zimmerma Date: Wed, 5 Sep 2018 10:20:36 +0000 Subject: [src/exp_2.c] fix for 16-bit limb [tests/texp.c] improve error message git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@13135 280ebfd0-de03-0410-8827-d642c229c3f4 --- src/exp_2.c | 40 +++++++++++++++++++++++++++++++++++++--- tests/texp.c | 6 ++++-- 2 files changed, 41 insertions(+), 5 deletions(-) diff --git a/src/exp_2.c b/src/exp_2.c index 06d2e9d52..316b3d555 100644 --- a/src/exp_2.c +++ b/src/exp_2.c @@ -33,6 +33,42 @@ mpz_normalize (mpz_t, mpz_t, mpfr_exp_t); static mpfr_exp_t mpz_normalize2 (mpz_t, mpz_t, mpfr_exp_t, mpfr_exp_t); +/* count the number of significant bits of n, i.e., + nbits(unsigned long) - count_leading_zeros (n) */ +static int +nbits_ulong (unsigned long n) +{ + int nbits = 0; + + MPFR_ASSERTD (n > 0); + while (n >= 0x10000) + { + n >>= 16; + nbits += 16; + } + MPFR_ASSERTD (n <= 0xffff); + if (n >= 0x100) + { + n >>= 8; + nbits += 8; + } + MPFR_ASSERTD (n <= 0xff); + if (n >= 0x10) + { + n >>= 4; + nbits += 4; + } + MPFR_ASSERTD (n <= 0xf); + if (n >= 4) + { + n >>= 2; + nbits += 2; + } + MPFR_ASSERTD (n <= 3); + /* now n = 1, 2, or 3 */ + return nbits + 1 + (n >= 2); +} + /* if k = the number of bits of z > q, divides z by 2^(k-q) and returns k-q. Otherwise do nothing and return 0. */ @@ -146,9 +182,7 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) error_r = 0; else { - count_leading_zeros (error_r, - (mp_limb_t) SAFE_ABS (unsigned long, n) + 1); - error_r = GMP_NUMB_BITS - error_r; + error_r = nbits_ulong (SAFE_ABS (unsigned long, n) + 1); /* we have |x| <= 2^error_r * log(2) */ } diff --git a/tests/texp.c b/tests/texp.c index 57fee2801..8755f019b 100644 --- a/tests/texp.c +++ b/tests/texp.c @@ -618,6 +618,8 @@ bug20080731 (void) mpfr_set_str (x, "-2.c5c85fdf473de6af278ece700fcbdabd03cd0cb9ca62d8b62c@7", 16, MPFR_RNDN); + /* exp(x) is just below 0xf.fffffffffffffffp-1073741828 */ + mpfr_init2 (y1, prec); mpfr_exp (y1, x, MPFR_RNDU); @@ -629,9 +631,9 @@ bug20080731 (void) if (mpfr_cmp0 (y1, y2) != 0) { printf ("Error in bug20080731\nExpected "); - mpfr_out_str (stdout, 16, 0, y2, MPFR_RNDN); + mpfr_dump (y2); printf ("\nGot "); - mpfr_out_str (stdout, 16, 0, y1, MPFR_RNDN); + mpfr_dump (y1); printf ("\n"); exit (1); } -- cgit v1.2.1