| // 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) |
| } |