This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix lgamma (negative) inaccuracy (bug 2542, bug 2543, bug 2558) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Thu, 10 Sep 2015 22:29:00 +0000
- Subject: Fix lgamma (negative) inaccuracy (bug 2542, bug 2543, bug 2558) [committed]
- Authentication-results: sourceware.org; auth=none
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.
Committed.
(auto-libm-test-out diffs omitted below.)
2015-09-10 Joseph Myers <joseph@codesourcery.com>
[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.
diff --git a/math/Makefile b/math/Makefile
index c98c3c4..48e7e4c 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -62,7 +62,7 @@ libm-calls = e_acos e_acosh e_asin e_atan2 e_atanh e_cosh e_exp e_fmod \
s_casinh s_cacosh s_catanh s_csqrt s_cpow s_cproj s_clog10 \
s_fma s_lrint s_llrint s_lround s_llround e_exp10 w_log2 \
s_isinf_ns s_issignaling $(calls:s_%=m_%) x2y2m1 k_casinh \
- gamma_product k_standard
+ gamma_product k_standard lgamma_neg lgamma_product
dbl-only-routines := branred doasin dosincos halfulp mpa mpatan2 \
mpatan mpexp mplog mpsqrt mptan sincos32 slowexp \
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 49f1c55..46c3e72 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1975,7 +1975,6 @@ lgamma max
lgamma 1
lgamma 3
lgamma 0.5
-lgamma -0.5
lgamma 0.7
lgamma 1.2
lgamma 0x3.8p56
@@ -2028,6 +2027,438 @@ lgamma -0x1p-16445
lgamma 0x1p-16494
lgamma -0x1p-16494
+lgamma -0x1.fa471547c2fe5p+1
+lgamma -0x1.9260dcp+1
+
+lgamma -0xffffffp-1
+lgamma -0x1fffffffffffffp-1
+lgamma -0xffffffffffffffffp-1
+lgamma -0x3ffffffffffffffffffffffffffp-1
+lgamma -0x1ffffffffffffffffffffffffffffp-1
+
+lgamma -0.25
+lgamma -0.5
+lgamma -0.75
+lgamma -1.25
+lgamma -1.5
+lgamma -1.75
+lgamma -0x2.08p0
+lgamma -0x2.1p0
+lgamma -0x2.18p0
+lgamma -0x2.2p0
+lgamma -0x2.28p0
+lgamma -0x2.3p0
+lgamma -0x2.38p0
+lgamma -0x2.4p0
+lgamma -0x2.48p0
+lgamma -0x2.5p0
+lgamma -0x2.58p0
+lgamma -0x2.6p0
+lgamma -0x2.68p0
+lgamma -0x2.7p0
+lgamma -0x2.78p0
+lgamma -0x2.8p0
+lgamma -0x2.88p0
+lgamma -0x2.9p0
+lgamma -0x2.98p0
+lgamma -0x2.ap0
+lgamma -0x2.a8p0
+lgamma -0x2.bp0
+lgamma -0x2.b8p0
+lgamma -0x2.cp0
+lgamma -0x2.c8p0
+lgamma -0x2.dp0
+lgamma -0x2.d8p0
+lgamma -0x2.ep0
+lgamma -0x2.e8p0
+lgamma -0x2.fp0
+lgamma -0x2.f8p0
+lgamma -0x3.08p0
+lgamma -0x3.1p0
+lgamma -0x3.18p0
+lgamma -0x3.2p0
+lgamma -0x3.28p0
+lgamma -0x3.3p0
+lgamma -0x3.38p0
+lgamma -0x3.4p0
+lgamma -0x3.48p0
+lgamma -0x3.5p0
+lgamma -0x3.58p0
+lgamma -0x3.6p0
+lgamma -0x3.68p0
+lgamma -0x3.7p0
+lgamma -0x3.78p0
+lgamma -0x3.8p0
+lgamma -0x3.88p0
+lgamma -0x3.9p0
+lgamma -0x3.98p0
+lgamma -0x3.ap0
+lgamma -0x3.a8p0
+lgamma -0x3.bp0
+lgamma -0x3.b8p0
+lgamma -0x3.cp0
+lgamma -0x3.c8p0
+lgamma -0x3.dp0
+lgamma -0x3.d8p0
+lgamma -0x3.ep0
+lgamma -0x3.e8p0
+lgamma -0x3.fp0
+lgamma -0x3.f8p0
+lgamma -4.25
+lgamma -4.5
+lgamma -4.75
+lgamma -5.25
+lgamma -5.5
+lgamma -5.75
+lgamma -6.25
+lgamma -6.5
+lgamma -6.75
+lgamma -7.25
+lgamma -7.5
+lgamma -7.75
+lgamma -8.25
+lgamma -8.5
+lgamma -8.75
+lgamma -9.25
+lgamma -9.5
+lgamma -9.75
+lgamma -10.25
+lgamma -10.5
+lgamma -10.75
+lgamma -11.25
+lgamma -11.5
+lgamma -11.75
+lgamma -12.25
+lgamma -12.5
+lgamma -12.75
+lgamma -13.25
+lgamma -13.5
+lgamma -13.75
+lgamma -14.25
+lgamma -14.5
+lgamma -14.75
+lgamma -15.25
+lgamma -15.5
+lgamma -15.75
+lgamma -16.25
+lgamma -16.5
+lgamma -16.75
+lgamma -17.25
+lgamma -17.5
+lgamma -17.75
+lgamma -18.25
+lgamma -18.5
+lgamma -18.75
+lgamma -19.25
+lgamma -19.5
+lgamma -19.75
+lgamma -20.25
+lgamma -20.5
+lgamma -20.75
+lgamma -21.25
+lgamma -21.5
+lgamma -21.75
+lgamma -22.25
+lgamma -22.5
+lgamma -22.75
+lgamma -23.25
+lgamma -23.5
+lgamma -23.75
+lgamma -24.25
+lgamma -24.5
+lgamma -24.75
+lgamma -25.25
+lgamma -25.5
+lgamma -25.75
+lgamma -26.25
+lgamma -26.5
+lgamma -26.75
+lgamma -27.25
+lgamma -27.5
+lgamma -27.75
+lgamma -28.25
+lgamma -28.5
+lgamma -28.75
+lgamma -29.25
+lgamma -29.5
+lgamma -29.75
+lgamma -30.25
+lgamma -30.5
+lgamma -30.75
+lgamma -31.25
+lgamma -31.5
+lgamma -31.75
+lgamma -32.25
+lgamma -32.5
+lgamma -32.75
+lgamma -33.25
+lgamma -33.5
+lgamma -33.75
+lgamma -34.25
+lgamma -34.5
+lgamma -34.75
+lgamma -35.25
+lgamma -35.5
+lgamma -35.75
+lgamma -36.25
+lgamma -36.5
+lgamma -36.75
+lgamma -37.25
+lgamma -37.5
+lgamma -37.75
+lgamma -38.25
+lgamma -38.5
+lgamma -38.75
+lgamma -39.25
+lgamma -39.5
+lgamma -39.75
+lgamma -40.25
+lgamma -40.5
+lgamma -40.75
+lgamma -41.25
+lgamma -41.5
+lgamma -41.75
+lgamma -42.25
+lgamma -42.5
+lgamma -42.75
+lgamma -43.25
+lgamma -43.5
+lgamma -43.75
+lgamma -44.25
+lgamma -44.5
+lgamma -44.75
+lgamma -45.25
+lgamma -45.5
+lgamma -45.75
+lgamma -46.25
+lgamma -46.5
+lgamma -46.75
+lgamma -47.25
+lgamma -47.5
+lgamma -47.75
+lgamma -48.25
+lgamma -48.5
+lgamma -48.75
+lgamma -49.25
+lgamma -49.5
+lgamma -49.75
+lgamma -50.25
+lgamma -50.5
+lgamma -50.75
+lgamma -51.25
+lgamma -51.5
+lgamma -51.75
+lgamma -52.25
+lgamma -52.5
+lgamma -52.75
+lgamma -53.25
+lgamma -53.5
+lgamma -53.75
+lgamma -54.25
+lgamma -54.5
+lgamma -54.75
+lgamma -55.25
+lgamma -55.5
+lgamma -55.75
+lgamma -56.25
+lgamma -56.5
+lgamma -56.75
+lgamma -57.25
+lgamma -57.5
+lgamma -57.75
+lgamma -58.25
+lgamma -58.5
+lgamma -58.75
+lgamma -59.25
+lgamma -59.5
+lgamma -59.75
+lgamma -60.25
+lgamma -60.5
+lgamma -60.75
+
+# Integers +/- 1ulp for ldbl-128 (gen-auto-libm-tests will round these
+# to produce integers +/- 1ulp for other formats).
+lgamma -0xf.fffffffffffffffffffffffffff8p-4
+lgamma -0x1.0000000000000000000000000001p+0
+lgamma -0x1.ffffffffffffffffffffffffffffp+0
+lgamma -0x2.0000000000000000000000000002p+0
+lgamma -0x2.fffffffffffffffffffffffffffep+0
+lgamma -0x3.0000000000000000000000000002p+0
+lgamma -0x3.fffffffffffffffffffffffffffep+0
+lgamma -0x4.0000000000000000000000000004p+0
+lgamma -0x4.fffffffffffffffffffffffffffcp+0
+lgamma -0x5.0000000000000000000000000004p+0
+lgamma -0x5.fffffffffffffffffffffffffffcp+0
+lgamma -0x6.0000000000000000000000000004p+0
+lgamma -0x6.fffffffffffffffffffffffffffcp+0
+lgamma -0x7.0000000000000000000000000004p+0
+lgamma -0x7.fffffffffffffffffffffffffffcp+0
+lgamma -0x8.0000000000000000000000000008p+0
+lgamma -0x8.fffffffffffffffffffffffffff8p+0
+lgamma -0x9.0000000000000000000000000008p+0
+lgamma -0x9.fffffffffffffffffffffffffff8p+0
+lgamma -0xa.0000000000000000000000000008p+0
+lgamma -0xa.fffffffffffffffffffffffffff8p+0
+lgamma -0xb.0000000000000000000000000008p+0
+lgamma -0xb.fffffffffffffffffffffffffff8p+0
+lgamma -0xc.0000000000000000000000000008p+0
+lgamma -0xc.fffffffffffffffffffffffffff8p+0
+lgamma -0xd.0000000000000000000000000008p+0
+lgamma -0xd.fffffffffffffffffffffffffff8p+0
+lgamma -0xe.0000000000000000000000000008p+0
+lgamma -0xe.fffffffffffffffffffffffffff8p+0
+lgamma -0xf.0000000000000000000000000008p+0
+lgamma -0xf.fffffffffffffffffffffffffff8p+0
+lgamma -0x1.0000000000000000000000000001p+4
+lgamma -0x1.0fffffffffffffffffffffffffffp+4
+lgamma -0x1.1000000000000000000000000001p+4
+lgamma -0x1.1fffffffffffffffffffffffffffp+4
+lgamma -0x1.2000000000000000000000000001p+4
+lgamma -0x1.2fffffffffffffffffffffffffffp+4
+lgamma -0x1.3000000000000000000000000001p+4
+lgamma -0x1.3fffffffffffffffffffffffffffp+4
+lgamma -0x1.4000000000000000000000000001p+4
+lgamma -0x1.4fffffffffffffffffffffffffffp+4
+lgamma -0x1.5000000000000000000000000001p+4
+lgamma -0x1.5fffffffffffffffffffffffffffp+4
+lgamma -0x1.6000000000000000000000000001p+4
+lgamma -0x1.6fffffffffffffffffffffffffffp+4
+lgamma -0x1.7000000000000000000000000001p+4
+lgamma -0x1.7fffffffffffffffffffffffffffp+4
+lgamma -0x1.8000000000000000000000000001p+4
+lgamma -0x1.8fffffffffffffffffffffffffffp+4
+lgamma -0x1.9000000000000000000000000001p+4
+lgamma -0x1.9fffffffffffffffffffffffffffp+4
+lgamma -0x1.a000000000000000000000000001p+4
+lgamma -0x1.afffffffffffffffffffffffffffp+4
+lgamma -0x1.b000000000000000000000000001p+4
+lgamma -0x1.bfffffffffffffffffffffffffffp+4
+lgamma -0x1.c000000000000000000000000001p+4
+lgamma -0x1.cfffffffffffffffffffffffffffp+4
+lgamma -0x1.d000000000000000000000000001p+4
+lgamma -0x1.dfffffffffffffffffffffffffffp+4
+lgamma -0x1.e000000000000000000000000001p+4
+lgamma -0x1.efffffffffffffffffffffffffffp+4
+lgamma -0x1.f000000000000000000000000001p+4
+lgamma -0x1.ffffffffffffffffffffffffffffp+4
+lgamma -0x2.0000000000000000000000000002p+4
+lgamma -0x2.0ffffffffffffffffffffffffffep+4
+lgamma -0x2.1000000000000000000000000002p+4
+lgamma -0x2.1ffffffffffffffffffffffffffep+4
+lgamma -0x2.2000000000000000000000000002p+4
+lgamma -0x2.2ffffffffffffffffffffffffffep+4
+lgamma -0x2.3000000000000000000000000002p+4
+lgamma -0x2.3ffffffffffffffffffffffffffep+4
+lgamma -0x2.4000000000000000000000000002p+4
+lgamma -0x2.4ffffffffffffffffffffffffffep+4
+lgamma -0x2.5000000000000000000000000002p+4
+lgamma -0x2.5ffffffffffffffffffffffffffep+4
+lgamma -0x2.6000000000000000000000000002p+4
+lgamma -0x2.6ffffffffffffffffffffffffffep+4
+lgamma -0x2.7000000000000000000000000002p+4
+lgamma -0x2.7ffffffffffffffffffffffffffep+4
+lgamma -0x2.8000000000000000000000000002p+4
+lgamma -0x2.8ffffffffffffffffffffffffffep+4
+lgamma -0x2.9000000000000000000000000002p+4
+lgamma -0x2.9ffffffffffffffffffffffffffep+4
+lgamma -0x2.a000000000000000000000000002p+4
+lgamma -0x2.affffffffffffffffffffffffffep+4
+lgamma -0x2.b000000000000000000000000002p+4
+lgamma -0x2.bffffffffffffffffffffffffffep+4
+lgamma -0x2.c000000000000000000000000002p+4
+lgamma -0x2.cffffffffffffffffffffffffffep+4
+lgamma -0x2.d000000000000000000000000002p+4
+lgamma -0x2.dffffffffffffffffffffffffffep+4
+lgamma -0x2.e000000000000000000000000002p+4
+lgamma -0x2.effffffffffffffffffffffffffep+4
+lgamma -0x2.f000000000000000000000000002p+4
+lgamma -0x2.fffffffffffffffffffffffffffep+4
+lgamma -0x3.0000000000000000000000000002p+4
+lgamma -0x3.0ffffffffffffffffffffffffffep+4
+lgamma -0x3.1000000000000000000000000002p+4
+lgamma -0x3.1ffffffffffffffffffffffffffep+4
+lgamma -0x3.2000000000000000000000000002p+4
+lgamma -0x3.2ffffffffffffffffffffffffffep+4
+lgamma -0x3.3000000000000000000000000002p+4
+lgamma -0x3.3ffffffffffffffffffffffffffep+4
+lgamma -0x3.4000000000000000000000000002p+4
+lgamma -0x3.4ffffffffffffffffffffffffffep+4
+lgamma -0x3.5000000000000000000000000002p+4
+lgamma -0x3.5ffffffffffffffffffffffffffep+4
+lgamma -0x3.6000000000000000000000000002p+4
+lgamma -0x3.6ffffffffffffffffffffffffffep+4
+lgamma -0x3.7000000000000000000000000002p+4
+lgamma -0x3.7ffffffffffffffffffffffffffep+4
+lgamma -0x3.8000000000000000000000000002p+4
+lgamma -0x3.8ffffffffffffffffffffffffffep+4
+lgamma -0x3.9000000000000000000000000002p+4
+lgamma -0x3.9ffffffffffffffffffffffffffep+4
+lgamma -0x3.a000000000000000000000000002p+4
+lgamma -0x3.affffffffffffffffffffffffffep+4
+lgamma -0x3.b000000000000000000000000002p+4
+lgamma -0x3.bffffffffffffffffffffffffffep+4
+lgamma -0x3.c000000000000000000000000002p+4
+
+# Zeroes of lgamma, until the point where they just duplicate integers
+# +/- 1ulp.
+lgamma -0x2.74ff92c01f0d82abec9f315f1a0712c334804d9cp+0
+lgamma -0x2.bf6821437b20197995a4b4641eaebf4b00b482ap+0
+lgamma -0x3.24c1b793cb35efb8be699ad3d9ba65454cb7fac8p+0
+lgamma -0x3.f48e2a8f85fca170d4561291236cc320a4887d1cp+0
+lgamma -0x4.0a139e16656030c39f0b0de18112ac17bfd6be9p+0
+lgamma -0x4.fdd5de9bbabf3510d0aa4076988501d7d7812528p+0
+lgamma -0x5.021a95fc2db6432a4c56e595394decc6af0430d8p+0
+lgamma -0x5.ffa4bd647d0357dd4ed62cbd31edf8e3f8e5deb8p+0
+lgamma -0x6.005ac9625f233b607c2d96d16385cb86ac56934p+0
+lgamma -0x6.fff2fddae1bbff3d626b65c23fd21f40300a3ba8p+0
+lgamma -0x7.000cff7b7f87adf4482dcdb98782ab2661ca58bp+0
+lgamma -0x7.fffe5fe05673c3ca9e82b522b0ca9d2e8837cd2p+0
+lgamma -0x8.0001a01459fc9f60cb3cec1cec8576677ca538ep+0
+lgamma -0x8.ffffd1c425e80ffc864e95749259e7e20210e8p+0
+lgamma -0x9.00002e3bb47d86d6d843fedc351deb7ad09ec5fp+0
+lgamma -0x9.fffffb606bdfdcd062ae77a50547c69d2eb6f34p+0
+lgamma -0xa.0000049f93bb9927b45d95e15441e03086db914p+0
+lgamma -0xa.ffffff9466e9f1b36dacd2adbd18d05a4e45806p+0
+lgamma -0xb.0000006b9915315d965a6ffea40e4bea39000ddp+0
+lgamma -0xb.fffffff7089387387de41acc3d3c978bd839c8cp+0
+lgamma -0xc.00000008f76c7731567c0f0250f387920df5676p+0
+lgamma -0xc.ffffffff4f6dcf617f97a5ffc757d548d2890cdp+0
+lgamma -0xd.00000000b092309c06683dd1b903e3700857a16p+0
+lgamma -0xd.fffffffff36345ab9e184a3e09d1176dc48e47fp+0
+lgamma -0xe.000000000c9cba545e94e75ec5718f753e2501ep+0
+lgamma -0xe.ffffffffff28c060c6604ef30371f89d37357cap+0
+lgamma -0xf.0000000000d73f9f399bd0e420f85e9ee31b0b9p+0
+lgamma -0xf.fffffffffff28c060c6621f512e72e4d113626ap+0
+lgamma -0x1.000000000000d73f9f399da1424bf93b91f177dp+4
+lgamma -0x1.0ffffffffffff3569c47e7a93e1c46a08a2e008ap+4
+lgamma -0x1.1000000000000ca963b8185688876ca5a3a64ec2p+4
+lgamma -0x1.1fffffffffffff4bec3ce234132d08b2b726187cp+4
+lgamma -0x1.20000000000000b413c31dcbeca4c3b2ffacbb4ap+4
+lgamma -0x1.2ffffffffffffff685b25cbf5f545ced932e3848p+4
+lgamma -0x1.30000000000000097a4da340a0ab81b7b1f1f002p+4
+lgamma -0x1.3fffffffffffffff86af516ff7f76bd67e720d58p+4
+lgamma -0x1.40000000000000007950ae9008089413ccc8a354p+4
+lgamma -0x1.4ffffffffffffffffa391c4248c2a39cfdd49d4ap+4
+lgamma -0x1.500000000000000005c6e3bdb73d5c62f55ed532p+4
+lgamma -0x1.5fffffffffffffffffbcc71a49201eb5aeb96c74p+4
+lgamma -0x1.6000000000000000004338e5b6dfe14a513fb4dp+4
+lgamma -0x1.6ffffffffffffffffffd13c97d9d38fcc4d08d7p+4
+lgamma -0x1.70000000000000000002ec368262c7033b2f6f32p+4
+lgamma -0x1.7fffffffffffffffffffe0d30fe68d0a88335b4cp+4
+lgamma -0x1.800000000000000000001f2cf01972f577cca4b4p+4
+lgamma -0x1.8ffffffffffffffffffffec0c3322e9a0572b1bcp+4
+lgamma -0x1.90000000000000000000013f3ccdd165fa8d4e44p+4
+lgamma -0x1.9ffffffffffffffffffffff3b8bd01cad8d32e38p+4
+lgamma -0x1.a0000000000000000000000c4742fe35272cd1c8p+4
+lgamma -0x1.afffffffffffffffffffffff8b9538f48cc5737ep+4
+lgamma -0x1.b00000000000000000000000746ac70b733a8c82p+4
+lgamma -0x1.bffffffffffffffffffffffffbd79d7672bde8b2p+4
+lgamma -0x1.c00000000000000000000000042862898d42174ep+4
+lgamma -0x1.cfffffffffffffffffffffffffdb4c0ce9794ea6p+4
+lgamma -0x1.d000000000000000000000000024b3f31686b15ap+4
+lgamma -0x1.dffffffffffffffffffffffffffec6cd3afb82ap+4
+lgamma -0x1.e0000000000000000000000000013932c5047d6p+4
+
lgamma 0x8.8d2d5p+0
lgamma 0x1.6a324ap+52
lgamma 0x9.62f59p+0
@@ -2038,6 +2469,10 @@ lgamma 0xb.01191p+0
lgamma 0xb.26fdap+0
lgamma 0xb.4ad0ap+0
lgamma 0xe.7a678p+20
+lgamma -0x2.dea4ccp-4
+lgamma -0x2.dd306p-4
+lgamma -0x1.bdc8bp+0
+lgamma -0x4.0a82e8p-4
log 1
log e
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index 0ab547d..6aea864 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -382,6 +382,22 @@ extern double __gamma_product (double x, double x_eps, int n, double *eps);
extern long double __gamma_productl (long double x, long double x_eps,
int n, long double *eps);
+/* Compute lgamma of a negative argument X, if it is in a range
+ (depending on the floating-point format) for which expansion around
+ zeros is used, setting *SIGNGAMP accordingly. */
+extern float __lgamma_negf (float x, int *signgamp);
+extern double __lgamma_neg (double x, int *signgamp);
+extern long double __lgamma_negl (long double x, int *signgamp);
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+extern double __lgamma_product (double t, double x, double x_eps, int n);
+extern long double __lgamma_productl (long double t, long double x,
+ long double x_eps, int n);
+
#ifndef math_opt_barrier
# define math_opt_barrier(x) \
({ __typeof (x) __x = (x); __asm ("" : "+m" (__x)); __x; })
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 1cbf0db..2b81704 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1558,36 +1558,36 @@ ildouble: 4
ldouble: 4
Function: "gamma":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 2
-ldouble: 2
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 3
+ldouble: 3
Function: "gamma_downward":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 6
-ldouble: 6
+double: 4
+float: 4
+idouble: 4
+ifloat: 4
+ildouble: 7
+ldouble: 7
Function: "gamma_towardzero":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 6
-ldouble: 6
+double: 4
+float: 4
+idouble: 4
+ifloat: 4
+ildouble: 7
+ldouble: 7
Function: "gamma_upward":
-double: 2
-float: 3
-idouble: 2
-ifloat: 3
-ildouble: 4
-ldouble: 4
+double: 3
+float: 4
+idouble: 3
+ifloat: 4
+ildouble: 5
+ldouble: 5
Function: "hypot":
double: 1
@@ -1710,36 +1710,36 @@ ildouble: 5
ldouble: 5
Function: "lgamma":
-double: 1
-float: 1
-idouble: 1
-ifloat: 1
-ildouble: 2
-ldouble: 2
+double: 3
+float: 3
+idouble: 3
+ifloat: 3
+ildouble: 3
+ldouble: 3
Function: "lgamma_downward":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 6
-ldouble: 6
+double: 4
+float: 4
+idouble: 4
+ifloat: 4
+ildouble: 7
+ldouble: 7
Function: "lgamma_towardzero":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 6
-ldouble: 6
+double: 4
+float: 4
+idouble: 4
+ifloat: 4
+ildouble: 7
+ldouble: 7
Function: "lgamma_upward":
-double: 2
-float: 3
-idouble: 2
-ifloat: 3
-ildouble: 4
-ldouble: 4
+double: 3
+float: 4
+idouble: 3
+ifloat: 4
+ildouble: 5
+ldouble: 5
Function: "log":
double: 1
diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
index fc6f594..ea8a9b4 100644
--- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c
+++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c
@@ -226,6 +226,8 @@ __ieee754_lgamma_r(double x, int *signgamp)
if(__builtin_expect(ix>=0x43300000, 0))
/* |x|>=2**52, must be -integer */
return x/zero;
+ if (x < -2.0 && x > -28.0)
+ return __lgamma_neg (x, signgamp);
t = sin_pi(x);
if(t==zero) return one/fabsf(t); /* -integer */
nadj = __ieee754_log(pi/fabs(t*x));
diff --git a/sysdeps/ieee754/dbl-64/lgamma_neg.c b/sysdeps/ieee754/dbl-64/lgamma_neg.c
new file mode 100644
index 0000000..8f54a0f
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/lgamma_neg.c
@@ -0,0 +1,399 @@
+/* lgamma expanding around zeros.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
+
+static const double lgamma_zeros[][2] =
+ {
+ { -0x2.74ff92c01f0d8p+0, -0x2.abec9f315f1ap-56 },
+ { -0x2.bf6821437b202p+0, 0x6.866a5b4b9be14p-56 },
+ { -0x3.24c1b793cb35ep+0, -0xf.b8be699ad3d98p-56 },
+ { -0x3.f48e2a8f85fcap+0, -0x1.70d4561291237p-56 },
+ { -0x4.0a139e1665604p+0, 0xf.3c60f4f21e7fp-56 },
+ { -0x4.fdd5de9bbabf4p+0, 0xa.ef2f55bf89678p-56 },
+ { -0x5.021a95fc2db64p+0, -0x3.2a4c56e595394p-56 },
+ { -0x5.ffa4bd647d034p+0, -0x1.7dd4ed62cbd32p-52 },
+ { -0x6.005ac9625f234p+0, 0x4.9f83d2692e9c8p-56 },
+ { -0x6.fff2fddae1bcp+0, 0xc.29d949a3dc03p-60 },
+ { -0x7.000cff7b7f87cp+0, 0x1.20bb7d2324678p-52 },
+ { -0x7.fffe5fe05673cp+0, -0x3.ca9e82b522b0cp-56 },
+ { -0x8.0001a01459fc8p+0, -0x1.f60cb3cec1cedp-52 },
+ { -0x8.ffffd1c425e8p+0, -0xf.fc864e9574928p-56 },
+ { -0x9.00002e3bb47d8p+0, -0x6.d6d843fedc35p-56 },
+ { -0x9.fffffb606bep+0, 0x2.32f9d51885afap-52 },
+ { -0xa.0000049f93bb8p+0, -0x1.927b45d95e154p-52 },
+ { -0xa.ffffff9466eap+0, 0xe.4c92532d5243p-56 },
+ { -0xb.0000006b9915p+0, -0x3.15d965a6ffea4p-52 },
+ { -0xb.fffffff708938p+0, -0x7.387de41acc3d4p-56 },
+ { -0xc.00000008f76c8p+0, 0x8.cea983f0fdafp-56 },
+ { -0xc.ffffffff4f6ep+0, 0x3.09e80685a0038p-52 },
+ { -0xd.00000000b092p+0, -0x3.09c06683dd1bap-52 },
+ { -0xd.fffffffff3638p+0, 0x3.a5461e7b5c1f6p-52 },
+ { -0xe.000000000c9c8p+0, -0x3.a545e94e75ec6p-52 },
+ { -0xe.ffffffffff29p+0, 0x3.f9f399fb10cfcp-52 },
+ { -0xf.0000000000d7p+0, -0x3.f9f399bd0e42p-52 },
+ { -0xf.fffffffffff28p+0, -0xc.060c6621f513p-56 },
+ { -0x1.000000000000dp+4, -0x7.3f9f399da1424p-52 },
+ { -0x1.0ffffffffffffp+4, -0x3.569c47e7a93e2p-52 },
+ { -0x1.1000000000001p+4, 0x3.569c47e7a9778p-52 },
+ { -0x1.2p+4, 0xb.413c31dcbecdp-56 },
+ { -0x1.2p+4, -0xb.413c31dcbeca8p-56 },
+ { -0x1.3p+4, 0x9.7a4da340a0ab8p-60 },
+ { -0x1.3p+4, -0x9.7a4da340a0ab8p-60 },
+ { -0x1.4p+4, 0x7.950ae90080894p-64 },
+ { -0x1.4p+4, -0x7.950ae90080894p-64 },
+ { -0x1.5p+4, 0x5.c6e3bdb73d5c8p-68 },
+ { -0x1.5p+4, -0x5.c6e3bdb73d5c8p-68 },
+ { -0x1.6p+4, 0x4.338e5b6dfe14cp-72 },
+ { -0x1.6p+4, -0x4.338e5b6dfe14cp-72 },
+ { -0x1.7p+4, 0x2.ec368262c7034p-76 },
+ { -0x1.7p+4, -0x2.ec368262c7034p-76 },
+ { -0x1.8p+4, 0x1.f2cf01972f578p-80 },
+ { -0x1.8p+4, -0x1.f2cf01972f578p-80 },
+ { -0x1.9p+4, 0x1.3f3ccdd165fa9p-84 },
+ { -0x1.9p+4, -0x1.3f3ccdd165fa9p-84 },
+ { -0x1.ap+4, 0xc.4742fe35272dp-92 },
+ { -0x1.ap+4, -0xc.4742fe35272dp-92 },
+ { -0x1.bp+4, 0x7.46ac70b733a8cp-96 },
+ { -0x1.bp+4, -0x7.46ac70b733a8cp-96 },
+ { -0x1.cp+4, 0x4.2862898d42174p-100 },
+ };
+
+static const double e_hi = 0x2.b7e151628aed2p+0, e_lo = 0xa.6abf7158809dp-56;
+
+/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
+ approximation to lgamma function. */
+
+static const double lgamma_coeff[] =
+ {
+ 0x1.5555555555555p-4,
+ -0xb.60b60b60b60b8p-12,
+ 0x3.4034034034034p-12,
+ -0x2.7027027027028p-12,
+ 0x3.72a3c5631fe46p-12,
+ -0x7.daac36664f1f4p-12,
+ 0x1.a41a41a41a41ap-8,
+ -0x7.90a1b2c3d4e6p-8,
+ 0x2.dfd2c703c0dp-4,
+ -0x1.6476701181f3ap+0,
+ 0xd.672219167003p+0,
+ -0x9.cd9292e6660d8p+4,
+ };
+
+#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
+
+/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
+ the integer end-point of the half-integer interval containing x and
+ x0 is the zero of lgamma in that half-integer interval. Each
+ polynomial is expressed in terms of x-xm, where xm is the midpoint
+ of the interval for which the polynomial applies. */
+
+static const double poly_coeff[] =
+ {
+ /* Interval [-2.125, -2] (polynomial degree 10). */
+ -0x1.0b71c5c54d42fp+0,
+ -0xc.73a1dc05f3758p-4,
+ -0x1.ec84140851911p-4,
+ -0xe.37c9da23847e8p-4,
+ -0x1.03cd87cdc0ac6p-4,
+ -0xe.ae9aedce12eep-4,
+ 0x9.b11a1780cfd48p-8,
+ -0xe.f25fc460bdebp-4,
+ 0x2.6e984c61ca912p-4,
+ -0xf.83fea1c6d35p-4,
+ 0x4.760c8c8909758p-4,
+ /* Interval [-2.25, -2.125] (polynomial degree 11). */
+ -0xf.2930890d7d678p-4,
+ -0xc.a5cfde054eaa8p-4,
+ 0x3.9c9e0fdebd99cp-4,
+ -0x1.02a5ad35619d9p+0,
+ 0x9.6e9b1167c164p-4,
+ -0x1.4d8332eba090ap+0,
+ 0x1.1c0c94b1b2b6p+0,
+ -0x1.c9a70d138c74ep+0,
+ 0x1.d7d9cf1d4c196p+0,
+ -0x2.91fbf4cd6abacp+0,
+ 0x2.f6751f74b8ff8p+0,
+ -0x3.e1bb7b09e3e76p+0,
+ /* Interval [-2.375, -2.25] (polynomial degree 12). */
+ -0xd.7d28d505d618p-4,
+ -0xe.69649a3040958p-4,
+ 0xb.0d74a2827cd6p-4,
+ -0x1.924b09228a86ep+0,
+ 0x1.d49b12bcf6175p+0,
+ -0x3.0898bb530d314p+0,
+ 0x4.207a6be8fda4cp+0,
+ -0x6.39eef56d4e9p+0,
+ 0x8.e2e42acbccec8p+0,
+ -0xd.0d91c1e596a68p+0,
+ 0x1.2e20d7099c585p+4,
+ -0x1.c4eb6691b4ca9p+4,
+ 0x2.96a1a11fd85fep+4,
+ /* Interval [-2.5, -2.375] (polynomial degree 13). */
+ -0xb.74ea1bcfff948p-4,
+ -0x1.2a82bd590c376p+0,
+ 0x1.88020f828b81p+0,
+ -0x3.32279f040d7aep+0,
+ 0x5.57ac8252ce868p+0,
+ -0x9.c2aedd093125p+0,
+ 0x1.12c132716e94cp+4,
+ -0x1.ea94dfa5c0a6dp+4,
+ 0x3.66b61abfe858cp+4,
+ -0x6.0cfceb62a26e4p+4,
+ 0xa.beeba09403bd8p+4,
+ -0x1.3188d9b1b288cp+8,
+ 0x2.37f774dd14c44p+8,
+ -0x3.fdf0a64cd7136p+8,
+ /* Interval [-2.625, -2.5] (polynomial degree 13). */
+ -0x3.d10108c27ebbp-4,
+ 0x1.cd557caff7d2fp+0,
+ 0x3.819b4856d36cep+0,
+ 0x6.8505cbacfc42p+0,
+ 0xb.c1b2e6567a4dp+0,
+ 0x1.50a53a3ce6c73p+4,
+ 0x2.57adffbb1ec0cp+4,
+ 0x4.2b15549cf400cp+4,
+ 0x7.698cfd82b3e18p+4,
+ 0xd.2decde217755p+4,
+ 0x1.7699a624d07b9p+8,
+ 0x2.98ecf617abbfcp+8,
+ 0x4.d5244d44d60b4p+8,
+ 0x8.e962bf7395988p+8,
+ /* Interval [-2.75, -2.625] (polynomial degree 12). */
+ -0x6.b5d252a56e8a8p-4,
+ 0x1.28d60383da3a6p+0,
+ 0x1.db6513ada89bep+0,
+ 0x2.e217118fa8c02p+0,
+ 0x4.450112c651348p+0,
+ 0x6.4af990f589b8cp+0,
+ 0x9.2db5963d7a238p+0,
+ 0xd.62c03647da19p+0,
+ 0x1.379f81f6416afp+4,
+ 0x1.c5618b4fdb96p+4,
+ 0x2.9342d0af2ac4ep+4,
+ 0x3.d9cdf56d2b186p+4,
+ 0x5.ab9f91d5a27a4p+4,
+ /* Interval [-2.875, -2.75] (polynomial degree 11). */
+ -0x8.a41b1e4f36ff8p-4,
+ 0xc.da87d3b69dbe8p-4,
+ 0x1.1474ad5c36709p+0,
+ 0x1.761ecb90c8c5cp+0,
+ 0x1.d279bff588826p+0,
+ 0x2.4e5d003fb36a8p+0,
+ 0x2.d575575566842p+0,
+ 0x3.85152b0d17756p+0,
+ 0x4.5213d921ca13p+0,
+ 0x5.55da7dfcf69c4p+0,
+ 0x6.acef729b9404p+0,
+ 0x8.483cc21dd0668p+0,
+ /* Interval [-3, -2.875] (polynomial degree 11). */
+ -0xa.046d667e468f8p-4,
+ 0x9.70b88dcc006cp-4,
+ 0xa.a8a39421c94dp-4,
+ 0xd.2f4d1363f98ep-4,
+ 0xd.ca9aa19975b7p-4,
+ 0xf.cf09c2f54404p-4,
+ 0x1.04b1365a9adfcp+0,
+ 0x1.22b54ef213798p+0,
+ 0x1.2c52c25206bf5p+0,
+ 0x1.4aa3d798aace4p+0,
+ 0x1.5c3f278b504e3p+0,
+ 0x1.7e08292cc347bp+0,
+ };
+
+static const size_t poly_deg[] =
+ {
+ 10,
+ 11,
+ 12,
+ 13,
+ 13,
+ 12,
+ 11,
+ 11,
+ };
+
+static const size_t poly_end[] =
+ {
+ 10,
+ 22,
+ 35,
+ 49,
+ 63,
+ 76,
+ 88,
+ 100,
+ };
+
+/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
+
+static double
+lg_sinpi (double x)
+{
+ if (x <= 0.25)
+ return __sin (M_PI * x);
+ else
+ return __cos (M_PI * (0.5 - x));
+}
+
+/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
+
+static double
+lg_cospi (double x)
+{
+ if (x <= 0.25)
+ return __cos (M_PI * x);
+ else
+ return __sin (M_PI * (0.5 - x));
+}
+
+/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
+
+static double
+lg_cotpi (double x)
+{
+ return lg_cospi (x) / lg_sinpi (x);
+}
+
+/* Compute lgamma of a negative argument -28 < X < -2, setting
+ *SIGNGAMP accordingly. */
+
+double
+__lgamma_neg (double x, int *signgamp)
+{
+ /* Determine the half-integer region X lies in, handle exact
+ integers and determine the sign of the result. */
+ int i = __floor (-2 * x);
+ if ((i & 1) == 0 && i == -2 * x)
+ return 1.0 / 0.0;
+ double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
+ i -= 4;
+ *signgamp = ((i & 2) == 0 ? -1 : 1);
+
+ SET_RESTORE_ROUND (FE_TONEAREST);
+
+ /* Expand around the zero X0 = X0_HI + X0_LO. */
+ double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
+ double xdiff = x - x0_hi - x0_lo;
+
+ /* For arguments in the range -3 to -2, use polynomial
+ approximations to an adjusted version of the gamma function. */
+ if (i < 2)
+ {
+ int j = __floor (-8 * x) - 16;
+ double xm = (-33 - 2 * j) * 0.0625;
+ double x_adj = x - xm;
+ size_t deg = poly_deg[j];
+ size_t end = poly_end[j];
+ double g = poly_coeff[end];
+ for (size_t j = 1; j <= deg; j++)
+ g = g * x_adj + poly_coeff[end - j];
+ return __log1p (g * xdiff / (x - xn));
+ }
+
+ /* The result we want is log (sinpi (X0) / sinpi (X))
+ + log (gamma (1 - X0) / gamma (1 - X)). */
+ double x_idiff = fabs (xn - x), x0_idiff = fabs (xn - x0_hi - x0_lo);
+ double log_sinpi_ratio;
+ if (x0_idiff < x_idiff * 0.5)
+ /* Use log not log1p to avoid inaccuracy from log1p of arguments
+ close to -1. */
+ log_sinpi_ratio = __ieee754_log (lg_sinpi (x0_idiff)
+ / lg_sinpi (x_idiff));
+ else
+ {
+ /* Use log1p not log to avoid inaccuracy from log of arguments
+ close to 1. X0DIFF2 has positive sign if X0 is further from
+ XN than X is from XN, negative sign otherwise. */
+ double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5;
+ double sx0d2 = lg_sinpi (x0diff2);
+ double cx0d2 = lg_cospi (x0diff2);
+ log_sinpi_ratio = __log1p (2 * sx0d2
+ * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
+ }
+
+ double log_gamma_ratio;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ double y0_tmp = 1 - x0_hi;
+ double y0 = y0_tmp;
+ double y0_eps = -x0_hi + (1 - y0) - x0_lo;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ double y_tmp = 1 - x;
+ double y = y_tmp;
+ double y_eps = -x + (1 - y);
+ /* We now wish to compute LOG_GAMMA_RATIO
+ = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
+ accurately approximates the difference Y0 + Y0_EPS - Y -
+ Y_EPS. Use Stirling's approximation. First, we may need to
+ adjust into the range where Stirling's approximation is
+ sufficiently accurate. */
+ double log_gamma_adj = 0;
+ if (i < 6)
+ {
+ int n_up = (7 - i) / 2;
+ double ny0, ny0_eps, ny, ny_eps;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ double y0_tmp = y0 + n_up;
+ ny0 = y0_tmp;
+ ny0_eps = y0 - (ny0 - n_up) + y0_eps;
+ y0 = ny0;
+ y0_eps = ny0_eps;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ double y_tmp = y + n_up;
+ ny = y_tmp;
+ ny_eps = y - (ny - n_up) + y_eps;
+ y = ny;
+ y_eps = ny_eps;
+ double prodm1 = __lgamma_product (xdiff, y - n_up, y_eps, n_up);
+ log_gamma_adj = -__log1p (prodm1);
+ }
+ double log_gamma_high
+ = (xdiff * __log1p ((y0 - e_hi - e_lo + y0_eps) / e_hi)
+ + (y - 0.5 + y_eps) * __log1p (xdiff / y) + log_gamma_adj);
+ /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
+ double y0r = 1 / y0, yr = 1 / y;
+ double y0r2 = y0r * y0r, yr2 = yr * yr;
+ double rdiff = -xdiff / (y * y0);
+ double bterm[NCOEFF];
+ double dlast = rdiff, elast = rdiff * yr * (yr + y0r);
+ bterm[0] = dlast * lgamma_coeff[0];
+ for (size_t j = 1; j < NCOEFF; j++)
+ {
+ double dnext = dlast * y0r2 + elast;
+ double enext = elast * yr2;
+ bterm[j] = dnext * lgamma_coeff[j];
+ dlast = dnext;
+ elast = enext;
+ }
+ double log_gamma_low = 0;
+ for (size_t j = 0; j < NCOEFF; j++)
+ log_gamma_low += bterm[NCOEFF - 1 - j];
+ log_gamma_ratio = log_gamma_high + log_gamma_low;
+
+ return log_sinpi_ratio + log_gamma_ratio;
+}
diff --git a/sysdeps/ieee754/dbl-64/lgamma_product.c b/sysdeps/ieee754/dbl-64/lgamma_product.c
new file mode 100644
index 0000000..8f877a8
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/lgamma_product.c
@@ -0,0 +1,82 @@
+/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Calculate X * Y exactly and store the result in *HI + *LO. It is
+ given that the values are small enough that no overflow occurs and
+ large enough (or zero) that no underflow occurs. */
+
+static void
+mul_split (double *hi, double *lo, double x, double y)
+{
+#ifdef __FP_FAST_FMA
+ /* Fast built-in fused multiply-add. */
+ *hi = x * y;
+ *lo = __builtin_fma (x, y, -*hi);
+#elif defined FP_FAST_FMA
+ /* Fast library fused multiply-add, compiler before GCC 4.6. */
+ *hi = x * y;
+ *lo = __fma (x, y, -*hi);
+#else
+ /* Apply Dekker's algorithm. */
+ *hi = x * y;
+# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
+ double x1 = x * C;
+ double y1 = y * C;
+# undef C
+ x1 = (x - x1) + x1;
+ y1 = (y - y1) + y1;
+ double x2 = x - x1;
+ double y2 = y - y1;
+ *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
+#endif
+}
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+
+double
+__lgamma_product (double t, double x, double x_eps, int n)
+{
+ double ret = 0, ret_eps = 0;
+ for (int i = 0; i < n; i++)
+ {
+ double xi = x + i;
+ double quot = t / xi;
+ double mhi, mlo;
+ mul_split (&mhi, &mlo, quot, xi);
+ double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi);
+ /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */
+ double rhi, rlo;
+ mul_split (&rhi, &rlo, ret, quot);
+ double rpq = ret + quot;
+ double rpq_eps = (ret - rpq) + quot;
+ double nret = rpq + rhi;
+ double nret_eps = (rpq - nret) + rhi;
+ ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot
+ + quot_lo + quot_lo * (ret + ret_eps));
+ ret = nret;
+ }
+ return ret + ret_eps;
+}
diff --git a/sysdeps/ieee754/flt-32/e_lgammaf_r.c b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
index c0bf415..424c4e7 100644
--- a/sysdeps/ieee754/flt-32/e_lgammaf_r.c
+++ b/sysdeps/ieee754/flt-32/e_lgammaf_r.c
@@ -161,6 +161,9 @@ __ieee754_lgammaf_r(float x, int *signgamp)
if(hx<0) {
if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */
return x/zero;
+ if (ix > 0x40000000 /* X < 2.0f. */
+ && ix < 0x41700000 /* X > -15.0f. */)
+ return __lgamma_negf (x, signgamp);
t = sin_pif(x);
if(t==zero) return one/fabsf(t); /* -integer */
nadj = __ieee754_logf(pi/fabsf(t*x));
diff --git a/sysdeps/ieee754/flt-32/lgamma_negf.c b/sysdeps/ieee754/flt-32/lgamma_negf.c
new file mode 100644
index 0000000..ed96598
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/lgamma_negf.c
@@ -0,0 +1,288 @@
+/* lgammaf expanding around zeros.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
+
+static const float lgamma_zeros[][2] =
+ {
+ { -0x2.74ff94p+0f, 0x1.3fe0f2p-24f },
+ { -0x2.bf682p+0f, -0x1.437b2p-24f },
+ { -0x3.24c1b8p+0f, 0x6.c34cap-28f },
+ { -0x3.f48e2cp+0f, 0x1.707a04p-24f },
+ { -0x4.0a13ap+0f, 0x1.e99aap-24f },
+ { -0x4.fdd5ep+0f, 0x1.64454p-24f },
+ { -0x5.021a98p+0f, 0x2.03d248p-24f },
+ { -0x5.ffa4cp+0f, 0x2.9b82fcp-24f },
+ { -0x6.005ac8p+0f, -0x1.625f24p-24f },
+ { -0x6.fff3p+0f, 0x2.251e44p-24f },
+ { -0x7.000dp+0f, 0x8.48078p-28f },
+ { -0x7.fffe6p+0f, 0x1.fa98c4p-28f },
+ { -0x8.0001ap+0f, -0x1.459fcap-28f },
+ { -0x8.ffffdp+0f, -0x1.c425e8p-24f },
+ { -0x9.00003p+0f, 0x1.c44b82p-24f },
+ { -0xap+0f, 0x4.9f942p-24f },
+ { -0xap+0f, -0x4.9f93b8p-24f },
+ { -0xbp+0f, 0x6.b9916p-28f },
+ { -0xbp+0f, -0x6.b9915p-28f },
+ { -0xcp+0f, 0x8.f76c8p-32f },
+ { -0xcp+0f, -0x8.f76c7p-32f },
+ { -0xdp+0f, 0xb.09231p-36f },
+ { -0xdp+0f, -0xb.09231p-36f },
+ { -0xep+0f, 0xc.9cba5p-40f },
+ { -0xep+0f, -0xc.9cba5p-40f },
+ { -0xfp+0f, 0xd.73f9fp-44f },
+ };
+
+static const float e_hi = 0x2.b7e15p+0f, e_lo = 0x1.628aeep-24f;
+
+/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
+ approximation to lgamma function. */
+
+static const float lgamma_coeff[] =
+ {
+ 0x1.555556p-4f,
+ -0xb.60b61p-12f,
+ 0x3.403404p-12f,
+ };
+
+#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
+
+/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
+ the integer end-point of the half-integer interval containing x and
+ x0 is the zero of lgamma in that half-integer interval. Each
+ polynomial is expressed in terms of x-xm, where xm is the midpoint
+ of the interval for which the polynomial applies. */
+
+static const float poly_coeff[] =
+ {
+ /* Interval [-2.125, -2] (polynomial degree 5). */
+ -0x1.0b71c6p+0f,
+ -0xc.73a1ep-4f,
+ -0x1.ec8462p-4f,
+ -0xe.37b93p-4f,
+ -0x1.02ed36p-4f,
+ -0xe.cbe26p-4f,
+ /* Interval [-2.25, -2.125] (polynomial degree 5). */
+ -0xf.29309p-4f,
+ -0xc.a5cfep-4f,
+ 0x3.9c93fcp-4f,
+ -0x1.02a2fp+0f,
+ 0x9.896bep-4f,
+ -0x1.519704p+0f,
+ /* Interval [-2.375, -2.25] (polynomial degree 5). */
+ -0xd.7d28dp-4f,
+ -0xe.6964cp-4f,
+ 0xb.0d4f1p-4f,
+ -0x1.9240aep+0f,
+ 0x1.dadabap+0f,
+ -0x3.1778c4p+0f,
+ /* Interval [-2.5, -2.375] (polynomial degree 6). */
+ -0xb.74ea2p-4f,
+ -0x1.2a82cp+0f,
+ 0x1.880234p+0f,
+ -0x3.320c4p+0f,
+ 0x5.572a38p+0f,
+ -0x9.f92bap+0f,
+ 0x1.1c347ep+4f,
+ /* Interval [-2.625, -2.5] (polynomial degree 6). */
+ -0x3.d10108p-4f,
+ 0x1.cd5584p+0f,
+ 0x3.819c24p+0f,
+ 0x6.84cbb8p+0f,
+ 0xb.bf269p+0f,
+ 0x1.57fb12p+4f,
+ 0x2.7b9854p+4f,
+ /* Interval [-2.75, -2.625] (polynomial degree 6). */
+ -0x6.b5d25p-4f,
+ 0x1.28d604p+0f,
+ 0x1.db6526p+0f,
+ 0x2.e20b38p+0f,
+ 0x4.44c378p+0f,
+ 0x6.62a08p+0f,
+ 0x9.6db3ap+0f,
+ /* Interval [-2.875, -2.75] (polynomial degree 5). */
+ -0x8.a41b2p-4f,
+ 0xc.da87fp-4f,
+ 0x1.147312p+0f,
+ 0x1.7617dap+0f,
+ 0x1.d6c13p+0f,
+ 0x2.57a358p+0f,
+ /* Interval [-3, -2.875] (polynomial degree 5). */
+ -0xa.046d6p-4f,
+ 0x9.70b89p-4f,
+ 0xa.a89a6p-4f,
+ 0xd.2f2d8p-4f,
+ 0xd.e32b4p-4f,
+ 0xf.fb741p-4f,
+ };
+
+static const size_t poly_deg[] =
+ {
+ 5,
+ 5,
+ 5,
+ 6,
+ 6,
+ 6,
+ 5,
+ 5,
+ };
+
+static const size_t poly_end[] =
+ {
+ 5,
+ 11,
+ 17,
+ 24,
+ 31,
+ 38,
+ 44,
+ 50,
+ };
+
+/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
+
+static float
+lg_sinpi (float x)
+{
+ if (x <= 0.25f)
+ return __sinf ((float) M_PI * x);
+ else
+ return __cosf ((float) M_PI * (0.5f - x));
+}
+
+/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
+
+static float
+lg_cospi (float x)
+{
+ if (x <= 0.25f)
+ return __cosf ((float) M_PI * x);
+ else
+ return __sinf ((float) M_PI * (0.5f - x));
+}
+
+/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
+
+static float
+lg_cotpi (float x)
+{
+ return lg_cospi (x) / lg_sinpi (x);
+}
+
+/* Compute lgamma of a negative argument -15 < X < -2, setting
+ *SIGNGAMP accordingly. */
+
+float
+__lgamma_negf (float x, int *signgamp)
+{
+ /* Determine the half-integer region X lies in, handle exact
+ integers and determine the sign of the result. */
+ int i = __floorf (-2 * x);
+ if ((i & 1) == 0 && i == -2 * x)
+ return 1.0f / 0.0f;
+ float xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
+ i -= 4;
+ *signgamp = ((i & 2) == 0 ? -1 : 1);
+
+ SET_RESTORE_ROUNDF (FE_TONEAREST);
+
+ /* Expand around the zero X0 = X0_HI + X0_LO. */
+ float x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
+ float xdiff = x - x0_hi - x0_lo;
+
+ /* For arguments in the range -3 to -2, use polynomial
+ approximations to an adjusted version of the gamma function. */
+ if (i < 2)
+ {
+ int j = __floorf (-8 * x) - 16;
+ float xm = (-33 - 2 * j) * 0.0625f;
+ float x_adj = x - xm;
+ size_t deg = poly_deg[j];
+ size_t end = poly_end[j];
+ float g = poly_coeff[end];
+ for (size_t j = 1; j <= deg; j++)
+ g = g * x_adj + poly_coeff[end - j];
+ return __log1pf (g * xdiff / (x - xn));
+ }
+
+ /* The result we want is log (sinpi (X0) / sinpi (X))
+ + log (gamma (1 - X0) / gamma (1 - X)). */
+ float x_idiff = fabsf (xn - x), x0_idiff = fabsf (xn - x0_hi - x0_lo);
+ float log_sinpi_ratio;
+ if (x0_idiff < x_idiff * 0.5f)
+ /* Use log not log1p to avoid inaccuracy from log1p of arguments
+ close to -1. */
+ log_sinpi_ratio = __ieee754_logf (lg_sinpi (x0_idiff)
+ / lg_sinpi (x_idiff));
+ else
+ {
+ /* Use log1p not log to avoid inaccuracy from log of arguments
+ close to 1. X0DIFF2 has positive sign if X0 is further from
+ XN than X is from XN, negative sign otherwise. */
+ float x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5f;
+ float sx0d2 = lg_sinpi (x0diff2);
+ float cx0d2 = lg_cospi (x0diff2);
+ log_sinpi_ratio = __log1pf (2 * sx0d2
+ * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
+ }
+
+ float log_gamma_ratio;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ float y0_tmp = 1 - x0_hi;
+ float y0 = y0_tmp;
+ float y0_eps = -x0_hi + (1 - y0) - x0_lo;
+#if FLT_EVAL_METHOD != 0
+ volatile
+#endif
+ float y_tmp = 1 - x;
+ float y = y_tmp;
+ float y_eps = -x + (1 - y);
+ /* We now wish to compute LOG_GAMMA_RATIO
+ = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
+ accurately approximates the difference Y0 + Y0_EPS - Y -
+ Y_EPS. Use Stirling's approximation. */
+ float log_gamma_high
+ = (xdiff * __log1pf ((y0 - e_hi - e_lo + y0_eps) / e_hi)
+ + (y - 0.5f + y_eps) * __log1pf (xdiff / y));
+ /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
+ float y0r = 1 / y0, yr = 1 / y;
+ float y0r2 = y0r * y0r, yr2 = yr * yr;
+ float rdiff = -xdiff / (y * y0);
+ float bterm[NCOEFF];
+ float dlast = rdiff, elast = rdiff * yr * (yr + y0r);
+ bterm[0] = dlast * lgamma_coeff[0];
+ for (size_t j = 1; j < NCOEFF; j++)
+ {
+ float dnext = dlast * y0r2 + elast;
+ float enext = elast * yr2;
+ bterm[j] = dnext * lgamma_coeff[j];
+ dlast = dnext;
+ elast = enext;
+ }
+ float log_gamma_low = 0;
+ for (size_t j = 0; j < NCOEFF; j++)
+ log_gamma_low += bterm[NCOEFF - 1 - j];
+ log_gamma_ratio = log_gamma_high + log_gamma_low;
+
+ return log_sinpi_ratio + log_gamma_ratio;
+}
diff --git a/sysdeps/ieee754/flt-32/lgamma_productf.c b/sysdeps/ieee754/flt-32/lgamma_productf.c
new file mode 100644
index 0000000..1cc8931
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/lgamma_productf.c
@@ -0,0 +1 @@
+/* Not needed. */
diff --git a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
index d8a5e5b..abf0f15 100644
--- a/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
+++ b/sysdeps/ieee754/ldbl-128/e_lgammal_r.c
@@ -781,6 +781,8 @@ __ieee754_lgammal_r (long double x, int *signgamp)
if (x < 0.0L)
{
+ if (x < -2.0L && x > (LDBL_MANT_DIG == 106 ? -48.0L : -50.0L))
+ return __lgamma_negl (x, signgamp);
q = -x;
p = __floorl (q);
if (p == q)
diff --git a/sysdeps/ieee754/ldbl-128/lgamma_negl.c b/sysdeps/ieee754/ldbl-128/lgamma_negl.c
new file mode 100644
index 0000000..d68b936
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/lgamma_negl.c
@@ -0,0 +1,551 @@
+/* lgammal expanding around zeros.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
+
+static const long double lgamma_zeros[][2] =
+ {
+ { -0x2.74ff92c01f0d82abec9f315f1a08p+0L, 0xe.d3ccb7fb2658634a2b9f6b2ba81p-116L },
+ { -0x2.bf6821437b20197995a4b4641eaep+0L, -0xb.f4b00b4829f961e428533e6ad048p-116L },
+ { -0x3.24c1b793cb35efb8be699ad3d9bap+0L, -0x6.5454cb7fac60e3f16d9d7840c2ep-116L },
+ { -0x3.f48e2a8f85fca170d4561291236cp+0L, -0xc.320a4887d1cb4c711828a75d5758p-116L },
+ { -0x4.0a139e16656030c39f0b0de18114p+0L, 0x1.53e84029416e1242006b2b3d1cfp-112L },
+ { -0x4.fdd5de9bbabf3510d0aa40769884p+0L, -0x1.01d7d78125286f78d1e501f14966p-112L },
+ { -0x5.021a95fc2db6432a4c56e595394cp+0L, -0x1.ecc6af0430d4fe5746fa7233356fp-112L },
+ { -0x5.ffa4bd647d0357dd4ed62cbd31ecp+0L, -0x1.f8e3f8e5deba2d67dbd70dd96ce1p-112L },
+ { -0x6.005ac9625f233b607c2d96d16384p+0L, -0x1.cb86ac569340cf1e5f24df7aab7bp-112L },
+ { -0x6.fff2fddae1bbff3d626b65c23fd4p+0L, 0x1.e0bfcff5c457ebcf4d3ad9674167p-112L },
+ { -0x7.000cff7b7f87adf4482dcdb98784p+0L, 0x1.54d99e35a74d6407b80292df199fp-112L },
+ { -0x7.fffe5fe05673c3ca9e82b522b0ccp+0L, 0x1.62d177c832e0eb42c2faffd1b145p-112L },
+ { -0x8.0001a01459fc9f60cb3cec1cec88p+0L, 0x2.8998835ac7277f7bcef67c47f188p-112L },
+ { -0x8.ffffd1c425e80ffc864e95749258p+0L, -0x1.e7e20210e7f81cf781b44e9d2b02p-112L },
+ { -0x9.00002e3bb47d86d6d843fedc352p+0L, 0x2.14852f613a16291751d2ab751f7ep-112L },
+ { -0x9.fffffb606bdfdcd062ae77a50548p+0L, 0x3.962d1490cc2e8f031c7007eaa1ap-116L },
+ { -0xa.0000049f93bb9927b45d95e1544p+0L, -0x1.e03086db9146a9287bd4f2172d5ap-112L },
+ { -0xa.ffffff9466e9f1b36dacd2adbd18p+0L, -0xd.05a4e458062f3f95345a4d9c9b6p-116L },
+ { -0xb.0000006b9915315d965a6ffea41p+0L, 0x1.b415c6fff233e7b7fdc3a094246fp-112L },
+ { -0xb.fffffff7089387387de41acc3d4p+0L, 0x3.687427c6373bd74a10306e10a28ep-112L },
+ { -0xc.00000008f76c7731567c0f0250fp+0L, -0x3.87920df5675833859190eb128ef6p-112L },
+ { -0xc.ffffffff4f6dcf617f97a5ffc758p+0L, 0x2.ab72d76f32eaee2d1a42ed515d3ap-116L },
+ { -0xd.00000000b092309c06683dd1b9p+0L, -0x3.e3700857a15c19ac5a611de9688ap-112L },
+ { -0xd.fffffffff36345ab9e184a3e09dp+0L, -0x1.176dc48e47f62d917973dd44e553p-112L },
+ { -0xe.000000000c9cba545e94e75ec57p+0L, -0x1.8f753e2501e757a17cf2ecbeeb89p-112L },
+ { -0xe.ffffffffff28c060c6604ef3037p+0L, -0x1.f89d37357c9e3dc17c6c6e63becap-112L },
+ { -0xf.0000000000d73f9f399bd0e420f8p+0L, -0x5.e9ee31b0b890744fc0e3fbc01048p-116L },
+ { -0xf.fffffffffff28c060c6621f512e8p+0L, 0xd.1b2eec9d960bd9adc5be5f5fa5p-116L },
+ { -0x1.000000000000d73f9f399da1424cp+4L, 0x6.c46e0e88305d2800f0e414c506a8p-116L },
+ { -0x1.0ffffffffffff3569c47e7a93e1cp+4L, -0x4.6a08a2e008a998ebabb8087efa2cp-112L },
+ { -0x1.1000000000000ca963b818568887p+4L, -0x6.ca5a3a64ec15db0a95caf2c9ffb4p-112L },
+ { -0x1.1fffffffffffff4bec3ce234132dp+4L, -0x8.b2b726187c841cb92cd5221e444p-116L },
+ { -0x1.20000000000000b413c31dcbeca5p+4L, 0x3.c4d005344b6cd0e7231120294abcp-112L },
+ { -0x1.2ffffffffffffff685b25cbf5f54p+4L, -0x5.ced932e38485f7dd296b8fa41448p-112L },
+ { -0x1.30000000000000097a4da340a0acp+4L, 0x7.e484e0e0ffe38d406ebebe112f88p-112L },
+ { -0x1.3fffffffffffffff86af516ff7f7p+4L, -0x6.bd67e720d57854502b7db75e1718p-112L },
+ { -0x1.40000000000000007950ae900809p+4L, 0x6.bec33375cac025d9c073168c5d9p-112L },
+ { -0x1.4ffffffffffffffffa391c4248c3p+4L, 0x5.c63022b62b5484ba346524db607p-112L },
+ { -0x1.500000000000000005c6e3bdb73dp+4L, -0x5.c62f55ed5322b2685c5e9a51e6a8p-112L },
+ { -0x1.5fffffffffffffffffbcc71a492p+4L, -0x1.eb5aeb96c74d7ad25e060528fb5p-112L },
+ { -0x1.6000000000000000004338e5b6ep+4L, 0x1.eb5aec04b2f2eb663e4e3d8a018cp-112L },
+ { -0x1.6ffffffffffffffffffd13c97d9dp+4L, -0x3.8fcc4d08d6fe5aa56ab04307ce7ep-112L },
+ { -0x1.70000000000000000002ec368263p+4L, 0x3.8fcc4d090cee2f5d0b69a99c353cp-112L },
+ { -0x1.7fffffffffffffffffffe0d30fe7p+4L, 0x7.2f577cca4b4c8cb1dc14001ac5ecp-112L },
+ { -0x1.800000000000000000001f2cf019p+4L, -0x7.2f577cca4b3442e35f0040b3b9e8p-112L },
+ { -0x1.8ffffffffffffffffffffec0c332p+4L, -0x2.e9a0572b1bb5b95f346a92d67a6p-112L },
+ { -0x1.90000000000000000000013f3ccep+4L, 0x2.e9a0572b1bb5c371ddb3561705ap-112L },
+ { -0x1.9ffffffffffffffffffffff3b8bdp+4L, -0x1.cad8d32e386fd783e97296d63dcbp-116L },
+ { -0x1.a0000000000000000000000c4743p+4L, 0x1.cad8d32e386fd7c1ab8c1fe34c0ep-116L },
+ { -0x1.afffffffffffffffffffffff8b95p+4L, -0x3.8f48cc5737d5979c39db806c5406p-112L },
+ { -0x1.b00000000000000000000000746bp+4L, 0x3.8f48cc5737d5979c3b3a6bda06f6p-112L },
+ { -0x1.bffffffffffffffffffffffffbd8p+4L, 0x6.2898d42174dcf171470d8c8c6028p-112L },
+ { -0x1.c000000000000000000000000428p+4L, -0x6.2898d42174dcf171470d18ba412cp-112L },
+ { -0x1.cfffffffffffffffffffffffffdbp+4L, -0x4.c0ce9794ea50a839e311320bde94p-112L },
+ { -0x1.d000000000000000000000000025p+4L, 0x4.c0ce9794ea50a839e311322f7cf8p-112L },
+ { -0x1.dfffffffffffffffffffffffffffp+4L, 0x3.932c5047d60e60caded4c298a174p-112L },
+ { -0x1.e000000000000000000000000001p+4L, -0x3.932c5047d60e60caded4c298973ap-112L },
+ { -0x1.fp+4L, 0xa.1a6973c1fade2170f7237d35fe3p-116L },
+ { -0x1.fp+4L, -0xa.1a6973c1fade2170f7237d35fe08p-116L },
+ { -0x2p+4L, 0x5.0d34b9e0fd6f10b87b91be9aff1p-120L },
+ { -0x2p+4L, -0x5.0d34b9e0fd6f10b87b91be9aff0cp-120L },
+ { -0x2.1p+4L, 0x2.73024a9ba1aa36a7059bff52e844p-124L },
+ { -0x2.1p+4L, -0x2.73024a9ba1aa36a7059bff52e844p-124L },
+ { -0x2.2p+4L, 0x1.2710231c0fd7a13f8a2b4af9d6b7p-128L },
+ { -0x2.2p+4L, -0x1.2710231c0fd7a13f8a2b4af9d6b7p-128L },
+ { -0x2.3p+4L, 0x8.6e2ce38b6c8f9419e3fad3f0312p-136L },
+ { -0x2.3p+4L, -0x8.6e2ce38b6c8f9419e3fad3f0312p-136L },
+ { -0x2.4p+4L, 0x3.bf30652185952560d71a254e4eb8p-140L },
+ { -0x2.4p+4L, -0x3.bf30652185952560d71a254e4eb8p-140L },
+ { -0x2.5p+4L, 0x1.9ec8d1c94e85af4c78b15c3d89d3p-144L },
+ { -0x2.5p+4L, -0x1.9ec8d1c94e85af4c78b15c3d89d3p-144L },
+ { -0x2.6p+4L, 0xa.ea565ce061d57489e9b85276274p-152L },
+ { -0x2.6p+4L, -0xa.ea565ce061d57489e9b85276274p-152L },
+ { -0x2.7p+4L, 0x4.7a6512692eb37804111dabad30ecp-156L },
+ { -0x2.7p+4L, -0x4.7a6512692eb37804111dabad30ecp-156L },
+ { -0x2.8p+4L, 0x1.ca8ed42a12ae3001a07244abad2bp-160L },
+ { -0x2.8p+4L, -0x1.ca8ed42a12ae3001a07244abad2bp-160L },
+ { -0x2.9p+4L, 0xb.2f30e1ce812063f12e7e8d8d96e8p-168L },
+ { -0x2.9p+4L, -0xb.2f30e1ce812063f12e7e8d8d96e8p-168L },
+ { -0x2.ap+4L, 0x4.42bd49d4c37a0db136489772e428p-172L },
+ { -0x2.ap+4L, -0x4.42bd49d4c37a0db136489772e428p-172L },
+ { -0x2.bp+4L, 0x1.95db45257e5122dcbae56def372p-176L },
+ { -0x2.bp+4L, -0x1.95db45257e5122dcbae56def372p-176L },
+ { -0x2.cp+4L, 0x9.3958d81ff63527ecf993f3fb6f48p-184L },
+ { -0x2.cp+4L, -0x9.3958d81ff63527ecf993f3fb6f48p-184L },
+ { -0x2.dp+4L, 0x3.47970e4440c8f1c058bd238c9958p-188L },
+ { -0x2.dp+4L, -0x3.47970e4440c8f1c058bd238c9958p-188L },
+ { -0x2.ep+4L, 0x1.240804f65951062ca46e4f25c608p-192L },
+ { -0x2.ep+4L, -0x1.240804f65951062ca46e4f25c608p-192L },
+ { -0x2.fp+4L, 0x6.36a382849fae6de2d15362d8a394p-200L },
+ { -0x2.fp+4L, -0x6.36a382849fae6de2d15362d8a394p-200L },
+ { -0x3p+4L, 0x2.123680d6dfe4cf4b9b1bcb9d8bdcp-204L },
+ { -0x3p+4L, -0x2.123680d6dfe4cf4b9b1bcb9d8bdcp-204L },
+ { -0x3.1p+4L, 0xa.d21786ff5842eca51fea0870919p-212L },
+ { -0x3.1p+4L, -0xa.d21786ff5842eca51fea0870919p-212L },
+ { -0x3.2p+4L, 0x3.766dedc259af040be140a68a6c04p-216L },
+ };
+
+static const long double e_hi = 0x2.b7e151628aed2a6abf7158809cf4p+0L;
+static const long double e_lo = 0xf.3c762e7160f38b4da56a784d9048p-116L;
+
+
+/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
+ approximation to lgamma function. */
+
+static const long double lgamma_coeff[] =
+ {
+ 0x1.5555555555555555555555555555p-4L,
+ -0xb.60b60b60b60b60b60b60b60b60b8p-12L,
+ 0x3.4034034034034034034034034034p-12L,
+ -0x2.7027027027027027027027027028p-12L,
+ 0x3.72a3c5631fe46ae1d4e700dca8f2p-12L,
+ -0x7.daac36664f1f207daac36664f1f4p-12L,
+ 0x1.a41a41a41a41a41a41a41a41a41ap-8L,
+ -0x7.90a1b2c3d4e5f708192a3b4c5d7p-8L,
+ 0x2.dfd2c703c0cfff430edfd2c703cp-4L,
+ -0x1.6476701181f39edbdb9ce625987dp+0L,
+ 0xd.672219167002d3a7a9c886459cp+0L,
+ -0x9.cd9292e6660d55b3f712eb9e07c8p+4L,
+ 0x8.911a740da740da740da740da741p+8L,
+ -0x8.d0cc570e255bf59ff6eec24b49p+12L,
+ 0xa.8d1044d3708d1c219ee4fdc446ap+16L,
+ -0xe.8844d8a169abbc406169abbc406p+20L,
+ 0x1.6d29a0f6433b79890cede62433b8p+28L,
+ -0x2.88a233b3c8cddaba9809357125d8p+32L,
+ 0x5.0dde6f27500939a85c40939a85c4p+36L,
+ -0xb.4005bde03d4642a243581714af68p+40L,
+ 0x1.bc8cd6f8f1f755c78753cdb5d5c9p+48L,
+ -0x4.bbebb143bb94de5a0284fa7ec424p+52L,
+ 0xe.2e1337f5af0bed90b6b0a352d4fp+56L,
+ -0x2.e78250162b62405ad3e4bfe61b38p+64L,
+ 0xa.5f7eef9e71ac7c80326ab4cc8bfp+68L,
+ -0x2.83be0395e550213369924971b21ap+76L,
+ 0xa.8ebfe48da17dd999790760b0cep+80L,
+ };
+
+#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
+
+/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
+ the integer end-point of the half-integer interval containing x and
+ x0 is the zero of lgamma in that half-integer interval. Each
+ polynomial is expressed in terms of x-xm, where xm is the midpoint
+ of the interval for which the polynomial applies. */
+
+static const long double poly_coeff[] =
+ {
+ /* Interval [-2.125, -2] (polynomial degree 23). */
+ -0x1.0b71c5c54d42eb6c17f30b7aa8f5p+0L,
+ -0xc.73a1dc05f34951602554c6d7506p-4L,
+ -0x1.ec841408528b51473e6c425ee5ffp-4L,
+ -0xe.37c9da26fc3c9a3c1844c8c7f1cp-4L,
+ -0x1.03cd87c519305703b021fa33f827p-4L,
+ -0xe.ae9ada65e09aa7f1c75216128f58p-4L,
+ 0x9.b11855a4864b5731cf85736015a8p-8L,
+ -0xe.f28c133e697a95c28607c9701dep-4L,
+ 0x2.6ec14a1c586a72a7cc33ee569d6ap-4L,
+ -0xf.57cab973e14464a262fc24723c38p-4L,
+ 0x4.5b0fc25f16e52997b2886bbae808p-4L,
+ -0xf.f50e59f1a9b56e76e988dac9ccf8p-4L,
+ 0x6.5f5eae15e9a93369e1d85146c6fcp-4L,
+ -0x1.0d2422daac459e33e0994325ed23p+0L,
+ 0x8.82000a0e7401fb1117a0e6606928p-4L,
+ -0x1.1f492f178a3f1b19f58a2ca68e55p+0L,
+ 0xa.cb545f949899a04c160b19389abp-4L,
+ -0x1.36165a1b155ba3db3d1b77caf498p+0L,
+ 0xd.44c5d5576f74302e5cf79e183eep-4L,
+ -0x1.51f22e0cdd33d3d481e326c02f3ep+0L,
+ 0xf.f73a349c08244ac389c007779bfp-4L,
+ -0x1.73317bf626156ba716747c4ca866p+0L,
+ 0x1.379c3c97b9bc71e1c1c4802dd657p+0L,
+ -0x1.a72a351c54f902d483052000f5dfp+0L,
+ /* Interval [-2.25, -2.125] (polynomial degree 24). */
+ -0xf.2930890d7d675a80c36afb0fd5e8p-4L,
+ -0xc.a5cfde054eab5c6770daeca577f8p-4L,
+ 0x3.9c9e0fdebb07cdf89c61d41c9238p-4L,
+ -0x1.02a5ad35605fcf4af65a6dbacb84p+0L,
+ 0x9.6e9b1185bb48be9de1918e00a2e8p-4L,
+ -0x1.4d8332f3cfbfa116fd611e9ce90dp+0L,
+ 0x1.1c0c8cb4d9f4b1d490e1a41fae4dp+0L,
+ -0x1.c9a6f5ae9130cd0299e293a42714p+0L,
+ 0x1.d7e9307fd58a2ea997f29573a112p+0L,
+ -0x2.921cb3473d96178ca2a11d2a8d46p+0L,
+ 0x2.e8d59113b6f3409ff8db226e9988p+0L,
+ -0x3.cbab931625a1ae2b26756817f264p+0L,
+ 0x4.7d9f0f05d5296d18663ca003912p+0L,
+ -0x5.ade9cba12a14ea485667b7135bbp+0L,
+ 0x6.dc983a5da74fb48e767b7fec0a3p+0L,
+ -0x8.8d9ed454ae31d9e138dd8ee0d1a8p+0L,
+ 0xa.6fa099d4e7c202e0c0fd6ed8492p+0L,
+ -0xc.ebc552a8090a0f0115e92d4ebbc8p+0L,
+ 0xf.d695e4772c0d829b53fba9ca5568p+0L,
+ -0x1.38c32ae38e5e9eb79b2a4c5570a9p+4L,
+ 0x1.8035145646cfab49306d0999a51bp+4L,
+ -0x1.d930adbb03dd342a4c2a8c4e1af6p+4L,
+ 0x2.45c2edb1b4943ddb3686cd9c6524p+4L,
+ -0x2.e818ebbfafe2f916fa21abf7756p+4L,
+ 0x3.9804ce51d0fb9a430a711fd7307p+4L,
+ /* Interval [-2.375, -2.25] (polynomial degree 25). */
+ -0xd.7d28d505d6181218a25f31d5e45p-4L,
+ -0xe.69649a3040985140cdf946829fap-4L,
+ 0xb.0d74a2827d053a8d44595012484p-4L,
+ -0x1.924b0922853617cac181afbc08ddp+0L,
+ 0x1.d49b12bccf0a568582e2d3c410f3p+0L,
+ -0x3.0898bb7d8c4093e636279c791244p+0L,
+ 0x4.207a6cac711cb53868e8a5057eep+0L,
+ -0x6.39ee63ea4fb1dcab0c9144bf3ddcp+0L,
+ 0x8.e2e2556a797b649bf3f53bd26718p+0L,
+ -0xd.0e83ac82552ef12af508589e7a8p+0L,
+ 0x1.2e4525e0ce6670563c6484a82b05p+4L,
+ -0x1.b8e350d6a8f2b222fa390a57c23dp+4L,
+ 0x2.805cd69b919087d8a80295892c2cp+4L,
+ -0x3.a42585424a1b7e64c71743ab014p+4L,
+ 0x5.4b4f409f98de49f7bfb03c05f984p+4L,
+ -0x7.b3c5827fbe934bc820d6832fb9fcp+4L,
+ 0xb.33b7b90cc96c425526e0d0866e7p+4L,
+ -0x1.04b77047ac4f59ee3775ca10df0dp+8L,
+ 0x1.7b366f5e94a34f41386eac086313p+8L,
+ -0x2.2797338429385c9849ca6355bfc2p+8L,
+ 0x3.225273cf92a27c9aac1b35511256p+8L,
+ -0x4.8f078aa48afe6cb3a4e89690f898p+8L,
+ 0x6.9f311d7b6654fc1d0b5195141d04p+8L,
+ -0x9.a0c297b6b4621619ca9bacc48ed8p+8L,
+ 0xe.ce1f06b6f90d92138232a76e4cap+8L,
+ -0x1.5b0e6806fa064daf011613e43b17p+12L,
+ /* Interval [-2.5, -2.375] (polynomial degree 27). */
+ -0xb.74ea1bcfff94b2c01afba9daa7d8p-4L,
+ -0x1.2a82bd590c37538cab143308de4dp+0L,
+ 0x1.88020f828b966fec66b8649fd6fcp+0L,
+ -0x3.32279f040eb694970e9db24863dcp+0L,
+ 0x5.57ac82517767e68a721005853864p+0L,
+ -0x9.c2aedcfe22833de43834a0a6cc4p+0L,
+ 0x1.12c132f1f5577f99e1a0ed3538e1p+4L,
+ -0x1.ea94e26628a3de3597f7bb55a948p+4L,
+ 0x3.66b4ac4fa582f58b59f96b2f7c7p+4L,
+ -0x6.0cf746a9cf4cba8c39afcc73fc84p+4L,
+ 0xa.c102ef2c20d75a342197df7fedf8p+4L,
+ -0x1.31ebff06e8f14626782df58db3b6p+8L,
+ 0x2.1fd6f0c0e710994e059b9dbdb1fep+8L,
+ -0x3.c6d76040407f447f8b5074f07706p+8L,
+ 0x6.b6d18e0d8feb4c2ef5af6a40ed18p+8L,
+ -0xb.efaf542c529f91e34217f24ae6a8p+8L,
+ 0x1.53852d873210e7070f5d9eb2296p+12L,
+ -0x2.5b977c0ddc6d540717173ac29fc8p+12L,
+ 0x4.310d452ae05100eff1e02343a724p+12L,
+ -0x7.73a5d8f20c4f986a7dd1912b2968p+12L,
+ 0xd.3f5ea2484f3fca15eab1f4d1a218p+12L,
+ -0x1.78d18aac156d1d93a2ffe7e08d3fp+16L,
+ 0x2.9df49ca75e5b567f5ea3e47106cp+16L,
+ -0x4.a7149af8961a08aa7c3233b5bb94p+16L,
+ 0x8.3db10ffa742c707c25197d989798p+16L,
+ -0xe.a26d6dd023cadd02041a049ec368p+16L,
+ 0x1.c825d90514e7c57c7fa5316f947cp+20L,
+ -0x3.34bb81e5a0952df8ca1abdc6684cp+20L,
+ /* Interval [-2.625, -2.5] (polynomial degree 28). */
+ -0x3.d10108c27ebafad533c20eac32bp-4L,
+ 0x1.cd557caff7d2b2085f41dbec5106p+0L,
+ 0x3.819b4856d399520dad9776ea2cacp+0L,
+ 0x6.8505cbad03dc34c5e42e8b12eb78p+0L,
+ 0xb.c1b2e653a9e38f82b399c94e7f08p+0L,
+ 0x1.50a53a38f148138105124df65419p+4L,
+ 0x2.57ae00cbe5232cbeeed34d89727ap+4L,
+ 0x4.2b156301b8604db85a601544bfp+4L,
+ 0x7.6989ed23ca3ca7579b3462592b5cp+4L,
+ 0xd.2dd2976557939517f831f5552cc8p+4L,
+ 0x1.76e1c3430eb860969bce40cd494p+8L,
+ 0x2.9a77bf5488742466db3a2c7c1ec6p+8L,
+ 0x4.a0d62ed7266e8eb36f725a8ebcep+8L,
+ 0x8.3a6184dd3021067df2f8b91e99c8p+8L,
+ 0xe.a0ade1538245bf55d39d7e436b1p+8L,
+ 0x1.a01359fae8617b5826dd74428e9p+12L,
+ 0x2.e3b0a32caae77251169acaca1ad4p+12L,
+ 0x5.2301257c81589f62b38fb5993ee8p+12L,
+ 0x9.21c9275db253d4e719b73b18cb9p+12L,
+ 0x1.03c104bc96141cda3f3fa4b112bcp+16L,
+ 0x1.cdc8ed65119196a08b0c78f1445p+16L,
+ 0x3.34f31d2eaacf34382cdb0073572ap+16L,
+ 0x5.b37628cadf12bf0000907d0ef294p+16L,
+ 0xa.22d8b332c0b1e6a616f425dfe5ap+16L,
+ 0x1.205b01444804c3ff922cd78b4c42p+20L,
+ 0x1.fe8f0cea9d1e0ff25be2470b4318p+20L,
+ 0x3.8872aebeb368399aee02b39340aep+20L,
+ 0x6.ebd560d351e84e26a4381f5b293cp+20L,
+ 0xc.c3644d094b0dae2fbcbf682cd428p+20L,
+ /* Interval [-2.75, -2.625] (polynomial degree 26). */
+ -0x6.b5d252a56e8a75458a27ed1c2dd4p-4L,
+ 0x1.28d60383da3ac721aed3c5794da9p+0L,
+ 0x1.db6513ada8a66ea77d87d9a8827bp+0L,
+ 0x2.e217118f9d348a27f7506a707e6ep+0L,
+ 0x4.450112c5cbf725a0fb9802396c9p+0L,
+ 0x6.4af99151eae7810a75df2a0303c4p+0L,
+ 0x9.2db598b4a97a7f69aeef32aec758p+0L,
+ 0xd.62bef9c22471f5ee47ea1b9c0b5p+0L,
+ 0x1.379f294e412bd62328326d4222f9p+4L,
+ 0x1.c5827349d8865f1e8825c37c31c6p+4L,
+ 0x2.93a7e7a75b7568cc8cbe8c016c12p+4L,
+ 0x3.bf9bb882afe57edb383d41879d3ap+4L,
+ 0x5.73c737828cee095c43a5566731c8p+4L,
+ 0x7.ee4653493a7f81e0442062b3823cp+4L,
+ 0xb.891c6b83fc8b55bd973b5d962d6p+4L,
+ 0x1.0c775d7de3bf9b246c0208e0207ep+8L,
+ 0x1.867ee43ec4bd4f4fd56abc05110ap+8L,
+ 0x2.37fe9ba6695821e9822d8c8af0a6p+8L,
+ 0x3.3a2c667e37c942f182cd3223a936p+8L,
+ 0x4.b1b500eb59f3f782c7ccec88754p+8L,
+ 0x6.d3efd3b65b3d0d8488d30b79fa4cp+8L,
+ 0x9.ee8224e65bed5ced8b75eaec609p+8L,
+ 0xe.72416e510cca77d53fc615c1f3dp+8L,
+ 0x1.4fb538b0a2dfe567a8904b7e0445p+12L,
+ 0x1.e7f56a9266cf525a5b8cf4cb76cep+12L,
+ 0x2.f0365c983f68c597ee49d099cce8p+12L,
+ 0x4.53aa229e1b9f5b5e59625265951p+12L,
+ /* Interval [-2.875, -2.75] (polynomial degree 24). */
+ -0x8.a41b1e4f36ff88dc820815607d68p-4L,
+ 0xc.da87d3b69dc0f2f9c6f368b8ca1p-4L,
+ 0x1.1474ad5c36158a7bea04fd2f98c6p+0L,
+ 0x1.761ecb90c555df6555b7dba955b6p+0L,
+ 0x1.d279bff9ae291caf6c4b4bcb3202p+0L,
+ 0x2.4e5d00559a6e2b9b5d7fe1f6689cp+0L,
+ 0x2.d57545a75cee8743ae2b17bc8d24p+0L,
+ 0x3.8514eee3aac88b89bec2307021bap+0L,
+ 0x4.5235e3b6e1891ffeb87fed9f8a24p+0L,
+ 0x5.562acdb10eef3c9a773b3e27a864p+0L,
+ 0x6.8ec8965c76efe03c26bff60b1194p+0L,
+ 0x8.15251aca144877af32658399f9b8p+0L,
+ 0x9.f08d56aba174d844138af782c0f8p+0L,
+ 0xc.3dbbeda2679e8a1346ccc3f6da88p+0L,
+ 0xf.0f5bfd5eacc26db308ffa0556fa8p+0L,
+ 0x1.28a6ccd84476fbc713d6bab49ac9p+4L,
+ 0x1.6d0a3ae2a3b1c8ff400641a3a21fp+4L,
+ 0x1.c15701b28637f87acfb6a91d33b5p+4L,
+ 0x2.28fbe0eccf472089b017651ca55ep+4L,
+ 0x2.a8a453004f6e8ffaacd1603bc3dp+4L,
+ 0x3.45ae4d9e1e7cd1a5dba0e4ec7f6cp+4L,
+ 0x4.065fbfacb7fad3e473cb577a61e8p+4L,
+ 0x4.f3d1473020927acac1944734a39p+4L,
+ 0x6.54bb091245815a36fb74e314dd18p+4L,
+ 0x7.d7f445129f7fb6c055e582d3f6ep+4L,
+ /* Interval [-3, -2.875] (polynomial degree 23). */
+ -0xa.046d667e468f3e44dcae1afcc648p-4L,
+ 0x9.70b88dcc006c214d8d996fdf5ccp-4L,
+ 0xa.a8a39421c86d3ff24931a0929fp-4L,
+ 0xd.2f4d1363f324da2b357c8b6ec94p-4L,
+ 0xd.ca9aa1a3a5c00de11bf60499a97p-4L,
+ 0xf.cf09c31eeb52a45dfa7ebe3778dp-4L,
+ 0x1.04b133a39ed8a09691205660468bp+0L,
+ 0x1.22b547a06edda944fcb12fd9b5ecp+0L,
+ 0x1.2c57fce7db86a91df09602d344b3p+0L,
+ 0x1.4aade4894708f84795212fe257eep+0L,
+ 0x1.579c8b7b67ec4afed5b28c8bf787p+0L,
+ 0x1.776820e7fc80ae5284239733078ap+0L,
+ 0x1.883ab28c7301fde4ca6b8ec26ec8p+0L,
+ 0x1.aa2ef6e1ae52eb42c9ee83b206e3p+0L,
+ 0x1.bf4ad50f0a9a9311300cf0c51ee7p+0L,
+ 0x1.e40206e0e96b1da463814dde0d09p+0L,
+ 0x1.fdcbcffef3a21b29719c2bd9feb1p+0L,
+ 0x2.25e2e8948939c4d42cf108fae4bep+0L,
+ 0x2.44ce14d2b59c1c0e6bf2cfa81018p+0L,
+ 0x2.70ee80bbd0387162be4861c43622p+0L,
+ 0x2.954b64d2c2ebf3489b949c74476p+0L,
+ 0x2.c616e133a811c1c9446105208656p+0L,
+ 0x3.05a69dfe1a9ba1079f90fcf26bd4p+0L,
+ 0x3.410d2ad16a0506de29736e6aafdap+0L,
+ };
+
+static const size_t poly_deg[] =
+ {
+ 23,
+ 24,
+ 25,
+ 27,
+ 28,
+ 26,
+ 24,
+ 23,
+ };
+
+static const size_t poly_end[] =
+ {
+ 23,
+ 48,
+ 74,
+ 102,
+ 131,
+ 158,
+ 183,
+ 207,
+ };
+
+/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_sinpi (long double x)
+{
+ if (x <= 0.25L)
+ return __sinl (M_PIl * x);
+ else
+ return __cosl (M_PIl * (0.5L - x));
+}
+
+/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cospi (long double x)
+{
+ if (x <= 0.25L)
+ return __cosl (M_PIl * x);
+ else
+ return __sinl (M_PIl * (0.5L - x));
+}
+
+/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cotpi (long double x)
+{
+ return lg_cospi (x) / lg_sinpi (x);
+}
+
+/* Compute lgamma of a negative argument -50 < X < -2, setting
+ *SIGNGAMP accordingly. */
+
+long double
+__lgamma_negl (long double x, int *signgamp)
+{
+ /* Determine the half-integer region X lies in, handle exact
+ integers and determine the sign of the result. */
+ int i = __floorl (-2 * x);
+ if ((i & 1) == 0 && i == -2 * x)
+ return 1.0L / 0.0L;
+ long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
+ i -= 4;
+ *signgamp = ((i & 2) == 0 ? -1 : 1);
+
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+
+ /* Expand around the zero X0 = X0_HI + X0_LO. */
+ long double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
+ long double xdiff = x - x0_hi - x0_lo;
+
+ /* For arguments in the range -3 to -2, use polynomial
+ approximations to an adjusted version of the gamma function. */
+ if (i < 2)
+ {
+ int j = __floorl (-8 * x) - 16;
+ long double xm = (-33 - 2 * j) * 0.0625L;
+ long double x_adj = x - xm;
+ size_t deg = poly_deg[j];
+ size_t end = poly_end[j];
+ long double g = poly_coeff[end];
+ for (size_t j = 1; j <= deg; j++)
+ g = g * x_adj + poly_coeff[end - j];
+ return __log1pl (g * xdiff / (x - xn));
+ }
+
+ /* The result we want is log (sinpi (X0) / sinpi (X))
+ + log (gamma (1 - X0) / gamma (1 - X)). */
+ long double x_idiff = fabsl (xn - x), x0_idiff = fabsl (xn - x0_hi - x0_lo);
+ long double log_sinpi_ratio;
+ if (x0_idiff < x_idiff * 0.5L)
+ /* Use log not log1p to avoid inaccuracy from log1p of arguments
+ close to -1. */
+ log_sinpi_ratio = __ieee754_logl (lg_sinpi (x0_idiff)
+ / lg_sinpi (x_idiff));
+ else
+ {
+ /* Use log1p not log to avoid inaccuracy from log of arguments
+ close to 1. X0DIFF2 has positive sign if X0 is further from
+ XN than X is from XN, negative sign otherwise. */
+ long double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5L;
+ long double sx0d2 = lg_sinpi (x0diff2);
+ long double cx0d2 = lg_cospi (x0diff2);
+ log_sinpi_ratio = __log1pl (2 * sx0d2
+ * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
+ }
+
+ long double log_gamma_ratio;
+ long double y0 = 1 - x0_hi;
+ long double y0_eps = -x0_hi + (1 - y0) - x0_lo;
+ long double y = 1 - x;
+ long double y_eps = -x + (1 - y);
+ /* We now wish to compute LOG_GAMMA_RATIO
+ = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
+ accurately approximates the difference Y0 + Y0_EPS - Y -
+ Y_EPS. Use Stirling's approximation. First, we may need to
+ adjust into the range where Stirling's approximation is
+ sufficiently accurate. */
+ long double log_gamma_adj = 0;
+ if (i < 20)
+ {
+ int n_up = (21 - i) / 2;
+ long double ny0, ny0_eps, ny, ny_eps;
+ ny0 = y0 + n_up;
+ ny0_eps = y0 - (ny0 - n_up) + y0_eps;
+ y0 = ny0;
+ y0_eps = ny0_eps;
+ ny = y + n_up;
+ ny_eps = y - (ny - n_up) + y_eps;
+ y = ny;
+ y_eps = ny_eps;
+ long double prodm1 = __lgamma_productl (xdiff, y - n_up, y_eps, n_up);
+ log_gamma_adj = -__log1pl (prodm1);
+ }
+ long double log_gamma_high
+ = (xdiff * __log1pl ((y0 - e_hi - e_lo + y0_eps) / e_hi)
+ + (y - 0.5L + y_eps) * __log1pl (xdiff / y) + log_gamma_adj);
+ /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
+ long double y0r = 1 / y0, yr = 1 / y;
+ long double y0r2 = y0r * y0r, yr2 = yr * yr;
+ long double rdiff = -xdiff / (y * y0);
+ long double bterm[NCOEFF];
+ long double dlast = rdiff, elast = rdiff * yr * (yr + y0r);
+ bterm[0] = dlast * lgamma_coeff[0];
+ for (size_t j = 1; j < NCOEFF; j++)
+ {
+ long double dnext = dlast * y0r2 + elast;
+ long double enext = elast * yr2;
+ bterm[j] = dnext * lgamma_coeff[j];
+ dlast = dnext;
+ elast = enext;
+ }
+ long double log_gamma_low = 0;
+ for (size_t j = 0; j < NCOEFF; j++)
+ log_gamma_low += bterm[NCOEFF - 1 - j];
+ log_gamma_ratio = log_gamma_high + log_gamma_low;
+
+ return log_sinpi_ratio + log_gamma_ratio;
+}
diff --git a/sysdeps/ieee754/ldbl-128/lgamma_productl.c b/sysdeps/ieee754/ldbl-128/lgamma_productl.c
new file mode 100644
index 0000000..cf0c778
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/lgamma_productl.c
@@ -0,0 +1,82 @@
+/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Calculate X * Y exactly and store the result in *HI + *LO. It is
+ given that the values are small enough that no overflow occurs and
+ large enough (or zero) that no underflow occurs. */
+
+static void
+mul_split (long double *hi, long double *lo, long double x, long double y)
+{
+#ifdef __FP_FAST_FMAL
+ /* Fast built-in fused multiply-add. */
+ *hi = x * y;
+ *lo = __builtin_fmal (x, y, -*hi);
+#elif defined FP_FAST_FMAL
+ /* Fast library fused multiply-add, compiler before GCC 4.6. */
+ *hi = x * y;
+ *lo = __fmal (x, y, -*hi);
+#else
+ /* Apply Dekker's algorithm. */
+ *hi = x * y;
+# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
+ long double x1 = x * C;
+ long double y1 = y * C;
+# undef C
+ x1 = (x - x1) + x1;
+ y1 = (y - y1) + y1;
+ long double x2 = x - x1;
+ long double y2 = y - y1;
+ *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
+#endif
+}
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+
+long double
+__lgamma_productl (long double t, long double x, long double x_eps, int n)
+{
+ long double ret = 0, ret_eps = 0;
+ for (int i = 0; i < n; i++)
+ {
+ long double xi = x + i;
+ long double quot = t / xi;
+ long double mhi, mlo;
+ mul_split (&mhi, &mlo, quot, xi);
+ long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi);
+ /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */
+ long double rhi, rlo;
+ mul_split (&rhi, &rlo, ret, quot);
+ long double rpq = ret + quot;
+ long double rpq_eps = (ret - rpq) + quot;
+ long double nret = rpq + rhi;
+ long double nret_eps = (rpq - nret) + rhi;
+ ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot
+ + quot_lo + quot_lo * (ret + ret_eps));
+ ret = nret;
+ }
+ return ret + ret_eps;
+}
diff --git a/sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c b/sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c
new file mode 100644
index 0000000..cd076ec
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c
@@ -0,0 +1,532 @@
+/* lgammal expanding around zeros.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
+
+static const long double lgamma_zeros[][2] =
+ {
+ { -0x2.74ff92c01f0d82abec9f315f1ap+0L, -0x7.12c334804d9a79cb5d46094d46p-112L },
+ { -0x2.bf6821437b20197995a4b4641fp+0L, 0x5.140b4ff4b7d6069e1bd7acc196p-108L },
+ { -0x3.24c1b793cb35efb8be699ad3dap+0L, 0x4.59abab3480539f1c0e926287cp-108L },
+ { -0x3.f48e2a8f85fca170d456129123p+0L, -0x6.cc320a4887d1cb4c711828a75ep-108L },
+ { -0x4.0a139e16656030c39f0b0de182p+0L, 0xe.d53e84029416e1242006b2b3dp-108L },
+ { -0x4.fdd5de9bbabf3510d0aa407698p+0L, -0x8.501d7d78125286f78d1e501f14p-108L },
+ { -0x5.021a95fc2db6432a4c56e5953ap+0L, 0xb.2133950fbcf2b01a8b9058dcccp-108L },
+ { -0x5.ffa4bd647d0357dd4ed62cbd32p+0L, 0x1.2071c071a2145d2982428f2269p-108L },
+ { -0x6.005ac9625f233b607c2d96d164p+0L, 0x7.a347953a96cbf30e1a0db20856p-108L },
+ { -0x6.fff2fddae1bbff3d626b65c24p+0L, 0x2.de0bfcff5c457ebcf4d3ad9674p-108L },
+ { -0x7.000cff7b7f87adf4482dcdb988p+0L, 0x7.d54d99e35a74d6407b80292df2p-108L },
+ { -0x7.fffe5fe05673c3ca9e82b522bp+0L, -0xc.a9d2e8837cd1f14bd3d05002e4p-108L },
+ { -0x8.0001a01459fc9f60cb3cec1cecp+0L, -0x8.576677ca538d88084310983b8p-108L },
+ { -0x8.ffffd1c425e80ffc864e957494p+0L, 0x1.a6181dfdef1807e3087e4bb163p-104L },
+ { -0x9.00002e3bb47d86d6d843fedc34p+0L, -0x1.1deb7ad09ec5e9d6e8ae2d548bp-104L },
+ { -0x9.fffffb606bdfdcd062ae77a504p+0L, -0x1.47c69d2eb6f33d170fce38ff818p-104L },
+ { -0xa.0000049f93bb9927b45d95e154p+0L, -0x4.1e03086db9146a9287bd4f2172p-108L },
+ { -0xa.ffffff9466e9f1b36dacd2adbcp+0L, -0x1.18d05a4e458062f3f95345a4dap-104L },
+ { -0xb.0000006b9915315d965a6ffea4p+0L, -0xe.4bea39000dcc1848023c5f6bdcp-112L },
+ { -0xb.fffffff7089387387de41acc3cp+0L, -0x1.3c978bd839c8c428b5efcf91ef8p-104L },
+ { -0xc.00000008f76c7731567c0f025p+0L, -0xf.387920df5675833859190eb128p-108L },
+ { -0xc.ffffffff4f6dcf617f97a5ffc8p+0L, 0xa.82ab72d76f32eaee2d1a42ed5p-108L },
+ { -0xd.00000000b092309c06683dd1b8p+0L, -0x1.03e3700857a15c19ac5a611de98p-104L },
+ { -0xd.fffffffff36345ab9e184a3e08p+0L, -0x1.d1176dc48e47f62d917973dd45p-104L },
+ { -0xe.000000000c9cba545e94e75ec4p+0L, -0x1.718f753e2501e757a17cf2ecbfp-104L },
+ { -0xe.ffffffffff28c060c6604ef304p+0L, 0x8.e0762c8ca8361c23e8393919c4p-108L },
+ { -0xf.0000000000d73f9f399bd0e42p+0L, -0xf.85e9ee31b0b890744fc0e3fbcp-108L },
+ { -0xf.fffffffffff28c060c6621f514p+0L, 0x1.18d1b2eec9d960bd9adc5be5f6p-104L },
+ { -0x1.000000000000d73f9f399da1428p+4L, 0x3.406c46e0e88305d2800f0e414cp-104L },
+ { -0x1.0ffffffffffff3569c47e7a93ep+4L, -0x1.c46a08a2e008a998ebabb8087fp-104L },
+ { -0x1.1000000000000ca963b81856888p+4L, -0x7.6ca5a3a64ec15db0a95caf2cap-108L },
+ { -0x1.1fffffffffffff4bec3ce23413p+4L, -0x2.d08b2b726187c841cb92cd5222p-104L },
+ { -0x1.20000000000000b413c31dcbec8p+4L, -0x2.4c3b2ffacbb4932f18dceedfd7p-104L },
+ { -0x1.2ffffffffffffff685b25cbf5f8p+4L, 0x2.ba3126cd1c7b7a0822d694705cp-104L },
+ { -0x1.30000000000000097a4da340a08p+4L, -0x2.b81b7b1f1f001c72bf914141efp-104L },
+ { -0x1.3fffffffffffffff86af516ff8p+4L, 0x8.9429818df2a87abafd48248a2p-108L },
+ { -0x1.40000000000000007950ae9008p+4L, -0x8.9413ccc8a353fda263f8ce973cp-108L },
+ { -0x1.4ffffffffffffffffa391c4249p+4L, 0x3.d5c63022b62b5484ba346524dbp-104L },
+ { -0x1.500000000000000005c6e3bdb7p+4L, -0x3.d5c62f55ed5322b2685c5e9a52p-104L },
+ { -0x1.5fffffffffffffffffbcc71a49p+4L, -0x2.01eb5aeb96c74d7ad25e060529p-104L },
+ { -0x1.6000000000000000004338e5b7p+4L, 0x2.01eb5aec04b2f2eb663e4e3d8ap-104L },
+ { -0x1.6ffffffffffffffffffd13c97d8p+4L, -0x1.d38fcc4d08d6fe5aa56ab04308p-104L },
+ { -0x1.70000000000000000002ec36828p+4L, 0x1.d38fcc4d090cee2f5d0b69a99cp-104L },
+ { -0x1.7fffffffffffffffffffe0d31p+4L, 0x1.972f577cca4b4c8cb1dc14001bp-104L },
+ { -0x1.800000000000000000001f2cfp+4L, -0x1.972f577cca4b3442e35f0040b38p-104L },
+ { -0x1.8ffffffffffffffffffffec0c3p+4L, -0x3.22e9a0572b1bb5b95f346a92d6p-104L },
+ { -0x1.90000000000000000000013f3dp+4L, 0x3.22e9a0572b1bb5c371ddb35617p-104L },
+ { -0x1.9ffffffffffffffffffffff3b88p+4L, -0x3.d01cad8d32e386fd783e97296dp-104L },
+ { -0x1.a0000000000000000000000c478p+4L, 0x3.d01cad8d32e386fd7c1ab8c1fep-104L },
+ { -0x1.afffffffffffffffffffffff8b8p+4L, -0x1.538f48cc5737d5979c39db806c8p-104L },
+ { -0x1.b00000000000000000000000748p+4L, 0x1.538f48cc5737d5979c3b3a6bdap-104L },
+ { -0x1.bffffffffffffffffffffffffcp+4L, 0x2.862898d42174dcf171470d8c8cp-104L },
+ { -0x1.c0000000000000000000000004p+4L, -0x2.862898d42174dcf171470d18bap-104L },
+ { -0x1.dp+4L, 0x2.4b3f31686b15af57c61ceecdf4p-104L },
+ { -0x1.dp+4L, -0x2.4b3f31686b15af57c61ceecdd1p-104L },
+ { -0x1.ep+4L, 0x1.3932c5047d60e60caded4c298ap-108L },
+ { -0x1.ep+4L, -0x1.3932c5047d60e60caded4c29898p-108L },
+ { -0x1.fp+4L, 0xa.1a6973c1fade2170f7237d36p-116L },
+ { -0x1.fp+4L, -0xa.1a6973c1fade2170f7237d36p-116L },
+ { -0x2p+4L, 0x5.0d34b9e0fd6f10b87b91be9bp-120L },
+ { -0x2p+4L, -0x5.0d34b9e0fd6f10b87b91be9bp-120L },
+ { -0x2.1p+4L, 0x2.73024a9ba1aa36a7059bff52e8p-124L },
+ { -0x2.1p+4L, -0x2.73024a9ba1aa36a7059bff52e8p-124L },
+ { -0x2.2p+4L, 0x1.2710231c0fd7a13f8a2b4af9d68p-128L },
+ { -0x2.2p+4L, -0x1.2710231c0fd7a13f8a2b4af9d68p-128L },
+ { -0x2.3p+4L, 0x8.6e2ce38b6c8f9419e3fad3f03p-136L },
+ { -0x2.3p+4L, -0x8.6e2ce38b6c8f9419e3fad3f03p-136L },
+ { -0x2.4p+4L, 0x3.bf30652185952560d71a254e4fp-140L },
+ { -0x2.4p+4L, -0x3.bf30652185952560d71a254e4fp-140L },
+ { -0x2.5p+4L, 0x1.9ec8d1c94e85af4c78b15c3d8ap-144L },
+ { -0x2.5p+4L, -0x1.9ec8d1c94e85af4c78b15c3d8ap-144L },
+ { -0x2.6p+4L, 0xa.ea565ce061d57489e9b8527628p-152L },
+ { -0x2.6p+4L, -0xa.ea565ce061d57489e9b8527628p-152L },
+ { -0x2.7p+4L, 0x4.7a6512692eb37804111dabad3p-156L },
+ { -0x2.7p+4L, -0x4.7a6512692eb37804111dabad3p-156L },
+ { -0x2.8p+4L, 0x1.ca8ed42a12ae3001a07244abadp-160L },
+ { -0x2.8p+4L, -0x1.ca8ed42a12ae3001a07244abadp-160L },
+ { -0x2.9p+4L, 0xb.2f30e1ce812063f12e7e8d8d98p-168L },
+ { -0x2.9p+4L, -0xb.2f30e1ce812063f12e7e8d8d98p-168L },
+ { -0x2.ap+4L, 0x4.42bd49d4c37a0db136489772e4p-172L },
+ { -0x2.ap+4L, -0x4.42bd49d4c37a0db136489772e4p-172L },
+ { -0x2.bp+4L, 0x1.95db45257e5122dcbae56def37p-176L },
+ { -0x2.bp+4L, -0x1.95db45257e5122dcbae56def37p-176L },
+ { -0x2.cp+4L, 0x9.3958d81ff63527ecf993f3fb7p-184L },
+ { -0x2.cp+4L, -0x9.3958d81ff63527ecf993f3fb7p-184L },
+ { -0x2.dp+4L, 0x3.47970e4440c8f1c058bd238c99p-188L },
+ { -0x2.dp+4L, -0x3.47970e4440c8f1c058bd238c99p-188L },
+ { -0x2.ep+4L, 0x1.240804f65951062ca46e4f25c6p-192L },
+ { -0x2.ep+4L, -0x1.240804f65951062ca46e4f25c6p-192L },
+ { -0x2.fp+4L, 0x6.36a382849fae6de2d15362d8a4p-200L },
+ { -0x2.fp+4L, -0x6.36a382849fae6de2d15362d8a4p-200L },
+ { -0x3p+4L, 0x2.123680d6dfe4cf4b9b1bcb9d8cp-204L },
+ };
+
+static const long double e_hi = 0x2.b7e151628aed2a6abf7158809dp+0L;
+static const long double e_lo = -0xb.0c389d18e9f0c74b25a9587b28p-112L;
+
+/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
+ approximation to lgamma function. */
+
+static const long double lgamma_coeff[] =
+ {
+ 0x1.555555555555555555555555558p-4L,
+ -0xb.60b60b60b60b60b60b60b60b6p-12L,
+ 0x3.4034034034034034034034034p-12L,
+ -0x2.7027027027027027027027027p-12L,
+ 0x3.72a3c5631fe46ae1d4e700dca9p-12L,
+ -0x7.daac36664f1f207daac36664f2p-12L,
+ 0x1.a41a41a41a41a41a41a41a41a4p-8L,
+ -0x7.90a1b2c3d4e5f708192a3b4c5ep-8L,
+ 0x2.dfd2c703c0cfff430edfd2c704p-4L,
+ -0x1.6476701181f39edbdb9ce625988p+0L,
+ 0xd.672219167002d3a7a9c886459cp+0L,
+ -0x9.cd9292e6660d55b3f712eb9e08p+4L,
+ 0x8.911a740da740da740da740da74p+8L,
+ -0x8.d0cc570e255bf59ff6eec24b48p+12L,
+ 0xa.8d1044d3708d1c219ee4fdc448p+16L,
+ -0xe.8844d8a169abbc406169abbc4p+20L,
+ 0x1.6d29a0f6433b79890cede624338p+28L,
+ -0x2.88a233b3c8cddaba9809357126p+32L,
+ 0x5.0dde6f27500939a85c40939a86p+36L,
+ -0xb.4005bde03d4642a243581714bp+40L,
+ 0x1.bc8cd6f8f1f755c78753cdb5d6p+48L,
+ -0x4.bbebb143bb94de5a0284fa7ec4p+52L,
+ 0xe.2e1337f5af0bed90b6b0a352d4p+56L,
+ -0x2.e78250162b62405ad3e4bfe61bp+64L,
+ 0xa.5f7eef9e71ac7c80326ab4cc8cp+68L,
+ -0x2.83be0395e550213369924971b2p+76L,
+ };
+
+#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
+
+/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
+ the integer end-point of the half-integer interval containing x and
+ x0 is the zero of lgamma in that half-integer interval. Each
+ polynomial is expressed in terms of x-xm, where xm is the midpoint
+ of the interval for which the polynomial applies. */
+
+static const long double poly_coeff[] =
+ {
+ /* Interval [-2.125, -2] (polynomial degree 21). */
+ -0x1.0b71c5c54d42eb6c17f30b7aa9p+0L,
+ -0xc.73a1dc05f34951602554c6d76cp-4L,
+ -0x1.ec841408528b51473e6c42f1c58p-4L,
+ -0xe.37c9da26fc3c9a3c1844c04b84p-4L,
+ -0x1.03cd87c519305703b00b046ce4p-4L,
+ -0xe.ae9ada65e09aa7f1c817c91048p-4L,
+ 0x9.b11855a4864b571b6a4f571c88p-8L,
+ -0xe.f28c133e697a95ba2dabb97584p-4L,
+ 0x2.6ec14a1c586a7ddb6c4be90fe1p-4L,
+ -0xf.57cab973e14496f0900851c0d4p-4L,
+ 0x4.5b0fc25f16b0df37175495c70cp-4L,
+ -0xf.f50e59f1a8fb8c402091e3cd3cp-4L,
+ 0x6.5f5eae1681d1e50e575c3d4d36p-4L,
+ -0x1.0d2422dac7ea8a52db6bf0d14fp+0L,
+ 0x8.820008f221eae5a36e15913bacp-4L,
+ -0x1.1f492eec53b9481ea23a7e944ep+0L,
+ 0xa.cb55b4d662945e8cf1f81ee5b4p-4L,
+ -0x1.3616863983e131d7935700ccd48p+0L,
+ 0xd.43c783ebab66074d18709d5cap-4L,
+ -0x1.51d5dbc56bc85976871c6e51f78p+0L,
+ 0x1.06253af656eb6b2ed998387aabp+0L,
+ -0x1.7d910a0aadc63d7a1ef7690dbb8p+0L,
+ /* Interval [-2.25, -2.125] (polynomial degree 22). */
+ -0xf.2930890d7d675a80c36afb0fd4p-4L,
+ -0xc.a5cfde054eab5c6770daeca684p-4L,
+ 0x3.9c9e0fdebb07cdf89c61d434adp-4L,
+ -0x1.02a5ad35605fcf4af65a67fe8a8p+0L,
+ 0x9.6e9b1185bb48be9de18d8bbeb8p-4L,
+ -0x1.4d8332f3cfbfa116fdf648372cp+0L,
+ 0x1.1c0c8cb4d9f4b1d495142b53ebp+0L,
+ -0x1.c9a6f5ae9130ccfb9b7e39136f8p+0L,
+ 0x1.d7e9307fd58a2e85209d0e83eap+0L,
+ -0x2.921cb3473d96462f22c171712fp+0L,
+ 0x2.e8d59113b6f3fc1ed3b556b62cp+0L,
+ -0x3.cbab931624e3b6cf299cea1213p+0L,
+ 0x4.7d9f0f05d2c4cf91e41ea1f048p+0L,
+ -0x5.ade9cba31affa276fe516135eep+0L,
+ 0x6.dc983a62cf6ddc935ae3c5b9ap+0L,
+ -0x8.8d9ed100b2a7813f82cbd83e3cp+0L,
+ 0xa.6fa0926892835a9a29c9b8db8p+0L,
+ -0xc.ebc90aff4ffe319d70bef0d61p+0L,
+ 0xf.d69cf50ab226bacece014c0b44p+0L,
+ -0x1.389964ac7cfef4578eec028e5c8p+4L,
+ 0x1.7ff0d2090164e25901f97cab3bp+4L,
+ -0x1.e9e6d282da6bd004619d073071p+4L,
+ 0x2.5d719ab6ad4be8b5c32b0fba2ap+4L,
+ /* Interval [-2.375, -2.25] (polynomial degree 24). */
+ -0xd.7d28d505d6181218a25f31d5e4p-4L,
+ -0xe.69649a3040985140cdf946827cp-4L,
+ 0xb.0d74a2827d053a8d4459500f88p-4L,
+ -0x1.924b0922853617cac181b097e48p+0L,
+ 0x1.d49b12bccf0a568582e2dbf8ep+0L,
+ -0x3.0898bb7d8c4093e6360d26bbc5p+0L,
+ 0x4.207a6cac711cb538684f74619ep+0L,
+ -0x6.39ee63ea4fb1dcac86ab337e3cp+0L,
+ 0x8.e2e2556a797b64a1b9328a3978p+0L,
+ -0xd.0e83ac82552ee5596df1706ff4p+0L,
+ 0x1.2e4525e0ce666e48fac68ddcdep+4L,
+ -0x1.b8e350d6a8f6597ed2eb3c2eff8p+4L,
+ 0x2.805cd69b9197ee0089dd1b1c46p+4L,
+ -0x3.a42585423e4d00db075f2d687ep+4L,
+ 0x5.4b4f409f874e2a7dcd8aa4a62ap+4L,
+ -0x7.b3c5829962ca1b95535db9cc4ep+4L,
+ 0xb.33b7b928986ec6b219e2e15a98p+4L,
+ -0x1.04b76dec4115106bb16316d9cd8p+8L,
+ 0x1.7b366d8d46f179d5c5302d6534p+8L,
+ -0x2.2799846ddc54813d40da622b99p+8L,
+ 0x3.2253a862c1078a3ccabac65bebp+8L,
+ -0x4.8d92cebc90a4a29816f4952f4ep+8L,
+ 0x6.9ebb8f9d72c66c80c4f4492e7ap+8L,
+ -0xa.2850a483f9ba0e43f5848b5cd8p+8L,
+ 0xe.e1b6bdce83b27944edab8c428p+8L,
+ /* Interval [-2.5, -2.375] (polynomial degree 25). */
+ -0xb.74ea1bcfff94b2c01afba9daa8p-4L,
+ -0x1.2a82bd590c37538cab143308e3p+0L,
+ 0x1.88020f828b966fec66b8648d16p+0L,
+ -0x3.32279f040eb694970e9db0308bp+0L,
+ 0x5.57ac82517767e68a72142041b4p+0L,
+ -0x9.c2aedcfe22833de438786dc658p+0L,
+ 0x1.12c132f1f5577f99dbfb7ecb408p+4L,
+ -0x1.ea94e26628a3de3557dc349db8p+4L,
+ 0x3.66b4ac4fa582f5cbe7e19d10c6p+4L,
+ -0x6.0cf746a9cf4cbcb0004cb01f66p+4L,
+ 0xa.c102ef2c20d5a313cbfd37f5b8p+4L,
+ -0x1.31ebff06e8f08f58d1c35eacfdp+8L,
+ 0x2.1fd6f0c0e788660ba1f1573722p+8L,
+ -0x3.c6d760404305e75356a86a11d6p+8L,
+ 0x6.b6d18e0c31a2ba4d5b5ac78676p+8L,
+ -0xb.efaf5426343e6b41a823ed6c44p+8L,
+ 0x1.53852db2fe01305b9f336d132d8p+12L,
+ -0x2.5b977cb2b568382e71ca93a36bp+12L,
+ 0x4.310d090a6119c7d85a2786a616p+12L,
+ -0x7.73a518387ef1d4d04917dfb25cp+12L,
+ 0xd.3f965798601aabd24bdaa6e68cp+12L,
+ -0x1.78db20b0b166480c93cf0031198p+16L,
+ 0x2.9be0068b65cf13bd1cf71f0eccp+16L,
+ -0x4.a221230466b9cd51d5b811d6b6p+16L,
+ 0x8.f6f8c13e2b52aa3e30a4ce6898p+16L,
+ -0x1.02145337ff16b44fa7c2adf7f28p+20L,
+ /* Interval [-2.625, -2.5] (polynomial degree 26). */
+ -0x3.d10108c27ebafad533c20eac33p-4L,
+ 0x1.cd557caff7d2b2085f41dbec538p+0L,
+ 0x3.819b4856d399520dad9776ebb9p+0L,
+ 0x6.8505cbad03dc34c5e42e89c4b4p+0L,
+ 0xb.c1b2e653a9e38f82b3997134a8p+0L,
+ 0x1.50a53a38f1481381051544750ep+4L,
+ 0x2.57ae00cbe5232cbeef4e94eb2cp+4L,
+ 0x4.2b156301b8604db82856d5767p+4L,
+ 0x7.6989ed23ca3ca751fc9c32eb88p+4L,
+ 0xd.2dd29765579396f3a456772c44p+4L,
+ 0x1.76e1c3430eb8630991d1aa8a248p+8L,
+ 0x2.9a77bf548873743fe65d025f56p+8L,
+ 0x4.a0d62ed7266389753842d7be74p+8L,
+ 0x8.3a6184dd32d31ec73fc6f2d37cp+8L,
+ 0xe.a0ade153a3bf0247db49e11ae8p+8L,
+ 0x1.a01359fa74d4eaf8858bbc35f68p+12L,
+ 0x2.e3b0a32845cbc135bae4a5216cp+12L,
+ 0x5.23012653815fe88456170a7dc6p+12L,
+ 0x9.21c92dcde748ec199bc9c65738p+12L,
+ 0x1.03c0f3621b4c67d2d86e5e813d8p+16L,
+ 0x1.cdc884edcc9f5404f2708551cb8p+16L,
+ 0x3.35025f0b1624d1ffc86688bf03p+16L,
+ 0x5.b3bd9562ebf2409c5ce99929ep+16L,
+ 0xa.1a229b1986d9f89cb80abccfdp+16L,
+ 0x1.1e69136ebd520146d51837f3308p+20L,
+ 0x2.2d2738c72449db2524171b9271p+20L,
+ 0x4.036e80cc6621b836f94f426834p+20L,
+ /* Interval [-2.75, -2.625] (polynomial degree 24). */
+ -0x6.b5d252a56e8a75458a27ed1c2ep-4L,
+ 0x1.28d60383da3ac721aed3c57949p+0L,
+ 0x1.db6513ada8a66ea77d87d9a796p+0L,
+ 0x2.e217118f9d348a27f7506c4b4fp+0L,
+ 0x4.450112c5cbf725a0fb982fc44cp+0L,
+ 0x6.4af99151eae7810a75a5fceac8p+0L,
+ 0x9.2db598b4a97a7f69ab7be31128p+0L,
+ 0xd.62bef9c22471f5f17955733c6p+0L,
+ 0x1.379f294e412bd6255506135f4a8p+4L,
+ 0x1.c5827349d8865d858d4f85f3c38p+4L,
+ 0x2.93a7e7a75b755bbea1785a1349p+4L,
+ 0x3.bf9bb882afed66a08b22ed7a45p+4L,
+ 0x5.73c737828d2044aca95fdef33ep+4L,
+ 0x7.ee46534920f1c81574db260f0ep+4L,
+ 0xb.891c6b837b513eaf1592fe78ccp+4L,
+ 0x1.0c775d815bf741526a3dd66ded8p+8L,
+ 0x1.867ee44cf11f26455a8924a56bp+8L,
+ 0x2.37fe968baa1018e55cae680f1dp+8L,
+ 0x3.3a2c557f686679eb5d8e960fd1p+8L,
+ 0x4.b1ba0539d4d80cc9174738b992p+8L,
+ 0x6.d3fd80155b6d2211956cb6bc5ap+8L,
+ 0x9.eb5a96b0ee3d9ca523f5fbc1fp+8L,
+ 0xe.6b37429c1acc7dc19ef312dda4p+8L,
+ 0x1.621132d6aa138b203a28e4792fp+12L,
+ 0x2.09610219270e2ce11a985d4d36p+12L,
+ /* Interval [-2.875, -2.75] (polynomial degree 23). */
+ -0x8.a41b1e4f36ff88dc820815607cp-4L,
+ 0xc.da87d3b69dc0f2f9c6f368b8c8p-4L,
+ 0x1.1474ad5c36158a7bea04fd30b28p+0L,
+ 0x1.761ecb90c555df6555b7dbb9ce8p+0L,
+ 0x1.d279bff9ae291caf6c4b17497f8p+0L,
+ 0x2.4e5d00559a6e2b9b5d7e35b575p+0L,
+ 0x2.d57545a75cee8743b1ff6e22b8p+0L,
+ 0x3.8514eee3aac88b89d2d4ddef4ep+0L,
+ 0x4.5235e3b6e1891fd9c975383318p+0L,
+ 0x5.562acdb10eef3c14a780490e3cp+0L,
+ 0x6.8ec8965c76f0b261bc41b5e532p+0L,
+ 0x8.15251aca144a98a1e1c0981388p+0L,
+ 0x9.f08d56ab9e7eee9515a457214cp+0L,
+ 0xc.3dbbeda2620d5be4fe8621ce6p+0L,
+ 0xf.0f5bfd65b3feb6d745a2cdbf9cp+0L,
+ 0x1.28a6ccd8dd27fb90fcaa31d37dp+4L,
+ 0x1.6d0a3a3091c3d64cfd1a3c5769p+4L,
+ 0x1.c1570107e02d5ab0b8bea6d6c98p+4L,
+ 0x2.28fc9b295b583fa469de7acceap+4L,
+ 0x2.a8a4cac0217026bbdbce34f4adp+4L,
+ 0x3.4532c98bce75262ac0ede53edep+4L,
+ 0x4.062fd9ba18e00e55c25a4f0688p+4L,
+ 0x5.22e00e6d9846a3451fad5587f8p+4L,
+ 0x6.5d0f7ce92a0bf928d4a30e92c6p+4L,
+ /* Interval [-3, -2.875] (polynomial degree 22). */
+ -0xa.046d667e468f3e44dcae1afcc8p-4L,
+ 0x9.70b88dcc006c214d8d996fdf7p-4L,
+ 0xa.a8a39421c86d3ff24931a093c4p-4L,
+ 0xd.2f4d1363f324da2b357c850124p-4L,
+ 0xd.ca9aa1a3a5c00de11bf5d7047p-4L,
+ 0xf.cf09c31eeb52a45dfb25e50ebcp-4L,
+ 0x1.04b133a39ed8a096914cc78812p+0L,
+ 0x1.22b547a06edda9447f516a2ee7p+0L,
+ 0x1.2c57fce7db86a91c8d0f12077b8p+0L,
+ 0x1.4aade4894708fb8b78365e9bf88p+0L,
+ 0x1.579c8b7b67ec5179ecc4e9c7dp+0L,
+ 0x1.776820e7fc7361c50e7ef40a88p+0L,
+ 0x1.883ab28c72ef238ada6c480ab18p+0L,
+ 0x1.aa2ef6e1d11b9fcea06a1dcab1p+0L,
+ 0x1.bf4ad50f2dd2aeb02395ea08648p+0L,
+ 0x1.e40206a5477615838e02279dfc8p+0L,
+ 0x1.fdcbcfd4b0777fb173b85d5b398p+0L,
+ 0x2.25e32b3b3c89e833029169a17bp+0L,
+ 0x2.44ce344ff0bda6570fe3d0a76dp+0L,
+ 0x2.70bfba6fa079faf2dbf31d2216p+0L,
+ 0x2.953e22a97725cc179ad21024fap+0L,
+ 0x2.d8ccc51524659a499eee0f267p+0L,
+ 0x3.080fbb09c14936c2171c8a51bcp+0L,
+ };
+
+static const size_t poly_deg[] =
+ {
+ 21,
+ 22,
+ 24,
+ 25,
+ 26,
+ 24,
+ 23,
+ 22,
+ };
+
+static const size_t poly_end[] =
+ {
+ 21,
+ 44,
+ 69,
+ 95,
+ 122,
+ 147,
+ 171,
+ 194,
+ };
+
+/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_sinpi (long double x)
+{
+ if (x <= 0.25L)
+ return __sinl (M_PIl * x);
+ else
+ return __cosl (M_PIl * (0.5L - x));
+}
+
+/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cospi (long double x)
+{
+ if (x <= 0.25L)
+ return __cosl (M_PIl * x);
+ else
+ return __sinl (M_PIl * (0.5L - x));
+}
+
+/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cotpi (long double x)
+{
+ return lg_cospi (x) / lg_sinpi (x);
+}
+
+/* Compute lgamma of a negative argument -48 < X < -2, setting
+ *SIGNGAMP accordingly. */
+
+long double
+__lgamma_negl (long double x, int *signgamp)
+{
+ /* Determine the half-integer region X lies in, handle exact
+ integers and determine the sign of the result. */
+ int i = __floorl (-2 * x);
+ if ((i & 1) == 0 && i == -2 * x)
+ return 1.0L / 0.0L;
+ long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
+ i -= 4;
+ *signgamp = ((i & 2) == 0 ? -1 : 1);
+
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+
+ /* Expand around the zero X0 = X0_HI + X0_LO. */
+ long double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
+ long double xdiff = x - x0_hi - x0_lo;
+
+ /* For arguments in the range -3 to -2, use polynomial
+ approximations to an adjusted version of the gamma function. */
+ if (i < 2)
+ {
+ int j = __floorl (-8 * x) - 16;
+ long double xm = (-33 - 2 * j) * 0.0625L;
+ long double x_adj = x - xm;
+ size_t deg = poly_deg[j];
+ size_t end = poly_end[j];
+ long double g = poly_coeff[end];
+ for (size_t j = 1; j <= deg; j++)
+ g = g * x_adj + poly_coeff[end - j];
+ return __log1pl (g * xdiff / (x - xn));
+ }
+
+ /* The result we want is log (sinpi (X0) / sinpi (X))
+ + log (gamma (1 - X0) / gamma (1 - X)). */
+ long double x_idiff = fabsl (xn - x), x0_idiff = fabsl (xn - x0_hi - x0_lo);
+ long double log_sinpi_ratio;
+ if (x0_idiff < x_idiff * 0.5L)
+ /* Use log not log1p to avoid inaccuracy from log1p of arguments
+ close to -1. */
+ log_sinpi_ratio = __ieee754_logl (lg_sinpi (x0_idiff)
+ / lg_sinpi (x_idiff));
+ else
+ {
+ /* Use log1p not log to avoid inaccuracy from log of arguments
+ close to 1. X0DIFF2 has positive sign if X0 is further from
+ XN than X is from XN, negative sign otherwise. */
+ long double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5L;
+ long double sx0d2 = lg_sinpi (x0diff2);
+ long double cx0d2 = lg_cospi (x0diff2);
+ log_sinpi_ratio = __log1pl (2 * sx0d2
+ * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
+ }
+
+ long double log_gamma_ratio;
+ long double y0 = 1 - x0_hi;
+ long double y0_eps = -x0_hi + (1 - y0) - x0_lo;
+ long double y = 1 - x;
+ long double y_eps = -x + (1 - y);
+ /* We now wish to compute LOG_GAMMA_RATIO
+ = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
+ accurately approximates the difference Y0 + Y0_EPS - Y -
+ Y_EPS. Use Stirling's approximation. First, we may need to
+ adjust into the range where Stirling's approximation is
+ sufficiently accurate. */
+ long double log_gamma_adj = 0;
+ if (i < 18)
+ {
+ int n_up = (19 - i) / 2;
+ long double ny0, ny0_eps, ny, ny_eps;
+ ny0 = y0 + n_up;
+ ny0_eps = y0 - (ny0 - n_up) + y0_eps;
+ y0 = ny0;
+ y0_eps = ny0_eps;
+ ny = y + n_up;
+ ny_eps = y - (ny - n_up) + y_eps;
+ y = ny;
+ y_eps = ny_eps;
+ long double prodm1 = __lgamma_productl (xdiff, y - n_up, y_eps, n_up);
+ log_gamma_adj = -__log1pl (prodm1);
+ }
+ long double log_gamma_high
+ = (xdiff * __log1pl ((y0 - e_hi - e_lo + y0_eps) / e_hi)
+ + (y - 0.5L + y_eps) * __log1pl (xdiff / y) + log_gamma_adj);
+ /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
+ long double y0r = 1 / y0, yr = 1 / y;
+ long double y0r2 = y0r * y0r, yr2 = yr * yr;
+ long double rdiff = -xdiff / (y * y0);
+ long double bterm[NCOEFF];
+ long double dlast = rdiff, elast = rdiff * yr * (yr + y0r);
+ bterm[0] = dlast * lgamma_coeff[0];
+ for (size_t j = 1; j < NCOEFF; j++)
+ {
+ long double dnext = dlast * y0r2 + elast;
+ long double enext = elast * yr2;
+ bterm[j] = dnext * lgamma_coeff[j];
+ dlast = dnext;
+ elast = enext;
+ }
+ long double log_gamma_low = 0;
+ for (size_t j = 0; j < NCOEFF; j++)
+ log_gamma_low += bterm[NCOEFF - 1 - j];
+ log_gamma_ratio = log_gamma_high + log_gamma_low;
+
+ return log_sinpi_ratio + log_gamma_ratio;
+}
diff --git a/sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c b/sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c
new file mode 100644
index 0000000..9d9c43a
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c
@@ -0,0 +1,38 @@
+/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+
+long double
+__lgamma_productl (long double t, long double x, long double x_eps, int n)
+{
+ long double x_full = x + x_eps;
+ long double ret = 0;
+ for (int i = 0; i < n; i++)
+ /* FIXME: no extra precision used. */
+ ret += (t / (x_full + i)) * (1 + ret);
+ return ret;
+}
diff --git a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
index 0cc35f9..a80002b 100644
--- a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
+++ b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c
@@ -306,6 +306,8 @@ __ieee754_lgammal_r (long double x, int *signgamp)
}
if (se & 0x8000)
{
+ if (x < -2.0L && x > -33.0L)
+ return __lgamma_negl (x, signgamp);
t = sin_pi (x);
if (t == zero)
return one / fabsl (t); /* -integer */
diff --git a/sysdeps/ieee754/ldbl-96/lgamma_negl.c b/sysdeps/ieee754/ldbl-96/lgamma_negl.c
new file mode 100644
index 0000000..28535a8
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/lgamma_negl.c
@@ -0,0 +1,418 @@
+/* lgammal expanding around zeros.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
+
+static const long double lgamma_zeros[][2] =
+ {
+ { -0x2.74ff92c01f0d82acp+0L, 0x1.360cea0e5f8ed3ccp-68L },
+ { -0x2.bf6821437b201978p+0L, -0x1.95a4b4641eaebf4cp-64L },
+ { -0x3.24c1b793cb35efb8p+0L, -0xb.e699ad3d9ba6545p-68L },
+ { -0x3.f48e2a8f85fca17p+0L, -0xd.4561291236cc321p-68L },
+ { -0x4.0a139e16656030cp+0L, -0x3.9f0b0de18112ac18p-64L },
+ { -0x4.fdd5de9bbabf351p+0L, -0xd.0aa4076988501d8p-68L },
+ { -0x5.021a95fc2db64328p+0L, -0x2.4c56e595394decc8p-64L },
+ { -0x5.ffa4bd647d0357ep+0L, 0x2.b129d342ce12071cp-64L },
+ { -0x6.005ac9625f233b6p+0L, -0x7.c2d96d16385cb868p-68L },
+ { -0x6.fff2fddae1bbff4p+0L, 0x2.9d949a3dc02de0cp-64L },
+ { -0x7.000cff7b7f87adf8p+0L, 0x3.b7d23246787d54d8p-64L },
+ { -0x7.fffe5fe05673c3c8p+0L, -0x2.9e82b522b0ca9d3p-64L },
+ { -0x8.0001a01459fc9f6p+0L, -0xc.b3cec1cec857667p-68L },
+ { -0x8.ffffd1c425e81p+0L, 0x3.79b16a8b6da6181cp-64L },
+ { -0x9.00002e3bb47d86dp+0L, -0x6.d843fedc351deb78p-64L },
+ { -0x9.fffffb606bdfdcdp+0L, -0x6.2ae77a50547c69dp-68L },
+ { -0xa.0000049f93bb992p+0L, -0x7.b45d95e15441e03p-64L },
+ { -0xa.ffffff9466e9f1bp+0L, -0x3.6dacd2adbd18d05cp-64L },
+ { -0xb.0000006b9915316p+0L, 0x2.69a590015bf1b414p-64L },
+ { -0xb.fffffff70893874p+0L, 0x7.821be533c2c36878p-64L },
+ { -0xc.00000008f76c773p+0L, -0x1.567c0f0250f38792p-64L },
+ { -0xc.ffffffff4f6dcf6p+0L, -0x1.7f97a5ffc757d548p-64L },
+ { -0xd.00000000b09230ap+0L, 0x3.f997c22e46fc1c9p-64L },
+ { -0xd.fffffffff36345bp+0L, 0x4.61e7b5c1f62ee89p-64L },
+ { -0xe.000000000c9cba5p+0L, -0x4.5e94e75ec5718f78p-64L },
+ { -0xe.ffffffffff28c06p+0L, -0xc.6604ef30371f89dp-68L },
+ { -0xf.0000000000d73fap+0L, 0xc.6642f1bdf07a161p-68L },
+ { -0xf.fffffffffff28cp+0L, -0x6.0c6621f512e72e5p-64L },
+ { -0x1.000000000000d74p+4L, 0x6.0c6625ebdb406c48p-64L },
+ { -0x1.0ffffffffffff356p+4L, -0x9.c47e7a93e1c46a1p-64L },
+ { -0x1.1000000000000caap+4L, 0x9.c47e7a97778935ap-64L },
+ { -0x1.1fffffffffffff4cp+4L, 0x1.3c31dcbecd2f74d4p-64L },
+ { -0x1.20000000000000b4p+4L, -0x1.3c31dcbeca4c3b3p-64L },
+ { -0x1.2ffffffffffffff6p+4L, -0x8.5b25cbf5f545ceep-64L },
+ { -0x1.300000000000000ap+4L, 0x8.5b25cbf5f547e48p-64L },
+ { -0x1.4p+4L, 0x7.950ae90080894298p-64L },
+ { -0x1.4p+4L, -0x7.950ae9008089414p-64L },
+ { -0x1.5p+4L, 0x5.c6e3bdb73d5c63p-68L },
+ { -0x1.5p+4L, -0x5.c6e3bdb73d5c62f8p-68L },
+ { -0x1.6p+4L, 0x4.338e5b6dfe14a518p-72L },
+ { -0x1.6p+4L, -0x4.338e5b6dfe14a51p-72L },
+ { -0x1.7p+4L, 0x2.ec368262c7033b3p-76L },
+ { -0x1.7p+4L, -0x2.ec368262c7033b3p-76L },
+ { -0x1.8p+4L, 0x1.f2cf01972f577ccap-80L },
+ { -0x1.8p+4L, -0x1.f2cf01972f577ccap-80L },
+ { -0x1.9p+4L, 0x1.3f3ccdd165fa8d4ep-84L },
+ { -0x1.9p+4L, -0x1.3f3ccdd165fa8d4ep-84L },
+ { -0x1.ap+4L, 0xc.4742fe35272cd1cp-92L },
+ { -0x1.ap+4L, -0xc.4742fe35272cd1cp-92L },
+ { -0x1.bp+4L, 0x7.46ac70b733a8c828p-96L },
+ { -0x1.bp+4L, -0x7.46ac70b733a8c828p-96L },
+ { -0x1.cp+4L, 0x4.2862898d42174ddp-100L },
+ { -0x1.cp+4L, -0x4.2862898d42174ddp-100L },
+ { -0x1.dp+4L, 0x2.4b3f31686b15af58p-104L },
+ { -0x1.dp+4L, -0x2.4b3f31686b15af58p-104L },
+ { -0x1.ep+4L, 0x1.3932c5047d60e60cp-108L },
+ { -0x1.ep+4L, -0x1.3932c5047d60e60cp-108L },
+ { -0x1.fp+4L, 0xa.1a6973c1fade217p-116L },
+ { -0x1.fp+4L, -0xa.1a6973c1fade217p-116L },
+ { -0x2p+4L, 0x5.0d34b9e0fd6f10b8p-120L },
+ { -0x2p+4L, -0x5.0d34b9e0fd6f10b8p-120L },
+ { -0x2.1p+4L, 0x2.73024a9ba1aa36a8p-124L },
+ };
+
+static const long double e_hi = 0x2.b7e151628aed2a6cp+0L;
+static const long double e_lo = -0x1.408ea77f630b0c38p-64L;
+
+/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) in Stirling's
+ approximation to lgamma function. */
+
+static const long double lgamma_coeff[] =
+ {
+ 0x1.5555555555555556p-4L,
+ -0xb.60b60b60b60b60bp-12L,
+ 0x3.4034034034034034p-12L,
+ -0x2.7027027027027028p-12L,
+ 0x3.72a3c5631fe46aep-12L,
+ -0x7.daac36664f1f208p-12L,
+ 0x1.a41a41a41a41a41ap-8L,
+ -0x7.90a1b2c3d4e5f708p-8L,
+ 0x2.dfd2c703c0cfff44p-4L,
+ -0x1.6476701181f39edcp+0L,
+ 0xd.672219167002d3ap+0L,
+ -0x9.cd9292e6660d55bp+4L,
+ 0x8.911a740da740da7p+8L,
+ -0x8.d0cc570e255bf5ap+12L,
+ 0xa.8d1044d3708d1c2p+16L,
+ -0xe.8844d8a169abbc4p+20L,
+ };
+
+#define NCOEFF (sizeof (lgamma_coeff) / sizeof (lgamma_coeff[0]))
+
+/* Polynomial approximations to (|gamma(x)|-1)(x-n)/(x-x0), where n is
+ the integer end-point of the half-integer interval containing x and
+ x0 is the zero of lgamma in that half-integer interval. Each
+ polynomial is expressed in terms of x-xm, where xm is the midpoint
+ of the interval for which the polynomial applies. */
+
+static const long double poly_coeff[] =
+ {
+ /* Interval [-2.125, -2] (polynomial degree 13). */
+ -0x1.0b71c5c54d42eb6cp+0L,
+ -0xc.73a1dc05f349517p-4L,
+ -0x1.ec841408528b6baep-4L,
+ -0xe.37c9da26fc3b492p-4L,
+ -0x1.03cd87c5178991ap-4L,
+ -0xe.ae9ada65ece2f39p-4L,
+ 0x9.b1185505edac18dp-8L,
+ -0xe.f28c130b54d3cb2p-4L,
+ 0x2.6ec1666cf44a63bp-4L,
+ -0xf.57cb2774193bbd5p-4L,
+ 0x4.5ae64671a41b1c4p-4L,
+ -0xf.f48ea8b5bd3a7cep-4L,
+ 0x6.7d73788a8d30ef58p-4L,
+ -0x1.11e0e4b506bd272ep+0L,
+ /* Interval [-2.25, -2.125] (polynomial degree 13). */
+ -0xf.2930890d7d675a8p-4L,
+ -0xc.a5cfde054eab5cdp-4L,
+ 0x3.9c9e0fdebb0676e4p-4L,
+ -0x1.02a5ad35605f0d8cp+0L,
+ 0x9.6e9b1185d0b92edp-4L,
+ -0x1.4d8332f3d6a3959p+0L,
+ 0x1.1c0c8cacd0ced3eap+0L,
+ -0x1.c9a6f592a67b1628p+0L,
+ 0x1.d7e9476f96aa4bd6p+0L,
+ -0x2.921cedb488bb3318p+0L,
+ 0x2.e8b3fd6ca193e4c8p+0L,
+ -0x3.cb69d9d6628e4a2p+0L,
+ 0x4.95f12c73b558638p+0L,
+ -0x5.d392d0b97c02ab6p+0L,
+ /* Interval [-2.375, -2.25] (polynomial degree 14). */
+ -0xd.7d28d505d618122p-4L,
+ -0xe.69649a304098532p-4L,
+ 0xb.0d74a2827d055c5p-4L,
+ -0x1.924b09228531c00ep+0L,
+ 0x1.d49b12bccee4f888p+0L,
+ -0x3.0898bb7dbb21e458p+0L,
+ 0x4.207a6cad6fa10a2p+0L,
+ -0x6.39ee630b46093ad8p+0L,
+ 0x8.e2e25211a3fb5ccp+0L,
+ -0xd.0e85ccd8e79c08p+0L,
+ 0x1.2e45882bc17f9e16p+4L,
+ -0x1.b8b6e841815ff314p+4L,
+ 0x2.7ff8bf7504fa04dcp+4L,
+ -0x3.c192e9c903352974p+4L,
+ 0x5.8040b75f4ef07f98p+4L,
+ /* Interval [-2.5, -2.375] (polynomial degree 15). */
+ -0xb.74ea1bcfff94b2cp-4L,
+ -0x1.2a82bd590c375384p+0L,
+ 0x1.88020f828b968634p+0L,
+ -0x3.32279f040eb80fa4p+0L,
+ 0x5.57ac825175943188p+0L,
+ -0x9.c2aedcfe10f129ep+0L,
+ 0x1.12c132f2df02881ep+4L,
+ -0x1.ea94e26c0b6ffa6p+4L,
+ 0x3.66b4a8bb0290013p+4L,
+ -0x6.0cf735e01f5990bp+4L,
+ 0xa.c10a8db7ae99343p+4L,
+ -0x1.31edb212b315feeap+8L,
+ 0x2.1f478592298b3ebp+8L,
+ -0x3.c546da5957ace6ccp+8L,
+ 0x7.0e3d2a02579ba4bp+8L,
+ -0xc.b1ea961c39302f8p+8L,
+ /* Interval [-2.625, -2.5] (polynomial degree 16). */
+ -0x3.d10108c27ebafad4p-4L,
+ 0x1.cd557caff7d2b202p+0L,
+ 0x3.819b4856d3995034p+0L,
+ 0x6.8505cbad03dd3bd8p+0L,
+ 0xb.c1b2e653aa0b924p+0L,
+ 0x1.50a53a38f05f72d6p+4L,
+ 0x2.57ae00cbd06efb34p+4L,
+ 0x4.2b1563077a577e9p+4L,
+ 0x7.6989ed790138a7f8p+4L,
+ 0xd.2dd28417b4f8406p+4L,
+ 0x1.76e1b71f0710803ap+8L,
+ 0x2.9a7a096254ac032p+8L,
+ 0x4.a0e6109e2a039788p+8L,
+ 0x8.37ea17a93c877b2p+8L,
+ 0xe.9506a641143612bp+8L,
+ 0x1.b680ed4ea386d52p+12L,
+ 0x3.28a2130c8de0ae84p+12L,
+ /* Interval [-2.75, -2.625] (polynomial degree 15). */
+ -0x6.b5d252a56e8a7548p-4L,
+ 0x1.28d60383da3ac72p+0L,
+ 0x1.db6513ada8a6703ap+0L,
+ 0x2.e217118f9d34aa7cp+0L,
+ 0x4.450112c5cbd6256p+0L,
+ 0x6.4af99151e972f92p+0L,
+ 0x9.2db598b5b183cd6p+0L,
+ 0xd.62bef9c9adcff6ap+0L,
+ 0x1.379f290d743d9774p+4L,
+ 0x1.c58271ff823caa26p+4L,
+ 0x2.93a871b87a06e73p+4L,
+ 0x3.bf9db66103d7ec98p+4L,
+ 0x5.73247c111fbf197p+4L,
+ 0x7.ec8b9973ba27d008p+4L,
+ 0xb.eca5f9619b39c03p+4L,
+ 0x1.18f2e46411c78b1cp+8L,
+ /* Interval [-2.875, -2.75] (polynomial degree 14). */
+ -0x8.a41b1e4f36ff88ep-4L,
+ 0xc.da87d3b69dc0f34p-4L,
+ 0x1.1474ad5c36158ad2p+0L,
+ 0x1.761ecb90c5553996p+0L,
+ 0x1.d279bff9ae234f8p+0L,
+ 0x2.4e5d0055a16c5414p+0L,
+ 0x2.d57545a783902f8cp+0L,
+ 0x3.8514eec263aa9f98p+0L,
+ 0x4.5235e338245f6fe8p+0L,
+ 0x5.562b1ef200b256c8p+0L,
+ 0x6.8ec9782b93bd565p+0L,
+ 0x8.14baf4836483508p+0L,
+ 0x9.efaf35dc712ea79p+0L,
+ 0xc.8431f6a226507a9p+0L,
+ 0xf.80358289a768401p+0L,
+ /* Interval [-3, -2.875] (polynomial degree 13). */
+ -0xa.046d667e468f3e4p-4L,
+ 0x9.70b88dcc006c216p-4L,
+ 0xa.a8a39421c86ce9p-4L,
+ 0xd.2f4d1363f321e89p-4L,
+ 0xd.ca9aa1a3ab2f438p-4L,
+ 0xf.cf09c31f05d02cbp-4L,
+ 0x1.04b133a195686a38p+0L,
+ 0x1.22b54799d0072024p+0L,
+ 0x1.2c5802b869a36ae8p+0L,
+ 0x1.4aadf23055d7105ep+0L,
+ 0x1.5794078dd45c55d6p+0L,
+ 0x1.7759069da18bcf0ap+0L,
+ 0x1.8e672cefa4623f34p+0L,
+ 0x1.b2acfa32c17145e6p+0L,
+ };
+
+static const size_t poly_deg[] =
+ {
+ 13,
+ 13,
+ 14,
+ 15,
+ 16,
+ 15,
+ 14,
+ 13,
+ };
+
+static const size_t poly_end[] =
+ {
+ 13,
+ 27,
+ 42,
+ 58,
+ 75,
+ 91,
+ 106,
+ 120,
+ };
+
+/* Compute sin (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_sinpi (long double x)
+{
+ if (x <= 0.25L)
+ return __sinl (M_PIl * x);
+ else
+ return __cosl (M_PIl * (0.5L - x));
+}
+
+/* Compute cos (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cospi (long double x)
+{
+ if (x <= 0.25L)
+ return __cosl (M_PIl * x);
+ else
+ return __sinl (M_PIl * (0.5L - x));
+}
+
+/* Compute cot (pi * X) for -0.25 <= X <= 0.5. */
+
+static long double
+lg_cotpi (long double x)
+{
+ return lg_cospi (x) / lg_sinpi (x);
+}
+
+/* Compute lgamma of a negative argument -33 < X < -2, setting
+ *SIGNGAMP accordingly. */
+
+long double
+__lgamma_negl (long double x, int *signgamp)
+{
+ /* Determine the half-integer region X lies in, handle exact
+ integers and determine the sign of the result. */
+ int i = __floorl (-2 * x);
+ if ((i & 1) == 0 && i == -2 * x)
+ return 1.0L / 0.0L;
+ long double xn = ((i & 1) == 0 ? -i / 2 : (-i - 1) / 2);
+ i -= 4;
+ *signgamp = ((i & 2) == 0 ? -1 : 1);
+
+ SET_RESTORE_ROUNDL (FE_TONEAREST);
+
+ /* Expand around the zero X0 = X0_HI + X0_LO. */
+ long double x0_hi = lgamma_zeros[i][0], x0_lo = lgamma_zeros[i][1];
+ long double xdiff = x - x0_hi - x0_lo;
+
+ /* For arguments in the range -3 to -2, use polynomial
+ approximations to an adjusted version of the gamma function. */
+ if (i < 2)
+ {
+ int j = __floorl (-8 * x) - 16;
+ long double xm = (-33 - 2 * j) * 0.0625L;
+ long double x_adj = x - xm;
+ size_t deg = poly_deg[j];
+ size_t end = poly_end[j];
+ long double g = poly_coeff[end];
+ for (size_t j = 1; j <= deg; j++)
+ g = g * x_adj + poly_coeff[end - j];
+ return __log1pl (g * xdiff / (x - xn));
+ }
+
+ /* The result we want is log (sinpi (X0) / sinpi (X))
+ + log (gamma (1 - X0) / gamma (1 - X)). */
+ long double x_idiff = fabsl (xn - x), x0_idiff = fabsl (xn - x0_hi - x0_lo);
+ long double log_sinpi_ratio;
+ if (x0_idiff < x_idiff * 0.5L)
+ /* Use log not log1p to avoid inaccuracy from log1p of arguments
+ close to -1. */
+ log_sinpi_ratio = __ieee754_logl (lg_sinpi (x0_idiff)
+ / lg_sinpi (x_idiff));
+ else
+ {
+ /* Use log1p not log to avoid inaccuracy from log of arguments
+ close to 1. X0DIFF2 has positive sign if X0 is further from
+ XN than X is from XN, negative sign otherwise. */
+ long double x0diff2 = ((i & 1) == 0 ? xdiff : -xdiff) * 0.5L;
+ long double sx0d2 = lg_sinpi (x0diff2);
+ long double cx0d2 = lg_cospi (x0diff2);
+ log_sinpi_ratio = __log1pl (2 * sx0d2
+ * (-sx0d2 + cx0d2 * lg_cotpi (x_idiff)));
+ }
+
+ long double log_gamma_ratio;
+ long double y0 = 1 - x0_hi;
+ long double y0_eps = -x0_hi + (1 - y0) - x0_lo;
+ long double y = 1 - x;
+ long double y_eps = -x + (1 - y);
+ /* We now wish to compute LOG_GAMMA_RATIO
+ = log (gamma (Y0 + Y0_EPS) / gamma (Y + Y_EPS)). XDIFF
+ accurately approximates the difference Y0 + Y0_EPS - Y -
+ Y_EPS. Use Stirling's approximation. First, we may need to
+ adjust into the range where Stirling's approximation is
+ sufficiently accurate. */
+ long double log_gamma_adj = 0;
+ if (i < 8)
+ {
+ int n_up = (9 - i) / 2;
+ long double ny0, ny0_eps, ny, ny_eps;
+ ny0 = y0 + n_up;
+ ny0_eps = y0 - (ny0 - n_up) + y0_eps;
+ y0 = ny0;
+ y0_eps = ny0_eps;
+ ny = y + n_up;
+ ny_eps = y - (ny - n_up) + y_eps;
+ y = ny;
+ y_eps = ny_eps;
+ long double prodm1 = __lgamma_productl (xdiff, y - n_up, y_eps, n_up);
+ log_gamma_adj = -__log1pl (prodm1);
+ }
+ long double log_gamma_high
+ = (xdiff * __log1pl ((y0 - e_hi - e_lo + y0_eps) / e_hi)
+ + (y - 0.5L + y_eps) * __log1pl (xdiff / y) + log_gamma_adj);
+ /* Compute the sum of (B_2k / 2k(2k-1))(Y0^-(2k-1) - Y^-(2k-1)). */
+ long double y0r = 1 / y0, yr = 1 / y;
+ long double y0r2 = y0r * y0r, yr2 = yr * yr;
+ long double rdiff = -xdiff / (y * y0);
+ long double bterm[NCOEFF];
+ long double dlast = rdiff, elast = rdiff * yr * (yr + y0r);
+ bterm[0] = dlast * lgamma_coeff[0];
+ for (size_t j = 1; j < NCOEFF; j++)
+ {
+ long double dnext = dlast * y0r2 + elast;
+ long double enext = elast * yr2;
+ bterm[j] = dnext * lgamma_coeff[j];
+ dlast = dnext;
+ elast = enext;
+ }
+ long double log_gamma_low = 0;
+ for (size_t j = 0; j < NCOEFF; j++)
+ log_gamma_low += bterm[NCOEFF - 1 - j];
+ log_gamma_ratio = log_gamma_high + log_gamma_low;
+
+ return log_sinpi_ratio + log_gamma_ratio;
+}
diff --git a/sysdeps/ieee754/ldbl-96/lgamma_product.c b/sysdeps/ieee754/ldbl-96/lgamma_product.c
new file mode 100644
index 0000000..ada1526
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/lgamma_product.c
@@ -0,0 +1,37 @@
+/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+
+double
+__lgamma_product (double t, double x, double x_eps, int n)
+{
+ long double x_full = (long double) x + (long double) x_eps;
+ long double ret = 0;
+ for (int i = 0; i < n; i++)
+ ret += (t / (x_full + i)) * (1 + ret);
+ return ret;
+}
diff --git a/sysdeps/ieee754/ldbl-96/lgamma_productl.c b/sysdeps/ieee754/ldbl-96/lgamma_productl.c
new file mode 100644
index 0000000..cf0c778
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-96/lgamma_productl.c
@@ -0,0 +1,82 @@
+/* Compute a product of 1 + (T/X), 1 + (T/(X+1)), ....
+ Copyright (C) 2015 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#include <math.h>
+#include <math_private.h>
+#include <float.h>
+
+/* Calculate X * Y exactly and store the result in *HI + *LO. It is
+ given that the values are small enough that no overflow occurs and
+ large enough (or zero) that no underflow occurs. */
+
+static void
+mul_split (long double *hi, long double *lo, long double x, long double y)
+{
+#ifdef __FP_FAST_FMAL
+ /* Fast built-in fused multiply-add. */
+ *hi = x * y;
+ *lo = __builtin_fmal (x, y, -*hi);
+#elif defined FP_FAST_FMAL
+ /* Fast library fused multiply-add, compiler before GCC 4.6. */
+ *hi = x * y;
+ *lo = __fmal (x, y, -*hi);
+#else
+ /* Apply Dekker's algorithm. */
+ *hi = x * y;
+# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1)
+ long double x1 = x * C;
+ long double y1 = y * C;
+# undef C
+ x1 = (x - x1) + x1;
+ y1 = (y - y1) + y1;
+ long double x2 = x - x1;
+ long double y2 = y - y1;
+ *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2;
+#endif
+}
+
+/* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS +
+ 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that
+ all the values X + 1, ..., X + N - 1 are exactly representable, and
+ X_EPS / X is small enough that factors quadratic in it can be
+ neglected. */
+
+long double
+__lgamma_productl (long double t, long double x, long double x_eps, int n)
+{
+ long double ret = 0, ret_eps = 0;
+ for (int i = 0; i < n; i++)
+ {
+ long double xi = x + i;
+ long double quot = t / xi;
+ long double mhi, mlo;
+ mul_split (&mhi, &mlo, quot, xi);
+ long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi);
+ /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */
+ long double rhi, rlo;
+ mul_split (&rhi, &rlo, ret, quot);
+ long double rpq = ret + quot;
+ long double rpq_eps = (ret - rpq) + quot;
+ long double nret = rpq + rhi;
+ long double nret_eps = (rpq - nret) + rhi;
+ ret_eps += (rpq_eps + nret_eps + rlo + ret_eps * quot
+ + quot_lo + quot_lo * (ret + ret_eps));
+ ret = nret;
+ }
+ return ret + ret_eps;
+}
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 19d263e..b4a65d0 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1642,36 +1642,36 @@ ildouble: 4
ldouble: 4
Function: "gamma":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 2
-ldouble: 2
-
-Function: "gamma_downward":
-double: 4
+double: 3
float: 3
-idouble: 4
+idouble: 3
ifloat: 3
-ildouble: 6
-ldouble: 6
+ildouble: 3
+ldouble: 3
-Function: "gamma_towardzero":
-double: 4
-float: 3
-idouble: 4
-ifloat: 3
-ildouble: 6
-ldouble: 6
+Function: "gamma_downward":
+double: 5
+float: 4
+idouble: 5
+ifloat: 4
+ildouble: 7
+ldouble: 7
-Function: "gamma_upward":
-double: 4
+Function: "gamma_towardzero":
+double: 5
float: 4
-idouble: 4
+idouble: 5
ifloat: 4
-ildouble: 4
-ldouble: 4
+ildouble: 7
+ldouble: 7
+
+Function: "gamma_upward":
+double: 5
+float: 5
+idouble: 5
+ifloat: 5
+ildouble: 5
+ldouble: 5
Function: "hypot":
double: 1
@@ -1794,36 +1794,36 @@ ildouble: 5
ldouble: 5
Function: "lgamma":
-double: 2
-float: 2
-idouble: 2
-ifloat: 2
-ildouble: 2
-ldouble: 2
-
-Function: "lgamma_downward":
-double: 4
+double: 3
float: 3
-idouble: 4
+idouble: 3
ifloat: 3
-ildouble: 6
-ldouble: 6
+ildouble: 3
+ldouble: 3
-Function: "lgamma_towardzero":
-double: 4
-float: 3
-idouble: 4
-ifloat: 3
-ildouble: 6
-ldouble: 6
+Function: "lgamma_downward":
+double: 5
+float: 4
+idouble: 5
+ifloat: 4
+ildouble: 7
+ldouble: 7
-Function: "lgamma_upward":
-double: 4
+Function: "lgamma_towardzero":
+double: 5
float: 4
-idouble: 4
+idouble: 5
ifloat: 4
-ildouble: 4
-ldouble: 4
+ildouble: 7
+ldouble: 7
+
+Function: "lgamma_upward":
+double: 5
+float: 5
+idouble: 5
+ifloat: 5
+ildouble: 5
+ldouble: 5
Function: "log":
float: 1
--
Joseph S. Myers
joseph@codesourcery.com