diff -Naurd mpfr-4.0.2-a/PATCHES mpfr-4.0.2-b/PATCHES
--- mpfr-4.0.2-a/PATCHES	2020-06-02 11:26:59.031114881 +0000
+++ mpfr-4.0.2-b/PATCHES	2020-06-02 11:26:59.175114040 +0000
@@ -0,0 +1 @@
+bernoulli-ziv
diff -Naurd mpfr-4.0.2-a/VERSION mpfr-4.0.2-b/VERSION
--- mpfr-4.0.2-a/VERSION	2020-04-03 13:54:03.891945830 +0000
+++ mpfr-4.0.2-b/VERSION	2020-06-02 11:26:59.175114040 +0000
@@ -1 +1 @@
-4.0.2-p7
+4.0.2-p8
diff -Naurd mpfr-4.0.2-a/doc/README.dev mpfr-4.0.2-b/doc/README.dev
--- mpfr-4.0.2-a/doc/README.dev	2019-01-07 14:06:05.000000000 +0000
+++ mpfr-4.0.2-b/doc/README.dev	2020-06-02 11:26:59.055114741 +0000
@@ -540,7 +540,9 @@
                         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 -Naurd mpfr-4.0.2-a/src/bernoulli.c mpfr-4.0.2-b/src/bernoulli.c
--- mpfr-4.0.2-a/src/bernoulli.c	2019-01-07 13:53:20.000000000 +0000
+++ mpfr-4.0.2-b/src/bernoulli.c	2020-06-02 11:26:59.055114741 +0000
@@ -39,7 +39,7 @@
    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_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 @@
   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 -Naurd mpfr-4.0.2-a/src/mpfr.h mpfr-4.0.2-b/src/mpfr.h
--- mpfr-4.0.2-a/src/mpfr.h	2020-04-03 13:54:03.891945830 +0000
+++ mpfr-4.0.2-b/src/mpfr.h	2020-06-02 11:26:59.175114040 +0000
@@ -27,7 +27,7 @@
 #define MPFR_VERSION_MAJOR 4
 #define MPFR_VERSION_MINOR 0
 #define MPFR_VERSION_PATCHLEVEL 2
-#define MPFR_VERSION_STRING "4.0.2-p7"
+#define MPFR_VERSION_STRING "4.0.2-p8"
 
 /* User macros:
    MPFR_USE_FILE:        Define it to make MPFR define functions dealing
diff -Naurd mpfr-4.0.2-a/src/version.c mpfr-4.0.2-b/src/version.c
--- mpfr-4.0.2-a/src/version.c	2020-04-03 13:54:03.891945830 +0000
+++ mpfr-4.0.2-b/src/version.c	2020-06-02 11:26:59.175114040 +0000
@@ -25,5 +25,5 @@
 const char *
 mpfr_get_version (void)
 {
-  return "4.0.2-p7";
+  return "4.0.2-p8";
 }
diff -Naurd mpfr-4.0.2-a/tests/tgamma.c mpfr-4.0.2-b/tests/tgamma.c
--- mpfr-4.0.2-a/tests/tgamma.c	2019-01-07 13:53:20.000000000 +0000
+++ mpfr-4.0.2-b/tests/tgamma.c	2020-06-02 11:26:59.055114741 +0000
@@ -1055,6 +1055,48 @@
   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 @@
 
   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;
 }