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>

diff --git a/src/math/all_test.go b/src/math/all_test.go
index 882400527b24d2aa52898e92049ed1c3546900b6..c77dcfea6e33b4615707b54659c780fb31f89f47 100644 (file)
--- 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 841ec114a0a84fc6831d681dcf921e090c507c51..c965f275e7601341a2aa779246a543cca5585221 100644 (file)
--- 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