1/*
2"A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
3"Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
4"An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
5
6approximation method:
7
8 (x - 0.5) S(x)
9Gamma(x) = (x + g - 0.5) * ----------------
10 exp(x + g - 0.5)
11
12with
13 a1 a2 a3 aN
14S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
15 x + 1 x + 2 x + 3 x + N
16
17with a0, a1, a2, a3,.. aN constants which depend on g.
18
19for x < 0 the following reflection formula is used:
20
21Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
22
23most ideas and constants are from boost and python
24*/
25#include "libm.h"
26
27static const double pi = 3.141592653589793238462643383279502884;
28
29/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
30static double sinpi(double x)
31{
32 int n;
33
34 /* argument reduction: x = |x| mod 2 */
35 /* spurious inexact when x is odd int */
36 x = x * 0.5;
37 x = 2 * (x - floor(x));
38
39 /* reduce x into [-.25,.25] */
40 n = 4 * x;
41 n = (n+1)/2;
42 x -= n * 0.5;
43
44 x *= pi;
45 switch (n) {
46 default: /* case 4 */
47 case 0:
48 return __sin(x, 0, 0);
49 case 1:
50 return __cos(x, 0);
51 case 2:
52 return __sin(-x, 0, 0);
53 case 3:
54 return -__cos(x, 0);
55 }
56}
57
58#define N 12
59//static const double g = 6.024680040776729583740234375;
60static const double gmhalf = 5.524680040776729583740234375;
61static const double Snum[N+1] = {
62 23531376880.410759688572007674451636754734846804940,
63 42919803642.649098768957899047001988850926355848959,
64 35711959237.355668049440185451547166705960488635843,
65 17921034426.037209699919755754458931112671403265390,
66 6039542586.3520280050642916443072979210699388420708,
67 1439720407.3117216736632230727949123939715485786772,
68 248874557.86205415651146038641322942321632125127801,
69 31426415.585400194380614231628318205362874684987640,
70 2876370.6289353724412254090516208496135991145378768,
71 186056.26539522349504029498971604569928220784236328,
72 8071.6720023658162106380029022722506138218516325024,
73 210.82427775157934587250973392071336271166969580291,
74 2.5066282746310002701649081771338373386264310793408,
75};
76static const double Sden[N+1] = {
77 0, 39916800, 120543840, 150917976, 105258076, 45995730, 13339535,
78 2637558, 357423, 32670, 1925, 66, 1,
79};
80/* n! for small integer n */
81static const double fact[] = {
82 1, 1, 2, 6, 24, 120, 720, 5040.0, 40320.0, 362880.0, 3628800.0, 39916800.0,
83 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0, 20922789888000.0,
84 355687428096000.0, 6402373705728000.0, 121645100408832000.0,
85 2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
86};
87
88/* S(x) rational function for positive x */
89static double S(double x)
90{
91 double_t num = 0, den = 0;
92 int i;
93
94 /* to avoid overflow handle large x differently */
95 if (x < 8)
96 for (i = N; i >= 0; i--) {
97 num = num * x + Snum[i];
98 den = den * x + Sden[i];
99 }
100 else
101 for (i = 0; i <= N; i++) {
102 num = num / x + Snum[i];
103 den = den / x + Sden[i];
104 }
105 return num/den;
106}
107
108double tgamma(double x)
109{
110 union {double f; uint64_t i;} u = {x};
111 double absx, y;
112 double_t dy, z, r;
113 uint32_t ix = u.i>>32 & 0x7fffffff;
114 int sign = u.i>>63;
115
116 /* special cases */
117 if (ix >= 0x7ff00000)
118 /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
119 return x + INFINITY;
120 if (ix < (0x3ff-54)<<20)
121 /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
122 return 1/x;
123
124 /* integer arguments */
125 /* raise inexact when non-integer */
126 if (x == floor(x)) {
127 if (sign)
128 return 0/0.0;
129 if (x <= sizeof fact/sizeof *fact)
130 return fact[(int)x - 1];
131 }
132
133 /* x >= 172: tgamma(x)=inf with overflow */
134 /* x =< -184: tgamma(x)=+-0 with underflow */
135 if (ix >= 0x40670000) { /* |x| >= 184 */
136 if (sign) {
137 FORCE_EVAL((float)(0x1p-126/x));
138 if (floor(x) * 0.5 == floor(x * 0.5))
139 return 0;
140 return -0.0;
141 }
142 x *= 0x1p1023;
143 return x;
144 }
145
146 absx = sign ? -x : x;
147
148 /* handle the error of x + g - 0.5 */
149 y = absx + gmhalf;
150 if (absx > gmhalf) {
151 dy = y - absx;
152 dy -= gmhalf;
153 } else {
154 dy = y - gmhalf;
155 dy -= absx;
156 }
157
158 z = absx - 0.5;
159 r = S(absx) * exp(-y);
160 if (x < 0) {
161 /* reflection formula for negative x */
162 /* sinpi(absx) is not 0, integers are already handled */
163 r = -pi / (sinpi(absx) * absx * r);
164 dy = -dy;
165 z = -z;
166 }
167 r += dy * (gmhalf+0.5) * r / y;
168 z = pow(y, 0.5*z);
169 y = r * z * z;
170 return y;
171}
172
173#if 0
174double __lgamma_r(double x, int *sign)
175{
176 double r, absx;
177
178 *sign = 1;
179
180 /* special cases */
181 if (!isfinite(x))
182 /* lgamma(nan)=nan, lgamma(+-inf)=inf */
183 return x*x;
184
185 /* integer arguments */
186 if (x == floor(x) && x <= 2) {
187 /* n <= 0: lgamma(n)=inf with divbyzero */
188 /* n == 1,2: lgamma(n)=0 */
189 if (x <= 0)
190 return 1/0.0;
191 return 0;
192 }
193
194 absx = fabs(x);
195
196 /* lgamma(x) ~ -log(|x|) for tiny |x| */
197 if (absx < 0x1p-54) {
198 *sign = 1 - 2*!!signbit(x);
199 return -log(absx);
200 }
201
202 /* use tgamma for smaller |x| */
203 if (absx < 128) {
204 x = tgamma(x);
205 *sign = 1 - 2*!!signbit(x);
206 return log(fabs(x));
207 }
208
209 /* second term (log(S)-g) could be more precise here.. */
210 /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
211 r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
212 if (x < 0) {
213 /* reflection formula for negative x */
214 x = sinpi(absx);
215 *sign = 2*!!signbit(x) - 1;
216 r = log(pi/(fabs(x)*absx)) - r;
217 }
218 return r;
219}
220
221weak_alias(__lgamma_r, lgamma_r);
222#endif