From e95cc9b70040edf2e3a83fd319cbac1c396b2a13 Mon Sep 17 00:00:00 2001
From: Bruce Evans <bde@FreeBSD.org>
Date: Wed, 2 May 2007 15:24:49 +0000
Subject: [PATCH] Fix tgamma() on some special args: (1) tgamma(-Inf) returned
 +Inf and failed to raise any exception, but     should always have raised an
 exception, and should behave like     tgamma(negative integer). (2)
 tgamma(negative integer) returned +Inf and raised divide-by-zero,     but
 should return NaN and raise "invalid" on any IEEEish system. (3) About half
 of the 2**52 negative intgers between -2**53 and -2**52     were
 misclassified as non-integers by using floor(x + 0.5) to round     to
 nearest, so tgamma(x) was wrong (+-0 instead of +Inf and now NaN)     on
 these args.  The floor() expression is hard to use since rounding     of (x +
 0.5) may give x or x + 1, depending on |x| and the current     rounding mode.
  The fixed version uses ceil(x) to classify x before     operating on x and
 ends up being more efficient since ceil(x) is     needed anyway. (4) On at
 least the problematic args in (3), tgamma() raised a spurious     inexact.
 (5) tgamma(large positive) raised divide-by-zero but should raise overflow.
 (6) tgamma(+Inf) raised divide-by-zero but should not raise any exception.
 (7) Raise inexact for tiny |x| in a way that has some chance of not being    
 optimized away.

The fix for (5) and (6), and probably for (2), also prevents -O optimizing
away the exception.

PR:		112180 (2)
Standards:	Annex F in C99 (IEC 60559 binding) requires (1), (2) and (6).
---
 lib/msun/bsdsrc/b_tgamma.c | 29 +++++++++++++++--------------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/lib/msun/bsdsrc/b_tgamma.c b/lib/msun/bsdsrc/b_tgamma.c
index 784e563cdba2..5593bbd95dc3 100644
--- a/lib/msun/bsdsrc/b_tgamma.c
+++ b/lib/msun/bsdsrc/b_tgamma.c
@@ -49,7 +49,7 @@ __FBSDID("$FreeBSD$");
 
 /* METHOD:
  * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
- * 	At negative integers, return +Inf and raise divide-by-zero.
+ * 	At negative integers, return NaN and raise invalid.
  *
  * x < 6.5:
  *	Use argument reduction G(x+1) = xG(x) to reach the
@@ -66,12 +66,12 @@ __FBSDID("$FreeBSD$");
  *	avoid premature round-off.
  *
  * Special values:
- *	-Inf:			return +Inf (without raising any exception!);
- *	negative integer:	return +Inf and raise divide-by-zero;
+ *	-Inf:			return NaN and raise invalid;
+ *	negative integer:	return NaN and raise invalid;
  *	other x ~< 177.79:	return +-0 and raise underflow;
  *	+-0:			return +-Inf and raise divide-by-zero;
- *	finite x ~> 171.63:	return +Inf and raise divide-by-zero(!);
- *	+Inf:			return +Inf and raise divide-by-zero(!);
+ *	finite x ~> 171.63:	return +Inf and raise overflow;
+ *	+Inf:			return +Inf;
  *	NaN: 			return NaN.
  *
  * Accuracy: tgamma(x) is accurate to within
@@ -135,7 +135,7 @@ tgamma(x)
 
 	if (x >= 6) {
 		if(x > 171.63)
-			return(one/zero);
+			return (x / zero);
 		u = large_gam(x);
 		return(__exp__D(u.a, u.b));
 	} else if (x >= 1.0 + LEFT + x0)
@@ -143,12 +143,11 @@ tgamma(x)
 	else if (x > 1.e-17)
 		return (smaller_gam(x));
 	else if (x > -1.e-17) {
-		if (x == 0.0)
-			return (one/x);
-		one+1e-20;		/* Raise inexact flag. */
+		if (x != 0.0)
+			u.a = one - tiny;	/* raise inexact */
 		return (one/x);
 	} else if (!finite(x))
-		return (x*x);		/* x = NaN, -Inf */
+		return (x - x);		/* x is NaN or -Inf */
 	else
 		return (neg_gam(x));
 }
@@ -282,11 +281,13 @@ neg_gam(x)
 	struct Double lg, lsine;
 	double y, z;
 
-	y = floor(x + .5);
+	y = ceil(x);
 	if (y == x)		/* Negative integer. */
-		return (one/zero);
-	z = fabs(x - y);
-	y = .5*ceil(x);
+		return ((x - x) / zero);
+	z = y - x;
+	if (z > 0.5)
+		z = one - z;
+	y = 0.5 * y;
 	if (y == ceil(y))
 		sgn = -1;
 	if (z < .25)