[src/digamma.c] Fixed bug in mpfr_digamma: if x is large, the

intermediate precision can be large. This issue may still occur,
but with a very low probability.
[tests/tdigamma.c] Testcase.
(merged changesets r14391-14398 from the trunk)

git-svn-id: https://scm.gforge.inria.fr/anonscm/svn/mpfr/branches/4.1@14424 280ebfd0-de03-0410-8827-d642c229c3f4

2 files changed, 76 insertions, 11 deletions

diff --git a/src/digamma.c b/src/digamma.c
index 9d9bae8d4..6dba748e6 100644
--- a/src/digamma.c
+++ b/src/digamma.c
@@ -214,19 +214,27 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
      ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y, inex));
 
-  /* compute a precision q such that x+1 is exact */
-  if (MPFR_PREC(x) < MPFR_GET_EXP(x))
-    q = MPFR_EXP(x);
-  else
-    q = MPFR_PREC(x) + 1;
-
-  /* for very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)) */
-  if (MPFR_PREC(y) + 10 < MPFR_EXP(x))
+  /* For very large x, use |digamma(x) - log(x)| < 1/x < 2^(1-EXP(x)).
+     However, for a fixed value of GUARD, MPFR_CAN_ROUND() might fail
+     with probability 1/2^GUARD, in which case the default code will
+     fail since it requires x+1 to be exact, thus a huge precision if
+     x is huge. There are two workarounds:
+     * either perform a Ziv's loop, by increasing GUARD at each step.
+       However, this might fail if x is moderately large, in which case
+       more terms of the asymptotic expansion would be needed.
+     * implement a full asymptotic expansion (with Ziv's loop). */
+#define GUARD 30
+  if (MPFR_PREC(y) + GUARD < MPFR_EXP(x))
     {
       /* this ensures EXP(x) >= 3, thus x >= 4, thus log(x) > 1 */
-      mpfr_init2 (t, MPFR_PREC(y) + 10);
-      mpfr_log (t, x, MPFR_RNDZ);
-      if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + 10, MPFR_PREC(y), rnd_mode))
+      mpfr_init2 (t, MPFR_PREC(y) + GUARD);
+      mpfr_log (t, x, MPFR_RNDN);
+      /* |t - digamma(x)| <= 1/2*ulp(t) + |digamma(x) - log(x)|
+                          <= 1/2*ulp(t) + 2^(1-EXP(x))
+                          <= 1/2*ulp(t) + 2^(-PREC(y)-GUARD)
+                          <= ulp(t)
+         since |t| >= 1 thus ulp(t) >= 2^(1-PREC(y)-GUARD) */
+      if (MPFR_CAN_ROUND (t, MPFR_PREC(y) + GUARD, MPFR_PREC(y), rnd_mode))
         {
           inex = mpfr_set (y, t, rnd_mode);
           mpfr_clear (t);
@@ -235,6 +243,21 @@ mpfr_digamma_positive (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
       mpfr_clear (t);
     }
 
+  /* compute a precision q such that x+1 is exact */
+  if (MPFR_PREC(x) < MPFR_GET_EXP(x))
+    {
+      /* The goal of the first assertion is to let the compiler ignore
+         the second one when MPFR_EMAX_MAX <= MPFR_PREC_MAX. */
+      MPFR_ASSERTD (MPFR_EXP(x) <= MPFR_EMAX_MAX);
+      MPFR_ASSERTN (MPFR_EXP(x) <= MPFR_PREC_MAX);
+      q = MPFR_EXP(x);
+    }
+  else
+    q = MPFR_PREC(x) + 1;
+
+  /* FIXME: q can be much too large, e.g. equal to the maximum exponent! */
+  MPFR_LOG_MSG (("q=%Pu\n", q));
+
   mpfr_init2 (x_plus_j, q);
 
   mpfr_init2 (t, p);
diff --git a/tests/tdigamma.c b/tests/tdigamma.c
index 83eb448de..6841220b0 100644
--- a/tests/tdigamma.c
+++ b/tests/tdigamma.c
@@ -49,12 +49,54 @@ special (void)
   mpfr_clear (y);
 }
 
+/* With some GMP_CHECK_RANDOMIZE values, test_generic triggers an error
+     tests_addsize(): too much memory (576460752303432776 bytes)
+  Each time on prec = 200, n = 3, xprec = 140.
+  The following test is a more general testcase.
+*/
+static void
+bug20210206 (void)
+{
+#define NPREC 4
+  mpfr_t x, y[NPREC], z;
+  mpfr_exp_t emin, emax;
+  int i, precx, precy[NPREC] = { 200, 400, 520, 1416 };
+
+  emin = mpfr_get_emin ();
+  emax = mpfr_get_emax ();
+  set_emin (MPFR_EMIN_MIN);
+  set_emax (MPFR_EMAX_MAX);
+
+  for (i = 0; i < NPREC; i++)
+    mpfr_init2 (y[i], precy[i]);
+  mpfr_init2 (z, precy[0]);
+
+  for (precx = MPFR_PREC_MIN; precx < 150; precx++)
+    {
+      mpfr_init2 (x, precx);
+      mpfr_setmax (x, __gmpfr_emax);
+      for (i = 0; i < NPREC; i++)
+        mpfr_digamma (y[i], x, MPFR_RNDA);
+      mpfr_set (z, y[1], MPFR_RNDA);
+      MPFR_ASSERTN (mpfr_equal_p (y[0], z));
+      mpfr_clear (x);
+    }
+
+  for (i = 0; i < NPREC; i++)
+    mpfr_clear (y[i]);
+  mpfr_clear (z);
+
+  set_emin (emin);
+  set_emax (emax);
+}
+
 int
 main (int argc, char *argv[])
 {
   tests_start_mpfr ();
 
   special ();
+  bug20210206 ();
 
   test_generic (MPFR_PREC_MIN, 200, 20);