-
Notifications
You must be signed in to change notification settings - Fork 5
/
gamma.cpp
167 lines (141 loc) · 4.71 KB
/
gamma.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.
#include <cmath>
#include <cfloat>
#include <sstream>
#include <stdexcept>
#include "gamma.h"
// Note that the functions Gamma and LogGamma are mutually dependent.
double Gamma(double x)
{
if (x <= 0.0)
{
std::stringstream os;
os << "Invalid Gamma() input argument " << x << ". Argument must be positive.";
throw std::invalid_argument(os.str());
}
// Split the function domain into three intervals:
// (0, 0.001), [0.001, 12), and (12, infinity)
///////////////////////////////////////////////////////////////////////////
// First interval: (0, 0.001)
//
// For small x, 1/Gamma(x) has power series x + gamma x^2 - ...
// So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
// The relative error over this interval is less than 6e-7.
const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
if (x < 0.001)
return 1.0/(x*(1.0 + gamma*x));
///////////////////////////////////////////////////////////////////////////
// Second interval: [0.001, 12)
if (x < 12.0)
{
// The algorithm directly approximates gamma over (1,2) and uses
// reduction identities to reduce other arguments to this interval.
double y = x;
int n = 0;
bool arg_was_less_than_one = (y < 1.0);
// Add or subtract integers as necessary to bring y into (1,2)
// Will correct for this below
if (arg_was_less_than_one)
y += 1.0;
else
{
n = static_cast<int> (floor(y)) - 1; // will use n later
y -= n;
}
// numerator coefficients for approximation over the interval (1,2)
static const double p[] =
{
-1.71618513886549492533811E+0,
2.47656508055759199108314E+1,
-3.79804256470945635097577E+2,
6.29331155312818442661052E+2,
8.66966202790413211295064E+2,
-3.14512729688483675254357E+4,
-3.61444134186911729807069E+4,
6.64561438202405440627855E+4
};
// denominator coefficients for approximation over the interval (1,2)
static const double q[] =
{
-3.08402300119738975254353E+1,
3.15350626979604161529144E+2,
-1.01515636749021914166146E+3,
-3.10777167157231109440444E+3,
2.25381184209801510330112E+4,
4.75584627752788110767815E+3,
-1.34659959864969306392456E+5,
-1.15132259675553483497211E+5
};
double num = 0.0;
double den = 1.0;
int i;
double z = y - 1;
for (i=0; i < 8; i++)
{
num = (num+p[i])*z;
den = den*z+q[i];
}
double result = num/den + 1.0;
// Apply correction if argument was not initially in (1,2)
if (arg_was_less_than_one)
{
// Use identity gamma(z) = gamma(z+1)/z
// The variable "result" now holds gamma of the original y + 1
// Thus we use y-1 to get back the orginal y.
result /= (y-1.0);
}
else
{
// Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
for (i=0; i < n; i++)
result *= y++;
}
return result;
}
///////////////////////////////////////////////////////////////////////////
// Third interval: [12, infinity)
if (x > 171.624)
{
// Correct answer too large to display. Force +infinity.
double temp = DBL_MAX;
return temp*2.0;
}
return exp(LogGamma(x));
}
double LogGamma(double x)
{
if (x <= 0.0)
{
std::stringstream os;
os << "Invalid LogGamma() input argument " << x << ". Argument must be positive.";
throw std::invalid_argument(os.str());
}
if (x < 12.0)
return log(fabs(Gamma(x)));
// Abramowitz and Stegun 6.1.41
// Asymptotic series should be good to at least 11 or 12 figures
// For error analysis, see Whittiker and Watson
// A Course in Modern Analysis (1927), page 252
static const double c[8] =
{
1.0/12.0,
-1.0/360.0,
1.0/1260.0,
-1.0/1680.0,
1.0/1188.0,
-691.0/360360.0,
1.0/156.0,
-3617.0/122400.0
};
double z = 1.0/(x*x);
double sum = c[7];
for (int i=6; i >= 0; i--)
{
sum *= z;
sum += c[i];
}
double series = sum/x;
static const double halfLogTwoPi = 0.91893853320467274178032973640562;
double logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;
return logGamma;
}