Black Lives Matter. Support the Equal Justice Initiative.

# Source file src/math/gamma.go

## Documentation: math

```     1  // Copyright 2010 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
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  //	Gamma(+Inf) = +Inf
125  //	Gamma(+0) = +Inf
126  //	Gamma(-0) = -Inf
127  //	Gamma(x) = NaN for integer x < 0
128  //	Gamma(-Inf) = NaN
129  //	Gamma(NaN) = NaN
130  func Gamma(x float64) float64 {
131  	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
132  	// special cases
133  	switch {
134  	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
135  		return NaN()
136  	case IsInf(x, 1):
137  		return Inf(1)
138  	case x == 0:
139  		if Signbit(x) {
140  			return Inf(-1)
141  		}
142  		return Inf(1)
143  	}
144  	q := Abs(x)
145  	p := Floor(q)
146  	if q > 33 {
147  		if x >= 0 {
148  			y1, y2 := stirling(x)
149  			return y1 * y2
150  		}
151  		// Note: x is negative but (checked above) not a negative integer,
152  		// so x must be small enough to be in range for conversion to int64.
153  		// If |x| were >= 2⁶³ it would have to be an integer.
154  		signgam := 1
155  		if ip := int64(p); ip&1 == 0 {
156  			signgam = -1
157  		}
158  		z := q - p
159  		if z > 0.5 {
160  			p = p + 1
161  			z = q - p
162  		}
163  		z = q * Sin(Pi*z)
164  		if z == 0 {
165  			return Inf(signgam)
166  		}
167  		sq1, sq2 := stirling(q)
168  		absz := Abs(z)
169  		d := absz * sq1 * sq2
170  		if IsInf(d, 0) {
171  			z = Pi / absz / sq1 / sq2
172  		} else {
173  			z = Pi / d
174  		}
175  		return float64(signgam) * z
176  	}
177
178  	// Reduce argument
179  	z := 1.0
180  	for x >= 3 {
181  		x = x - 1
182  		z = z * x
183  	}
184  	for x < 0 {
185  		if x > -1e-09 {
186  			goto small
187  		}
188  		z = z / x
189  		x = x + 1
190  	}
191  	for x < 2 {
192  		if x < 1e-09 {
193  			goto small
194  		}
195  		z = z / x
196  		x = x + 1
197  	}
198
199  	if x == 2 {
200  		return z
201  	}
202
203  	x = x - 2
204  	p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
205  	q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
206  	return z * p / q
207
208  small:
209  	if x == 0 {
210  		return Inf(1)
211  	}
212  	return z / ((1 + Euler*x) * x)
213  }
214
215  func isNegInt(x float64) bool {
216  	if x < 0 {
217  		_, xf := Modf(x)
218  		return xf == 0
219  	}
220  	return false
221  }
222
```

View as plain text