diff options
-rw-r--r-- | doc/README.dev | 4 | ||||
-rw-r--r-- | src/bernoulli.c | 8 | ||||
-rw-r--r-- | tests/tgamma.c | 47 |
3 files changed, 54 insertions, 5 deletions
diff --git a/doc/README.dev b/doc/README.dev index ae97da440..331ca57d5 100644 --- a/doc/README.dev +++ b/doc/README.dev @@ -543,7 +543,9 @@ Environment variables that affect the tests: By default, a fixed seed is used. Only developers and testers should change the seed. -+ MPFR_CHECK_LARGEMEM: Define to enable expensive tests. ++ MPFR_CHECK_LARGEMEM: Define to enable tests that take a lot of memory. + ++ MPFR_CHECK_EXPENSIVE: Define to enable tests that take a lot of time. + MPFR_CHECK_LIBC_PRINTF: Define to enable comparisons with the printf diff --git a/src/bernoulli.c b/src/bernoulli.c index ded7651b0..c52e45f10 100644 --- a/src/bernoulli.c +++ b/src/bernoulli.c @@ -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; + prec += p + mpz_sizeinbase (den, 2); + /* the +2 term ensures no rounding failure up to n=10000 */ + prec += __gmpfr_ceil_log2 (prec) + 2; } try_again: @@ -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); diff --git a/tests/tgamma.c b/tests/tgamma.c index da59dce51..b8c929e09 100644 --- a/tests/tgamma.c +++ b/tests/tgamma.c @@ -1055,6 +1055,48 @@ exp_lgamma_tests (void) set_emax (emax); } +/* Bug reported by Frithjof Blomquist on May 19, 2020. + For the record, this bug was present since r8981 + (in mpfr_bernoulli_internal, den was wrongly reset to 1 in case + of failure in Ziv's loop). The bug only occurred up from r8986 + where the initial precision was reduced, but was potentially + present in any case of failure of Ziv's loop. */ +static void +bug20200519 (void) +{ + mpfr_prec_t prec = 25093; + mpfr_t x, y, z, d; + double dd; + size_t min_memory_limit, old_memory_limit; + + old_memory_limit = tests_memory_limit; + min_memory_limit = 24000000; + if (tests_memory_limit > 0 && tests_memory_limit < min_memory_limit) + tests_memory_limit = min_memory_limit; + + mpfr_init2 (x, prec); + mpfr_init2 (y, prec); + mpfr_init2 (z, prec + 100); + mpfr_init2 (d, 24); + mpfr_set_d (x, 2.5, MPFR_RNDN); + mpfr_gamma (y, x, MPFR_RNDN); + mpfr_gamma (z, x, MPFR_RNDN); + mpfr_sub (d, y, z, MPFR_RNDN); + mpfr_mul_2si (d, d, prec - mpfr_get_exp (y), MPFR_RNDN); + dd = mpfr_get_d (d, MPFR_RNDN); + if (dd < -0.5 || 0.5 < dd) + { + printf ("Error in bug20200519: dd=%f\n", dd); + exit (1); + } + mpfr_clear (x); + mpfr_clear (y); + mpfr_clear (z); + mpfr_clear (d); + + tests_memory_limit = old_memory_limit; +} + int main (int argc, char *argv[]) { @@ -1086,6 +1128,11 @@ main (int argc, char *argv[]) data_check ("data/gamma", mpfr_gamma, "mpfr_gamma"); + /* this test takes about one minute */ + if (getenv ("MPFR_CHECK_EXPENSIVE") != NULL && + getenv ("MPFR_CHECK_LARGEMEM") != NULL) + bug20200519 (); + tests_end_mpfr (); return 0; } |