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 | |
6 | approximation method: |
7 | |
8 | (x - 0.5) S(x) |
9 | Gamma(x) = (x + g - 0.5) * ---------------- |
10 | exp(x + g - 0.5) |
11 | |
12 | with |
13 | a1 a2 a3 aN |
14 | S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ] |
15 | x + 1 x + 2 x + 3 x + N |
16 | |
17 | with a0, a1, a2, a3,.. aN constants which depend on g. |
18 | |
19 | for x < 0 the following reflection formula is used: |
20 | |
21 | Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) |
22 | |
23 | most ideas and constants are from boost and python |
24 | */ |
25 | extern crate core; |
26 | use super::{exp, floor, k_cos, k_sin, pow}; |
27 | |
28 | const PI: f64 = 3.141592653589793238462643383279502884; |
29 | |
30 | /* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */ |
31 | fn sinpi(mut x: f64) -> f64 { |
32 | let mut n: isize; |
33 | |
34 | /* argument reduction: x = |x| mod 2 */ |
35 | /* spurious inexact when x is odd int */ |
36 | x = x * 0.5; |
37 | x = 2.0 * (x - floor(x)); |
38 | |
39 | /* reduce x into [-.25,.25] */ |
40 | n = (4.0 * x) as isize; |
41 | n = div!(n + 1, 2); |
42 | x -= (n as f64) * 0.5; |
43 | |
44 | x *= PI; |
45 | match n { |
46 | 1 => k_cos(x, y:0.0), |
47 | 2 => k_sin(-x, y:0.0, iy:0), |
48 | 3 => -k_cos(x, y:0.0), |
49 | 0 | _ => k_sin(x, y:0.0, iy:0), |
50 | } |
51 | } |
52 | |
53 | const N: usize = 12; |
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, |
70 | ]; |
71 | const SDEN: [f64; N + 1] = [ |
72 | 0.0, |
73 | 39916800.0, |
74 | 120543840.0, |
75 | 150917976.0, |
76 | 105258076.0, |
77 | 45995730.0, |
78 | 13339535.0, |
79 | 2637558.0, |
80 | 357423.0, |
81 | 32670.0, |
82 | 1925.0, |
83 | 66.0, |
84 | 1.0, |
85 | ]; |
86 | /* n! for small integer n */ |
87 | const FACT: [f64; 23] = [ |
88 | 1.0, |
89 | 1.0, |
90 | 2.0, |
91 | 6.0, |
92 | 24.0, |
93 | 120.0, |
94 | 720.0, |
95 | 5040.0, |
96 | 40320.0, |
97 | 362880.0, |
98 | 3628800.0, |
99 | 39916800.0, |
100 | 479001600.0, |
101 | 6227020800.0, |
102 | 87178291200.0, |
103 | 1307674368000.0, |
104 | 20922789888000.0, |
105 | 355687428096000.0, |
106 | 6402373705728000.0, |
107 | 121645100408832000.0, |
108 | 2432902008176640000.0, |
109 | 51090942171709440000.0, |
110 | 1124000727777607680000.0, |
111 | ]; |
112 | |
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; |
117 | |
118 | /* to avoid overflow handle large x differently */ |
119 | if x < 8.0 { |
120 | for i: usize in (0..=N).rev() { |
121 | num = num * x + i!(SNUM, i); |
122 | den = den * x + i!(SDEN, i); |
123 | } |
124 | } else { |
125 | for i: usize in 0..=N { |
126 | num = num / x + i!(SNUM, i); |
127 | den = den / x + i!(SDEN, i); |
128 | } |
129 | } |
130 | return num / den; |
131 | } |
132 | |
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(); |
136 | let absx: f64; |
137 | let mut y: f64; |
138 | let mut dy: f64; |
139 | let mut z: f64; |
140 | let mut r: f64; |
141 | let ix: u32 = ((u >> 32) as u32) & 0x7fffffff; |
142 | let sign: bool = (u >> 63) != 0; |
143 | |
144 | /* special cases */ |
145 | if ix >= 0x7ff00000 { |
146 | /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */ |
147 | return x + core::f64::INFINITY; |
148 | } |
149 | if ix < ((0x3ff - 54) << 20) { |
150 | /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */ |
151 | return 1.0 / x; |
152 | } |
153 | |
154 | /* integer arguments */ |
155 | /* raise inexact when non-integer */ |
156 | if x == floor(x) { |
157 | if sign { |
158 | return 0.0 / 0.0; |
159 | } |
160 | if x <= FACT.len() as f64 { |
161 | return i!(FACT, (x as usize) - 1); |
162 | } |
163 | } |
164 | |
165 | /* x >= 172: tgamma(x)=inf with overflow */ |
166 | /* x =< -184: tgamma(x)=+-0 with underflow */ |
167 | if ix >= 0x40670000 { |
168 | /* |x| >= 184 */ |
169 | if sign { |
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) { |
173 | return 0.0; |
174 | } else { |
175 | return -0.0; |
176 | } |
177 | } |
178 | let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 == 2^1023 |
179 | x *= x1p1023; |
180 | return x; |
181 | } |
182 | |
183 | absx = if sign { -x } else { x }; |
184 | |
185 | /* handle the error of x + g - 0.5 */ |
186 | y = absx + GMHALF; |
187 | if absx > GMHALF { |
188 | dy = y - absx; |
189 | dy -= GMHALF; |
190 | } else { |
191 | dy = y - GMHALF; |
192 | dy -= absx; |
193 | } |
194 | |
195 | z = absx - 0.5; |
196 | r = s(absx) * exp(-y); |
197 | if x < 0.0 { |
198 | /* reflection formula for negative x */ |
199 | /* sinpi(absx) is not 0, integers are already handled */ |
200 | r = -PI / (sinpi(absx) * absx * r); |
201 | dy = -dy; |
202 | z = -z; |
203 | } |
204 | r += dy * (GMHALF + 0.5) * r / y; |
205 | z = pow(y, 0.5 * z); |
206 | y = r * z * z; |
207 | return y; |
208 | } |
209 | |