From 4f40e4b30754c678b7f93cfb9a545deb534f7e4e Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Sun, 22 Dec 2013 20:50:16 +0000 Subject: [PATCH] Fix ldbl-128 lgammal for small negative arguments (bug 16337). This patch fixes bug 16337, ldbl-128 lgammal spurious overflows for small negative arguments (the arguments in question are already in the testsuite). The implementation uses the reflection formula to compute lgamma of negative x from lgamma of -x, effectively resulting in a calculation -log(x^2) + log(-x); cancellation isn't problematic in this case (bugs for problematic cancellation in lgamma are 2542, 2543, 2558), but the x^2 calculation can underflow (in which case there is spurious logic to return an overflowing value - lgamma can only ever correctly overflow for large positive arguments, though tgamma can overflow for small arguments of either sign as well as large positive arguments). The fix is simply to calculate the result directly with logl when the argument is a small enough negative number. Tested mips64. * sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r): Calculate results for small negative arguments directly rather than using reflection formula with special underflow handling. --- ChangeLog | 5 +++++ NEWS | 2 +- sysdeps/ieee754/ldbl-128/e_lgammal_r.c | 4 ++-- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/ChangeLog b/ChangeLog index 3e806a20f8..bfca05b4f8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,10 @@ 2013-12-22 Joseph Myers + [BZ #16337] + * sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r): + Calculate results for small negative arguments directly rather + than using reflection formula with special underflow handling. + * sysdeps/mach/hurd/Implies: Change unix/bsd/bsd4.4 to unix/bsd. * sysdeps/unix/bsd/syscalls.list (chflags): Add entry from sysdeps/unix/bsd/bsd4.4/syscalls.list. diff --git a/NEWS b/NEWS index 75c69ecdf3..a07df8cd8a 100644 --- a/NEWS +++ b/NEWS @@ -22,7 +22,7 @@ Version 2.19 15966, 15985, 15988, 15997, 16032, 16034, 16036, 16037, 16038, 16041, 16055, 16071, 16072, 16074, 16077, 16078, 16103, 16112, 16143, 16144, 16146, 16150, 16151, 16153, 16167, 16172, 16195, 16214, 16245, 16271, - 16274, 16283, 16289, 16293, 16314, 16316, 16330, 16338, 16356. + 16274, 16283, 16289, 16293, 16314, 16316, 16330, 16337, 16338, 16356. * The public headers no longer use __unused nor __block. This change is to support compiling programs that are derived from BSD sources and use diff --git a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c index 2b44afb759..23ab9b9a43 100644 --- a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c +++ b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c @@ -782,6 +782,8 @@ __ieee754_lgammal_r (long double x, int *signgamp) *signgamp = -1; else *signgamp = 1; + if (q < 0x1p-120L) + return -__logl (q); z = q - p; if (z > 0.5L) { @@ -789,8 +791,6 @@ __ieee754_lgammal_r (long double x, int *signgamp) z = p - q; } z = q * __sinl (PIL * z); - if (z == 0.0L) - return (*signgamp * huge * huge); w = __ieee754_lgammal_r (q, &i); z = __logl (PIL / z) - w; return (z);