diff options
author | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
---|---|---|
committer | Akira Kakuto <kakuto@fuk.kindai.ac.jp> | 2017-12-26 11:50:03 +0000 |
commit | 6ddc203e5368cf1d4108e64c18b0bf6c7d53d176 (patch) | |
tree | dda1d9b65325699c60418cfcd02fa5d800cde1f8 /Build/source/libs/mpfr/mpfr-src/src/li2.c | |
parent | bbfe0bab9260ac7a06a69f710b053dcc73cecad2 (diff) |
mpfr-4.0.0
git-svn-id: svn://tug.org/texlive/trunk@46140 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/li2.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/li2.c | 35 |
1 files changed, 16 insertions, 19 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/li2.c b/Build/source/libs/mpfr/mpfr-src/src/li2.c index d8191f7e16a..cfa2d382309 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/li2.c +++ b/Build/source/libs/mpfr/mpfr-src/src/li2.c @@ -33,18 +33,17 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., static int li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) { - int i, Bm, Bmax; + int i; mpfr_t s, u, v, w; mpfr_prec_t sump, p; mpfr_exp_t se, err; - mpz_t *B; MPFR_ZIV_DECL (loop); /* The series converges for |z| < 2 pi, but in mpfr_li2 the argument is reduced so that 0 < z <= log(2). Here is additionnal check that z is (nearly) correct */ MPFR_ASSERTD (MPFR_IS_STRICTPOS (z)); - MPFR_ASSERTD (mpfr_cmp_d (z, 0.6953125) <= 0); + MPFR_ASSERTD (mpfr_cmp_ui_2exp (z, 89, -7) <= 0); /* z <= 0.6953125 */ sump = MPFR_PREC (sum); /* target precision */ p = sump + MPFR_INT_CEIL_LOG2 (sump) + 4; /* the working precision */ @@ -53,9 +52,6 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) mpfr_init2 (v, p); mpfr_init2 (w, p); - B = mpfr_bernoulli_internal ((mpz_t *) 0, 0); - Bm = Bmax = 1; - MPFR_ZIV_INIT (loop, p); for (;;) { @@ -67,9 +63,6 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) for (i = 1;; i++) { - if (i >= Bmax) - B = mpfr_bernoulli_internal (B, Bmax++); /* B_2i*(2i+1)!, exact */ - mpfr_mul (v, u, v, MPFR_RNDU); mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); @@ -77,7 +70,7 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) mpfr_div_ui (v, v, 2 * i + 1, MPFR_RNDU); /* here, v_2i = v_{2i-2} / (2i * (2i+1))^2 */ - mpfr_mul_z (w, v, B[i], MPFR_RNDN); + mpfr_mul_z (w, v, mpfr_bernoulli_cache(i), MPFR_RNDN); /* here, w_2i = v_2i * B_2i * (2i+1)! with error(w_2i) < 2^(5 * i + 8) ulp(w_2i) (see algorithms.tex) */ @@ -107,10 +100,6 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) MPFR_ZIV_FREE (loop); mpfr_set (sum, s, rnd_mode); - Bm = Bmax; - while (Bm--) - mpz_clear (B[Bm]); - (*__gmp_free_func) (B, Bmax * sizeof (mpz_t)); mpfr_clears (s, u, v, w, (mpfr_ptr) 0); /* Let K be the returned value. @@ -129,7 +118,7 @@ li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) thus |Li2(x) - g(x)| <= 2/x. Assumes x >= 38, which ensures log(x)^2/2 >= 2*Pi^2/3, and g(x) <= -3.3. Return 0 if asymptotic expansion failed (unable to round), otherwise - returns correct ternary value. + returns 1 for RNDF, and correct ternary value otherwise. */ static int mpfr_li2_asympt_pos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) @@ -161,7 +150,11 @@ mpfr_li2_asympt_pos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) bounded by 16 ulp(g). */ if ((MPFR_EXP (x) >= (mpfr_exp_t) w - MPFR_EXP (g)) && MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) - inex = mpfr_set (y, g, rnd_mode); + { + inex = mpfr_set (y, g, rnd_mode); + if (rnd_mode == MPFR_RNDF) + inex = 1; + } mpfr_clear (g); mpfr_clear (h); @@ -175,7 +168,7 @@ mpfr_li2_asympt_pos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) |Li2(x) - g(x)| <= 1/|x|. Assumes x <= -7, which ensures |log(-x)^2/2| >= Pi^2/6, and g(x) <= -3.5. Return 0 if asymptotic expansion failed (unable to round), otherwise - returns correct ternary value. + returns 1 for RNDF, and correct ternary value otherwise. */ static int mpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) @@ -204,7 +197,11 @@ mpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) total error is bounded by 16 ulp(g). */ if ((MPFR_EXP (x) >= (mpfr_exp_t) (w - 2) - MPFR_EXP (g)) && MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) - inex = mpfr_neg (y, g, rnd_mode); + { + inex = mpfr_neg (y, g, rnd_mode); + if (rnd_mode == MPFR_RNDF) + inex = 1; + } mpfr_clear (g); mpfr_clear (h); @@ -264,7 +261,7 @@ mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) yp = MPFR_PREC (y); m = yp + MPFR_INT_CEIL_LOG2 (yp) + 13; - if (MPFR_LIKELY ((mpfr_cmp_ui (x, 0) > 0) && (mpfr_cmp_d (x, 0.5) <= 0))) + if (MPFR_LIKELY ((mpfr_cmp_ui (x, 0) > 0) && (mpfr_cmp_ui_2exp (x, 1, -1) <= 0))) /* 0 < x <= 1/2: Li2(x) = S(-log(1-x))-log^2(1-x)/4 */ { mpfr_t s, u; |