diff options
Diffstat (limited to 'src/pkg/math/gamma.go')
| -rw-r--r-- | src/pkg/math/gamma.go | 188 |
1 files changed, 0 insertions, 188 deletions
diff --git a/src/pkg/math/gamma.go b/src/pkg/math/gamma.go deleted file mode 100644 index 73ca0e53a..000000000 --- a/src/pkg/math/gamma.go +++ /dev/null @@ -1,188 +0,0 @@ -// 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 := 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*_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) -} |