Improve tgamma accuracy (bug 18613).

In non-default rounding modes, tgamma can be slightly less accurate
than permitted by glibc's accuracy goals.

Part of the problem is error accumulation, addressed in this patch by
setting round-to-nearest for internal computations.  However, there
was also a bug in the code dealing with computing pow (x + n, x + n)
where x + n is not exactly representable, providing another source of
error even in round-to-nearest mode; it was necessary to address both
bugs to get errors for all testcases within glibc's accuracy goals.
Given this second fix, accuracy in round-to-nearest mode is also
improved (hence regeneration of ulps for tgamma should be from scratch
- truncate libm-test-ulps or at least remove existing tgamma entries -
so that the expected ulps can be reduced).

Some additional complications also arose.  Certain tgamma tests should
strictly, according to IEEE semantics, overflow or not depending on
the rounding mode; this is beyond the scope of glibc's accuracy goals
for any function without exactly-determined results, but
gen-auto-libm-tests doesn't handle being lax there as it does for
underflow.  (libm-test.inc also doesn't handle being lax about whether
the result in cases very close to the overflow threshold is infinity
or a finite value close to overflow, but that doesn't cause problems
in this case though I've seen it cause problems with random test
generation for some functions.)  Thus, spurious-overflow markings,
with a comment, are added to auto-libm-test-in (no bug in Bugzilla
because the issue is with the testsuite, not a user-visible bug in
glibc).  And on x86, after the patch I saw ERANGE issues as previously
reported by Carlos (see my commentary in
<https://sourceware.org/ml/libc-alpha/2015-01/msg00485.html>), which
needed addressing by ensuring excess range and precision were
eliminated at various points if FLT_EVAL_METHOD != 0.

I also noticed and fixed a cosmetic issue where 1.0f was used in long
double functions and should have been 1.0L.

This completes the move of all functions to testing in all rounding
modes with ALL_RM_TEST, so gen-libm-have-vector-test.sh is updated to
remove the workaround for some functions not using ALL_RM_TEST.

Tested for x86_64, x86, mips64 and powerpc.

	[BZ #18613]
	* sysdeps/ieee754/dbl-64/e_gamma_r.c (gamma_positive): Take log of
	X_ADJ not X when adjusting exponent.
	(__ieee754_gamma_r): Do intermediate computations in
	round-to-nearest then adjust overflowing and underflowing results
	as needed.
	* sysdeps/ieee754/flt-32/e_gammaf_r.c (gammaf_positive): Take log
	of X_ADJ not X when adjusting exponent.
	(__ieee754_gammaf_r): Do intermediate computations in
	round-to-nearest then adjust overflowing and underflowing results
	as needed.
	* sysdeps/ieee754/ldbl-128/e_gammal_r.c (gammal_positive): Take
	log of X_ADJ not X when adjusting exponent.
	(__ieee754_gammal_r): Do intermediate computations in
	round-to-nearest then adjust overflowing and underflowing results
	as needed.  Use 1.0L not 1.0f as numerator of division.
	* sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c (gammal_positive): Take
	log of X_ADJ not X when adjusting exponent.
	(__ieee754_gammal_r): Do intermediate computations in
	round-to-nearest then adjust overflowing and underflowing results
	as needed.  Use 1.0L not 1.0f as numerator of division.
	* sysdeps/ieee754/ldbl-96/e_gammal_r.c (gammal_positive): Take log
	of X_ADJ not X when adjusting exponent.
	(__ieee754_gammal_r): Do intermediate computations in
	round-to-nearest then adjust overflowing and underflowing results
	as needed.  Use 1.0L not 1.0f as numerator of division.
	* math/libm-test.inc (tgamma_test_data): Remove one test.  Moved
	to auto-libm-test-in.
	(tgamma_test): Use ALL_RM_TEST.
	* math/auto-libm-test-in: Add one test of tgamma.  Mark some other
	tests of tgamma with spurious-overflow.
	* math/auto-libm-test-out: Regenerated.
	* math/gen-libm-have-vector-test.sh: Do not check for START.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

1 files changed, 49 insertions, 25 deletions

diff --git a/sysdeps/ieee754/ldbl-128/e_gammal_r.c b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
index 39c0638ac3..c51b050e0e 100644
--- a/sysdeps/ieee754/ldbl-128/e_gammal_r.c
+++ b/sysdeps/ieee754/ldbl-128/e_gammal_r.c
@@ -109,7 +109,7 @@ gammal_positive (long double x, int *exp2_adj)
 			 * __ieee754_expl (-x_adj)
 			 * __ieee754_sqrtl (2 * M_PIl / x_adj)
 			 / prod);
-      exp_adj += x_eps * __ieee754_logl (x);
+      exp_adj += x_eps * __ieee754_logl (x_adj);
       long double bsum = gamma_coeff[NCOEFF - 1];
       long double x_adj2 = x_adj * x_adj;
       for (size_t i = 1; i <= NCOEFF - 1; i++)
@@ -124,6 +124,7 @@ __ieee754_gammal_r (long double x, int *signgamp)
 {
   int64_t hx;
   u_int64_t lx;
+  long double ret;
 
   GET_LDOUBLE_WORDS64 (hx, lx, x);
 
@@ -159,35 +160,58 @@ __ieee754_gammal_r (long double x, int *signgamp)
       *signgamp = 0;
       return LDBL_MAX * LDBL_MAX;
     }
-  else if (x > 0.0L)
+  else
     {
-      *signgamp = 0;
-      int exp2_adj;
-      long double ret = gammal_positive (x, &exp2_adj);
-      return __scalbnl (ret, exp2_adj);
+      SET_RESTORE_ROUNDL (FE_TONEAREST);
+      if (x > 0.0L)
+	{
+	  *signgamp = 0;
+	  int exp2_adj;
+	  ret = gammal_positive (x, &exp2_adj);
+	  ret = __scalbnl (ret, exp2_adj);
+	}
+      else if (x >= -LDBL_EPSILON / 4.0L)
+	{
+	  *signgamp = 0;
+	  ret = 1.0L / x;
+	}
+      else
+	{
+	  long double tx = __truncl (x);
+	  *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
+	  if (x <= -1775.0L)
+	    /* Underflow.  */
+	    ret = LDBL_MIN * LDBL_MIN;
+	  else
+	    {
+	      long double frac = tx - x;
+	      if (frac > 0.5L)
+		frac = 1.0L - frac;
+	      long double sinpix = (frac <= 0.25L
+				    ? __sinl (M_PIl * frac)
+				    : __cosl (M_PIl * (0.5L - frac)));
+	      int exp2_adj;
+	      ret = M_PIl / (-x * sinpix
+			     * gammal_positive (-x, &exp2_adj));
+	      ret = __scalbnl (ret, -exp2_adj);
+	    }
+	}
     }
-  else if (x >= -LDBL_EPSILON / 4.0L)
+  if (isinf (ret) && x != 0)
     {
-      *signgamp = 0;
-      return 1.0f / x;
+      if (*signgamp < 0)
+	return -(-__copysignl (LDBL_MAX, ret) * LDBL_MAX);
+      else
+	return __copysignl (LDBL_MAX, ret) * LDBL_MAX;
     }
-  else
+  else if (ret == 0)
     {
-      long double tx = __truncl (x);
-      *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
-      if (x <= -1775.0L)
-	/* Underflow.  */
-	return LDBL_MIN * LDBL_MIN;
-      long double frac = tx - x;
-      if (frac > 0.5L)
-	frac = 1.0L - frac;
-      long double sinpix = (frac <= 0.25L
-			    ? __sinl (M_PIl * frac)
-			    : __cosl (M_PIl * (0.5L - frac)));
-      int exp2_adj;
-      long double ret = M_PIl / (-x * sinpix
-				 * gammal_positive (-x, &exp2_adj));
-      return __scalbnl (ret, -exp2_adj);
+      if (*signgamp < 0)
+	return -(-__copysignl (LDBL_MIN, ret) * LDBL_MIN);
+      else
+	return __copysignl (LDBL_MIN, ret) * LDBL_MIN;
     }
+  else
+    return ret;
 }
 strong_alias (__ieee754_gammal_r, __gammal_r_finite)