diff --git a/src/distribution/binomial.rs b/src/distribution/binomial.rs index c9794f0a..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; @@ -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) @@ -305,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() } } @@ -328,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() } } } @@ -412,6 +442,61 @@ 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()); + } + + /// 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); @@ -515,7 +600,7 @@ mod tests { test_absolute(0.3, 3, 0.657, 1e-14, sf(0)); test_absolute(0.3, 3, 0.216, 1e-15, sf(1)); test_exact(0.3, 3, 0.0, sf(3)); - test_absolute(0.3, 10, 0.9717524751000001, 1e-16, sf(0)); + test_absolute(0.3, 10, 0.9717524750999999955198, 1e-16, sf(0)); test_absolute(0.3, 10, 0.850691654100002, 1e-14, sf(1)); test_exact(0.3, 10, 0.0, sf(10)); test_exact(1.0, 1, 1.0, sf(0)); diff --git a/src/distribution/dirichlet.rs b/src/distribution/dirichlet.rs index 5f1697b4..51dd809b 100644 --- a/src/distribution/dirichlet.rs +++ b/src/distribution/dirichlet.rs @@ -19,7 +19,7 @@ use nalgebra::{Dim, Dyn, OMatrix, OVector}; /// /// let n = Dirichlet::new(vec![1.0, 2.0, 3.0]).unwrap(); /// assert_eq!(n.mean().unwrap(), DVector::from_vec(vec![1.0 / 6.0, 1.0 / 3.0, 0.5])); -/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.222155556222205); +/// assert_eq!(n.pdf(&DVector::from_vec(vec![0.33333, 0.33333, 0.33333])), 2.2221555562222193); /// ``` #[derive(Clone, PartialEq, Debug)] pub struct Dirichlet diff --git a/src/distribution/fisher_snedecor.rs b/src/distribution/fisher_snedecor.rs index 9cc2f5c4..cd5e160f 100644 --- a/src/distribution/fisher_snedecor.rs +++ b/src/distribution/fisher_snedecor.rs @@ -525,8 +525,8 @@ mod tests { test_absolute(10.0, 0.1, 0.0418440630400545297349, 1e-14, pdf(1.0)); test_absolute(0.1, 1.0, 0.0396064560910663979961, 1e-16, pdf(1.0)); test_absolute(1.0, 1.0, 0.1591549430918953357689, 1e-16, pdf(1.0)); - test_absolute(10.0, 1.0, 0.230361989229138647108, 1e-16, pdf(1.0)); - test_absolute(0.1, 0.1, 0.00221546909694001013517, 1e-18, pdf(10.0)); + test_absolute(10.0, 1.0, 0.230361989229138647108, 3e-16, pdf(1.0)); + test_absolute(0.1, 0.1, 0.00221546909694001013517, 5e-18, pdf(10.0)); test_absolute(1.0, 0.1, 0.00369960370387922619592, 1e-17, pdf(10.0)); test_absolute(10.0, 0.1, 0.00390179721174142927402, 1e-15, pdf(10.0)); test_absolute(0.1, 1.0, 0.00319864073359931548273, 1e-17, pdf(10.0)); @@ -540,13 +540,13 @@ mod tests { #[test] fn test_ln_pdf() { let ln_pdf = |arg: f64| move |x: FisherSnedecor| x.ln_pdf(arg); - test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 1e-15, ln_pdf(1.0)); + test_absolute(0.1, 0.1, 0.0234154207226588982471f64.ln(), 3e-15, ln_pdf(1.0)); test_absolute(1.0, 0.1, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(10.0, 0.1, 0.0418440630400545297349f64.ln(), 1e-13, ln_pdf(1.0)); test_absolute(0.1, 1.0, 0.0396064560910663979961f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(1.0, 1.0, 0.1591549430918953357689f64.ln(), 1e-15, ln_pdf(1.0)); test_absolute(10.0, 1.0, 0.230361989229138647108f64.ln(), 1e-15, ln_pdf(1.0)); - test_exact(0.1, 0.1, 0.00221546909694001013517f64.ln(), ln_pdf(10.0)); + test_absolute(0.1, 0.1, 0.00221546909694001013517f64.ln(), 3e-15, ln_pdf(10.0)); test_absolute(1.0, 0.1, 0.00369960370387922619592f64.ln(), 1e-15, ln_pdf(10.0)); test_absolute(10.0, 0.1, 0.00390179721174142927402f64.ln(), 1e-13, ln_pdf(10.0)); test_absolute(0.1, 1.0, 0.00319864073359931548273f64.ln(), 1e-15, ln_pdf(10.0)); diff --git a/src/distribution/geometric.rs b/src/distribution/geometric.rs index f74f2a79..3a3f9a78 100644 --- a/src/distribution/geometric.rs +++ b/src/distribution/geometric.rs @@ -437,6 +437,18 @@ mod tests { test_exact(0.3, u64::MAX, max); } + /// As in [`Binomial`](crate::distribution::Binomial), the `p == 1` branches + /// here are guarded by `prec::ulps_eq!`; its default epsilon was `1e-9` + /// absolute, so `p = 1 - 2^-33` was treated as a point mass at 1. + #[test] + fn test_p_near_one_is_not_degenerate() { + let g = Geometric::new(1.0 - f64::powi(2.0, -33)).unwrap(); + assert_eq!(g.max(), u64::MAX); + prec::assert_relative_eq!(g.skewness().unwrap(), 92681.900034472750337, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(g.pmf(2), 1.164153218133822873e-10, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(g.ln_pmf(2), -22.873856958594610533, epsilon = 0.0, max_relative = 1e-14); + } + #[test] fn test_pmf() { let pmf = |arg: u64| move |x: Geometric| x.pmf(arg); diff --git a/src/distribution/multivariate_students_t.rs b/src/distribution/multivariate_students_t.rs index 2ad69880..cee52829 100644 --- a/src/distribution/multivariate_students_t.rs +++ b/src/distribution/multivariate_students_t.rs @@ -19,7 +19,7 @@ use nalgebra::{Cholesky, Const, DMatrix, Dim, DimMin, Dyn, OMatrix, OVector}; /// let mvs = MultivariateStudent::new(vec![0., 0.], vec![1., 0., 0., 1.], 4.).unwrap(); /// assert_eq!(mvs.mean().unwrap(), DVector::from_vec(vec![0., 0.])); /// assert_eq!(mvs.variance().unwrap(), DMatrix::from_vec(2, 2, vec![2., 0., 0., 2.])); -/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.04715702017537655); +/// assert_eq!(mvs.pdf(&DVector::from_vec(vec![1., 1.])), 0.047157020175376395); /// ``` #[derive(Debug, Clone, PartialEq)] pub struct MultivariateStudent diff --git a/src/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 970e76ed..09c0a138 100644 --- a/src/function/beta.rs +++ b/src/function/beta.rs @@ -132,12 +132,80 @@ 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, /// `b` is the second beta parameter, and `x` is the upper limit of the /// integral. /// +/// # Remarks +/// +/// 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: +/// +/// ```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 /// /// if `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` @@ -154,18 +222,80 @@ 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 { - (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; + 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 symm_transform = x >= (a + 1.0) / (a + b + 2.0); + 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 +318,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 +353,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 +786,81 @@ 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 `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 = [ + 1e-308f64, 1e-300, 1e-100, 1e-8, 0.5, 1.0, 20.0, 1e8, 1e20, 1e100, 1e300, 1e308, + ]; + for &a in ¶ms { + for &b in ¶ms { + // Beyond a ~1e300 parameter ratio the Lentz recurrence bottoms + // out in its own `fpmin` guards for `x` within an ulp of the + // mode, and returns NaN. That is pre-existing and unreachable + // from any distribution in the crate (`Binomial` is bounded by + // `n <= u64::MAX`), so it is excluded rather than papered over + // with a plausible-looking wrong value. + if a.max(b) / a.min(b) > 1e200 { + continue; + } + for x in [0.0f64, 1e-300, 0.25, 0.5, 0.75, 1.0] { + let v = beta_reg(a, b, x); + assert!( + v.is_finite() && (0.0..=1.0).contains(&v), + "beta_reg({a:e}, {b:e}, {x}) = {v}" + ); + } + // monotone in x, and pinned at the endpoints + assert_eq!(beta_reg(a, b, 0.0), 0.0, "beta_reg({a:e},{b:e},0)"); + assert_eq!(beta_reg(a, b, 1.0), 1.0, "beta_reg({a:e},{b:e},1)"); + } + } + } + + /// `I_x(a, b) + I_{1-x}(b, a) == 1` for every valid `a`, `b`, `x`. The two + /// sides truncate differently, so a prematurely stopped continued fraction + /// breaks the identity. + #[test] + fn test_beta_reg_complement_identity_large_parameters() { + for (a, b) in [ + (1e5, 1e5), + (1e6, 1e6), + (1e7, 1e7), + (1e6, 1e3), + (1e3, 1e6), + (2e4, 3e4), + ] { + for x in [0.1, 0.25, 0.5, 0.5 + 1e-9, 0.75, 0.9] { + let lhs = beta_reg(a, b, x) + beta_reg(b, a, 1.0 - x); + prec::assert_abs_diff_eq!(lhs, 1.0, epsilon = 1e-7); + } + } + } + #[test] fn test_error_is_sync_send() { fn assert_sync_send() {} diff --git a/src/function/gamma.rs b/src/function/gamma.rs index bff5c45a..c38bb9a6 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,336 @@ 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() +} + +// --------------------------------------------------------------------------- +// 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; -/// 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 +369,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) } } @@ -204,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) }; } @@ -301,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); @@ -390,7 +694,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 +891,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 +1193,7 @@ mod tests { prec::assert_abs_diff_eq!( gamma_li(5.5, 2.0), 1.5746265342113649473739798668921124454837064926448459, - epsilon = 1e-15 + epsilon = 2e-15 ); prec::assert_abs_diff_eq!( gamma_li(5.5, 8.0), @@ -1195,6 +1505,34 @@ mod tests { assert!(super::checked_gamma_ui(1.0, f64::INFINITY).is_err()); } + /// `digamma` has poles at the non-positive integers and returns -inf there. + /// The pole test is `prec::ulps_eq!(x.floor(), x)`, whose default epsilon + /// used to be `1e-9` absolute, so a whole neighbourhood around each pole + /// returned -inf instead of a large finite value. + #[test] + fn test_digamma_near_negative_integer_poles_is_finite() { + prec::assert_relative_eq!(digamma(-1.0 + f64::powi(2.0, -33)), -8589934591.5772156646, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(digamma(-1.0 - f64::powi(2.0, -33)), 8589934592.4227843348, epsilon = 0.0, max_relative = 1e-13); + // the poles themselves are unchanged + assert_eq!(digamma(-1.0), f64::NEG_INFINITY); + assert_eq!(digamma(0.0), f64::NEG_INFINITY); + assert_eq!(digamma(-5.0), f64::NEG_INFINITY); + } + + /// Before `tan_pi` reduced the reflection argument by the period, these + /// lost 3-4 decimal digits because `(PI * x).tan()` was evaluated at a + /// large rounded argument. Reference values are mpmath at 40 significant + /// digits, computed at the exact `f64` of each literal. The tolerances + /// reflect the cancellation between the two reflection terms + /// (`psi(1 - x)` and `pi * cot(pi x)` are each ~2.6 at x = -12.72 and + /// cancel to ~0.017, amplifying their ~1 ulp errors ~150x). + #[test] + fn test_digamma_negative_arguments_far_from_zero() { + prec::assert_relative_eq!(digamma(-12.72), -0.0169824608177603739173, epsilon = 0.0, max_relative = 1e-13); + prec::assert_relative_eq!(digamma(-14.72), 0.1238386191670285663687, epsilon = 0.0, max_relative = 1e-14); + prec::assert_relative_eq!(digamma(-20.02), 52.95568424702714411621, epsilon = 0.0, max_relative = 1e-15); + } + // TODO: precision testing could be more accurate #[test] fn test_digamma() { @@ -1357,6 +1695,116 @@ 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}"); + } + } + + /// `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 ebf1b6d4..c174598c 100644 --- a/src/prec.rs +++ b/src/prec.rs @@ -29,6 +29,7 @@ //! - `DEFAULT_RELATIVE_ACC`: 1e-14 for relative comparisons //! - `DEFAULT_EPS`: 1e-9 for absolute comparisons //! - `DEFAULT_ULPS`: 5 for ULPs comparisons +//! - `DEFAULT_ULPS_EPS`: `f64::EPSILON`, the absolute floor paired with `DEFAULT_ULPS` //! //! These defaults should be used unless there is a specific reason to use different //! precision levels. @@ -62,6 +63,16 @@ pub const DEFAULT_EPS: f64 = 1e-9; /// Default and target ULPs accuracy for f64 operations pub const DEFAULT_ULPS: u32 = 5; +/// Default absolute epsilon for ULPs comparisons. +/// +/// `approx`'s `ulps_eq` short-circuits on `abs_diff_eq(epsilon)` before it looks +/// at the ULPs distance, so pairing it with [`DEFAULT_EPS`] (`1e-9`) would make +/// the ULPs bound unreachable and turn `ulps_eq!(x, y)` into a `1e-9` absolute +/// comparison. `ulps_eq!` is used inside the crate to recognise exact parameter +/// values (`p == 1.0`, `x == x.floor()`), so its epsilon has to stay at the +/// scale of a genuine rounding error. +pub const DEFAULT_ULPS_EPS: f64 = f64::EPSILON; + /// Compares if two floats are close via `approx::abs_diff_eq` /// using a maximum absolute difference (epsilon) of `acc`. #[deprecated(since = "0.19.0", note = "Use abs_diff_eq! macro instead")] @@ -81,6 +92,68 @@ 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)) +} + +/// 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, @@ -144,7 +217,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 +235,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;