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//   moshier@na-net.ornl.gov
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 polynomial is valid for 33 <= x <= 172.
95func stirling(x float64) float64 {
96	const (
97		SqrtTwoPi   = 2.506628274631000502417
98		MaxStirling = 143.01608
99	)
100	w := 1 / x
101	w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
102	y := Exp(x)
103	if x > MaxStirling { // avoid Pow() overflow
104		v := Pow(x, 0.5*x-0.25)
105		y = v * (v / y)
106	} else {
107		y = Pow(x, x-0.5) / y
108	}
109	y = SqrtTwoPi * y * w
110	return y
111}
112
113// Gamma returns the Gamma function of x.
114//
115// Special cases are:
116//	Gamma(+Inf) = +Inf
117//	Gamma(+0) = +Inf
118//	Gamma(-0) = -Inf
119//	Gamma(x) = NaN for integer x < 0
120//	Gamma(-Inf) = NaN
121//	Gamma(NaN) = NaN
122func Gamma(x float64) float64 {
123	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
124	// special cases
125	switch {
126	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
127		return NaN()
128	case x == 0:
129		if Signbit(x) {
130			return Inf(-1)
131		}
132		return Inf(1)
133	case x < -170.5674972726612 || x > 171.61447887182298:
134		return Inf(1)
135	}
136	q := Abs(x)
137	p := Floor(q)
138	if q > 33 {
139		if x >= 0 {
140			return stirling(x)
141		}
142		signgam := 1
143		if ip := int(p); ip&1 == 0 {
144			signgam = -1
145		}
146		z := q - p
147		if z > 0.5 {
148			p = p + 1
149			z = q - p
150		}
151		z = q * Sin(Pi*z)
152		if z == 0 {
153			return Inf(signgam)
154		}
155		z = Pi / (Abs(z) * stirling(q))
156		return float64(signgam) * z
157	}
158
159	// Reduce argument
160	z := 1.0
161	for x >= 3 {
162		x = x - 1
163		z = z * x
164	}
165	for x < 0 {
166		if x > -1e-09 {
167			goto small
168		}
169		z = z / x
170		x = x + 1
171	}
172	for x < 2 {
173		if x < 1e-09 {
174			goto small
175		}
176		z = z / x
177		x = x + 1
178	}
179
180	if x == 2 {
181		return z
182	}
183
184	x = x - 2
185	p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
186	q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
187	return z * p / q
188
189small:
190	if x == 0 {
191		return Inf(1)
192	}
193	return z / ((1 + Euler*x) * x)
194}
195
196func isNegInt(x float64) bool {
197	if x < 0 {
198		_, xf := Modf(x)
199		return xf == 0
200	}
201	return false
202}
203