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)
9 Gamma(x) = (x + g - 0.5) * ----------------
14 S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
15 x + 1 x + 2 x + 3 x + N
17 with a0, a1, a2, a3,.. aN constants which depend on g.
19 for x < 0 the following reflection formula is used:
21 Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
23 most ideas and constants are from boost and python
26 use super::{exp, floor, k_cos, k_sin, pow}
;
28 const PI
: f64 = 3.141592653589793238462643383279502884;
30 /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
31 fn sinpi(mut x
: f64) -> f64 {
34 /* argument reduction: x = |x| mod 2 */
35 /* spurious inexact when x is odd int */
37 x
= 2.0 * (x
- floor(x
));
39 /* reduce x into [-.25,.25] */
40 n
= (4.0 * x
) as isize;
42 x
-= (n
as f64) * 0.5;
47 2 => k_sin(-x
, 0.0, 0),
49 0 | _
=> k_sin(x
, 0.0, 0),
54 //static const double g = 6.024680040776729583740234375;
55 const GMHALF
: f64 = 5.524680040776729583740234375;
56 const SNUM
: [f64; N
+ 1] = [
57 23531376880.410759688572007674451636754734846804940,
58 42919803642.649098768957899047001988850926355848959,
59 35711959237.355668049440185451547166705960488635843,
60 17921034426.037209699919755754458931112671403265390,
61 6039542586.3520280050642916443072979210699388420708,
62 1439720407.3117216736632230727949123939715485786772,
63 248874557.86205415651146038641322942321632125127801,
64 31426415.585400194380614231628318205362874684987640,
65 2876370.6289353724412254090516208496135991145378768,
66 186056.26539522349504029498971604569928220784236328,
67 8071.6720023658162106380029022722506138218516325024,
68 210.82427775157934587250973392071336271166969580291,
69 2.5066282746310002701649081771338373386264310793408,
71 const SDEN
: [f64; N
+ 1] = [
86 /* n! for small integer n */
87 const FACT
: [f64; 23] = [
107 121645100408832000.0,
108 2432902008176640000.0,
109 51090942171709440000.0,
110 1124000727777607680000.0,
113 /* S(x) rational function for positive x */
114 fn s(x
: f64) -> f64 {
115 let mut num
: f64 = 0.0;
116 let mut den
: f64 = 0.0;
118 /* to avoid overflow handle large x differently */
120 for i
in (0..=N
).rev() {
121 num
= num
* x
+ i
!(SNUM
, i
);
122 den
= den
* x
+ i
!(SDEN
, i
);
126 num
= num
/ x
+ i
!(SNUM
, i
);
127 den
= den
/ x
+ i
!(SDEN
, i
);
133 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
134 pub fn tgamma(mut x
: f64) -> f64 {
135 let u
: u64 = x
.to_bits();
141 let ix
: u32 = ((u
>> 32) as u32) & 0x7fffffff;
142 let sign
: bool
= (u
>> 63) != 0;
145 if ix
>= 0x7ff00000 {
146 /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
147 return x
+ core
::f64::INFINITY
;
149 if ix
< ((0x3ff - 54) << 20) {
150 /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
154 /* integer arguments */
155 /* raise inexact when non-integer */
160 if x
<= FACT
.len() as f64 {
161 return i
!(FACT
, (x
as usize) - 1);
165 /* x >= 172: tgamma(x)=inf with overflow */
166 /* x =< -184: tgamma(x)=+-0 with underflow */
167 if ix
>= 0x40670000 {
170 let x1p_126
= f64::from_bits(0x3810000000000000); // 0x1p-126 == 2^-126
171 force_eval
!((x1p_126
/ x
) as f32);
172 if floor(x
) * 0.5 == floor(x
* 0.5) {
178 let x1p1023
= f64::from_bits(0x7fe0000000000000); // 0x1p1023 == 2^1023
183 absx
= if sign { -x }
else { x }
;
185 /* handle the error of x + g - 0.5 */
196 r
= s(absx
) * exp(-y
);
198 /* reflection formula for negative x */
199 /* sinpi(absx) is not 0, integers are already handled */
200 r
= -PI
/ (sinpi(absx
) * absx
* r
);
204 r
+= dy
* (GMHALF
+ 0.5) * r
/ y
;