qbit/go
1
0
Fork 0
mirror of https://github.com/golang/go synced 2024-11-18 13:44:48 -07:00
Code Releases Activity

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

Browse Source 
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>
This commit is contained in:
Russ Cox 2016-10-03 22:13:57 -04:00
parent aab849e429
commit a39920fdbb
2 changed files with 101 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
src/math/all_test.go
View File

@ -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])
}
}
}

10
src/math/gamma.go
View File

@ -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