libm/math/
tgamma.rs

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
6approximation method:
7
8                        (x - 0.5)         S(x)
9Gamma(x) = (x + g - 0.5)         *  ----------------
10                                    exp(x + g - 0.5)
11
12with
13                 a1      a2      a3            aN
14S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
15               x + 1   x + 2   x + 3         x + N
16
17with a0, a1, a2, a3,.. aN constants which depend on g.
18
19for x < 0 the following reflection formula is used:
20
21Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
22
23most ideas and constants are from boost and python
24*/
25use super::{exp, floor, k_cos, k_sin, pow};
26
27const PI: f64 = 3.141592653589793238462643383279502884;
28
29/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
30fn sinpi(mut x: f64) -> f64 {
31    let mut n: isize;
32
33    /* argument reduction: x = |x| mod 2 */
34    /* spurious inexact when x is odd int */
35    x = x * 0.5;
36    x = 2.0 * (x - floor(x));
37
38    /* reduce x into [-.25,.25] */
39    n = (4.0 * x) as isize;
40    n = div!(n + 1, 2);
41    x -= (n as f64) * 0.5;
42
43    x *= PI;
44    match n {
45        1 => k_cos(x, 0.0),
46        2 => k_sin(-x, 0.0, 0),
47        3 => -k_cos(x, 0.0),
48        0 | _ => k_sin(x, 0.0, 0),
49    }
50}
51
52const N: usize = 12;
53//static const double g = 6.024680040776729583740234375;
54const GMHALF: f64 = 5.524680040776729583740234375;
55const SNUM: [f64; N + 1] = [
56    23531376880.410759688572007674451636754734846804940,
57    42919803642.649098768957899047001988850926355848959,
58    35711959237.355668049440185451547166705960488635843,
59    17921034426.037209699919755754458931112671403265390,
60    6039542586.3520280050642916443072979210699388420708,
61    1439720407.3117216736632230727949123939715485786772,
62    248874557.86205415651146038641322942321632125127801,
63    31426415.585400194380614231628318205362874684987640,
64    2876370.6289353724412254090516208496135991145378768,
65    186056.26539522349504029498971604569928220784236328,
66    8071.6720023658162106380029022722506138218516325024,
67    210.82427775157934587250973392071336271166969580291,
68    2.5066282746310002701649081771338373386264310793408,
69];
70const SDEN: [f64; N + 1] = [
71    0.0,
72    39916800.0,
73    120543840.0,
74    150917976.0,
75    105258076.0,
76    45995730.0,
77    13339535.0,
78    2637558.0,
79    357423.0,
80    32670.0,
81    1925.0,
82    66.0,
83    1.0,
84];
85/* n! for small integer n */
86const FACT: [f64; 23] = [
87    1.0,
88    1.0,
89    2.0,
90    6.0,
91    24.0,
92    120.0,
93    720.0,
94    5040.0,
95    40320.0,
96    362880.0,
97    3628800.0,
98    39916800.0,
99    479001600.0,
100    6227020800.0,
101    87178291200.0,
102    1307674368000.0,
103    20922789888000.0,
104    355687428096000.0,
105    6402373705728000.0,
106    121645100408832000.0,
107    2432902008176640000.0,
108    51090942171709440000.0,
109    1124000727777607680000.0,
110];
111
112/* S(x) rational function for positive x */
113fn s(x: f64) -> f64 {
114    let mut num: f64 = 0.0;
115    let mut den: f64 = 0.0;
116
117    /* to avoid overflow handle large x differently */
118    if x < 8.0 {
119        for i in (0..=N).rev() {
120            num = num * x + i!(SNUM, i);
121            den = den * x + i!(SDEN, i);
122        }
123    } else {
124        for i in 0..=N {
125            num = num / x + i!(SNUM, i);
126            den = den / x + i!(SDEN, i);
127        }
128    }
129    return num / den;
130}
131
132#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
133pub 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 i!(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}
libm/math/tgamma.rs

libm/math/
tgamma.rs
