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 */
sinpi(mut x: f64) -> f6431 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 = (n + 1) / 2;
42 x -= (n as f64) * 0.5;
43
44 x *= PI;
45 match n {
46 1 => k_cos(x, 0.0),
47 2 => k_sin(-x, 0.0, 0),
48 3 => -k_cos(x, 0.0),
49 0 | _ => k_sin(x, 0.0, 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 */
s(x: f64) -> f64114 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 in (0..=N).rev() {
121 num = num * x + SNUM[i];
122 den = den * x + SDEN[i];
123 }
124 } else {
125 for i in 0..=N {
126 num = num / x + SNUM[i];
127 den = den / x + SDEN[i];
128 }
129 }
130 return num / den;
131 }
132
tgamma(mut x: f64) -> f64133 pub fn tgamma(mut x: f64) -> f64 {
134 let u: u64 = x.to_bits();
135 let absx: f64;
136 let mut y: f64;
137 let mut dy: f64;
138 let mut z: f64;
139 let mut r: f64;
140 let ix: u32 = ((u >> 32) as u32) & 0x7fffffff;
141 let sign: bool = (u >> 63) != 0;
142
143 /* special cases */
144 if ix >= 0x7ff00000 {
145 /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
146 return x + core::f64::INFINITY;
147 }
148 if ix < ((0x3ff - 54) << 20) {
149 /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
150 return 1.0 / x;
151 }
152
153 /* integer arguments */
154 /* raise inexact when non-integer */
155 if x == floor(x) {
156 if sign {
157 return 0.0 / 0.0;
158 }
159 if x <= FACT.len() as f64 {
160 return FACT[(x as usize) - 1];
161 }
162 }
163
164 /* x >= 172: tgamma(x)=inf with overflow */
165 /* x =< -184: tgamma(x)=+-0 with underflow */
166 if ix >= 0x40670000 {
167 /* |x| >= 184 */
168 if sign {
169 let x1p_126 = f64::from_bits(0x3810000000000000); // 0x1p-126 == 2^-126
170 force_eval!((x1p_126 / x) as f32);
171 if floor(x) * 0.5 == floor(x * 0.5) {
172 return 0.0;
173 } else {
174 return -0.0;
175 }
176 }
177 let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 == 2^1023
178 x *= x1p1023;
179 return x;
180 }
181
182 absx = if sign { -x } else { x };
183
184 /* handle the error of x + g - 0.5 */
185 y = absx + GMHALF;
186 if absx > GMHALF {
187 dy = y - absx;
188 dy -= GMHALF;
189 } else {
190 dy = y - GMHALF;
191 dy -= absx;
192 }
193
194 z = absx - 0.5;
195 r = s(absx) * exp(-y);
196 if x < 0.0 {
197 /* reflection formula for negative x */
198 /* sinpi(absx) is not 0, integers are already handled */
199 r = -PI / (sinpi(absx) * absx * r);
200 dy = -dy;
201 z = -z;
202 }
203 r += dy * (GMHALF + 0.5) * r / y;
204 z = pow(y, 0.5 * z);
205 y = r * z * z;
206 return y;
207 }
208