diff options
author | Joseph Myers <joseph@codesourcery.com> | 2015-05-15 17:21:08 +0000 |
---|---|---|
committer | Joseph Myers <joseph@codesourcery.com> | 2015-05-15 17:21:08 +0000 |
commit | ff069f024ae8cf15d53429e034d67ddcece0f67a (patch) | |
tree | d6ad47bf62f47e532cadb862d427eaa2f610c12c /sysdeps/ieee754/flt-32 | |
parent | b2fb25240813266e4f62f2dd039573f1a042ea2f (diff) | |
download | glibc-ff069f024ae8cf15d53429e034d67ddcece0f67a.tar.gz glibc-ff069f024ae8cf15d53429e034d67ddcece0f67a.tar.xz glibc-ff069f024ae8cf15d53429e034d67ddcece0f67a.zip |
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.
Diffstat (limited to 'sysdeps/ieee754/flt-32')
-rw-r--r-- | sysdeps/ieee754/flt-32/e_lgammaf_r.c | 6 |
1 files changed, 3 insertions, 3 deletions
diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c index 0dba9af8d7..4743bee438 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; |