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