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.
4 5 package math
6 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.
10 //
11 // tgamma.c
12 //
13 // Gamma function
14 //
15 // SYNOPSIS:
16 //
17 // double x, y, tgamma();
18 // extern int signgam;
19 //
20 // y = tgamma( x );
21 //
22 // DESCRIPTION:
23 //
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
28 // function lgamma().
29 //
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.
35 //
36 // ACCURACY:
37 //
38 // Relative error:
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
44 //
45 // Error for arguments outside the test range will be larger
46 // owing to error amplification by the exponential function.
47 //
48 // Cephes Math Library Release 2.8: June, 2000
49 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
50 //
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
57 // guarantee.
58 //
59 // The two known misprints in the book are repaired here in the
60 // source listings for the gamma function and the incomplete beta
61 // integral.
62 //
63 // Stephen L. Moshier
64 // moshier@na-net.ornl.gov
65 66 var _gamP = [...]float64{
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,
74 }
75 var _gamQ = [...]float64{
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,
84 }
85 var _gamS = [...]float64{
86 7.87311395793093628397e-04,
87 -2.29549961613378126380e-04,
88 -2.68132617805781232825e-03,
89 3.47222221605458667310e-03,
90 8.33333333333482257126e-02,
91 }
92 93 // Gamma function computed by Stirling's formula.
94 // The pair of results must be multiplied together to get the actual answer.
95 // The multiplication is left to the caller so that, if careful, the caller can avoid
96 // infinity for 172 <= x <= 180.
97 // The polynomial is valid for 33 <= x <= 172; larger values are only used
98 // in reciprocal and produce denormalized floats. The lower precision there
99 // masks any imprecision in the polynomial.
100 func stirling(x float64) (float64, float64) {
101 if x > 200 {
102 return Inf(1), 1
103 }
104 const (
105 SqrtTwoPi = 2.506628274631000502417
106 MaxStirling = 143.01608
107 )
108 w := 1 / x
109 w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
110 y1 := Exp(x)
111 y2 := 1.0
112 if x > MaxStirling { // avoid Pow() overflow
113 v := Pow(x, 0.5*x-0.25)
114 y1, y2 = v, v/y1
115 } else {
116 y1 = Pow(x, x-0.5) / y1
117 }
118 return y1, SqrtTwoPi * w * y2
119 }
120 121 // Gamma returns the Gamma function of x.
122 //
123 // Special cases are:
124 //
125 // Gamma(+Inf) = +Inf
126 // Gamma(+0) = +Inf
127 // Gamma(-0) = -Inf
128 // Gamma(x) = NaN for integer x < 0
129 // Gamma(-Inf) = NaN
130 // Gamma(NaN) = NaN
131 func Gamma(x float64) float64 {
132 const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
133 // special cases
134 switch {
135 case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
136 return NaN()
137 case IsInf(x, 1):
138 return Inf(1)
139 case x == 0:
140 if Signbit(x) {
141 return Inf(-1)
142 }
143 return Inf(1)
144 }
145 q := Abs(x)
146 p := Floor(q)
147 if q > 33 {
148 if x >= 0 {
149 y1, y2 := stirling(x)
150 return y1 * y2
151 }
152 // Note: x is negative but (checked above) not a negative integer,
153 // so x must be small enough to be in range for conversion to int64.
154 // If |x| were >= 2⁶³ it would have to be an integer.
155 signgam := 1
156 if ip := int64(p); ip&1 == 0 {
157 signgam = -1
158 }
159 z := q - p
160 if z > 0.5 {
161 p = p + 1
162 z = q - p
163 }
164 z = q * Sin(Pi*z)
165 if z == 0 {
166 return Inf(signgam)
167 }
168 sq1, sq2 := stirling(q)
169 absz := Abs(z)
170 d := absz * sq1 * sq2
171 if IsInf(d, 0) {
172 z = Pi / absz / sq1 / sq2
173 } else {
174 z = Pi / d
175 }
176 return float64(signgam) * z
177 }
178 179 // Reduce argument
180 z := 1.0
181 for x >= 3 {
182 x = x - 1
183 z = z * x
184 }
185 for x < 0 {
186 if x > -1e-09 {
187 goto small
188 }
189 z = z / x
190 x = x + 1
191 }
192 for x < 2 {
193 if x < 1e-09 {
194 goto small
195 }
196 z = z / x
197 x = x + 1
198 }
199 200 if x == 2 {
201 return z
202 }
203 204 x = x - 2
205 p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
206 q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
207 return z * p / q
208 209 small:
210 if x == 0 {
211 return Inf(1)
212 }
213 return z / ((1 + Euler*x) * x)
214 }
215 216 func isNegInt(x float64) bool {
217 if x < 0 {
218 _, xf := Modf(x)
219 return xf == 0
220 }
221 return false
222 }
223