32
32
*
33
33
****************************************************************************/
34
34
35
- /* "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
36
- * "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
37
- * "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
35
+ /* "A Precision Approximation of the Gamma Function"
36
+ * - Cornelius Lanczos (1964)
37
+ * "Lanczos Implementation of the Gamma Function"
38
+ * - Paul Godfrey (2001)
39
+ * "An Analysis of the Lanczos Gamma Approximation"
40
+ * - Glendon Ralph Pugh (2004)
38
41
*
39
42
* Approximation method:
40
43
*
@@ -133,9 +136,10 @@ static const double g_sden[N + 1] =
133
136
static const double g_fact [] =
134
137
{
135
138
1 , 1 , 2 , 6 , 24 , 120 , 720 , 5040.0 , 40320.0 , 362880.0 , 3628800.0 , 39916800.0 ,
136
- 479001600.0 , 6227020800.0 , 87178291200.0 , 1307674368000.0 , 20922789888000.0 ,
137
- 355687428096000.0 , 6402373705728000.0 , 121645100408832000.0 ,
138
- 2432902008176640000.0 , 51090942171709440000.0 , 1124000727777607680000.0 ,
139
+ 479001600.0 , 6227020800.0 , 87178291200.0 , 1307674368000.0 ,
140
+ 20922789888000.0 , 355687428096000.0 , 6402373705728000.0 ,
141
+ 121645100408832000.0 , 2432902008176640000.0 , 51090942171709440000.0 ,
142
+ 1124000727777607680000.0 ,
139
143
};
140
144
141
145
/* S(x) rational function for positive x */
@@ -151,6 +155,7 @@ static double sinpi(double x)
151
155
int n ;
152
156
153
157
/* argument reduction: x = |x| mod 2 */
158
+
154
159
/* spurious inexact when x is odd int */
155
160
156
161
x = x * 0.5 ;
@@ -205,7 +210,7 @@ static double s(double x)
205
210
}
206
211
}
207
212
208
- return num / den ;
213
+ return num / den ;
209
214
}
210
215
211
216
/****************************************************************************
@@ -219,6 +224,7 @@ double tgamma(double x)
219
224
double f ;
220
225
uint64_t i ;
221
226
} u ;
227
+
222
228
u .f = x ;
223
229
224
230
double absx ;
@@ -241,17 +247,19 @@ double tgamma(double x)
241
247
if (ix < (0x3ff - 54 ) << 20 )
242
248
{
243
249
/* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
250
+
244
251
return 1 / x ;
245
252
}
246
253
247
254
/* integer arguments */
255
+
248
256
/* raise inexact when non-integer */
249
257
250
258
if (x == floor (x ))
251
259
{
252
260
if (sign )
253
261
{
254
- return 0 / 0.0 ;
262
+ return NAN ;
255
263
}
256
264
257
265
if (x <= sizeof g_fact / sizeof * g_fact )
@@ -261,6 +269,7 @@ double tgamma(double x)
261
269
}
262
270
263
271
/* x >= 172: tgamma(x)=inf with overflow */
272
+
264
273
/* x =< -184: tgamma(x)=+-0 with underflow */
265
274
266
275
if (ix >= 0x40670000 )
@@ -269,11 +278,13 @@ double tgamma(double x)
269
278
270
279
if (sign )
271
280
{
272
- FORCE_EVAL ((float )(0x1p-126 / x ));
281
+ FORCE_EVAL ((float )(ldexp (1.0 , -126 ) / x ));
282
+
273
283
if (floor (x ) * 0.5 == floor (x * 0.5 ))
274
284
{
275
285
return 0 ;
276
286
}
287
+
277
288
return -0.0 ;
278
289
}
279
290
@@ -302,6 +313,7 @@ double tgamma(double x)
302
313
if (x < 0 )
303
314
{
304
315
/* reflection formula for negative x */
316
+
305
317
/* sinpi(absx) is not 0, integers are already handled */
306
318
307
319
r = - pi / (sinpi (absx ) * absx * r );
0 commit comments