blob: cc9e869496b114a2bbc6c20ba7d1010d4e0f0d68 [file] [log] [blame]
// 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 _gamP = [...]float64{
1.60119522476751861407e-04,
1.19135147006586384913e-03,
1.04213797561761569935e-02,
4.76367800457137231464e-02,
2.07448227648435975150e-01,
4.94214826801497100753e-01,
9.99999999999999996796e-01,
}
var _gamQ = [...]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 _gamS = [...]float64{
7.87311395793093628397e-04,
-2.29549961613378126380e-04,
-2.68132617805781232825e-03,
3.47222221605458667310e-03,
8.33333333333482257126e-02,
}
// Gamma function computed by Stirling's formula.
// 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
MaxStirling = 143.01608
)
w := 1 / x
w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
y1 := Exp(x)
y2 := 1.0
if x > MaxStirling { // avoid Pow() overflow
v := Pow(x, 0.5*x-0.25)
y1, y2 = v, v/y1
} else {
y1 = Pow(x, x-0.5) / y1
}
return y1, SqrtTwoPi * w * y2
}
// Gamma returns the Gamma function of x.
//
// Special cases are:
// Gamma(+Inf) = +Inf
// Gamma(+0) = +Inf
// Gamma(-0) = -Inf
// Gamma(x) = NaN for integer x < 0
// Gamma(-Inf) = NaN
// Gamma(NaN) = NaN
func Gamma(x float64) float64 {
const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
// special cases
switch {
case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
return NaN()
case IsInf(x, 1):
return Inf(1)
case x == 0:
if Signbit(x) {
return Inf(-1)
}
return Inf(1)
}
q := Abs(x)
p := Floor(q)
if q > 33 {
if x >= 0 {
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.
// If |x| were >= 2⁶³ it would have to be an integer.
signgam := 1
if ip := int64(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)
}
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
}
// Reduce argument
z := 1.0
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*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
return z * p / q
small:
if x == 0 {
return Inf(1)
}
return z / ((1 + Euler*x) * x)
}
func isNegInt(x float64) bool {
if x < 0 {
_, xf := Modf(x)
return xf == 0
}
return false
}