math: fix Gamma(x) for x < -170.5 and other corner cases

Fixes #11441.

Test tables generated by

	package main

	import (
		"bytes"
		"fmt"
		"log"
		"os/exec"
		"strconv"
		"strings"
	)

	var inputs = []float64{
		0.5,
		1.5,
		2.5,
		3.5,
		-0.5,
		-1.5,
		-2.5,
		-3.5,
		0.1,
		0.01,
		1e-8,
		1e-16,
		1e-3,
		1e-16,
		1e-308,
		5.6e-309,
		5.5e-309,
		1e-309,
		1e-323,
		5e-324,
		-0.1,
		-0.01,
		-1e-8,
		-1e-16,
		-1e-3,
		-1e-16,
		-1e-308,
		-5.6e-309,
		-5.5e-309,
		-1e-300 / 1e9,
		-1e-300 / 1e23,
		-5e-300 / 1e24,
		-0.9999999999999999,
		-1.0000000000000002,
		-1.9999999999999998,
		-2.0000000000000004,
		-100.00000000000001,
		-99.999999999999986,
		17,
		171,
		171.6,
		171.624,
		171.625,
		172,
		2000,
		-100.5,
		-160.5,
		-170.5,
		-171.5,
		-176.5,
		-177.5,
		-178.5,
		-179.5,
		-201.0001,
		-202.9999,
		-1000.5,
		-1000000000.3,
		-4503599627370495.5,
		-63.349078729022985,
		-127.45117632943295,
	}

	func main() {
		var buf bytes.Buffer
		for _, v := range inputs {
			fmt.Fprintf(&buf, "gamma(%.1000g)\n", v)
		}
		cmd := exec.Command("gp", "-q")
		cmd.Stdin = &buf
		out, err := cmd.CombinedOutput()
		if err != nil {
			log.Fatalf("gp: %v", err)
		}
		f := strings.Split(string(out), "\n")
		if len(f) > 0 && f[len(f)-1] == "" {
			f = f[:len(f)-1]
		}
		if len(f) != len(inputs) {
			log.Fatalf("gp: wrong output count\n%s\n", out)
		}
		for i, g := range f {
			gf, err := strconv.ParseFloat(strings.Replace(g, " E", "e", -1), 64)
			if err != nil {
				if strings.Contains(err.Error(), "value out of range") {
					if strings.HasPrefix(g, "-") {
						fmt.Printf("\t{%g, Inf(-1)},\n", inputs[i])
					} else {
						fmt.Printf("\t{%g, Inf(1)},\n", inputs[i])
					}
					continue
				}
				log.Fatal(err)
			}
			if gf == 0 && strings.HasPrefix(g, "-") {
				fmt.Printf("\t{%g, Copysign(0, -1)},\n", inputs[i])
				continue
			}
			fmt.Printf("\t{%g, %g},\n", inputs[i], gf)
		}
	}

Change-Id: Ie98c7751d92b8ffb40e8313f5ea10df0890e2feb
Reviewed-on: https://go-review.googlesource.com/30146
Run-TryBot: Russ Cox <rsc@golang.org>
Reviewed-by: Quentin Smith <quentin@golang.org>

2 files changed, 101 insertions, 21 deletions

diff --git a/src/math/all_test.go b/src/math/all_test.go
index 882400527b..c77dcfea6e 100644
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -1165,21 +1165,88 @@ var frexpSC = []fi{
 	{NaN(), 0},
 }
 
