[src/lngamma.c] Fixed handling of reduced exponent range in mpfr_lgamma

  (also removed a useless cast).
[tests/tlgamma.c] Added a corresponding test case.
(merged changesets r12089-12092 from the trunk)

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

2 files changed, 58 insertions, 2 deletions

diff --git a/src/lngamma.c b/src/lngamma.c
index cc63672af..aaff525e3 100644
--- a/src/lngamma.c
+++ b/src/lngamma.c
@@ -818,6 +818,9 @@ mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mpfr_rnd_t rnd)
           int ok, inex2;
           mpfr_prec_t w = MPFR_PREC (y) + 14;
           mpfr_exp_t expl;
+          MPFR_SAVE_EXPO_DECL (expo);
+
+          MPFR_SAVE_EXPO_MARK (expo);
 
           while (1)
             {
@@ -847,13 +850,18 @@ mpfr_lgamma (mpfr_ptr y, int *signp, mpfr_srcptr x, mpfr_rnd_t rnd)
               mpfr_clear (l);
               mpfr_clear (h);
               if (ok)
-                return inex;
+                {
+                  MPFR_SAVE_EXPO_FREE (expo);
+                  return mpfr_check_range (y, inex, rnd);
+                }
               /* if ulp(log(-x)) <= |x| there is no reason to loop,
                  since the width of [l, h] will be at least |x| */
-              if (expl < MPFR_EXP(x) + (mpfr_exp_t) w)
+              if (expl < MPFR_EXP (x) + w)
                 break;
               w += MPFR_INT_CEIL_LOG2(w) + 3;
             }
+
+          MPFR_SAVE_EXPO_FREE (expo);
         }
     }
 
diff --git a/tests/tlgamma.c b/tests/tlgamma.c
index 20cb37ab6..048d3c7a9 100644
--- a/tests/tlgamma.c
+++ b/tests/tlgamma.c
@@ -405,12 +405,60 @@ tcc_bug20160606 (void)
   mpfr_clear (y);
 }
 
+/* With r12088, mpfr_lgamma is much too slow with a reduced emax that
+   yields an overflow, even though this case is easier. In practice,
+   this test will hang. */
+static void
+bug20180110 (void)
+{
+  mpfr_exp_t emax, e;
+  mpfr_t x, y, z;
+  mpfr_flags_t flags, eflags;
+  int i, inex, sign;
+
+  emax = mpfr_get_emax ();
+
+  mpfr_init2 (x, 2);
+  mpfr_inits2 (64, y, z, (mpfr_ptr) 0);
+  eflags = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
+
+  for (i = 10; i <= 30; i++)
+    {
+      mpfr_set_si_2exp (x, -1, -(1L << i), MPFR_RNDN);  /* -2^(-2^i) */
+      mpfr_lgamma (y, &sign, x, MPFR_RNDZ);
+      e = mpfr_get_exp (y);
+      mpfr_set_emax (e - 1);
+      mpfr_clear_flags ();
+      inex = mpfr_lgamma (y, &sign, x, MPFR_RNDZ);
+      flags = __gmpfr_flags;
+      mpfr_set_inf (z, 1);
+      mpfr_nextbelow (z);
+      mpfr_set_emax (emax);
+      if (! (mpfr_equal_p (y, z) && SAME_SIGN (inex, -1) && flags == eflags))
+        {
+          printf ("Error in bug20180110 for i = %d:\n", i);
+          printf ("Expected ");
+          mpfr_dump (z);
+          printf ("with inex = %d and flags =", -1);
+          flags_out (eflags);
+          printf ("Got      ");
+          mpfr_dump (y);
+          printf ("with inex = %d and flags =", inex);
+          flags_out (flags);
+          exit (1);
+        }
+    }
+
+  mpfr_clears (x, y, z, (mpfr_ptr) 0);
+}
+
 int
 main (void)
 {
   tests_start_mpfr ();
 
   tcc_bug20160606 ();
+  bug20180110 ();
 
   special ();
   test_generic (MPFR_PREC_MIN, 100, 2);