From a6a48229aff72e129ae0d4f2cf9443eac41bf09d Mon Sep 17 00:00:00 2001 From: Felix Agene Date: Sun, 26 Jul 2026 09:53:28 -0500 Subject: [PATCH 1/3] 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/3] 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/3] 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() {}