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>

diff --git a/src/math/all_test.go b/src/math/all_test.go
index c77dcfea6e33b4615707b54659c780fb31f89f47..3d8cd7223de0681cc0609c54da3c1f4979826d95 100644 (file)
--- 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 c965f275e7601341a2aa779246a543cca5585221..514260be05b3a5848c10829673332da7f27f9527 100644 (file)
--- 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
        }