[src/digamma.c] Fixed a bug in mpfr_digamma_reflection.

[tests/tdigamma.c] Added testcase.
(merged changesets r14430,14432-14433 from the trunk)

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

2 files changed, 35 insertions, 13 deletions

diff --git a/src/digamma.c b/src/digamma.c
index b971930fa..0e8129b60 100644
--- a/src/digamma.c
+++ b/src/digamma.c
@@ -173,16 +173,19 @@ mpfr_digamma_reflection (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
       mpfr_digamma (v, u, MPFR_RNDN);   /* error <= 1/2 ulp */
       expv = MPFR_GET_EXP (v);
       mpfr_sub (v, v, t, MPFR_RNDN);
-      if (MPFR_GET_EXP (v) < MPFR_GET_EXP (t))
-        e1 += MPFR_EXP(t) - MPFR_EXP(v); /* scale error for t wrt new v */
-      /* now take into account the 1/2 ulp error for v */
-      if (expv - MPFR_EXP(v) - 1 > e1)
-        e1 = expv - MPFR_EXP(v) - 1;
-      else
-        e1 ++;
-      e1 ++; /* rounding error for mpfr_sub */
-      if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode))
-        break;
+      if (MPFR_NOTZERO(v))
+        {
+          if (MPFR_GET_EXP (v) < MPFR_GET_EXP (t))
+            e1 += MPFR_EXP(t) - MPFR_EXP(v); /* scale error for t wrt new v */
+          /* now take into account the 1/2 ulp error for v */
+          if (expv - MPFR_EXP(v) - 1 > e1)
+            e1 = expv - MPFR_EXP(v) - 1;
+          else
+            e1 ++;
+          e1 ++; /* rounding error for mpfr_sub */
+          if (MPFR_CAN_ROUND (v, p - e1, MPFR_PREC(y), rnd_mode))
+            break;
+        }
       MPFR_ZIV_NEXT (loop, p);
       mpfr_set_prec (t, p);
       mpfr_set_prec (v, p);
@@ -416,10 +419,8 @@ mpfr_digamma (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
         }
     }
 
-  if (MPFR_IS_NEG(x))
-    inex = mpfr_digamma_reflection (y, x, rnd_mode);
   /* if x < 1/2 we use the reflection formula */
-  else if (MPFR_EXP(x) < 0)
+  if (MPFR_IS_NEG(x) || MPFR_EXP(x) < 0)
     inex = mpfr_digamma_reflection (y, x, rnd_mode);
   else
     inex = mpfr_digamma_positive (y, x, rnd_mode);
diff --git a/tests/tdigamma.c b/tests/tdigamma.c
index 612e5e098..4c9fbd5c7 100644
--- a/tests/tdigamma.c
+++ b/tests/tdigamma.c
@@ -110,6 +110,26 @@ bug20210208 (void)
   mpfr_clear (y);
 }
 
+/* another test that fails with GMP_CHECK_RANDOMIZE=1613197421465830
+   on revision 14429 */
+static void
+bug20210215 (void)
+{
+  mpfr_t x, y;
+  int inex;
+
+  mpfr_init2 (x, 510);
+  mpfr_init2 (y, 4);
+  mpfr_set_str (x, "-8.2923051438433494998166335341807999322052669984208422481227138906096000469898717007386115912802685588348601663465077353194268894939972221117314512518182580e+35", 10, MPFR_RNDN);
+  mpfr_clear_flags ();
+  inex = mpfr_digamma (y, x, MPFR_RNDU);
+  MPFR_ASSERTN (mpfr_cmp_ui0 (y, 88) == 0);
+  MPFR_ASSERTN (inex > 0);
+  MPFR_ASSERTN (__gmpfr_flags == MPFR_FLAGS_INEXACT);
+  mpfr_clear (x);
+  mpfr_clear (y);
+}
+
 int
 main (int argc, char *argv[])
 {
@@ -118,6 +138,7 @@ main (int argc, char *argv[])
   special ();
   bug20210206 ();
   bug20210208 ();
+  bug20210215 ();
 
   test_generic (MPFR_PREC_MIN, 200, 20);