Source code
Revision control
Copy as Markdown
Other Tools
/*
"A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
"Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
"An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
approximation method:
(x - 0.5) S(x)
Gamma(x) = (x + g - 0.5) * ----------------
exp(x + g - 0.5)
with
a1 a2 a3 aN
S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
x + 1 x + 2 x + 3 x + N
with a0, a1, a2, a3,.. aN constants which depend on g.
for x < 0 the following reflection formula is used:
Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
most ideas and constants are from boost and python
*/
extern crate core;
use super::{exp, floor, k_cos, k_sin, pow};
const PI: f64 = 3.141592653589793238462643383279502884;
/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
fn sinpi(mut x: f64) -> f64 {
let mut n: isize;
/* argument reduction: x = |x| mod 2 */
/* spurious inexact when x is odd int */
x = x * 0.5;
x = 2.0 * (x - floor(x));
/* reduce x into [-.25,.25] */
n = (4.0 * x) as isize;
n = div!(n + 1, 2);
x -= (n as f64) * 0.5;
x *= PI;
match n {
1 => k_cos(x, 0.0),
2 => k_sin(-x, 0.0, 0),
3 => -k_cos(x, 0.0),
0 | _ => k_sin(x, 0.0, 0),
}
}
const N: usize = 12;
//static const double g = 6.024680040776729583740234375;
const GMHALF: f64 = 5.524680040776729583740234375;
const SNUM: [f64; N + 1] = [
23531376880.410759688572007674451636754734846804940,
42919803642.649098768957899047001988850926355848959,
35711959237.355668049440185451547166705960488635843,
17921034426.037209699919755754458931112671403265390,
6039542586.3520280050642916443072979210699388420708,
1439720407.3117216736632230727949123939715485786772,
248874557.86205415651146038641322942321632125127801,
31426415.585400194380614231628318205362874684987640,
2876370.6289353724412254090516208496135991145378768,
186056.26539522349504029498971604569928220784236328,
8071.6720023658162106380029022722506138218516325024,
210.82427775157934587250973392071336271166969580291,
2.5066282746310002701649081771338373386264310793408,
];
const SDEN: [f64; N + 1] = [
0.0,
39916800.0,
120543840.0,
150917976.0,
105258076.0,
45995730.0,
13339535.0,
2637558.0,
357423.0,
32670.0,
1925.0,
66.0,
1.0,
];
/* n! for small integer n */
const FACT: [f64; 23] = [
1.0,
1.0,
2.0,
6.0,
24.0,
120.0,
720.0,
5040.0,
40320.0,
362880.0,
3628800.0,
39916800.0,
479001600.0,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
121645100408832000.0,
2432902008176640000.0,
51090942171709440000.0,
1124000727777607680000.0,
];
/* S(x) rational function for positive x */
fn s(x: f64) -> f64 {
let mut num: f64 = 0.0;
let mut den: f64 = 0.0;
/* to avoid overflow handle large x differently */
if x < 8.0 {
for i in (0..=N).rev() {
num = num * x + i!(SNUM, i);
den = den * x + i!(SDEN, i);
}
} else {
for i in 0..=N {
num = num / x + i!(SNUM, i);
den = den / x + i!(SDEN, i);
}
}
return num / den;
}
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
pub fn tgamma(mut x: f64) -> f64 {
let u: u64 = x.to_bits();
let absx: f64;
let mut y: f64;
let mut dy: f64;
let mut z: f64;
let mut r: f64;
let ix: u32 = ((u >> 32) as u32) & 0x7fffffff;
let sign: bool = (u >> 63) != 0;
/* special cases */
if ix >= 0x7ff00000 {
/* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
return x + core::f64::INFINITY;
}
if ix < ((0x3ff - 54) << 20) {
/* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
return 1.0 / x;
}
/* integer arguments */
/* raise inexact when non-integer */
if x == floor(x) {
if sign {
return 0.0 / 0.0;
}
if x <= FACT.len() as f64 {
return i!(FACT, (x as usize) - 1);
}
}
/* x >= 172: tgamma(x)=inf with overflow */
/* x =< -184: tgamma(x)=+-0 with underflow */
if ix >= 0x40670000 {
/* |x| >= 184 */
if sign {
let x1p_126 = f64::from_bits(0x3810000000000000); // 0x1p-126 == 2^-126
force_eval!((x1p_126 / x) as f32);
if floor(x) * 0.5 == floor(x * 0.5) {
return 0.0;
} else {
return -0.0;
}
}
let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 == 2^1023
x *= x1p1023;
return x;
}
absx = if sign { -x } else { x };
/* handle the error of x + g - 0.5 */
y = absx + GMHALF;
if absx > GMHALF {
dy = y - absx;
dy -= GMHALF;
} else {
dy = y - GMHALF;
dy -= absx;
}
z = absx - 0.5;
r = s(absx) * exp(-y);
if x < 0.0 {
/* reflection formula for negative x */
/* sinpi(absx) is not 0, integers are already handled */
r = -PI / (sinpi(absx) * absx * r);
dy = -dy;
z = -z;
}
r += dy * (GMHALF + 0.5) * r / y;
z = pow(y, 0.5 * z);
y = r * z * z;
return y;
}