[Bug math/2542] Incorrect return from float gamma (-0X1.FA471547C2FE5P+1)

["cvs-commit at gcc dot gnu.org" <sourceware-bugzilla at sourceware dot org>]

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.