Fix lgammaf spurious underflows (bug 18220).

The flt-32 implementation of lgammaf produces spurious underflow
exceptions for some large arguments, because of calculations involving
x^-2 multiplied by small constants.  This patch fixes this by
adjusting the threshold for a simpler computation to 2**26 (the error
in the simpler computation is on the order of 0.5 * log (x), for a
result on the order of x * log (x)).

Tested for x86_64 and x86.

	[BZ #18220]
	* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Use
	2**26 not 2**58 as threshold for returning x * (log (x) - 1).
	* math/auto-libm-test-in: Add another test of lgamma.
	* math/auto-libm-test-out: Regenerated.

diff --git a/ChangeLog b/ChangeLog
index 7749398ad51de94915f677b13d65d3b7c0f27480..2b29214422bac94993bf35b0685eefe1075f3563 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2015-05-15  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #18220]
+       * sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Use
+       2**26 not 2**58 as threshold for returning x * (log (x) - 1).
+       * math/auto-libm-test-in: Add another test of lgamma.
+       * math/auto-libm-test-out: Regenerated.
+
 2015-05-15  Wilco Dijkstra  <wdijkstr@arm.com>
 
        * stdio-common/printf_fp.c (___printf_fp): Use abs.

diff --git a/NEWS b/NEWS
index bc6c0cda4f62504c9ad8ad3addf00432d09de22d..b85ba6017d5d7fa2c383b83b023f286b04810922 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -17,7 +17,7 @@ Version 2.22
   17999, 18007, 18019, 18020, 18029, 18030, 18032, 18036, 18038, 18039,
   18042, 18043, 18046, 18047, 18068, 18080, 18093, 18100, 18104, 18110,
   18111, 18125, 18128, 18138, 18185, 18196, 18197, 18206, 18210, 18211,
-  18217, 18247, 18287, 18319, 18333, 18346, 18397, 18409.
+  18217, 18220, 18247, 18287, 18319, 18333, 18346, 18397, 18409.
 
 * Cache information can be queried via sysconf() function on s390 e.g. with
   _SC_LEVEL1_ICACHE_SIZE as argument.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index c4bfe7439c1fe030052af29717535a80c4035963..2a88403aa72cc784cd7bafa2ba6d1ee0d82937bb 100644 (file)
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1757,6 +1757,7 @@ lgamma 0.5
 lgamma -0.5
 lgamma 0.7
 lgamma 1.2
+lgamma 0x3.8p56
 lgamma 0x1p-5
 lgamma -0x1p-5
 lgamma 0x1p-10

diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
index d717ea4431dcd8290abbc585beba334c39ec083c..0aa7b75c3b88dece4003dd19acf075bdf01deb26 100644 (file)
--- a/math/auto-libm-test-out
+++ b/math/auto-libm-test-out
@@ -139612,6 +139612,31 @@ lgamma 1.2
 = lgamma tonearest ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83d8p-4L 1 : inexact-ok
 = lgamma towardzero ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83dp-4L 1 : inexact-ok
 = lgamma upward ldbl-128ibm 0x1.33333333333333333333333333p+0L : -0x1.5db138c7d70c72cb0b57c9e83dp-4L 1 : inexact-ok
+lgamma 0x3.8p56
+= lgamma downward flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
+= lgamma tonearest flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
+= lgamma towardzero flt-32 0x3.8p+56f : 0x8.8bdd4p+60f 1 : inexact-ok
+= lgamma upward flt-32 0x3.8p+56f : 0x8.8bdd5p+60f 1 : inexact-ok
+= lgamma downward dbl-64 0x3.8p+56 : 0x8.8bdd41bf4484p+60 1 : inexact-ok
+= lgamma tonearest dbl-64 0x3.8p+56 : 0x8.8bdd41bf44848p+60 1 : inexact-ok
+= lgamma towardzero dbl-64 0x3.8p+56 : 0x8.8bdd41bf4484p+60 1 : inexact-ok
+= lgamma upward dbl-64 0x3.8p+56 : 0x8.8bdd41bf44848p+60 1 : inexact-ok
+= lgamma downward ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma tonearest ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma towardzero ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma upward ldbl-96-intel 0x3.8p+56L : 0x8.8bdd41bf4484606p+60L 1 : inexact-ok
+= lgamma downward ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma tonearest ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma towardzero ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484605p+60L 1 : inexact-ok
+= lgamma upward ldbl-96-m68k 0x3.8p+56L : 0x8.8bdd41bf4484606p+60L 1 : inexact-ok
+= lgamma downward ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d568p+60L 1 : inexact-ok
+= lgamma tonearest ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d57p+60L 1 : inexact-ok
+= lgamma towardzero ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d568p+60L 1 : inexact-ok
+= lgamma upward ldbl-128 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d57p+60L 1 : inexact-ok
+= lgamma downward ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
+= lgamma tonearest ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
+= lgamma towardzero ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d4p+60L 1 : inexact-ok
+= lgamma upward ldbl-128ibm 0x3.8p+56L : 0x8.8bdd41bf44846050819264e2d8p+60L 1 : inexact-ok
 lgamma 0x1p-5
 = lgamma downward flt-32 0x8p-8f : 0x3.72d02cp+0f 1 : inexact-ok
 = lgamma tonearest flt-32 0x8p-8f : 0x3.72d03p+0f 1 : inexact-ok

diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index 0dba9af8d795d59dc9af3f6f8acb74d03cb77867..4743bee438c014003a65ae2199bfd2a53bb59ada 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -219,15 +219,15 @@ __ieee754_lgammaf_r(float x, int *signgamp)
            case 3: z *= (y+(float)2.0);        /* FALLTHRU */
                    r += __ieee754_logf(z); break;
            }
-    /* 8.0 <= x < 2**58 */
-       } else if (ix < 0x5c800000) {
+    /* 8.0 <= x < 2**26 */
+       } else if (ix < 0x4c800000) {
            t = __ieee754_logf(x);
            z = one/x;
            y = z*z;
            w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
            r = (x-half)*(t-one)+w;
        } else
-    /* 2**58 <= x <= inf */
+    /* 2**26 <= x <= inf */
            r =  x*(__ieee754_logf(x)-one);
        if(hx<0) r = nadj - r;
        return r;