This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix lgammaf spurious underflows (bug 18220) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Fri, 15 May 2015 17:21:51 +0000
- Subject: Fix lgammaf spurious underflows (bug 18220) [committed]
- Authentication-results: sourceware.org; auth=none
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. Committed.
(auto-libm-test-out diffs omitted below.)
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.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index c4bfe74..2a88403 100644
--- 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/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index 0dba9af..4743bee 100644
--- 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;
--
Joseph S. Myers
joseph@codesourcery.com