-var vfgammaSC = []float64{
-	Inf(-1),
-	-3,
-	Copysign(0, -1),
-	0,
-	Inf(1),
-	NaN(),
-}
-var gammaSC = []float64{
-	NaN(),
-	NaN(),
-	Inf(-1),
-	Inf(1),
-	Inf(1),
-	NaN(),
+var vfgamma = [][2]float64{
+	{Inf(1), Inf(1)},
+	{Inf(-1), NaN()},
+	{0, Inf(1)},
+	{Copysign(0, -1), Inf(-1)},
+	{NaN(), NaN()},
+	{-1, NaN()},
+	{-2, NaN()},
+	{-3, NaN()},
+	{-1e16, NaN()},
+	{-1e300, NaN()},
+	{1.7e308, Inf(1)},
+
+	// Test inputs inspired by Python test suite.
+	// Outputs computed at high precision by PARI/GP.
+	// If recomputing table entries, be careful to use
+	// high-precision (%.1000g) formatting of the float64 inputs.
+	// For example, -2.0000000000000004 is the float64 with exact value
+	// -2.00000000000000044408920985626161695, and
+	// gamma(-2.0000000000000004) = -1249999999999999.5386078562728167651513, while
+	// gamma(-2.00000000000000044408920985626161695) = -1125899906826907.2044875028130093136826.
+	// Thus the table lists -1.1258999068426235e+15 as the answer.
+	{0.5, 1.772453850905516},
+	{1.5, 0.886226925452758},
+	{2.5, 1.329340388179137},
+	{3.5, 3.3233509704478426},
+	{-0.5, -3.544907701811032},
+	{-1.5, 2.363271801207355},
+	{-2.5, -0.9453087204829419},
+	{-3.5, 0.2700882058522691},
+	{0.1, 9.51350769866873},
+	{0.01, 99.4325851191506},
+	{1e-08, 9.999999942278434e+07},
+	{1e-16, 1e+16},
+	{0.001, 999.4237724845955},
+	{1e-16, 1e+16},
+	{1e-308, 1e+308},
+	{5.6e-309, 1.7857142857142864e+308},
+	{5.5e-309, Inf(1)},
+	{1e-309, Inf(1)},
+	{1e-323, Inf(1)},
+	{5e-324, Inf(1)},
+	{-0.1, -10.686287021193193},
+	{-0.01, -100.58719796441078},
+	{-1e-08, -1.0000000057721567e+08},
+	{-1e-16, -1e+16},
+	{-0.001, -1000.5782056293586},
+	{-1e-16, -1e+16},
+	{-1e-308, -1e+308},
+	{-5.6e-309, -1.7857142857142864e+308},
+	{-5.5e-309, Inf(-1)},
+	{-1e-309, Inf(-1)},
+	{-1e-323, Inf(-1)},
+	{-5e-324, Inf(-1)},
+	{-0.9999999999999999, -9.007199254740992e+15},
+	{-1.0000000000000002, 4.5035996273704955e+15},
+	{-1.9999999999999998, 2.2517998136852485e+15},
+	{-2.0000000000000004, -1.1258999068426235e+15},
+	{-100.00000000000001, -7.540083334883109e-145},
+	{-99.99999999999999, 7.540083334884096e-145},
+	{17, 2.0922789888e+13},
+	{171, 7.257415615307999e+306},
+	{171.6, 1.5858969096672565e+308},
+	{171.624, 1.7942117599248104e+308},
+	{171.625, Inf(1)},
+	{172, Inf(1)},
+	{2000, Inf(1)},
+	{-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
+	{-178.5, Copysign(0, -1)},
+	{-179.5, 0},
+	{-201.0001, 0},
+	{-202.9999, Copysign(0, -1)},
+	{-1000.5, Copysign(0, -1)},
+	{-1.0000000003e+09, Copysign(0, -1)},
+	{-4.5035996273704955e+15, 0},
+	{-63.349078729022985, 4.177797167776188e-88},
+	{-127.45117632943295, 1.183111089623681e-214},
 }
 
 var vfhypotSC = [][2]float64{
@@ -2147,9 +2214,18 @@ func TestGamma(t *testing.T) {
 			t.Errorf("Gamma(%g) = %g, want %g", vf[i], f, gamma[i])
 		}
 	}
-	for i := 0; i < len(vfgammaSC); i++ {
-		if f := Gamma(vfgammaSC[i]); !alike(gammaSC[i], f) {
-			t.Errorf("Gamma(%g) = %g, want %g", vfgammaSC[i], f, gammaSC[i])
+	for _, g := range vfgamma {
+		f := Gamma(g[0])
+		var ok bool
+		if IsNaN(g[1]) || IsInf(g[1], 0) || g[1] == 0 || f == 0 {
+			ok = alike(g[1], f)
+		} else if g[0] > -50 && g[0] <= 171 {
+			ok = veryclose(g[1], f)
+		} else {
+			ok = close(g[1], f)
+		}
+		if !ok {
+			t.Errorf("Gamma(%g) = %g, want %g", g[0], f, g[1])
 		}
 	}
 }
diff --git a/src/math/gamma.go b/src/math/gamma.go
index 841ec114a0..c965f275e7 100644
--- a/src/math/gamma.go
+++ b/src/math/gamma.go
@@ -93,6 +93,9 @@ 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)
+	}
 	const (
 		SqrtTwoPi   = 2.506628274631000502417
 		MaxStirling = 143.01608
@@ -130,8 +133,6 @@ func Gamma(x float64) float64 {
 			return Inf(-1)
 		}
 		return Inf(1)
-	case x < -170.5674972726612 || x > 171.61447887182298:
-		return Inf(1)
 	}
 	q := Abs(x)
 	p := Floor(q)
@@ -139,8 +140,11 @@ func Gamma(x float64) float64 {
 		if x >= 0 {
 			return stirling(x)
 		}
+		// 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.
+		// If |x| were >= 2⁶³ it would have to be an integer.
 		signgam := 1
-		if ip := int(p); ip&1 == 0 {
+		if ip := int64(p); ip&1 == 0 {
 			signgam = -1
 		}
 		z := q - p