diff options
author | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-06-07 05:48:48 +0000 |
---|---|---|
committer | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-06-07 05:48:48 +0000 |
commit | ef24de24a82fe681581cc130f342363c47c0969a (patch) | |
tree | 0d494f7e1a38b95c92426f58fe6eaa877303a86c /vendor/libm-0.1.4/src/math/tgamma.rs | |
parent | Releasing progress-linux version 1.74.1+dfsg1-1~progress7.99u1. (diff) | |
download | rustc-ef24de24a82fe681581cc130f342363c47c0969a.tar.xz rustc-ef24de24a82fe681581cc130f342363c47c0969a.zip |
Merging upstream version 1.75.0+dfsg1.
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'vendor/libm-0.1.4/src/math/tgamma.rs')
-rw-r--r-- | vendor/libm-0.1.4/src/math/tgamma.rs | 207 |
1 files changed, 0 insertions, 207 deletions
diff --git a/vendor/libm-0.1.4/src/math/tgamma.rs b/vendor/libm-0.1.4/src/math/tgamma.rs deleted file mode 100644 index f8ccf669a..000000000 --- a/vendor/libm-0.1.4/src/math/tgamma.rs +++ /dev/null @@ -1,207 +0,0 @@ -/* -"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 = (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 + SNUM[i]; - den = den * x + SDEN[i]; - } - } else { - for i in 0..=N { - num = num / x + SNUM[i]; - den = den / x + SDEN[i]; - } - } - return num / den; -} - -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 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; -} |