mirror of
https://github.com/golang/go
synced 2024-11-21 22:24:40 -07:00
math: add Gamma function
R=rsc CC=golang-dev https://golang.org/cl/649041
This commit is contained in:
Charles L. Dorian
Russ Cox
parent
64f33880e5
commit
26f0c83eb8
@ -53,6 +53,7 @@ ALLGOFILES=\
|
||||
floor.go\
|
||||
fmod.go\
|
||||
frexp.go\
|
||||
gamma.go\
|
||||
hypot.go\
|
||||
hypot_port.go\
|
||||
logb.go\
|
||||
|
||||
@ -286,6 +286,18 @@ var frexp = []fi{
|
||||
fi{9.1265404584042750000e-01, 1},
|
||||
fi{-5.4287029803597508250e-01, 4},
|
||||
}
|
||||
var gamma = []float64{
|
||||
2.3254348370739963835386613898e+01,
|
||||
2.991153837155317076427529816e+03,
|
||||
-4.561154336726758060575129109e+00,
|
||||
7.719403468842639065959210984e-01,
|
||||
1.6111876618855418534325755566e+05,
|
||||
1.8706575145216421164173224946e+00,
|
||||
3.4082787447257502836734201635e+01,
|
||||
1.579733951448952054898583387e+00,
|
||||
9.3834586598354592860187267089e-01,
|
||||
-2.093995902923148389186189429e-05,
|
||||
}
|
||||
var lgamma = []fi{
|
||||
fi{3.146492141244545774319734e+00, 1},
|
||||
fi{8.003414490659126375852113e+00, 1},
|
||||
@ -736,6 +748,21 @@ var frexpSC = []fi{
|
||||
fi{NaN(), 0},
|
||||
}
|
||||
|
||||
var vfgammaSC = []float64{
|
||||
Inf(-1),
|
||||
-3,
|
||||
0,
|
||||
Inf(1),
|
||||
NaN(),
|
||||
}
|
||||
var gammaSC = []float64{
|
||||
Inf(-1),
|
||||
Inf(1),
|
||||
Inf(1),
|
||||
Inf(1),
|
||||
NaN(),
|
||||
}
|
||||
|
||||
var vfhypotSC = [][2]float64{
|
||||
[2]float64{Inf(-1), Inf(-1)},
|
||||
[2]float64{Inf(-1), 0},
|
||||
@ -1278,6 +1305,19 @@ func TestFrexp(t *testing.T) {
|
||||
}
|
||||
}
|
||||
|
||||
func TestGamma(t *testing.T) {
|
||||
for i := 0; i < len(vf); i++ {
|
||||
if f := Gamma(vf[i]); !close(gamma[i], f) {
|
||||
t.Errorf("Gamma(%g) = %g, want %g\n", 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\n", vfgammaSC[i], f, gammaSC[i])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
func TestHypot(t *testing.T) {
|
||||
for i := 0; i < len(vf); i++ {
|
||||
a := Fabs(1e200 * tanh[i] * Sqrt(2))
|
||||
@ -1748,6 +1788,12 @@ func BenchmarkFrexp(b *testing.B) {
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkGamma(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
Gamma(2.5)
|
||||
}
|
||||
}
|
||||
|
||||
func BenchmarkHypot(b *testing.B) {
|
||||
for i := 0; i < b.N; i++ {
|
||||
Hypot(3, 4)
|
||||
|
||||
188
src/pkg/math/gamma.go
Normal file
188
src/pkg/math/gamma.go
Normal file
@ -0,0 +1,188 @@
|
||||
// Copyright 2010 The Go Authors. All rights reserved.
|
||||
// Use of this source code is governed by a BSD-style
|
||||
// license that can be found in the LICENSE file.
|
||||
|
||||
package math
|
||||
|
||||
// The original C code, the long comment, and the constants
|
||||
// below are from http://netlib.sandia.gov/cephes/cprob/gamma.c.
|
||||
// The go code is a simplified version of the original C.
|
||||
//
|
||||
// tgamma.c
|
||||
//
|
||||
// Gamma function
|
||||
//
|
||||
// SYNOPSIS:
|
||||
//
|
||||
// double x, y, tgamma();
|
||||
// extern int signgam;
|
||||
//
|
||||
// y = tgamma( x );
|
||||
//
|
||||
// DESCRIPTION:
|
||||
//
|
||||
// Returns gamma function of the argument. The result is
|
||||
// correctly signed, and the sign (+1 or -1) is also
|
||||
// returned in a global (extern) variable named signgam.
|
||||
// This variable is also filled in by the logarithmic gamma
|
||||
// function lgamma().
|
||||
//
|
||||
// Arguments |x| <= 34 are reduced by recurrence and the function
|
||||
// approximated by a rational function of degree 6/7 in the
|
||||
// interval (2,3). Large arguments are handled by Stirling's
|
||||
// formula. Large negative arguments are made positive using
|
||||
// a reflection formula.
|
||||
//
|
||||
// ACCURACY:
|
||||
//
|
||||
// Relative error:
|
||||
// arithmetic domain # trials peak rms
|
||||
// DEC -34, 34 10000 1.3e-16 2.5e-17
|
||||
// IEEE -170,-33 20000 2.3e-15 3.3e-16
|
||||
// IEEE -33, 33 20000 9.4e-16 2.2e-16
|
||||
// IEEE 33, 171.6 20000 2.3e-15 3.2e-16
|
||||
//
|
||||
// Error for arguments outside the test range will be larger
|
||||
// owing to error amplification by the exponential function.
|
||||
//
|
||||
// Cephes Math Library Release 2.8: June, 2000
|
||||
// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
|
||||
//
|
||||
// The readme file at http://netlib.sandia.gov/cephes/ says:
|
||||
// Some software in this archive may be from the book _Methods and
|
||||
// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
|
||||
// International, 1989) or from the Cephes Mathematical Library, a
|
||||
// commercial product. In either event, it is copyrighted by the author.
|
||||
// What you see here may be used freely but it comes with no support or
|
||||
// guarantee.
|
||||
//
|
||||
// The two known misprints in the book are repaired here in the
|
||||
// source listings for the gamma function and the incomplete beta
|
||||
// integral.
|
||||
//
|
||||
// Stephen L. Moshier
|
||||
// moshier@na-net.ornl.gov
|
||||
|
||||
var _P = []float64{
|
||||
1.60119522476751861407e-04,
|
||||
1.19135147006586384913e-03,
|
||||
1.04213797561761569935e-02,
|
||||
4.76367800457137231464e-02,
|
||||
2.07448227648435975150e-01,
|
||||
4.94214826801497100753e-01,
|
||||
9.99999999999999996796e-01,
|
||||
}
|
||||
var _Q = []float64{
|
||||
-2.31581873324120129819e-05,
|
||||
5.39605580493303397842e-04,
|
||||
-4.45641913851797240494e-03,
|
||||
1.18139785222060435552e-02,
|
||||
3.58236398605498653373e-02,
|
||||
-2.34591795718243348568e-01,
|
||||
7.14304917030273074085e-02,
|
||||
1.00000000000000000320e+00,
|
||||
}
|
||||
var _S = []float64{
|
||||
7.87311395793093628397e-04,
|
||||
-2.29549961613378126380e-04,
|
||||
-2.68132617805781232825e-03,
|
||||
3.47222221605458667310e-03,
|
||||
8.33333333333482257126e-02,
|
||||
}
|
||||
|
||||
// Gamma function computed by Stirling's formula.
|
||||
// The polynomial is valid for 33 <= x <= 172.
|
||||
func stirling(x float64) float64 {
|
||||
const (
|
||||
SqrtTwoPi = 2.506628274631000502417
|
||||
MaxStirling = 143.01608
|
||||
)
|
||||
w := 1 / x
|
||||
w = 1 + w*((((_S[0]*w+_S[1])*w+_S[2])*w+_S[3])*w+_S[4])
|
||||
y := Exp(x)
|
||||
if x > MaxStirling { // avoid Pow() overflow
|
||||
v := Pow(x, 0.5*x-0.25)
|
||||
y = v * (v / y)
|
||||
} else {
|
||||
y = Pow(x, x-0.5) / y
|
||||
}
|
||||
y = SqrtTwoPi * y * w
|
||||
return y
|
||||
}
|
||||
|
||||
// Gamma(x) returns the Gamma function of x.
|
||||
//
|
||||
// Special cases are:
|
||||
// Gamma(Inf) = Inf
|
||||
// Gamma(-Inf) = -Inf
|
||||
// Gamma(NaN) = NaN
|
||||
// Large values overflow to +Inf.
|
||||
// Negative integer values equal ±Inf.
|
||||
func Gamma(x float64) float64 {
|
||||
const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
|
||||
// special cases
|
||||
switch {
|
||||
case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x):
|
||||
return x
|
||||
case x < -170.5674972726612 || x > 171.61447887182298:
|
||||
return Inf(1)
|
||||
}
|
||||
q := Fabs(x)
|
||||
p := Floor(q)
|
||||
if q > 33 {
|
||||
if x >= 0 {
|
||||
return stirling(x)
|
||||
}
|
||||
signgam := 1
|
||||
if ip := int(p); ip&1 == 0 {
|
||||
signgam = -1
|
||||
}
|
||||
z := q - p
|
||||
if z > 0.5 {
|
||||
p = p + 1
|
||||
z = q - p
|
||||
}
|
||||
z = q * Sin(Pi*z)
|
||||
if z == 0 {
|
||||
return Inf(signgam)
|
||||
}
|
||||
z = Pi / (Fabs(z) * stirling(q))
|
||||
return float64(signgam) * z
|
||||
}
|
||||
|
||||
// Reduce argument
|
||||
z := float64(1)
|
||||
for x >= 3 {
|
||||
x = x - 1
|
||||
z = z * x
|
||||
}
|
||||
for x < 0 {
|
||||
if x > -1e-09 {
|
||||
goto small
|
||||
}
|
||||
z = z / x
|
||||
x = x + 1
|
||||
}
|
||||
for x < 2 {
|
||||
if x < 1e-09 {
|
||||
goto small
|
||||
}
|
||||
z = z / x
|
||||
x = x + 1
|
||||
}
|
||||
|
||||
if x == 2 {
|
||||
return z
|
||||
}
|
||||
|
||||
x = x - 2
|
||||
p = (((((x*_P[0]+_P[1])*x+_P[2])*x+_P[3])*x+_P[4])*x+_P[5])*x + _P[6]
|
||||
q = ((((((x*_Q[0]+_Q[1])*x+_Q[2])*x+_Q[3])*x+_Q[4])*x+_Q[5])*x+_Q[6])*x + _Q[7]
|
||||
return z * p / q
|
||||
|
||||
small:
|
||||
if x == 0 {
|
||||
return Inf(1)
|
||||
}
|
||||
return z / ((1 + Euler*x) * x)
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user