1 // Copyright 2010 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
7 // The original C code, the long comment, and the constants
8 // below are from http://netlib.sandia.gov/cephes/cprob/gamma.c.
9 // The go code is a simplified version of the original C.
17 // double x, y, tgamma();
18 // extern int signgam;
24 // Returns gamma function of the argument. The result is
25 // correctly signed, and the sign (+1 or -1) is also
26 // returned in a global (extern) variable named signgam.
27 // This variable is also filled in by the logarithmic gamma
30 // Arguments |x| <= 34 are reduced by recurrence and the function
31 // approximated by a rational function of degree 6/7 in the
32 // interval (2,3). Large arguments are handled by Stirling's
33 // formula. Large negative arguments are made positive using
34 // a reflection formula.
39 // arithmetic domain # trials peak rms
40 // DEC -34, 34 10000 1.3e-16 2.5e-17
41 // IEEE -170,-33 20000 2.3e-15 3.3e-16
42 // IEEE -33, 33 20000 9.4e-16 2.2e-16
43 // IEEE 33, 171.6 20000 2.3e-15 3.2e-16
45 // Error for arguments outside the test range will be larger
46 // owing to error amplification by the exponential function.
48 // Cephes Math Library Release 2.8: June, 2000
49 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
51 // The readme file at http://netlib.sandia.gov/cephes/ says:
52 // Some software in this archive may be from the book _Methods and
53 // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
54 // International, 1989) or from the Cephes Mathematical Library, a
55 // commercial product. In either event, it is copyrighted by the author.
56 // What you see here may be used freely but it comes with no support or
59 // The two known misprints in the book are repaired here in the
60 // source listings for the gamma function and the incomplete beta
64 // moshier@na-net.ornl.gov
67 1.60119522476751861407e-04,
68 1.19135147006586384913e-03,
69 1.04213797561761569935e-02,
70 4.76367800457137231464e-02,
71 2.07448227648435975150e-01,
72 4.94214826801497100753e-01,
73 9.99999999999999996796e-01,
76 -2.31581873324120129819e-05,
77 5.39605580493303397842e-04,
78 -4.45641913851797240494e-03,
79 1.18139785222060435552e-02,
80 3.58236398605498653373e-02,
81 -2.34591795718243348568e-01,
82 7.14304917030273074085e-02,
83 1.00000000000000000320e+00,
86 7.87311395793093628397e-04,
87 -2.29549961613378126380e-04,
88 -2.68132617805781232825e-03,
89 3.47222221605458667310e-03,
90 8.33333333333482257126e-02,
93 // Gamma function computed by Stirling's formula.
94 // The polynomial is valid for 33 <= x <= 172.
95 func stirling(x float64) float64 {
97 SqrtTwoPi = 2.506628274631000502417
98 MaxStirling = 143.01608
101 w = 1 + w*((((_S[0]*w+_S[1])*w+_S[2])*w+_S[3])*w+_S[4])
103 if x > MaxStirling { // avoid Pow() overflow
104 v := Pow(x, 0.5*x-0.25)
107 y = Pow(x, x-0.5) / y
109 y = SqrtTwoPi * y * w
113 // Gamma(x) returns the Gamma function of x.
115 // Special cases are:
117 // Gamma(-Inf) = -Inf
119 // Large values overflow to +Inf.
120 // Negative integer values equal ±Inf.
121 func Gamma(x float64) float64 {
122 const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
125 case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x):
127 case x < -170.5674972726612 || x > 171.61447887182298:
137 if ip := int(p); ip&1 == 0 {
149 z = Pi / (Abs(z) * stirling(q))
150 return float64(signgam) * z
179 p = (((((x*_P[0]+_P[1])*x+_P[2])*x+_P[3])*x+_P[4])*x+_P[5])*x + _P[6]
180 q = ((((((x*_Q[0]+_Q[1])*x+_Q[2])*x+_Q[3])*x+_Q[4])*x+_Q[5])*x+_Q[6])*x + _Q[7]
187 return z / ((1 + Euler*x) * x)