* [Bug math/2542] Incorrect return from float gamma (-0X1.FA471547C2FE5P+1)
[not found] <bug-2542-131@http.sourceware.org/bugzilla/>
2012-02-23 2:05 ` [Bug math/2542] Incorrect return from float gamma (-0X1.FA471547C2FE5P+1) jsm28 at gcc dot gnu.org
2012-02-23 2:08 ` jsm28 at gcc dot gnu.org
@ 2015-09-10 22:29 ` cvs-commit at gcc dot gnu.org
2015-09-10 22:30 ` jsm28 at gcc dot gnu.org
3 siblings, 0 replies; 5+ messages in thread
From: cvs-commit at gcc dot gnu.org @ 2015-09-10 22:29 UTC (permalink / raw)
To: glibc-bugs
https://sourceware.org/bugzilla/show_bug.cgi?id=2542
--- Comment #4 from cvs-commit at gcc dot gnu.org <cvs-commit at gcc dot gnu.org> ---
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU C Library master sources".
The branch, master has been updated
via 050f29c18873ec05ba04a4034bed8cb3f6ae4463 (commit)
from d18c36e6007b03533a38c890c68544daa78d301a (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=050f29c18873ec05ba04a4034bed8cb3f6ae4463
commit 050f29c18873ec05ba04a4034bed8cb3f6ae4463
Author: Joseph Myers <joseph@codesourcery.com>
Date: Thu Sep 10 22:27:58 2015 +0000
Fix lgamma (negative) inaccuracy (bug 2542, bug 2543, bug 2558).
The existing implementations of lgamma functions (except for the ia64
versions) use the reflection formula for negative arguments. This
suffers large inaccuracy from cancellation near zeros of lgamma (near
where the gamma function is +/- 1).
This patch fixes this inaccuracy. For arguments above -2, there are
no zeros and no large cancellation, while for sufficiently large
negative arguments the zeros are so close to integers that even for
integers +/- 1ulp the log(gamma(1-x)) term dominates and cancellation
is not significant. Thus, it is only necessary to take special care
about cancellation for arguments around a limited number of zeros.
Accordingly, this patch uses precomputed tables of relevant zeros,
expressed as the sum of two floating-point values. The log of the
ratio of two sines can be computed accurately using log1p in cases
where log would lose accuracy. The log of the ratio of two gamma(1-x)
values can be computed using Stirling's approximation (the difference
between two values of that approximation to lgamma being computable
without computing the two values and then subtracting), with
appropriate adjustments (which don't reduce accuracy too much) in
cases where 1-x is too small to use Stirling's approximation directly.
In the interval from -3 to -2, using the ratios of sines and of
gamma(1-x) can still produce too much cancellation between those two
parts of the computation (and that interval is also the worst interval
for computing the ratio between gamma(1-x) values, which computation
becomes more accurate, while being less critical for the final result,
for larger 1-x). Because this can result in errors slightly above
those accepted in glibc, this interval is instead dealt with by
polynomial approximations. Separate polynomial approximations to
(|gamma(x)|-1)(x-n)/(x-x0) are used for each interval of length 1/8
from -3 to -2, where n (-3 or -2) is the nearest integer to the
1/8-interval and x0 is the zero of lgamma in the relevant half-integer
interval (-3 to -2.5 or -2.5 to -2).
Together, the two approaches are intended to give sufficient accuracy
for all negative arguments in the problem range. Outside that range,
the previous implementation continues to be used.
Tested for x86_64, x86, mips64 and powerpc. The mips64 and powerpc
testing shows up pre-existing problems for ldbl-128 and ldbl-128ibm
with large negative arguments giving spurious "invalid" exceptions
(exposed by newly added tests for cases this patch doesn't affect the
logic for); I'll address those problems separately.
[BZ #2542]
[BZ #2543]
[BZ #2558]
* sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r): Call
__lgamma_neg for arguments from -28.0 to -2.0.
* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Call
__lgamma_negf for arguments from -15.0 to -2.0.
* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
Call __lgamma_negl for arguments from -48.0 or -50.0 to -2.0.
* sysdeps/ieee754/ldbl-96/e_lgammal_r.c (__ieee754_lgammal_r):
Call __lgamma_negl for arguments from -33.0 to -2.0.
* sysdeps/ieee754/dbl-64/lgamma_neg.c: New file.
* sysdeps/ieee754/dbl-64/lgamma_product.c: Likewise.
* sysdeps/ieee754/flt-32/lgamma_negf.c: Likewise.
* sysdeps/ieee754/flt-32/lgamma_productf.c: Likewise.
* sysdeps/ieee754/ldbl-128/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128/lgamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_negl.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_product.c: Likewise.
* sysdeps/ieee754/ldbl-96/lgamma_productl.c: Likewise.
* sysdeps/generic/math_private.h (__lgamma_negf): New prototype.
(__lgamma_neg): Likewise.
(__lgamma_negl): Likewise.
(__lgamma_product): Likewise.
(__lgamma_productl): Likewise.
* math/Makefile (libm-calls): Add lgamma_neg and lgamma_product.
* math/auto-libm-test-in: Add more tests of lgamma.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
-----------------------------------------------------------------------
Summary of changes:
ChangeLog | 35 +
NEWS | 10 +-
math/Makefile | 2 +-
math/auto-libm-test-in | 437 +-
math/auto-libm-test-out |21344 +++++++++++++++++++-
sysdeps/generic/math_private.h | 16 +
sysdeps/i386/fpu/libm-test-ulps | 96 +-
sysdeps/ieee754/dbl-64/e_lgamma_r.c | 2 +
sysdeps/ieee754/dbl-64/lgamma_neg.c | 399 +
sysdeps/ieee754/dbl-64/lgamma_product.c | 82 +
sysdeps/ieee754/flt-32/e_lgammaf_r.c | 3 +
sysdeps/ieee754/flt-32/lgamma_negf.c | 288 +
.../doasin.c => ieee754/flt-32/lgamma_productf.c} | 0
sysdeps/ieee754/ldbl-128/e_lgammal_r.c | 2 +
sysdeps/ieee754/ldbl-128/lgamma_negl.c | 551 +
sysdeps/ieee754/ldbl-128/lgamma_productl.c | 82 +
sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c | 532 +
sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c | 38 +
sysdeps/ieee754/ldbl-96/e_lgammal_r.c | 2 +
sysdeps/ieee754/ldbl-96/lgamma_negl.c | 418 +
sysdeps/ieee754/ldbl-96/lgamma_product.c | 37 +
sysdeps/ieee754/ldbl-96/lgamma_productl.c | 82 +
sysdeps/x86_64/fpu/libm-test-ulps | 96 +-
23 files changed, 24426 insertions(+), 128 deletions(-)
create mode 100644 sysdeps/ieee754/dbl-64/lgamma_neg.c
create mode 100644 sysdeps/ieee754/dbl-64/lgamma_product.c
create mode 100644 sysdeps/ieee754/flt-32/lgamma_negf.c
copy sysdeps/{i386/fpu/doasin.c => ieee754/flt-32/lgamma_productf.c} (100%)
create mode 100644 sysdeps/ieee754/ldbl-128/lgamma_negl.c
create mode 100644 sysdeps/ieee754/ldbl-128/lgamma_productl.c
create mode 100644 sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c
create mode 100644 sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c
create mode 100644 sysdeps/ieee754/ldbl-96/lgamma_negl.c
create mode 100644 sysdeps/ieee754/ldbl-96/lgamma_product.c
create mode 100644 sysdeps/ieee754/ldbl-96/lgamma_productl.c
--
You are receiving this mail because:
You are on the CC list for the bug.
^ permalink raw reply [flat|nested] 5+ messages in thread