Fix lgamma (negative) inaccuracy (bug 2542, bug 2543, bug 2558).

The existing implementations of lgamma functions (except for the ia64
versions) use the reflection formula for negative arguments.  This
suffers large inaccuracy from cancellation near zeros of lgamma (near
where the gamma function is +/- 1).

This patch fixes this inaccuracy.  For arguments above -2, there are
no zeros and no large cancellation, while for sufficiently large
negative arguments the zeros are so close to integers that even for
integers +/- 1ulp the log(gamma(1-x)) term dominates and cancellation
is not significant.  Thus, it is only necessary to take special care
about cancellation for arguments around a limited number of zeros.

Accordingly, this patch uses precomputed tables of relevant zeros,
expressed as the sum of two floating-point values.  The log of the
ratio of two sines can be computed accurately using log1p in cases
where log would lose accuracy.  The log of the ratio of two gamma(1-x)
values can be computed using Stirling's approximation (the difference
between two values of that approximation to lgamma being computable
without computing the two values and then subtracting), with
appropriate adjustments (which don't reduce accuracy too much) in
cases where 1-x is too small to use Stirling's approximation directly.

In the interval from -3 to -2, using the ratios of sines and of
gamma(1-x) can still produce too much cancellation between those two
parts of the computation (and that interval is also the worst interval
for computing the ratio between gamma(1-x) values, which computation
becomes more accurate, while being less critical for the final result,
for larger 1-x).  Because this can result in errors slightly above
those accepted in glibc, this interval is instead dealt with by
polynomial approximations.  Separate polynomial approximations to
(|gamma(x)|-1)(x-n)/(x-x0) are used for each interval of length 1/8
from -3 to -2, where n (-3 or -2) is the nearest integer to the
1/8-interval and x0 is the zero of lgamma in the relevant half-integer
interval (-3 to -2.5 or -2.5 to -2).

Together, the two approaches are intended to give sufficient accuracy
for all negative arguments in the problem range.  Outside that range,
the previous implementation continues to be used.

Tested for x86_64, x86, mips64 and powerpc.  The mips64 and powerpc
testing shows up pre-existing problems for ldbl-128 and ldbl-128ibm
with large negative arguments giving spurious "invalid" exceptions
(exposed by newly added tests for cases this patch doesn't affect the
logic for); I'll address those problems separately.

	[BZ #2542]
	[BZ #2543]
	[BZ #2558]
	* sysdeps/ieee754/dbl-64/e_lgamma_r.c (__ieee754_lgamma_r): Call
	__lgamma_neg for arguments from -28.0 to -2.0.
	* sysdeps/ieee754/flt-32/e_lgammaf_r.c (__ieee754_lgammaf_r): Call
	__lgamma_negf for arguments from -15.0 to -2.0.
	* sysdeps/ieee754/ldbl-128/e_lgammal_r.c (__ieee754_lgammal_r):
	Call __lgamma_negl for arguments from -48.0 or -50.0 to -2.0.
	* sysdeps/ieee754/ldbl-96/e_lgammal_r.c (__ieee754_lgammal_r):
	Call __lgamma_negl for arguments from -33.0 to -2.0.
	* sysdeps/ieee754/dbl-64/lgamma_neg.c: New file.
	* sysdeps/ieee754/dbl-64/lgamma_product.c: Likewise.
	* sysdeps/ieee754/flt-32/lgamma_negf.c: Likewise.
	* sysdeps/ieee754/flt-32/lgamma_productf.c: Likewise.
	* sysdeps/ieee754/ldbl-128/lgamma_negl.c: Likewise.
	* sysdeps/ieee754/ldbl-128/lgamma_productl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/lgamma_negl.c: Likewise.
	* sysdeps/ieee754/ldbl-128ibm/lgamma_productl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/lgamma_negl.c: Likewise.
	* sysdeps/ieee754/ldbl-96/lgamma_product.c: Likewise.
	* sysdeps/ieee754/ldbl-96/lgamma_productl.c: Likewise.
	* sysdeps/generic/math_private.h (__lgamma_negf): New prototype.
	(__lgamma_neg): Likewise.
	(__lgamma_negl): Likewise.
	(__lgamma_product): Likewise.
	(__lgamma_productl): Likewise.
	* math/Makefile (libm-calls): Add lgamma_neg and lgamma_product.
	* math/auto-libm-test-in: Add more tests of lgamma.
	* math/auto-libm-test-out: Regenerated.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

1 files changed, 82 insertions, 0 deletions

diff --git a/sysdeps/ieee754/dbl-64/lgamma_product.c b/sysdeps/ieee754/dbl-64/lgamma_product.c
new file mode 100644
index 0000000000..8f877a8069
--- /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;
+}