diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..021f3eba 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -237,7 +237,11 @@ impl Distribution for Binomial { } else { (0..self.n + 1).fold(0.0, |acc, x| { let p = self.pmf(x); - acc - p * p.ln() + // A mass that underflows to zero contributes nothing, taking + // `0 * ln 0 == 0` as usual for entropy. Evaluating it would give + // `0 * -inf == NaN` and poison the whole sum, which made + // `entropy()` return `NaN` for any `n` past roughly 1100. + if p > 0.0 { acc - p * p.ln() } else { acc } }) }; Some(entr) @@ -412,6 +416,39 @@ mod tests { test_exact(0.3, 10, 10, max); } + /// `p` very close to but not equal to 1 must take the general path, not the + /// degenerate `p == 1` branch. The branch is guarded by `prec::ulps_eq!`, + /// whose default epsilon was `1e-9` *absolute*, so any `p` within `1e-9` of + /// 1 collapsed to a point mass at `n`. + /// + /// `1 - 2^-33` is used rather than `1 - 1e-10` so that `1 - p` is exact and + /// the reference values are not limited by the representation of `p`. + #[test] + fn test_pmf_p_near_one_is_not_degenerate() { + let n = Binomial::new(1.0 - f64::powi(2.0, -33), 100).unwrap(); + prec::assert_relative_eq!(n.pmf(99), 1.1641532048523463366e-8, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(n.pmf(100), 0.99999998835846788439, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(n.ln_pmf(99), -18.268686784015220704, epsilon = 0.0, max_relative = 5e-15); + // entropy of a non-degenerate distribution is strictly positive + assert!(n.entropy().unwrap() > 0.0); + } + + /// `entropy()` sums `-p * ln(p)` over the whole support, so it used to + /// return `NaN` as soon as one mass underflowed to zero (`0 * -inf`). For + /// `p = 0.5` that starts at around `n = 1100`. + #[test] + fn test_entropy_with_underflowing_masses() { + for n in [1_100, 2_000, 5_000, 20_000] { + let d = Binomial::new(0.5, n).unwrap(); + let h = d.entropy().unwrap(); + assert!(h.is_finite(), "entropy for n={n} was {h}"); + // 0.5 * ln(2 * pi * e * n * p * q) is the asymptotic form + let approx = 0.5 * (2.0 * f64::consts::PI * f64::consts::E * n as f64 * 0.25).ln(); + prec::assert_relative_eq!(h, approx, epsilon = 0.0, max_relative = 1e-3); + } + assert!(Binomial::new(0.1, 5_000).unwrap().entropy().unwrap().is_finite()); + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Binomial| x.pmf(arg); @@ -515,7 +552,7 @@ mod tests { test_absolute(0.3, 3, 0.657, 1e-14, sf(0)); test_absolute(0.3, 3, 0.216, 1e-15, sf(1)); test_exact(0.3, 3, 0.0, sf(3)); - test_absolute(0.3, 10, 0.9717524751000001, 1e-16, sf(0)); + test_absolute(0.3, 10, 0.9717524750999999955198, 1e-16, sf(0)); test_absolute(0.3, 10, 0.850691654100002, 1e-14, sf(1)); test_exact(0.3, 10, 0.0, sf(10)); test_exact(1.0, 1, 1.0, sf(0)); diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index 5f1697b4..51dd809b 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -19,7 +19,7 @@ use nalgebra::{Dim, Dyn, OMatrix, OVector}; /// /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap(); /// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5])); -/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205); +/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.2221555562222193); /// ``` #[derive(Clone, PartialEq, Debug)] pub struct Dirichlet diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs index 9cc2f5c4..cd5e160f 100644 --- a/src/distribution/fisher_snedecor.rs +++ b/src/distribution/fisher_snedecor.rs @@ -525,8 +525,8 @@ mod tests { test_absolute(10.0, 0.1, 0.0418440630400545297349, 1e-14, pdf(1.0)); test_absolute(0.1, 1.0, 0.0396064560910663979961, 1e-16, pdf(1.0)); test_absolute(1.0, 1.0, 0.1591549430918953357689, 1e-16, pdf(1.0)); - test_absolute(10.0, 1.0, 0.230361989229138647108, 1e-16, pdf(1.0)); - test_absolute(0.1, 0.1, 0.00221546909694001013517, 1e-18, pdf(10.0)); + test_absolute(10.0, 1.0, 0.230361989229138647108, 3e-16, pdf(1.0)); + test_absolute(0.1, 0.1, 0.00221546909694001013517, 5e-18, pdf(10.0)); test_absolute(1.0, 0.1, 0.00369960370387922619592, 1e-17, pdf(10.0)); test_absolute(10.0, 0.1, 0.00390179721174142927402, 1e-15, pdf(10.0)); test_absolute(0.1, 1.0, 0.00319864073359931548273, 1e-17, pdf(10.0)); @@ -540,13 +540,13 @@ mod tests { #[test] fn test_ln_pdf() { let ln_pdf = |arg: f64| move |x: FisherSnedecor| x.ln_pdf(arg); - test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 1e-15, ln_pdf(1.0)); + test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 3e-15, ln_pdf(1.0)); test_absolute(1.0, 0.1, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(10.0, 0.1, 0.0418440630400545297349f64.ln(), 1e-13, ln_pdf(1.0)); test_absolute(0.1, 1.0, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(1.0, 1.0, 0.1591549430918953357689f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(10.0, 1.0, 0.230361989229138647108f64.ln(), 1e-15, ln_pdf(1.0)); - test_exact(0.1, 0.1, 0.00221546909694001013517f64.ln(), ln_pdf(10.0)); + test_absolute(0.1, 0.1, 0.00221546909694001013517f64.ln(), 3e-15, ln_pdf(10.0)); test_absolute(1.0, 0.1, 0.00369960370387922619592f64.ln(), 1e-15, ln_pdf(10.0)); test_absolute(10.0, 0.1, 0.00390179721174142927402f64.ln(), 1e-13, ln_pdf(10.0)); test_absolute(0.1, 1.0, 0.00319864073359931548273f64.ln(), 1e-15, ln_pdf(10.0)); diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index f74f2a79..3a3f9a78 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -437,6 +437,18 @@ mod tests { test_exact(0.3, u64::MAX, max); } + /// As in [`Binomial`](crate::distribution::Binomial), the `p == 1` branches + /// here are guarded by `prec::ulps_eq!`; its default epsilon was `1e-9` + /// absolute, so `p = 1 - 2^-33` was treated as a point mass at 1. + #[test] + fn test_p_near_one_is_not_degenerate() { + let g = Geometric::new(1.0 - f64::powi(2.0, -33)).unwrap(); + assert_eq!(g.max(), u64::MAX); + prec::assert_relative_eq!(g.skewness().unwrap(), 92681.900034472750337, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(g.pmf(2), 1.164153218133822873e-10, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(g.ln_pmf(2), -22.873856958594610533, epsilon = 0.0, max_relative = 1e-14); + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Geometric| x.pmf(arg); diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 2ad69880..cee52829 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -19,7 +19,7 @@ use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; /// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.).unwrap(); /// assert_eq!(mvs.mean().unwrap(), DVector::from_vec(vec![0., 0.])); /// assert_eq!(mvs.variance().unwrap(), DMatrix::from_vec(2, 2, vec![2., 0., 0., 2.])); -/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.04715702017537655); +/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.047157020175376395); /// ``` #[derive(Debug, Clone, PartialEq)] pub struct MultivariateStudent diff --git a/src/function/beta.rs b/src/function/beta.rs index 970e76ed..3c1ff6ed 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -138,6 +138,18 @@ pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 { /// `b` is the second beta parameter, and `x` is the upper limit of the /// integral. /// +/// # Remarks +/// +/// Relative accuracy degrades as `a + b` grows, because the leading factor is +/// evaluated as `exp(ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b) + ...)` and the +/// cancellation in that exponent grows with `ln_gamma(a + b)`. Measured against +/// the exact identity `I_{1/2}(a, a) == 1/2`, the relative error is about +/// `6e-11` at `a = b = 1e4` and `3e-9` at `1e6`. +/// +/// Past `min(a, b) ~ 1e16` the recurrence is truncated by its iteration bound and +/// the result is unreliable, though still clamped to `[0, 1]`. Evaluation is +/// bounded at roughly 10 ms in the worst case. +/// /// # Errors /// /// if `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` @@ -154,6 +166,16 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { return Err(BetaFuncError::XOutOfRange); } + // `I_0(a, b) == 0` and `I_1(a, b) == 1` for every `a` and `b`. Handling the + // endpoints here keeps them independent of the symmetry test below, which + // otherwise mapped `x == 0` to `1.0` once `a + b` overflowed. + if x == 0.0 { + return Ok(0.0); + } + if x == 1.0 { + return Ok(1.0); + } + let bt = if x == 0.0 || crate::prec::ulps_eq!(x, 1.0, epsilon = MODULE_EPS) { 0.0 } else { @@ -162,10 +184,65 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { + b * (1.0 - x).ln()) .exp() }; - let symm_transform = x >= (a + 1.0) / (a + b + 2.0); + let symm_transform = { + let denom = a + b + 2.0; + if denom.is_finite() { + x >= (a + 1.0) / denom + } else { + // `a + b` overflowed, which would collapse the threshold to zero and + // send every `x` down the transformed branch. Scaling numerator and + // denominator by `max(a, b)` keeps the ratio exact enough to compare. + let m = a.max(b); + x >= (a / m + 1.0 / m) / (a / m + b / m + 2.0 / m) + } + }; + + // The result is `bt * h / a` with the continued fraction `h` of order one, + // so a prefix that has underflowed pins the answer at the corresponding + // endpoint. Returning here also avoids forming `0.0 * h`, which is NaN + // whenever the recurrence overflowed (e.g. `beta_reg(1e300, 1e-300, 0.5)`), + // and skips the recurrence entirely in the regime where the distribution + // has concentrated to a step function. + if bt == 0.0 { + return Ok(if symm_transform { 1.0 } else { 0.0 }); + } + + // Fallback for the regime where the recurrence below is truncated by + // `max_iters` and can degenerate to a non-finite value: the distribution has + // concentrated around its mean, so this is the limiting step function. It is + // only consulted when the recurrence produced something unusable - see + // `finish`. + let saturated = { + let mean = a / (a + b); + if x < mean { + 0.0 + } else if x > mean { + 1.0 + } else { + 0.5 + } + }; + let eps = prec::F64_PREC; let fpmin = f64::MIN_POSITIVE / eps; + // Iterations the Lentz recurrence below needs before `del` settles. It is + // slowest at the centre of the distribution (`x ~ a / (a + b)`), where the + // worst case over `x` grows like `5 * min(a, b).cbrt()`; the bound here + // carries headroom on top of that. A fixed bound of 140 used to be applied + // regardless of `a` and `b`, which silently truncated the recurrence and + // returned a badly wrong value once `min(a, b)` passed ~1.5e4 (for example + // `I_0.5(1e6, 1e6)` came back as 0.491 instead of 0.5). The loop still + // stops as soon as it converges, so the typical few-dozen-iteration case is + // unchanged. + // + // The upper clamp bounds the work at roughly 10 ms. Past `min(a, b) ~ 1e15` + // the recurrence needs more iterations than that, but it has also stopped + // being able to deliver an accurate answer (its own rounding accumulates + // over millions of steps), so spending longer buys nothing - see the + // accuracy note on `checked_beta_reg`. + let max_iters = ((8.0 * a.min(b).cbrt()) as u32).clamp(140, 1_000_000); + let mut a = a; let mut b = b; let mut x = x; @@ -188,7 +265,7 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { d = 1.0 / d; let mut h = d; - for m in 1..141 { + for m in 1..=max_iters { let m = f64::from(m); let m2 = m * 2.0; let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2)); @@ -223,18 +300,28 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { h *= del; if (del - 1.0).abs() <= eps { - return if symm_transform { - Ok(1.0 - bt * h / a) - } else { - Ok(bt * h / a) - }; + return Ok(finish(symm_transform, bt, h, a, saturated)); } } - if symm_transform { - Ok(1.0 - bt * h / a) + Ok(finish(symm_transform, bt, h, a, saturated)) +} + +/// Assembles `I_x(a, b)` from the prefix and continued fraction, keeping the +/// result inside `[0, 1]`. +/// +/// Neither guard engages while the recurrence converges. Once it is truncated by +/// `max_iters` (only for `min(a, b)` past ~1e16) the raw value can drift outside +/// the unit interval - `-2.56` at `a = b = 1e20` - or become non-finite +/// entirely, and callers such as `Binomial::cdf` are contractually +/// probabilities. +fn finish(symm_transform: bool, bt: f64, h: f64, a: f64, saturated: f64) -> f64 { + let v = bt * h / a; + let v = if symm_transform { 1.0 - v } else { v }; + if v.is_finite() { + v.clamp(0.0, 1.0) } else { - Ok(bt * h / a) + saturated } } @@ -646,6 +733,61 @@ mod tests { assert!(checked_beta_reg(1.0, 1.0, 2.0).is_err()); } + /// `beta_reg` is a probability and must stay in `[0, 1]` and finite for + /// every valid input, including parameter ratios extreme enough to + /// over/underflow the intermediate quantities. Before the short-circuit on an + /// underflowed prefix and the `[0, 1]` clamp, this grid produced NaN (from + /// `0.0 * inf` once the recurrence overflowed) and `-2.56` at `a = b = 1e20`. + #[test] + fn test_beta_reg_extreme_parameters_stay_a_probability() { + let params = [ + 1e-308f64, 1e-300, 1e-100, 1e-8, 0.5, 1.0, 20.0, 1e8, 1e20, 1e100, 1e300, 1e308, + ]; + for &a in ¶ms { + for &b in ¶ms { + // Beyond a ~1e300 parameter ratio the Lentz recurrence bottoms + // out in its own `fpmin` guards for `x` within an ulp of the + // mode, and returns NaN. That is pre-existing and unreachable + // from any distribution in the crate (`Binomial` is bounded by + // `n <= u64::MAX`), so it is excluded rather than papered over + // with a plausible-looking wrong value. + if a.max(b) / a.min(b) > 1e200 { + continue; + } + for x in [0.0f64, 1e-300, 0.25, 0.5, 0.75, 1.0] { + let v = beta_reg(a, b, x); + assert!( + v.is_finite() && (0.0..=1.0).contains(&v), + "beta_reg({a:e}, {b:e}, {x}) = {v}" + ); + } + // monotone in x, and pinned at the endpoints + assert_eq!(beta_reg(a, b, 0.0), 0.0, "beta_reg({a:e},{b:e},0)"); + assert_eq!(beta_reg(a, b, 1.0), 1.0, "beta_reg({a:e},{b:e},1)"); + } + } + } + + /// `I_x(a, b) + I_{1-x}(b, a) == 1` for every valid `a`, `b`, `x`. The two + /// sides truncate differently, so a prematurely stopped continued fraction + /// breaks the identity. + #[test] + fn test_beta_reg_complement_identity_large_parameters() { + for (a, b) in [ + (1e5, 1e5), + (1e6, 1e6), + (1e7, 1e7), + (1e6, 1e3), + (1e3, 1e6), + (2e4, 3e4), + ] { + for x in [0.1, 0.25, 0.5, 0.5 + 1e-9, 0.75, 0.9] { + let lhs = beta_reg(a, b, x) + beta_reg(b, a, 1.0 - x); + prec::assert_abs_diff_eq!(lhs, 1.0, epsilon = 1e-7); + } + } + } + #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/function/gamma.rs b/src/function/gamma.rs index bff5c45a..13d2f126 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -2,7 +2,9 @@ //! related functions use crate::consts; +use crate::function::evaluate; use crate::prec; +use crate::prec::{dekker_product_err, two_diff}; use core::f64; #[cfg(not(feature = "std"))] use num_traits::Float as _; @@ -30,48 +32,183 @@ impl core::fmt::Display for GammaFuncError { impl core::error::Error for GammaFuncError {} +/// Computes `sin(PI * x)` for use in reflection formulas. +/// +/// Evaluating `(PI * x).sin()` directly loses the fractional part of `x` to +/// rounding once `PI * x` is large, and near the zeros of the sine (integer +/// `x`) the representation error of `PI` alone costs `~1.2e-16 / dist_to_pole` +/// of relative accuracy. Reducing with the period first (`x - round(x)` is +/// exact) keeps the error relative to the fractional part instead. +#[inline] +fn sin_pi(x: f64) -> f64 { + let m = x.round(); + let r = x - m; // exact, |r| <= 0.5 + // sin(PI * (m + r)) = (-1)^m * sin(PI * r); `m * 0.5` is exact, so its + // fractional part is 0.5 exactly when `m` is odd. (For |m| >= 2^53 every + // representable float is even.) + let sin_r = (f64::consts::PI * r).sin(); + if (m * 0.5).fract() == 0.0 { + sin_r + } else { + -sin_r + } +} + +/// Computes `tan(PI * x)` with the same period reduction as [`sin_pi`] +/// (`tan` has period `PI`, so no sign correction is needed). +#[inline] +fn tan_pi(x: f64) -> f64 { + let r = x - x.round(); // exact, |r| <= 0.5 + (f64::consts::PI * r).tan() +} + /// Auxiliary variable when evaluating the `gamma_ln` function const GAMMA_R: f64 = 10.900511; -/// Polynomial coefficients for approximating the `gamma_ln` function -const GAMMA_DK: &[f64] = &[ - 2.48574089138753565546e-5, - 1.05142378581721974210, - -3.45687097222016235469, - 4.51227709466894823700, - -2.98285225323576655721, - 1.05639711577126713077, - -1.95428773191645869583e-1, - 1.70970543404441224307e-2, - -5.71926117404305781283e-4, - 4.63399473359905636708e-6, - -2.71994908488607703910e-9, +// The Lanczos sum used by `ln_gamma` and `gamma` is Pugh's 10-term +// partial-fraction approximation with g = GAMMA_R ("An Analysis of the Lanczos +// Gamma Approximation", Glendon Ralph Pugh, 2004 p. 116): +// +// S(z) = d_0 + sum_{k=1..10} d_k / (z + k - 1) +// +// with residues d_k = [2.48574089138753565546e-5, 1.05142378581721974210, +// -3.45687097222016235469, 4.51227709466894823700, -2.98285225323576655721, +// 1.05639711577126713077, -1.95428773191645869583e-1, 1.70970543404441224307e-2, +// -5.71926117404305781283e-4, 4.63399473359905636708e-6, +// -2.71994908488607703910e-9]. +// +// The residues alternate in sign, so evaluating the sum directly cancels +// catastrophically (condition number ~3600 around z = 50, i.e. ~440 eps of +// relative error in the sum). The constants below are the mathematically +// identical single-fraction form S(z) = N(z) / D(z), derived from the d_k in +// exact rational arithmetic: D(z) = z (z+1) ... (z+9) (whose expanded +// coefficients are unsigned Stirling numbers of the first kind, exact in f64) +// and N = d_0 D + sum_k d_k D / (z + k - 1). Every coefficient of both +// polynomials is positive, so for z > 0 both Horner evaluations have condition +// number 1 and the sum is accurate to a few eps. + +/// Numerator `N(z)` of the Lanczos sum in single-fraction form (ascending). +const LANCZOS_NUM: &[f64] = &[ + 381540.6633973527, + 365505.352696257, + 157567.99949360118, + 40253.83538142639, + 6748.767525934567, + 775.8779405455635, + 61.94528891422096, + 3.391366244015308, + 0.12184807036444657, + 0.002594340508809067, + 2.4857408913875355e-5, +]; + +/// Denominator `D(z) = z (z+1) ... (z+9)` expanded (ascending; unsigned +/// Stirling numbers of the first kind, exact integers). +const LANCZOS_DENOM: &[f64] = &[ + 0.0, 362880.0, 1026576.0, 1172700.0, 723680.0, 269325.0, 63273.0, 9450.0, 870.0, 45.0, 1.0, +]; + +/// `LANCZOS_NUM` reversed: `N(z) / z^10` as a polynomial in `w = 1/z`. +const LANCZOS_NUM_REV: &[f64] = &[ + 2.4857408913875355e-5, + 0.002594340508809067, + 0.12184807036444657, + 3.391366244015308, + 61.94528891422096, + 775.8779405455635, + 6748.767525934567, + 40253.83538142639, + 157567.99949360118, + 365505.352696257, + 381540.6633973527, ]; +/// `LANCZOS_DENOM` reversed: `D(z) / z^10` as a polynomial in `w = 1/z`. +const LANCZOS_DENOM_REV: &[f64] = &[ + 1.0, 45.0, 870.0, 9450.0, 63273.0, 269325.0, 723680.0, 1172700.0, 1026576.0, 362880.0, 0.0, +]; + +/// Low half of Euler's number: `e == f64::consts::E + E_LO` to double-double +/// precision. +const E_LO: f64 = 1.4456468917292502e-16; + +/// Computes `((p + GAMMA_R) / e)^exponent`, compensating the rounding of the base. +/// +/// `powf` amplifies a relative error `eps` in its base by the exponent, so the +/// two roundings in `(p + GAMMA_R) / e` cost `~2 |p|` ulps - the dominant +/// error of the direct evaluation (~190 ulps in `gamma` by x = 122). The +/// rounding residuals of the addition (two-sum) and the division (Dekker +/// product, plus the low word of `e`) are recovered exactly and applied as the +/// first-order correction `(b (1 + delta))^p ~= b^p (1 + p delta)`. +/// `p_err` is the rounding residual of the exponent itself (the true exponent +/// is `p + p_err`), folded in as `b^(p_err) ~= 1 + p_err * ln(b)`; it must be +/// zero unless `exponent == p`. +/// +/// `exponent` is normally `p`. The overflow-avoiding path in [`gamma`] passes +/// `p / 2` and squares the result, which needs the *base* to keep coming from +/// the full `p` - halving `p` in the base as well would compute an entirely +/// different quantity. +#[inline] +fn lanczos_power(p: f64, p_err: f64, exponent: f64) -> f64 { + let (zgh, add_err) = two_diff(p, -GAMMA_R); + let b = zgh / f64::consts::E; + let pw = b * f64::consts::E; + // exact residual of the division against the double-double e + let div_err = (zgh - pw) - dekker_product_err(b, f64::consts::E, pw); + let delta = (div_err + add_err + p_err - b * E_LO) / zgh; + let mut corr = exponent * delta; + if p_err != 0.0 { + corr += p_err * b.ln(); + } + let base = b.powf(exponent); + // `(b (1 + delta))^p ~= b^p (1 + p delta)` only holds while `|p delta| << 1`. + // Since `delta` is of order `eps`, that fails once `p` reaches ~1e13 - far + // past where `b^p` overflows, so the uncorrected value is the right answer + // there. Applying it anyway would let a negative `corr` flip the sign and + // turn an overflowing `gamma` into `-inf`. + if corr.abs() < 0.25 { + base * (1.0 + corr) + } else { + base + } +} + +/// Evaluates the Lanczos sum `S(z)` for `z >= 0.5` via the well-conditioned +/// single-fraction form (see the derivation note above). +#[inline] +fn lanczos_sum(z: f64) -> f64 { + if z < 1e29 { + evaluate::polynomial(z, LANCZOS_NUM) / evaluate::polynomial(z, LANCZOS_DENOM) + } else { + // z^10 would overflow; divide both polynomials by z^10 and evaluate + // in w = 1/z instead. + let w = 1.0 / z; + evaluate::polynomial(w, LANCZOS_NUM_REV) / evaluate::polynomial(w, LANCZOS_DENOM_REV) + } +} + /// Computes the logarithm of the gamma function /// with an accuracy of 16 floating point digits. /// The implementation is derived from /// "An Analysis of the Lanczos Gamma Approximation", /// Glendon Ralph Pugh, 2004 p. 116 pub fn ln_gamma(x: f64) -> f64 { + // ln Gamma(n) = ln((n - 1)!) via the exact factorial table (~0.4 ulp, and + // exactly 0 at n = 1, 2, where the Lanczos formula's terms only cancel + // approximately). + if x.fract() == 0.0 && (1.0..=171.0).contains(&x) { + return crate::function::factorial::ln_factorial(x as u64 - 1); + } if x < 0.5 { - let s = GAMMA_DK - .iter() - .enumerate() - .skip(1) - .fold(GAMMA_DK[0], |s, t| s + t.1 / (t.0 as f64 - x)); + let s = lanczos_sum(1.0 - x); consts::LN_PI - - (f64::consts::PI * x).sin().ln() + - sin_pi(x).ln() - s.ln() - consts::LN_2_SQRT_E_OVER_PI - (0.5 - x) * ((0.5 - x + GAMMA_R) / f64::consts::E).ln() } else { - let s = GAMMA_DK - .iter() - .enumerate() - .skip(1) - .fold(GAMMA_DK[0], |s, t| s + t.1 / (x + t.0 as f64 - 1.0)); + let s = lanczos_sum(x); s.ln() + consts::LN_2_SQRT_E_OVER_PI @@ -79,31 +216,43 @@ pub fn ln_gamma(x: f64) -> f64 { } } -/// Computes the gamma function with an accuracy -/// of 16 floating point digits. The implementation -/// is derived from "An Analysis of the Lanczos Gamma Approximation", -/// Glendon Ralph Pugh, 2004 p. 116 +/// Computes the gamma function. The implementation is derived from +/// "An Analysis of the Lanczos Gamma Approximation", +/// Glendon Ralph Pugh, 2004 p. 116. +/// +/// Exact at the positive integers up to 171; elsewhere accurate to about +/// `4e-15` relative, which is the approximation floor of the (f64-rounded) +/// Pugh coefficient set itself - the evaluation is compensated to well below +/// that. pub fn gamma(x: f64) -> f64 { + // Gamma(n) = (n - 1)! exactly at the positive integers where the factorial + // is representable; the Lanczos path is only good to ~1 ulp there. + if x.fract() == 0.0 && (1.0..=171.0).contains(&x) { + return crate::function::factorial::factorial(x as u64 - 1); + } if x < 0.5 { - let s = GAMMA_DK - .iter() - .enumerate() - .skip(1) - .fold(GAMMA_DK[0], |s, t| s + t.1 / (t.0 as f64 - x)); + let s = lanczos_sum(1.0 - x); + // 0.5 - x rounds for x below the binade of 0.5; keep its residual + let (pw, pw_err) = two_diff(0.5, x); f64::consts::PI - / ((f64::consts::PI * x).sin() - * s - * consts::TWO_SQRT_E_OVER_PI - * ((0.5 - x + GAMMA_R) / f64::consts::E).powf(0.5 - x)) + / (sin_pi(x) * s * consts::TWO_SQRT_E_OVER_PI * lanczos_power(pw, pw_err, pw)) } else { - let s = GAMMA_DK - .iter() - .enumerate() - .skip(1) - .fold(GAMMA_DK[0], |s, t| s + t.1 / (x + t.0 as f64 - 1.0)); - - s * consts::TWO_SQRT_E_OVER_PI * ((x - 0.5 + GAMMA_R) / f64::consts::E).powf(x - 0.5) + let s = lanczos_sum(x); + + // x - 0.5 is exact for 0.5 <= x < 2^52 (same or finer grid) + let p = x - 0.5; + if p > 168.0 { + // `lanczos_power` alone overflows from about x = 169.7, while the + // full product stays representable up to x ~ 171.61 (the true + // overflow point of Gamma). Halve the exponent and square at the + // end so the intermediate stays in range; this costs a couple of + // ulp, well inside the approximation floor of the coefficient set. + let half = + s.sqrt() * consts::TWO_SQRT_E_OVER_PI.sqrt() * lanczos_power(p, 0.0, 0.5 * p); + return half * half; + } + s * consts::TWO_SQRT_E_OVER_PI * lanczos_power(p, 0.0, p) } } @@ -390,7 +539,11 @@ pub fn digamma(x: f64) -> f64 { return f64::NEG_INFINITY; } if x < 0.0 { - return digamma(1.0 - x) + f64::consts::PI / (-f64::consts::PI * x).tan(); + // Reflection formula `psi(x) = psi(1 - x) - PI / tan(PI * x)`, with the + // period reduction done by `tan_pi`. `(PI * x).tan()` evaluated directly + // lost up to ~6 decimal digits near the poles (5586 ulp at x = -12.72, + // 3e-7 relative at 1e-10 from a pole). + return digamma(1.0 - x) - f64::consts::PI / tan_pi(x); } if x <= s { return d1 - 1.0 / x + d2 * x; @@ -583,9 +736,11 @@ mod tests { 11.51291869289055371493077240324332039045238086972508869965363, epsilon = 1e-14 ); - assert_eq!( + prec::assert_relative_eq!( super::ln_gamma(1.000001e-2), - 4.599478872433667224554543378460164306444416156144779542513592 + 4.599478872433667224554543378460164306444416156144779542513592, + epsilon = 0.0, + max_relative = 1e-15 ); prec::assert_abs_diff_eq!( super::ln_gamma(0.1), @@ -883,7 +1038,7 @@ mod tests { prec::assert_abs_diff_eq!( gamma_li(5.5, 2.0), 1.5746265342113649473739798668921124454837064926448459, - epsilon = 1e-15 + epsilon = 2e-15 ); prec::assert_abs_diff_eq!( gamma_li(5.5, 8.0), @@ -1195,6 +1350,34 @@ mod tests { assert!(super::checked_gamma_ui(1.0, f64::INFINITY).is_err()); } + /// `digamma` has poles at the non-positive integers and returns -inf there. + /// The pole test is `prec::ulps_eq!(x.floor(), x)`, whose default epsilon + /// used to be `1e-9` absolute, so a whole neighbourhood around each pole + /// returned -inf instead of a large finite value. + #[test] + fn test_digamma_near_negative_integer_poles_is_finite() { + prec::assert_relative_eq!(digamma(-1.0 + f64::powi(2.0, -33)), -8589934591.5772156646, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)), 8589934592.4227843348, epsilon = 0.0, max_relative = 1e-13); + // the poles themselves are unchanged + assert_eq!(digamma(-1.0), f64::NEG_INFINITY); + assert_eq!(digamma(0.0), f64::NEG_INFINITY); + assert_eq!(digamma(-5.0), f64::NEG_INFINITY); + } + + /// Before `tan_pi` reduced the reflection argument by the period, these + /// lost 3-4 decimal digits because `(PI * x).tan()` was evaluated at a + /// large rounded argument. Reference values are mpmath at 40 significant + /// digits, computed at the exact `f64` of each literal. The tolerances + /// reflect the cancellation between the two reflection terms + /// (`psi(1 - x)` and `pi * cot(pi x)` are each ~2.6 at x = -12.72 and + /// cancel to ~0.017, amplifying their ~1 ulp errors ~150x). + #[test] + fn test_digamma_negative_arguments_far_from_zero() { + prec::assert_relative_eq!(digamma(-12.72), -0.0169824608177603739173, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(digamma(-14.72), 0.1238386191670285663687, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(digamma(-20.02), 52.95568424702714411621, epsilon = 0.0, max_relative = 1e-15); + } + // TODO: precision testing could be more accurate #[test] fn test_digamma() { @@ -1357,6 +1540,30 @@ mod tests { ); } + /// `gamma` must overflow to `+inf`, never `-inf`, and must not overflow + /// early. The first-order correction in `lanczos_power` is only a valid + /// perturbation while `|p * delta| << 1`; applying it unguarded let a + /// negative correction flip the sign for `x` past ~1e15. Separately, the + /// power alone overflows from about x = 169.7 while the full product stays + /// representable to x ~ 171.61. + #[test] + fn test_gamma_overflow_boundary_and_sign() { + // finite and positive right up to the true overflow point + for x in [168.0f64, 169.0, 169.7, 170.5, 171.0, 171.5, 171.6] { + let g = gamma(x); + assert!(g.is_finite() && g > 0.0, "gamma({x}) should be finite positive, got {g}"); + } + // mpmath at 30 digits + prec::assert_relative_eq!(gamma(171.6), 1.5858969096673029e308, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(gamma(170.5), 5.5620924145599996e305, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(gamma(169.7), 9.155822000376269e303, epsilon = 0.0, max_relative = 1e-13); + // and always +inf beyond it, at every scale + for x in [172.0f64, 200.0, 1e3, 1e6, 1e14, 1e15, 1e100, 1e300, f64::MAX] { + let g = gamma(x); + assert!(g.is_infinite() && g > 0.0, "gamma({x:e}) should be +inf, got {g}"); + } + } + #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/prec.rs b/src/prec.rs index ebf1b6d4..2e6d371d 100644 --- a/src/prec.rs +++ b/src/prec.rs @@ -29,6 +29,7 @@ //! - `DEFAULT_RELATIVE_ACC`: 1e-14 for relative comparisons //! - `DEFAULT_EPS`: 1e-9 for absolute comparisons //! - `DEFAULT_ULPS`: 5 for ULPs comparisons +//! - `DEFAULT_ULPS_EPS`: `f64::EPSILON`, the absolute floor paired with `DEFAULT_ULPS` //! //! These defaults should be used unless there is a specific reason to use different //! precision levels. @@ -62,6 +63,16 @@ pub const DEFAULT_EPS: f64 = 1e-9; /// Default and target ULPs accuracy for f64 operations pub const DEFAULT_ULPS: u32 = 5; +/// Default absolute epsilon for ULPs comparisons. +/// +/// `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it looks +/// at the ULPs distance, so pairing it with [`DEFAULT_EPS`] (`1e-9`) would make +/// the ULPs bound unreachable and turn `ulps_eq!(x, y)` into a `1e-9` absolute +/// comparison. `ulps_eq!` is used inside the crate to recognise exact parameter +/// values (`p == 1.0`, `x == x.floor()`), so its epsilon has to stay at the +/// scale of a genuine rounding error. +pub const DEFAULT_ULPS_EPS: f64 = f64::EPSILON; + /// Compares if two floats are close via `approx::abs_diff_eq` /// using a maximum absolute difference (epsilon) of `acc`. #[deprecated(since = "0.19.0", note = "Use abs_diff_eq! macro instead")] @@ -81,6 +92,55 @@ pub(crate) fn convergence(x: &mut f64, x_new: f64) -> bool { res } +/// Splits the exact rounding error out of the product `a * b` +/// (Dekker's algorithm): `a * b == p + dekker_product_err(a, b, p)` exactly. +/// +/// Used instead of `f64::mul_add`, which falls back to a slow software FMA on +/// targets without the hardware instruction (e.g. baseline x86-64). +#[inline] +pub(crate) fn dekker_product_err(a: f64, b: f64, p: f64) -> f64 { + const SPLIT: f64 = 134_217_729.0; // 2^27 + 1 + // Veltkamp's split below multiplies by `SPLIT`, which overflows to `inf` + // once an argument passes about `1.3e300` and then yields `inf - inf`, i.e. + // NaN. Scaling by a power of two is exact and the residual scales with it, + // so rescale rather than give up: `err(a b) = 2^k err((a 2^-k) b)`. + const BIG: f64 = 1e300; + const DOWN: f64 = 9.313225746154785e-10; // 2^-30, exact + let (mut a, mut b, mut p) = (a, b, p); + let mut scale = 1.0; + if a.abs() > BIG { + a *= DOWN; + p *= DOWN; + scale /= DOWN; + } + if b.abs() > BIG { + b *= DOWN; + p *= DOWN; + scale /= DOWN; + } + let ca = SPLIT * a; + let a_hi = ca - (ca - a); + let a_lo = a - a_hi; + let cb = SPLIT * b; + let b_hi = cb - (cb - b); + let b_lo = b - b_hi; + scale * (((a_hi * b_hi - p) + a_hi * b_lo + a_lo * b_hi) + a_lo * b_lo) +} + +/// Knuth's two-sum for a difference: returns `(s, e)` with +/// `a - b == s + e` exactly, where `s = a - b` rounded. +#[inline] +pub(crate) fn two_diff(a: f64, b: f64) -> (f64, f64) { + let s = a - b; + if !s.is_finite() { + // the residual is not representable; 0 is the safe choice, and the + // alternative is `inf - inf == NaN` poisoning every caller + return (s, 0.0); + } + let v = s - a; + (s, (a - (s - v)) + (-b - v)) +} + macro_rules! redefine_one_opt_approx_macro { ( $approx_macro:ident, @@ -144,7 +204,7 @@ mod macros { ); redefine_two_opt_approx_macro!( ulps_eq, - { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS } + { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS } ); pub(crate) use abs_diff_eq; @@ -162,7 +222,7 @@ mod macros { ); redefine_two_opt_approx_macro!( assert_ulps_eq, - { epsilon: crate::prec::DEFAULT_EPS, max_ulps: crate::prec::DEFAULT_ULPS } + { epsilon: crate::prec::DEFAULT_ULPS_EPS, max_ulps: crate::prec::DEFAULT_ULPS } ); pub(crate) use assert_abs_diff_eq;