diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c | 96 |
1 files changed, 63 insertions, 33 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c b/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c index ce397785234..983dc77a4a6 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c +++ b/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c @@ -1,6 +1,6 @@ /* mpfr_mpn_exp -- auxiliary function for mpfr_get_str and mpfr_set_str -Copyright 1999-2019 Free Software Foundation, Inc. +Copyright 1999-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -24,7 +24,43 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., #define MPFR_NEED_LONGLONG_H #include "mpfr-impl.h" -/* this function computes an approximation of b^e in {a, n}, with exponent +/* count the number of significant bits of e, i.e., + nbits(mpfr_exp_t) - count_leading_zeros (e) */ +static int +nbits_mpfr_exp_t (mpfr_exp_t e) +{ + int nbits = 0; + + MPFR_ASSERTD (e > 0); + while (e >= 0x10000) + { + e >>= 16; + nbits += 16; + } + MPFR_ASSERTD (e <= 0xffff); + if (e >= 0x100) + { + e >>= 8; + nbits += 8; + } + MPFR_ASSERTD (e <= 0xff); + if (e >= 0x10) + { + e >>= 4; + nbits += 4; + } + MPFR_ASSERTD (e <= 0xf); + if (e >= 4) + { + e >>= 2; + nbits += 2; + } + MPFR_ASSERTD (e <= 3); + /* now e = 1, 2, or 3 */ + return nbits + 1 + (e >= 2); +} + +/* this function computes an approximation to b^e in {a, n}, with exponent stored in exp_r. The computed value is rounded toward zero (truncated). It returns an integer f such that the final error is bounded by 2^f ulps, that is: @@ -36,23 +72,23 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., Return -2 if an overflow occurred in the computation of exp_r. */ -long +int mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) { mp_limb_t *c, B; mpfr_exp_t f, h; int i; - unsigned long t; /* number of bits in e */ - unsigned long bits; - size_t n1; - unsigned int error; /* (number - 1) of loop a^2b inexact */ + int t; /* number of bits in e */ + size_t bits, n1; + unsigned int error; /* (number - 1) of loops a^2b inexact */ /* error == t means no error */ int err_s_a2 = 0; int err_s_ab = 0; /* number of error when shift A^2, AB */ MPFR_TMP_DECL(marker); - MPFR_ASSERTN(e > 0); - MPFR_ASSERTN((2 <= b) && (b <= 62)); + MPFR_ASSERTN (n > 0 && n <= ((size_t) -1) / GMP_NUMB_BITS); + MPFR_ASSERTN (e > 0); + MPFR_ASSERTN (2 <= b && b <= 62); MPFR_TMP_MARK(marker); @@ -75,17 +111,16 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) f = h - (n - 1) * GMP_NUMB_BITS; /* determine number of bits in e */ - count_leading_zeros (t, (mp_limb_t) e); + t = nbits_mpfr_exp_t (e); - t = GMP_NUMB_BITS - t; /* number of bits of exponent e */ - - error = t; /* error <= GMP_NUMB_BITS */ + error = t; + /* t, error <= bitsize(mpfr_exp_t) */ + MPFR_ASSERTD (error >= 0); MPN_ZERO (c, 2 * n); for (i = t - 2; i >= 0; i--) { - /* determine precision needed */ bits = n * GMP_NUMB_BITS - mpn_scan1 (a, 0); n1 = (n * GMP_NUMB_BITS - bits) / GMP_NUMB_BITS; @@ -98,19 +133,19 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) /* set {c+n, 2n1-n} to 0 : {c, n} = {a, n}^2*K^n */ /* check overflow on f */ - if (MPFR_UNLIKELY(f < MPFR_EXP_MIN/2 || f > MPFR_EXP_MAX/2)) + if (MPFR_UNLIKELY (f < MPFR_EXP_MIN / 2 || f > MPFR_EXP_MAX / 2)) { overflow: MPFR_TMP_FREE(marker); return -2; } - /* FIXME: Could f = 2*f + n * GMP_NUMB_BITS be used? */ - f = 2*f; + /* FIXME: Could f = 2 * f + n * GMP_NUMB_BITS be used? */ + f = 2 * f; MPFR_SADD_OVERFLOW (f, f, n * GMP_NUMB_BITS, mpfr_exp_t, mpfr_uexp_t, MPFR_EXP_MIN, MPFR_EXP_MAX, goto overflow, goto overflow); - if ((c[2*n - 1] & MPFR_LIMB_HIGHBIT) == 0) + if (MPFR_LIMB_MSB (c[2*n - 1]) == 0) { /* shift A by one bit to the left */ mpn_lshift (a, c + n, n, 1); @@ -122,11 +157,11 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) else MPN_COPY (a, c + n, n); - if ((error == t) && (2 * n1 <= n) && - (mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS)) + if (error == t && 2 * n1 <= n && + mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS) error = i; - if (e & ((mpfr_exp_t) 1 << i)) + if ((e >> i) & 1) { /* multiply A by B */ c[2 * n - 1] = mpn_mul_1 (c + n - 1, a, n, B); @@ -143,18 +178,22 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) if (error != t) err_s_ab ++; } - if ((error == t) && (c[n - 1] != 0)) + if (error == t && c[n - 1] != 0) error = i; } } + MPFR_ASSERTD (error >= 0); + MPFR_ASSERTD (err_s_a2 >= 0); + MPFR_ASSERTD (err_s_ab >= 0); + MPFR_TMP_FREE(marker); *exp_r = f; if (error == t) return -1; /* result is exact */ - else /* error <= t-2 <= GMP_NUMB_BITS-2 + else /* error <= t-2 <= bitsize(mpfr_exp_t)-2 err_s_ab, err_s_a2 <= t-1 */ { /* if there are p loops after the first inexact result, with @@ -162,15 +201,6 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n) at most 2^(p+ceil((j+1)/2)+l+1)*ulp(res). This is bounded by 2^(5/2*t-1/2) where t is the number of bits of e. */ - error = error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */ -#if 0 - if ((error - 1) >= ((n * GMP_NUMB_BITS - 1) / 2)) - error = n * GMP_NUMB_BITS; /* result is completely wrong: - this is very unlikely since error is - at most 5/2*log_2(e), and - n * GMP_NUMB_BITS is at least - 3*log_2(e) */ -#endif - return error; + return error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */ } } |