From 4f3a641e6e44fab414f7e384ac55e5b9e8d6fc7e Mon Sep 17 00:00:00 2001
From: Russ Cox <rsc@golang.org>
Date: Thu, 6 Oct 2016 10:39:50 -0400
Subject: [PATCH] math: fix Gamma(-171.5) on all platforms

Using 387 mode was computing it without underflow to zero,
apparently due to an 80-bit intermediate. Avoid underflow even
with 64-bit floats.

This eliminates the TODOs in the test suite.

Fixes linux-386-387 build and fixes #11441.

Change-Id: I8abaa63bfdf040438a95625d1cb61042f0302473
Reviewed-on: https://go-review.googlesource.com/30540
Run-TryBot: Russ Cox <rsc@golang.org>
Reviewed-by: Brad Fitzpatrick <bradfitz@golang.org>
---
 src/math/all_test.go | 12 +++++++++---
 src/math/gamma.go    | 35 ++++++++++++++++++++++++-----------
 2 files changed, 33 insertions(+), 14 deletions(-)

diff --git a/src/math/all_test.go b/src/math/all_test.go
index c77dcfea6e..3d8cd7223d 100644
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -1235,9 +1235,9 @@ var vfgamma = [][2]float64{
 	{-100.5, -3.3536908198076787e-159},
 	{-160.5, -5.255546447007829e-286},
 	{-170.5, -3.3127395215386074e-308},
-	{-171.5, 0},               // TODO: 1.9316265431712e-310
-	{-176.5, Copysign(0, -1)}, // TODO: -1.196e-321
-	{-177.5, 0},               // TODO: 5e-324
+	{-171.5, 1.9316265431712e-310},
+	{-176.5, -1.196e-321},
+	{-177.5, 5e-324},
 	{-178.5, Copysign(0, -1)},
 	{-179.5, 0},
 	{-201.0001, 0},
@@ -1802,6 +1802,12 @@ var logbBC = []float64{
 }
 
 func tolerance(a, b, e float64) bool {
+	// Multiplying by e here can underflow denormal values to zero.
+	// Check a==b so that at least if a and b are small and identical
+	// we say they match.
+	if a == b {
+		return true
+	}
 	d := a - b
 	if d < 0 {
 		d = -d
diff --git a/src/math/gamma.go b/src/math/gamma.go
index c965f275e7..514260be05 100644
--- a/src/math/gamma.go
+++ b/src/math/gamma.go
@@ -91,10 +91,15 @@ var _gamS = [...]float64{
 }
 
 // Gamma function computed by Stirling's formula.
-// The polynomial is valid for 33 <= x <= 172.
-func stirling(x float64) float64 {
-	if x > 171.625 {
-		return Inf(1)
+// The pair of results must be multiplied together to get the actual answer.
+// The multiplication is left to the caller so that, if careful, the caller can avoid
+// infinity for 172 <= x <= 180.
+// The polynomial is valid for 33 <= x <= 172; larger values are only used
+// in reciprocal and produce denormalized floats. The lower precision there
+// masks any imprecision in the polynomial.
+func stirling(x float64) (float64, float64) {
+	if x > 200 {
+		return Inf(1), 1
 	}
 	const (
 		SqrtTwoPi   = 2.506628274631000502417
@@ -102,15 +107,15 @@ func stirling(x float64) float64 {
 	)
 	w := 1 / x
 	w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
-	y := Exp(x)
+	y1 := Exp(x)
+	y2 := 1.0
 	if x > MaxStirling { // avoid Pow() overflow
 		v := Pow(x, 0.5*x-0.25)
-		y = v * (v / y)
+		y1, y2 = v, v/y1
 	} else {
-		y = Pow(x, x-0.5) / y
+		y1 = Pow(x, x-0.5) / y1
 	}
-	y = SqrtTwoPi * y * w
-	return y
+	return y1, SqrtTwoPi * w * y2
 }
 
 // Gamma returns the Gamma function of x.
@@ -138,7 +143,8 @@ func Gamma(x float64) float64 {
 	p := Floor(q)
 	if q > 33 {
 		if x >= 0 {
-			return stirling(x)
+			y1, y2 := stirling(x)
+			return y1 * y2
 		}
 		// Note: x is negative but (checked above) not a negative integer,
 		// so x must be small enough to be in range for conversion to int64.
@@ -156,7 +162,14 @@ func Gamma(x float64) float64 {
 		if z == 0 {
 			return Inf(signgam)
 		}
-		z = Pi / (Abs(z) * stirling(q))
+		sq1, sq2 := stirling(q)
+		absz := Abs(z)
+		d := absz * sq1 * sq2
+		if IsInf(d, 0) {
+			z = Pi / absz / sq1 / sq2
+		} else {
+			z = Pi / d
+		}
 		return float64(signgam) * z
 	}