qbit/go
1
0
Fork 0
mirror of https://github.com/golang/go synced 2024-11-05 17:36:15 -07:00
Code Releases Activity

math: fix Gamma(-171.5) on all platforms

Browse Source 
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>
This commit is contained in:
Russ Cox 2016-10-06 10:39:50 -04:00
parent 20c48c9557
commit 4f3a641e6e
2 changed files with 33 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/math/all_test.go
View File

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

35
src/math/gamma.go
View File

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