diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/eint.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/eint.c | 44 |
1 files changed, 34 insertions, 10 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/eint.c b/Build/source/libs/mpfr/mpfr-src/src/eint.c index d14180e3344..4df65821d6f 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/eint.c +++ b/Build/source/libs/mpfr/mpfr-src/src/eint.c @@ -1,6 +1,6 @@ /* mpfr_eint, mpfr_eint1 -- the exponential integral -Copyright 2005-2019 Free Software Foundation, Inc. +Copyright 2005-2020 Free Software Foundation, Inc. Contributed by the AriC and Caramba projects, INRIA. This file is part of the GNU MPFR Library. @@ -31,7 +31,8 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., */ /* Compute in y an approximation of sum(x^k/k/k!, k=1..infinity), - and return e such that the absolute error is bound by 2^e ulp(y). + assuming x != 0, and return e such that the absolute error is + bounded by 2^e ulp(y). Return PREC(y) when the truncated series does not converge. */ static mpfr_exp_t @@ -45,6 +46,10 @@ mpfr_eint_aux (mpfr_t y, mpfr_srcptr x) unsigned long k; MPFR_GROUP_DECL (group); + MPFR_LOG_FUNC ( + ("x[%Pu]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x), + ("y[%Pu]=%.*Rg", mpfr_get_prec (y), mpfr_log_prec, y)); + /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x) where |R(x)| <= (x/2)^2/(1-|x|/2) <= 2*(x/2)^2 thus |R(x)/x| <= |x|/2 @@ -61,18 +66,32 @@ mpfr_eint_aux (mpfr_t y, mpfr_srcptr x) mpz_init (u); mpz_init (m); MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs); - e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */ - MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x)); + e = mpfr_get_z_2exp (m, x); /* x = m * 2^e with m != 0 */ + MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e)); + MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x)); /* since m != 0 */ if (MPFR_PREC (x) > w) { e += MPFR_PREC (x) - w; - mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w); + mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w); /* one still has m != 0 */ + MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e)); } - /* remove trailing zeroes from m: this will speed up much cases where - x is a small integer divided by a power of 2 */ - k = mpz_scan1 (m, 0); - mpz_tdiv_q_2exp (m, m, k); - e += k; + /* Remove trailing zeroes from m: this will speed up much cases where + x is a small integer divided by a power of 2. + Note: As shown above, m != 0. This is needed for the "e += ..." below, + otherwise n would take the largest value of mp_bitcnt_t and could be + too large. */ + { + mp_bitcnt_t n = mpz_scan1 (m, 0); + mpz_tdiv_q_2exp (m, m, n); + /* Since one initially has mpz_sizeinbase (m, 2) == MPFR_PREC (x) + and m has not increased, one can deduce that n <= MPFR_PREC (x), + so that the cast to mpfr_prec_t is valid. This cast is needed to + ensure that the operand e of the addition below is not converted + to an unsigned integer type, which could yield incorrect results + with some C implementations. */ + MPFR_ASSERTD (n <= MPFR_PREC (x)); + e += (mpfr_prec_t) n; + } /* initialize t to 2^w */ mpz_set_ui (t, 1); mpz_mul_2exp (t, t, w); @@ -145,6 +164,7 @@ mpfr_eint_aux (mpfr_t y, mpfr_srcptr x) mpz_clear (t); mpz_clear (u); mpz_clear (m); + MPFR_LOG_MSG (("e=%" MPFR_EXP_FSPEC "d\n", (mpfr_eexp_t) e)); return e; } @@ -161,6 +181,10 @@ mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x) unsigned long k; mpfr_exp_t err_exp; + MPFR_LOG_FUNC ( + ("x[%Pu]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x), + ("err_exp=%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) err_exp)); + mpfr_init2 (t, p); mpfr_init2 (invx, p); mpfr_init2 (err, 31); /* error in ulps on y */ |