From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 127408 invoked by alias); 10 Sep 2015 22:29:08 -0000 Mailing-List: contact glibc-bugs-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Post: List-Help: , Sender: glibc-bugs-owner@sourceware.org Received: (qmail 127306 invoked by uid 55); 10 Sep 2015 22:29:04 -0000 From: "cvs-commit at gcc dot gnu.org" To: glibc-bugs@sourceware.org Subject: [Bug math/2558] Incorrect return from double gamma (-0X1.FA471547C2FE5P+1) Date: Thu, 10 Sep 2015 22:29:00 -0000 X-Bugzilla-Reason: CC X-Bugzilla-Type: changed X-Bugzilla-Watch-Reason: None X-Bugzilla-Product: glibc X-Bugzilla-Component: math X-Bugzilla-Version: unspecified X-Bugzilla-Keywords: X-Bugzilla-Severity: normal X-Bugzilla-Who: cvs-commit at gcc dot gnu.org X-Bugzilla-Status: NEW X-Bugzilla-Resolution: X-Bugzilla-Priority: P2 X-Bugzilla-Assigned-To: unassigned at sourceware dot org X-Bugzilla-Target-Milestone: --- X-Bugzilla-Flags: X-Bugzilla-Changed-Fields: Message-ID: In-Reply-To: References: Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: 7bit X-Bugzilla-URL: http://sourceware.org/bugzilla/ Auto-Submitted: auto-generated MIME-Version: 1.0 X-SW-Source: 2015-09/txt/msg00126.txt.bz2 https://sourceware.org/bugzilla/show_bug.cgi?id=2558 --- Comment #4 from 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 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.