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
5package 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//   [email protected]
65
66var _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}
75var _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}
85var _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.
100func 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
131func 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
209small:
210	if x == 0 {
211		return Inf(1)
212	}
213	return z / ((1 + Euler*x) * x)
214}
215
216func isNegInt(x float64) bool {
217	if x < 0 {
218		_, xf := Modf(x)
219		return xf == 0
220	}
221	return false
222}
223