diff options
Diffstat (limited to 'src/math/gamma.go')
| -rw-r--r-- | src/math/gamma.go | 222 |
1 files changed, 222 insertions, 0 deletions
diff --git a/src/math/gamma.go b/src/math/gamma.go new file mode 100644 index 0000000..86c6723 --- /dev/null +++ b/src/math/gamma.go @@ -0,0 +1,222 @@ +// 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 +} |