Fix ldbl-128 / ldbl-128ibm lgammal spurious underflow (bug 16400).

This patch fixes bug 16400, spurious underflow exceptions for ldbl-128
/ ldbl-128ibm lgammal with small positive arguments, by just using
-__logl (x) as the result in the problem cases (similar to the
previous fix for problems with small negative arguments).

Tested powerpc32, and also tested on mips64 that this does not require
ulps regeneration for the ldbl-128 case.

	* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
	Return -__logl (x) for small positive arguments without evaluating
	a polynomial.

diff --git a/ChangeLog b/ChangeLog
index 1984691a78855b19820271cf93abfefec9d825c6..8bdea171790408de1b5c4619adbf11f63ebebcb8 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2014-01-06  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #16400]
+       * sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
+       Return -__logl (x) for small positive arguments without evaluating
+       a polynomial.
+
 2014-01-06  Mike Frysinger  <vapier@gentoo.org>
 
        * sysdeps/unix/sysv/linux/sys/ptrace.h (ptrace_peeksiginfo_args):

diff --git a/NEWS b/NEWS
index e6c5020a079012014c1be38c75b040d8dd2218d5..ab3b882631f46d08106a7e3d999e763a7c3ed9df 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -24,7 +24,7 @@ Version 2.19
   16112, 16143, 16144, 16146, 16150, 16151, 16153, 16167, 16172, 16195,
   16214, 16245, 16271, 16274, 16283, 16289, 16293, 16314, 16316, 16330,
   16337, 16338, 16356, 16365, 16366, 16369, 16372, 16375, 16379, 16384,
-  16385, 16386, 16390.
+  16385, 16386, 16390, 16400.
 
 * Slovenian translations for glibc messages have been contributed by the
   Translation Project's Slovenian team of translators.

diff --git a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
index 23ab9b9a43977259c8f1c6f002a8c8fd347b7c38..1961355a73e4016db104a1b172a0e81db2d19930 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
+++ b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
@@ -805,7 +805,9 @@ __ieee754_lgammal_r (long double x, int *signgamp)
        {
        case 0:
          /* log gamma (x + 1) = log(x) + log gamma(x) */
-         if (x <= 0.125)
+         if (x < 0x1p-120L)
+           return -__logl (x);
+         else if (x <= 0.125)
            {
              p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
            }