summaryrefslogtreecommitdiffstats
path: root/src/math/gamma.go
diff options
context:
space:
mode:
Diffstat (limited to 'src/math/gamma.go')
-rw-r--r--src/math/gamma.go222
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
+}