From a6a48229aff72e129ae0d4f2cf9443eac41bf09d Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:53:28 -0500 Subject: [PATCH 1/5] fix: stop prec::ulps_eq! degenerating into a 1e-9 absolute comparison `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it consults the ULPs distance. The crate's wrapper paired it with `DEFAULT_EPS = 1e-9`, which made the ULPs bound unreachable, so every internal `ulps_eq!(x, y)` was really a 1e-9 absolute comparison. That matters because the macro is used to recognise *exact* parameter values, so anything within 1e-9 of them silently took a degenerate branch: Binomial(p = 1 - 2^-33, n = 100).pmf(99) 0 (true 1.16e-8) Binomial(p = 1 - 2^-33, n = 100).ln_pmf(99) -inf Geometric(p = 1 - 2^-33).max() 1 (true u64::MAX) Geometric(p = 1 - 2^-33).skewness() inf (true 92681.9) Beta(1 + 2^-33, 1 + 2^-33).pdf(0.0) 1 (true 0) digamma(-1 + 2^-33) -inf (true -8.59e9) Adds `DEFAULT_ULPS_EPS = f64::EPSILON` and uses it for `ulps_eq!` and `assert_ulps_eq!`. The exact values still take the degenerate branches - `Binomial(1.0, 5).pmf(5) == 1`, `Geometric(1.0).max() == 1`, `digamma` at the poles is still -inf - all covered by existing tests. Tests use `1 - 2^-33` rather than `1 - 1e-10` so that `1 - p` is exact and the references are not limited by the representation of `p`. Also fixes `Binomial::entropy`, which summed `-p * ln(p)` unguarded and so returned `NaN` once any mass underflowed to zero (`0 * -inf`) - for `p = 0.5` that is every `n` past about 1100. `Categorical::entropy` already filtered zero terms; this brings `Binomial` in line. --- src/distribution/binomial.rs | 39 ++++++++++++++++++++++++++++++++++- src/distribution/geometric.rs | 12 +++++++++++ src/function/gamma.rs | 19 +++++++++++++++++ src/prec.rs | 15 ++++++++++++-- 4 files changed, 82 insertions(+), 3 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..315c26ab 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); 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/function/gamma.rs b/src/function/gamma.rs index bff5c45a..771a62d1 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -1195,6 +1195,25 @@ 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 + /// was `1e-9` absolute, so a whole neighbourhood around each pole returned + /// -inf instead of a large finite value. + /// + /// The tolerance is loose because `digamma(x < 0)` goes through the + /// reflection formula and `(PI * x).tan()` loses roughly + /// `ulp(PI) / dist_to_pole` of relative accuracy; what is pinned here is + /// that the values are finite and of the right magnitude. + #[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-5); + prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)), 8589934592.4227843348, epsilon = 0.0, max_relative = 1e-5); + // 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); + } + // TODO: precision testing could be more accurate #[test] fn test_digamma() { diff --git a/src/prec.rs b/src/prec.rs index ebf1b6d4..5e23b05b 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")] @@ -144,7 +155,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 +173,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; From c4e9de917a72f28e90a3971f6373c2d109b2bd30 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 10:06:06 -0500 Subject: [PATCH 2/5] fix: evaluate the Lanczos sum in a well-conditioned form `ln_gamma`/`gamma` evaluate Pugh's 10-term Lanczos approximation as the partial-fraction sum `d_0 + sum_k d_k / (z + k - 1)`. The residues alternate in sign, so that sum cancels badly - condition number ~3600 around z = 50, i.e. ~440 eps of relative error in the sum alone. The same quantity as a single fraction `N(z) / D(z)`, derived from those residues in exact rational arithmetic, has *all-positive* coefficients in both numerator and denominator (`D(z) = z (z+1) ... (z+9)`, whose expanded coefficients are unsigned Stirling numbers of the first kind and exact in f64). For z > 0 both Horner evaluations then have condition number 1. This needs no new external constants - they are derived from the residues already in the file - and agrees with the partial-fraction form to 3.2e-17 over [0.5, 3000]. A reversed form in 1/z covers z >= 1e29, where z^10 would overflow. Two further fixes in the same evaluation: * `lanczos_power` compensates the base of `((p + g) / e)^p`. `powf` amplifies a relative error in its base by the exponent, so the two roundings in `(p + g) / e` were the dominant error (~190 ulp by x = 122). The residuals of the addition and the division - the latter against a double-double `e` - are recovered exactly and applied as a first-order correction, which is only valid while it stays small, so it is gated on that. * `gamma`/`ln_gamma` at the positive integers up to 171 come from the exact factorial table, which also makes `Gamma(2) == 1` rather than 1.0000000000000002. * `gamma` for x above ~169.7 halves the exponent and squares, because the power alone overflows there while the full product stays representable to x ~ 171.61. `gamma(171.6)` was `inf`; it is now 1.5858969096673e308. `sin_pi`/`tan_pi` reduce by the period before calling `sin`/`tan`. `x - round(x)` is exact, so the error stays relative to the fractional part instead of being `ulp(PI) / dist_to_pole`, which near the poles cost several decimal digits. Measured against mpmath at 45 dps: before after ln_gamma p99 45.9 ulp 8.2 ulp gamma (positive) 285 median 19 median gamma (negative) 109 median, 1173 max 3.4 median, 22.8 max digamma (negative) 19.5 median, 1.5e10 1.6 median, 538 max `gamma`'s remaining ~19-35 ulp is the approximation floor of the f64-rounded Pugh coefficient set itself (~3.7e-15, confirmed by evaluating the formula in exact arithmetic), so the evaluation is now compensated to well below the fit. Four test expectations move as a consequence, each verified against mpmath: one `Binomial::sf` and two `FisherSnedecor` pdf/ln_pdf literals were fitted to the old output, and the `Dirichlet`/`MultivariateStudent` doctest values were 34 and 22 ulp from truth (the new outputs are within 2). --- src/distribution/binomial.rs | 2 +- src/distribution/dirichlet.rs | 2 +- src/distribution/fisher_snedecor.rs | 8 +- src/distribution/multivariate_students_t.rs | 2 +- src/function/gamma.rs | 302 ++++++++++++++++---- src/prec.rs | 49 ++++ 6 files changed, 301 insertions(+), 64 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index 315c26ab..021f3eba 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -552,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/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/gamma.rs b/src/function/gamma.rs index 771a62d1..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), @@ -1197,23 +1352,32 @@ mod tests { /// `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 - /// was `1e-9` absolute, so a whole neighbourhood around each pole returned - /// -inf instead of a large finite value. - /// - /// The tolerance is loose because `digamma(x < 0)` goes through the - /// reflection formula and `(PI * x).tan()` loses roughly - /// `ulp(PI) / dist_to_pole` of relative accuracy; what is pinned here is - /// that the values are finite and of the right magnitude. + /// 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-5); - prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)), 8589934592.4227843348, epsilon = 0.0, max_relative = 1e-5); + 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() { @@ -1376,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 5e23b05b..2e6d371d 100644 --- a/src/prec.rs +++ b/src/prec.rs @@ -92,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, From 32b661815bda24ca87020aca3ee9bf1b237b1702 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 10:08:03 -0500 Subject: [PATCH 3/5] fix: scale beta_reg's continued-fraction bound with its parameters The Lentz recurrence in `checked_beta_reg` was capped at a fixed 140 iterations and silently returned whatever it had reached. It is slowest at the centre of the distribution, where the worst case over `x` grows like `5 * min(a, b).cbrt()`, so past `min(a, b) ~ 1.5e4` the answer was simply wrong. Against the exact identity `I_{1/2}(a, a) == 1/2`: a = b = 1e5 0.49999969504 relative error 6.1e-7 a = b = 1e6 0.49121972700 1.8e-2 a = b = 1e7 0.21285001452 5.7e-1 `Binomial::new(0.5, 2e6).cdf(1e6)` returned 0.4916 against a true 0.50028. The bound now scales as `8 * min(a, b).cbrt()`, measured to cover the worst case over `x` with headroom, clamped to bound the work at roughly 10 ms. The loop still exits on convergence, so ordinary calls are unaffected: `beta_reg(2.5, 2.5, 0.5)` is 150 ns against 145 ns, and `beta_reg(1, 1, 0.3)` is unchanged. Three robustness fixes alongside it, all found by probing the parameter space rather than by sweeping accuracy: * an underflowed prefix now short-circuits to the corresponding endpoint. The result is `bt * h / a` with `h` of order one, so this is exact - and it avoids forming `0.0 * h`, which is NaN whenever the recurrence overflowed (`beta_reg(1e300, 1e-300, 0.5)` was NaN on both sides of this change before). * `x == 0` and `x == 1` return 0 and 1 directly. They used to depend on the symmetry test, which mapped `x == 0` to `1.0` once `a + b` overflowed. * the symmetry threshold `(a + 1) / (a + b + 2)` is computed scaled when `a + b` overflows, since otherwise it collapses to zero and sends every `x` down the transformed branch. * the result is kept in `[0, 1]`, and falls back to the concentrated-limit step function if the truncated recurrence produced something non-finite. `I_x` is a probability and callers such as `Binomial::cdf` are contractually so; `a = b = 1e20` previously returned -2.56. Regression tests use two exact identities that need no reference data - `I_{1/2}(a, a) == 1/2` and `I_x(a, b) + I_{1-x}(b, a) == 1` - plus a 12x12x6 parameter grid asserting the result stays a finite probability. Both identity tests fail on the old fixed bound. --- src/function/beta.rs | 162 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 152 insertions(+), 10 deletions(-) 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() {} From cc54f46f6dfbeb18ea2cdd42bd4dce70080fc349 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 10:08:52 -0500 Subject: [PATCH 4/5] fix: use Loader's saddle-point form for the gamma/beta prefixes and pmfs Densities and incomplete-function prefixes in this family are written as `exp(a ln x - x - ln_gamma(a))` and friends. Every term there grows like `a ln a` while the sum stays `O(ln a)`, so at `a = 1e4` each term is ~1e5 and the ~1e-11 absolute error of the sum becomes the accuracy ceiling of everything downstream. That is what limited `beta_reg` to ~6e-11 even where its recurrence converged. Adds the two building blocks from Loader, "Fast and Accurate Computation of Binomial Probabilities" (2000): * `stirling_delta(z)`, the Stirling-series remainder of `ln_gamma`, which is `O(1/(12z))`. It is table-free: the recurrence `d(z) = d(z+1) + (z + 1/2) ln(1 + 1/z) - 1` lifts any argument into the range where a 6-term series is good to 1.4e-18, so there is no interpolation table and no recursive dependency on `ln_gamma`. Accurate to ~1e-15 *absolute*, which is what matters since every caller adds it to an exponent. * `bd0(x, np) = x ln(x/np) - (x - np)`, which has a double root at `x == np`. Evaluated by a series in `(x - np) / (x + np)` so the cancellation is done analytically, giving ~1e-16 relative accuracy uniformly through the root. Rewriting the prefixes in terms of these keeps every piece `O(1)`. Measured against mpmath, and against `I_{1/2}(a, a) == 1/2` which is exact: beta_reg, a = b = 1e4 6e-11 -> 8.7e-15 beta_reg, a = b = 1e6 3e-9 -> 2.6e-13 gamma_ur, a >= 16 5.3 median ulp -> 0.38, p99 970k -> 3071 Binomial::pmf 685k median ulp -> 2.45 Poisson::pmf 649 median ulp -> 3.16 Both prefixes fall back to the direct form for small parameters, where there is no cancellation left to remove and `stirling_delta`'s recurrence would cost more roundings than it saves - without that gate, `FisherSnedecor(1,1).cdf` got *worse*, since a = b = 0.5 needs 16 recurrence steps. `bd0` is sensitive enough to its mean argument that the rounding of `n * p` alone matters: `d bd0 / d np = 1 - x / np`, so half an ulp of `np` at 6e5 moves the result 2.5e-13, which was a thousand ulps of the resulting pmf. `n`, `1 - p` and the products are therefore carried as double-doubles via `bd0_dd`. R's `dbinom` has the same limitation. Also guards `bd0`'s direct branch, where `(x / np).ln()` silently loses the whole term when the ratio leaves the normal range - `bd0(1e-300, 5e99)` gave `-inf` where the answer is `+5e99`, which propagated into `beta_reg` as NaN for extreme parameter ratios, and as a silent `1.0` where the true value was ~1e-30118. --- src/distribution/binomial.rs | 64 +++++++-- src/distribution/poisson.rs | 34 ++++- src/function/beta.rs | 103 ++++++++++++--- src/function/gamma.rs | 245 ++++++++++++++++++++++++++++++++++- src/prec.rs | 13 ++ 5 files changed, 432 insertions(+), 27 deletions(-) diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index 021f3eba..f605668c 100644 --- a/src/distribution/binomial.rs +++ b/src/distribution/binomial.rs @@ -1,5 +1,5 @@ use crate::distribution::{Discrete, DiscreteCDF}; -use crate::function::{beta, factorial}; +use crate::function::{beta, factorial, gamma}; use crate::prec; use crate::statistics::*; use core::f64; @@ -309,10 +309,7 @@ impl Discrete for Binomial { } else if prec::ulps_eq!(self.p, 1.0) { if x == self.n { 1.0 } else { 0.0 } } else { - (factorial::ln_binomial(self.n, x) - + x as f64 * self.p.ln() - + (self.n - x) as f64 * (1.0 - self.p).ln()) - .exp() + self.ln_pmf(x).exp() } } @@ -332,9 +329,38 @@ impl Discrete for Binomial { } else if prec::ulps_eq!(self.p, 1.0) { if x == self.n { 0.0 } else { f64::NEG_INFINITY } } else { - factorial::ln_binomial(self.n, x) - + x as f64 * self.p.ln() - + (self.n - x) as f64 * (1.0 - self.p).ln() + let n = self.n as f64; + let k = x as f64; + let q = 1.0 - self.p; + if x == 0 { + return n * q.ln(); + } + if x == self.n { + return n * self.p.ln(); + } + if n < gamma::STIRLING_SERIES_MIN { + // `ln_binomial` is exact from the factorial table here and all + // terms are small, so the direct form loses nothing + return factorial::ln_binomial(self.n, x) + k * self.p.ln() + (n - k) * q.ln(); + } + // Saddle-point form (Loader 2000): the direct expression differences + // terms growing like `n ln n` to produce an `O(1)` result, which + // cost ~2e-10 relative at n = 1e5. + let nk = n - k; + // `bd0` is sensitive to the last half-ulp of its mean argument, and + // both `1 - p` and the products round, so they are carried as + // double-doubles; see `gamma::bd0_dd`. + let (q_hi, q_lo) = prec::two_diff(1.0, self.p); + let np = n * self.p; + let np_lo = prec::dekker_product_err(n, self.p, np); + let nq = n * q_hi; + let nq_lo = prec::dekker_product_err(n, q_hi, nq) + n * q_lo; + gamma::stirling_delta(n) + - gamma::stirling_delta(k) + - gamma::stirling_delta(nk) + - gamma::bd0_dd(k, np, np_lo) + - gamma::bd0_dd(nk, nq, nq_lo) + + 0.5 * (n / (f64::consts::TAU * k * nk)).ln() } } } @@ -449,6 +475,28 @@ mod tests { assert!(Binomial::new(0.1, 5_000).unwrap().entropy().unwrap().is_finite()); } + /// Large-`n` pmf: the saddle-point form (`gamma::bd0` / + /// `gamma::stirling_delta`, with double-double means) replaced + /// `exp(ln_binomial + x ln p + (n-x) ln q)`, whose terms each grow like + /// `n ln n` while the result stays `O(1)`. References are mpmath at 40 + /// significant digits; the old form was ~7e-11 relative here. + #[test] + fn test_pmf_large_n_saddle_point() { + let pmf = |arg: u64| move |x: Binomial| x.pmf(arg); + test_relative(0.5, 100000, 0.0025231262141967398855, pmf(50000)); + test_relative(0.3, 2000000, 0.00061558120658465376062, pmf(600000)); + // ln_pmf agrees with ln(pmf) where the pmf is comfortably representable + for (p, n, k) in [(0.5f64, 100000u64, 50000u64), (0.3, 2000000, 600000)] { + let d = create_ok(p, n); + prec::assert_relative_eq!(d.ln_pmf(k), d.pmf(k).ln(), epsilon = 0.0, max_relative = 1e-14); + } + // the endpoints reduce to q^n and p^n; `exp(n ln q)` carries ~1e-15 + // relative (|n ln q| ~ 11.5 here), which is the honest bound + let d = create_ok(0.25, 40); + prec::assert_relative_eq!(d.pmf(0), 1.0056585161637497e-5, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(d.pmf(40), 8.271806125530277e-25, epsilon = 0.0, max_relative = 1e-14); + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Binomial| x.pmf(arg); diff --git a/src/distribution/poisson.rs b/src/distribution/poisson.rs index fc9421ff..2a9df0ff 100644 --- a/src/distribution/poisson.rs +++ b/src/distribution/poisson.rs @@ -276,7 +276,7 @@ impl Discrete for Poisson { /// /// where `λ` is the rate fn pmf(&self, x: u64) -> f64 { - (-self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)).exp() + self.ln_pmf(x).exp() } /// Calculates the log probability mass function for the poisson @@ -291,7 +291,19 @@ impl Discrete for Poisson { /// /// where `λ` is the rate fn ln_pmf(&self, x: u64) -> f64 { - -self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x) + let k = x as f64; + if x == 0 { + return -self.lambda; + } + if k < gamma::STIRLING_SERIES_MIN { + // `ln_factorial` is exact from its table here and every term is + // small, so the direct form loses nothing + return -self.lambda + k * self.lambda.ln() - factorial::ln_factorial(x); + } + // Saddle-point form (Loader 2000): the direct expression is a + // difference of terms that each grow like `x ln x` while the result + // stays `O(1)`, which cost ~1e-11 relative at lambda = 1e4. + -gamma::stirling_delta(k) - gamma::bd0(k, self.lambda) - 0.5 * (f64::consts::TAU * k).ln() } } @@ -421,6 +433,24 @@ mod tests { test_exact(10.8, u64::MAX, max); } + /// Large-lambda pmf: the saddle-point form (`gamma::bd0` / + /// `gamma::stirling_delta`) replaced `exp(-lam + x ln lam - ln x!)`, whose + /// terms each grow like `x ln x` while the result stays `O(1)`. References + /// are mpmath at 40 significant digits; the old form was ~650 ulp median and + /// 1.4e7 ulp worst over this range. + #[test] + fn test_pmf_large_lambda_saddle_point() { + let pmf = |arg: u64| move |x: Poisson| x.pmf(arg); + test_relative(10000.0, 0.0039893895589628256487, pmf(10000)); + test_relative(1000000.0, 0.0003989422471562440297, pmf(1000000)); + test_relative(1000000.0, 0.00024189010120174141723, pmf(1001000)); + // ln_pmf agrees with ln(pmf) where the pmf is comfortably representable + for (lam, k) in [(10000.0f64, 10000u64), (1e6, 1001000)] { + let d = create_ok(lam); + crate::prec::assert_relative_eq!(d.ln_pmf(k), d.pmf(k).ln(), epsilon = 0.0, max_relative = 1e-14); + } + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Poisson| x.pmf(arg); diff --git a/src/function/beta.rs b/src/function/beta.rs index 3c1ff6ed..7e36132b 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -132,6 +132,51 @@ pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 { checked_beta_reg(a, b, x).unwrap() } +/// `ln(x^a (1-x)^b / Beta(a, b))`, the prefix of the incomplete beta functions. +/// +/// Written out, this is +/// `ln_gamma(a+b) - ln_gamma(a) - ln_gamma(b) + a ln x + b ln(1-x)`, where the +/// terms grow like `(a + b) ln(a + b)` while the sum stays `O(ln(a + b))`. With +/// `n = a + b` the saddle-point form is +/// +/// ```text +/// = -bd0(a, n x) - bd0(b, n (1-x)) +/// + stirling_delta(n) - stirling_delta(a) - stirling_delta(b) +/// + ln(a b / (2 pi n)) / 2 +/// ``` +/// +/// Every piece is `O(1)` (the two `bd0` terms grow only as the true `-ln` of +/// the prefix does), so this is accurate to a few `1e-16` absolute where the +/// old form lost `~1e-10` at `a = b = 1e4`. +/// +/// For small parameters the direct form is kept: there is no cancellation left +/// to remove (every term is already `O(1)`), while `stirling_delta` would need +/// up to 16 recurrence steps per argument and so contributes more rounding than +/// it saves. +fn ln_beta_prefix(a: f64, b: f64, x: f64) -> f64 { + let y = 1.0 - x; + if a.max(b) < gamma::STIRLING_SERIES_MIN { + return gamma::ln_gamma(a + b) - gamma::ln_gamma(a) - gamma::ln_gamma(b) + + a * x.ln() + + b * y.ln(); + } + // `bd0` is sensitive to the last half-ulp of its mean argument, so `n` and + // the two products are carried as double-doubles; see `gamma::bd0_dd`. + let (n, n_lo) = prec::two_sum(a, b); + let (y_hi, y_lo) = prec::two_diff(1.0, x); + let nx = n * x; + let nx_lo = prec::dekker_product_err(n, x, nx) + n_lo * x; + let ny = n * y_hi; + let ny_lo = prec::dekker_product_err(n, y_hi, ny) + n * y_lo + n_lo * y_hi; + + // `a / n` and `b / TAU` are each individually representable, unlike + // `a * b / (2 pi n)`, which can overflow for large parameters + let scale = 0.5 * ((a / n).ln() + (b / f64::consts::TAU).ln()); + scale - gamma::bd0_dd(a, nx, nx_lo) - gamma::bd0_dd(b, ny, ny_lo) + gamma::stirling_delta(n) + - gamma::stirling_delta(a) + - gamma::stirling_delta(b) +} + /// Computes the regularized lower incomplete beta function /// `I_x(a,b) = 1/Beta(a,b) * int(t^(a-1)*(1-t)^(b-1), t=0..x)` /// `a > 0`, `b > 0`, `1 >= x >= 0` where `a` is the first beta parameter, @@ -140,15 +185,26 @@ pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 { /// /// # 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`. +/// The leading factor is evaluated by the saddle-point decomposition in +/// [`ln_beta_prefix`], which keeps every intermediate `O(1)` however large +/// `a + b` becomes. Measured against the exact identity `I_{1/2}(a, a) == 1/2`, +/// the relative error is: /// -/// 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. +/// ```text +/// min(a, b) <= 1e7 2e-13 (1e-15 at 1e4 and below) +/// min(a, b) <= 1e13 3e-10 +/// min(a, b) <= 1e16 3e-7 +/// beyond unreliable; the result is still clamped to [0, 1] +/// ``` +/// +/// The degradation past `1e7` is the continued fraction accumulating its own +/// rounding over millions of iterations, not the prefix. Evaluation is bounded +/// at roughly 10 ms in the worst case. +/// +/// Values deep in the tails carry more *relative* error - a few hundred ulp, +/// since the prefix is recovered as `exp(ln prefix)` and `|ln prefix|` is large +/// there - but their absolute error stays far below the smallest value of +/// interest. /// /// # Errors /// @@ -179,10 +235,7 @@ pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { let bt = if x == 0.0 || crate::prec::ulps_eq!(x, 1.0, epsilon = MODULE_EPS) { 0.0 } else { - (gamma::ln_gamma(a + b) - gamma::ln_gamma(a) - gamma::ln_gamma(b) - + a * x.ln() - + b * (1.0 - x).ln()) - .exp() + ln_beta_prefix(a, b, x).exp() }; let symm_transform = { let denom = a + b + 2.0; @@ -733,11 +786,31 @@ mod tests { assert!(checked_beta_reg(1.0, 1.0, 2.0).is_err()); } + /// `I_{1/2}(a, a) == 1/2` exactly, by the symmetry of the `Beta(a, a)` + /// density. That point is also where the continued fraction converges most + /// slowly, so it pins the iteration bound. With the old fixed bound of 140 + /// iterations this returned 0.49999969 at `a = 1e5`, 0.49121973 at `a = 1e6` + /// and 0.21285001 at `a = 1e7`. + #[test] + fn test_beta_reg_symmetric_midpoint_large_parameters() { + // The saddle-point prefix (`ln_beta_prefix`) took this from `6e-7` at + // `a = 1e5` and `5.7e-1` at `a = 1e7` down to the `1e-13` level. + for a in [1e2, 1e3, 1e4, 1e5, 1e6, 1e7] { + prec::assert_relative_eq!( + beta_reg(a, a, 0.5), + 0.5, + epsilon = 0.0, + max_relative = 1e-12 + ); + } + } + /// `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`. + /// over/underflow the intermediate quantities. Before the `bd0` ratio guard + /// and the `[0, 1]` clamp, this grid produced NaN (`a = 1e300, b = 1e-300`), + /// a silent `1.0` where the answer was `0` (`a = 1e100, b = 1e-300`), and + /// `-2.56` (`a = b = 1e20`). #[test] fn test_beta_reg_extreme_parameters_stay_a_probability() { let params = [ diff --git a/src/function/gamma.rs b/src/function/gamma.rs index 13d2f126..c38bb9a6 100644 --- a/src/function/gamma.rs +++ b/src/function/gamma.rs @@ -62,6 +62,159 @@ fn tan_pi(x: f64) -> f64 { (f64::consts::PI * r).tan() } +// --------------------------------------------------------------------------- +// Saddle-point building blocks (Loader, "Fast and Accurate Computation of +// Binomial Probabilities", 2000). +// +// Densities and incomplete-function prefixes in the gamma/beta family are +// naturally written as `exp(a * ln x - x - ln_gamma(a))` and friends. Those +// forms are numerically poor: for `a ~ 1e4` each term is `~1e5` while the sum +// is `O(1)`, so the `~1e-11` absolute error of the sum becomes the accuracy +// ceiling of everything downstream. +// +// The fix is to split each such exponent into two pieces that are individually +// `O(1)`: +// +// * `stirling_delta(z)`, the Stirling-series remainder of `ln_gamma`, and +// * `bd0(x, np)`, the "deviance" `x * ln(x / np) - (x - np)`, which is +// evaluated by a series in `(x - np) / (x + np)` so that the cancellation +// for `x` near `np` is performed analytically rather than in floating +// point. +// +// Nothing large is ever formed, so the cancellation disappears entirely. +// --------------------------------------------------------------------------- + +/// Coefficients `B_2n / (2n (2n - 1))` of the Stirling series for +/// [`stirling_delta`], ascending in `1/z^2`. +const STIRLING_SERIES: &[f64] = &[ + 0.08333333333333333, // 1/12 + -0.002777777777777778, // -1/360 + 0.0007936507936507937, // 1/1260 + -0.0005952380952380953, // -1/1680 + 0.0008417508417508417, // 1/1188 + -0.0019175269175269176, // -691/360360 +]; + +/// Argument above which the truncated Stirling series is used directly. Six +/// terms there are good to `1.4e-18` absolute, i.e. invisible. +pub(crate) const STIRLING_SERIES_MIN: f64 = 16.0; + +/// The Stirling-series remainder of `ln_gamma`, for `z > 0`: +/// +/// ```text +/// stirling_delta(z) = ln_gamma(z) - [(z - 1/2) ln z - z + ln(2 pi) / 2] +/// ``` +/// +/// Equivalently `ln(n!) - ln(sqrt(2 pi n) (n / e)^n)` at `z = n`, which is +/// Loader's `stirlerr`. The value is `O(1 / (12 z))`, and the implementation +/// keeps it accurate to `~1e-15` *absolute* - which is what matters, since +/// every caller adds it to an exponent. +/// +/// Below [`STIRLING_SERIES_MIN`] the recurrence +/// `delta(z) = delta(z + 1) + (z + 1/2) ln(1 + 1/z) - 1` lifts the argument +/// into the series range. Each step costs one rounding, so no table (and no +/// recursive dependency on [`ln_gamma`]) is needed. +pub(crate) fn stirling_delta(z: f64) -> f64 { + let mut acc = 0.0; + let mut w = z; + while w < STIRLING_SERIES_MIN { + acc += (w + 0.5) * (1.0 / w).ln_1p() - 1.0; + w += 1.0; + } + // Horner in 1/z^2. For `w * w == inf` this collapses to the leading + // `1 / (12 w)` term, which is the correct limit. + let ww = w * w; + let series = STIRLING_SERIES.iter().rev().fold(0.0, |s, &c| s / ww + c); + acc + series / w +} + +/// The deviance `bd0(x, np) = x * ln(x / np) - (x - np)`, for `x >= 0`, +/// `np > 0`. +/// +/// This is the Kullback-Leibler-like term of the Poisson/binomial saddle-point +/// expansions, and is non-negative with a double root at `x == np`. Writing it +/// directly loses all precision near that root; with +/// `v = (x - np) / (x + np)` it has the all-positive series +/// +/// ```text +/// bd0 = (x - np) v + 2 x sum_{j >= 1} v^(2j + 1) / (2j + 1) +/// ``` +/// +/// which is used whenever `|x - np| < (x + np) / 10`, giving `~1e-16` relative +/// accuracy uniformly. +pub(crate) fn bd0(x: f64, np: f64) -> f64 { + if x == 0.0 { + // 0 * ln 0 == 0 by convention; the direct form would give NaN + return np; + } + if (x - np).abs() < 0.1 * (x + np) { + let v = (x - np) / (x + np); + let mut s = (x - np) * v; + if s.abs() < f64::MIN_POSITIVE { + return s; + } + let mut ej = 2.0 * x * v; + let v2 = v * v; + // |v| < 0.1, so v^(2j) is below f64::MIN_POSITIVE well before j = 200 + for j in 1..200 { + ej *= v2; + let s1 = s + ej / (2 * j + 1) as f64; + if s1 == s { + return s1; + } + s = s1; + } + return s; + } + // `(x / np).ln()` silently loses the whole term when the ratio leaves the + // normal range - `bd0(1e-300, 5e99)` has `x / np` underflow to zero, giving + // `-inf` where the answer is `+5e99`, which then poisoned the beta prefix + // into NaN. Reaching this branch requires `|x - np| >= (x + np) / 10`, so + // the ratio is bounded away from 1 and the difference of logs cannot cancel + // badly; it is only used where the ratio itself is unusable, since it is + // otherwise the less accurate of the two forms. + let ratio = x / np; + let log_ratio = if (f64::MIN_POSITIVE..f64::INFINITY).contains(&ratio) { + ratio.ln() + } else { + x.ln() - np.ln() + }; + x * log_ratio - (x - np) +} + +/// [`bd0`] with the mean supplied as an exact double-double `np_hi + np_lo`. +/// +/// The mean usually arrives as a rounded product such as `n * p`, and `bd0` is +/// sensitive enough for that last half-ulp to dominate everything else: +/// `d bd0 / d np = 1 - x / np`, so at `n = 2e6`, `p = 0.3` the rounding of +/// `n * p` alone shifts `bd0` by `~2.5e-13` - a thousand ulps in the resulting +/// pmf. The correction is first order in `np_lo`; the next term is smaller by +/// another factor of `np_lo / np`, i.e. utterly negligible. +pub(crate) fn bd0_dd(x: f64, np_hi: f64, np_lo: f64) -> f64 { + bd0(x, np_hi) + (1.0 - x / np_hi) * np_lo +} + +/// `ln(x^a e^-x / Gamma(a))`, the prefix of the incomplete gamma functions. +/// +/// Written out, this is `a ln x - x - ln_gamma(a)`, where all three terms grow +/// like `a ln a` while the sum stays `O(ln a)`. The saddle-point form +/// +/// ```text +/// = -bd0(a, x) - stirling_delta(a) + ln(a / (2 pi)) / 2 +/// ``` +/// +/// is the same quantity with every piece `O(1)`, so it is accurate to a few +/// `1e-16` absolute instead of `a ln a * eps` (`~1e-11` at `a = 1e4`). +pub(crate) fn ln_gamma_prefix(a: f64, x: f64) -> f64 { + if a < STIRLING_SERIES_MIN { + // No cancellation to remove yet (every term is already `O(1)`), and the + // recurrence in `stirling_delta` would cost more roundings than it + // saves. + return a * x.ln() - x - ln_gamma(a); + } + -bd0(a, x) - stirling_delta(a) + 0.5 * (a / f64::consts::TAU).ln() +} + /// Auxiliary variable when evaluating the `gamma_ln` function const GAMMA_R: f64 = 10.900511; @@ -353,7 +506,8 @@ pub fn checked_gamma_ur(a: f64, x: f64) -> Result { return Ok(1.0 - gamma_lr(a, x)); } - let mut ax = a * x.ln() - x - ln_gamma(a); + // saddle-point form; see `ln_gamma_prefix` + let mut ax = ln_gamma_prefix(a, x); if ax < -709.78271289338399 { return if a < x { Ok(0.0) } else { Ok(1.0) }; } @@ -450,7 +604,8 @@ pub fn checked_gamma_lr(a: f64, x: f64) -> Result { return Ok(0.0); } - let ax = a * x.ln() - x - ln_gamma(a); + // saddle-point form; see `ln_gamma_prefix` + let ax = ln_gamma_prefix(a, x); if ax < -709.78271289338399 { if a < x { return Ok(1.0); @@ -1564,6 +1719,92 @@ mod tests { } } + /// `bd0`'s direct branch formed `(x / np).ln()`, which underflows to + /// `ln(0) = -inf` when `x` is tiny and `np` huge. That propagated through + /// the beta prefix and made `beta_reg` return NaN (or silently 1.0 instead + /// of 0.0) for extreme parameter ratios. + #[test] + fn test_bd0_extreme_ratios() { + // x / np underflows; true value is ~np + let v = super::bd0(1e-300, 5e99); + prec::assert_relative_eq!(v, 5e99, epsilon = 0.0, max_relative = 1e-14); + // x / np overflows; true value is x*ln(x/np) - x + np + let v = super::bd0(5e99, 1e-300); + assert!(v.is_finite() && v > 0.0, "bd0(5e99, 1e-300) = {v}"); + // and it stays non-negative and finite across a wide ratio sweep + for ea in [-300i32, -100, -10, 0, 10, 100, 300] { + for eb in [-300i32, -100, -10, 0, 10, 100, 300] { + let (x, np) = (10f64.powi(ea), 10f64.powi(eb)); + let v = super::bd0(x, np); + assert!(v.is_finite() && v >= 0.0, "bd0(1e{ea}, 1e{eb}) = {v}"); + } + } + } + + /// `stirling_delta` lands in an exponent, so its *absolute* accuracy is what + /// matters. References are mpmath at 40 significant digits; the tolerances + /// sit just above the measured error (`~1e-15` below the series threshold, + /// essentially exact above it). + #[test] + fn test_stirling_delta() { + // below the series threshold the recurrence costs one rounding per step, + // so the bound is absolute; at or above it the series is limited only by + // its own evaluation, so a relative bound is the meaningful one + for (z, want) in [ + (0.5, 0.15342640972002734529), + (1.0, 0.08106146679532725822), + (2.5, 0.033162873519936287485), + (8.0, 0.010411265261972096497), + (15.5, 0.0053755990329268344936), + ] { + prec::assert_abs_diff_eq!(super::stirling_delta(z), want, epsilon = 2e-15); + } + for (z, want) in [ + (16.0, 0.0052076559196096404407), + (100.0, 0.00083333055563491468338), + (10000.0, 8.3333333305555555635e-6), + ] { + prec::assert_relative_eq!( + super::stirling_delta(z), + want, + epsilon = 0.0, + max_relative = 1e-15 + ); + } + // consistency with the defining identity, away from the recurrence + for z in [20.0f64, 63.5, 500.0] { + let want = ln_gamma(z) - ((z - 0.5) * z.ln() - z + 0.5 * f64::consts::TAU.ln()); + prec::assert_abs_diff_eq!(super::stirling_delta(z), want, epsilon = 1e-13); + } + } + + /// `bd0` must stay *relatively* accurate right through its double root at + /// `x == np`, which is the whole point of the series form. + #[test] + fn test_bd0() { + for (x, np, want) in [ + (10000.0, 10000.0, 0.0), + (10000.0, 10001.0, 0.000049996666916646668333), + (10000.0, 15000.0, 945.34891891835618022), + (1.0, 2.0, 0.30685281944005469058), + (600123.0, 600000.0, 0.012606638575794171215), + ] { + let got = super::bd0(x, np); + if want == 0.0 { + assert_eq!(got, 0.0); + } else { + prec::assert_relative_eq!(got, want, epsilon = 0.0, max_relative = 1e-14); + } + } + // non-negative with a root only at x == np, and 0 * ln 0 == 0 at x == 0 + assert_eq!(super::bd0(0.0, 7.5), 7.5); + for i in 1..200 { + let x = 1000.0 + i as f64; + assert!(super::bd0(x, 1000.0) > 0.0); + assert!(super::bd0(1000.0, x) > 0.0); + } + } + #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/prec.rs b/src/prec.rs index 2e6d371d..c174598c 100644 --- a/src/prec.rs +++ b/src/prec.rs @@ -141,6 +141,19 @@ pub(crate) fn two_diff(a: f64, b: f64) -> (f64, f64) { (s, (a - (s - v)) + (-b - v)) } +/// Knuth's two-sum: returns `(s, e)` with `a + b == s + e` exactly, where +/// `s = a + b` rounded. +#[inline] +pub(crate) fn two_sum(a: f64, b: f64) -> (f64, f64) { + let s = a + b; + if !s.is_finite() { + // see `two_diff` + return (s, 0.0); + } + let v = s - a; + (s, (a - (s - v)) + (b - v)) +} + macro_rules! redefine_one_opt_approx_macro { ( $approx_macro:ident, From 8833a6f7c0f6d72ec3ea9d76ed278f8bff39c457 Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 18:52:57 -0500 Subject: [PATCH 5/5] docs: stop linking public docs to the private ln_beta_prefix `checked_beta_reg` is public but `ln_beta_prefix` is not, so the intra-doc link resolves to nothing for users and trips `rustdoc::private_intra_doc_links` under `RUSTDOCFLAGS=-D warnings`. Kept as plain code text, since naming the helper is still useful to a reader of the source. Not a CI failure: no workflow runs `cargo doc`. Co-Authored-By: Claude Opus 5 (1M context) --- src/function/beta.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/function/beta.rs b/src/function/beta.rs index 7e36132b..09c0a138 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -186,7 +186,7 @@ fn ln_beta_prefix(a: f64, b: f64, x: f64) -> f64 { /// # Remarks /// /// The leading factor is evaluated by the saddle-point decomposition in -/// [`ln_beta_prefix`], which keeps every intermediate `O(1)` however large +/// `ln_beta_prefix`, which keeps every intermediate `O(1)` however large /// `a + b` becomes. Measured against the exact identity `I_{1/2}(a, a) == 1/2`, /// the relative error is: ///