Fixed bug in the computation of the Bernoulli numbers, affecting4.0

mpfr_gamma (bug introduced in r8981). Added non-regression test
with mpfr_gamma.
  * doc/README.dev: document the new macro MPFR_CHECK_EXPENSIVE.
  * src/bernoulli.c: bug fix; initial precision improved.
  * tests/tgamma.c: non-regression test, performed only with
    MPFR_CHECK_EXPENSIVE and MPFR_CHECK_LARGEMEM set.
(merged changesets r13907-13912,13922 from the trunk)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/4.0@13924 280ebfd0-de03-0410-8827-d642c229c3f4

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;
 }