diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/bernoulli.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/bernoulli.c | 16 |
1 files changed, 8 insertions, 8 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c b/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c index d32069103bc..3656d80fdc8 100644 --- a/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c +++ b/Build/source/libs/mpfr/mpfr-src/src/bernoulli.c @@ -1,6 +1,6 @@ /* bernoulli -- internal function to compute Bernoulli numbers. -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. @@ -39,7 +39,7 @@ is_prime (unsigned long p) using Von Staudt–Clausen theorem, which says that the denominator of B[n] divides the product of all primes p such that p-1 divides n. Since B[n] = zeta(n) * 2*n!/(2pi)^n, we compute an approximation of - d * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ + (2n+1)! * zeta(n) * 2*n!/(2pi)^n and round it to the nearest integer. */ static void mpfr_bernoulli_internal (mpz_t *b, unsigned long n) { @@ -86,8 +86,9 @@ mpfr_bernoulli_internal (mpz_t *b, unsigned long n) mpfr_mul_ui (z, z, n, MPFR_RNDU); p = mpfr_get_ui (z, MPFR_RNDU); /* (n/e/2/pi)^n <= 2^p */ mpfr_clear (z); - /* the +14 term ensures no rounding failure up to n=10000 */ - prec += p + mpz_sizeinbase (den, 2) + 14; + MPFR_INC_PREC (prec, p + mpz_sizeinbase (den, 2)); + /* the +2 term ensures no rounding failure up to n=10000 */ + MPFR_INC_PREC (prec, __gmpfr_ceil_log2 (prec) + 2); } try_again: @@ -139,13 +140,13 @@ mpfr_bernoulli_internal (mpz_t *b, unsigned long n) Since z <= 2^prec * zeta(n) * 2*den*n!, ulp(z) <= 2*zeta(n) * 2*den*n!, thus (2^prec * zeta(n)-(p+1)) * 2*den*n! < z <= 2^prec * zeta(n) * 2*den*n! */ - mpfr_div_2exp (z, z, prec, MPFR_RNDZ); + mpfr_div_2ui (z, z, prec, MPFR_RNDZ); /* now (zeta(n) - (p+1)/2^prec) * 2*den*n! < z <= zeta(n) * 2*den*n! */ /* divide by (2pi)^n */ mpfr_init2 (y, prec); mpfr_const_pi (y, MPFR_RNDU); /* pi <= y <= pi * (1 + 2^(1-prec)) */ - mpfr_mul_2exp (y, y, 1, MPFR_RNDU); + mpfr_mul_2ui (y, y, 1, MPFR_RNDU); /* 2pi <= y <= 2pi * (1 + 2^(1-prec)) */ mpfr_pow_ui (y, y, n, MPFR_RNDU); /* (2pi)^n <= y <= (2pi)^n * (1 + 2^(1-prec))^(n+1) */ @@ -184,7 +185,6 @@ mpfr_bernoulli_internal (mpz_t *b, unsigned long n) mpz_mul_ui (t, t, n + 1); mpz_divexact (t, t, den); /* t was still n! */ mpz_mul (num, num, t); - mpz_set_ui (den, 1); mpfr_clear (y); mpfr_clear (z); @@ -194,7 +194,7 @@ mpfr_bernoulli_internal (mpz_t *b, unsigned long n) if (!ok) { - prec += prec / 10; + MPFR_INC_PREC (prec, prec / 10); goto try_again; } |