diff --git a/Cargo.toml b/Cargo.toml index 120e3bfa..3bed0943 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,14 +3,22 @@ name = "statrs" version = "0.19.1" authors = ["Michael Ma"] description = "Statistical computing library for Rust" -license = "MIT" +license = "MIT AND BSL-1.0" keywords = ["probability", "statistics", "stats", "distribution", "math"] categories = ["science"] homepage = "https://github.com/statrs-dev/statrs" repository = "https://github.com/statrs-dev/statrs" edition = "2024" -include = ["CHANGELOG.md", "LICENSE.md", "src/", "tests/"] +include = [ + "CHANGELOG.md", + "LICENSE.md", + "LICENSE-BOOST.md", + "THIRD_PARTY_NOTICES.md", + "src/", + "tests/", + "benches/", +] # When changing MSRV: Also update the README rust-version = "1.89.0" @@ -29,6 +37,11 @@ name = "density" harness = false required-features = ["rand", "std", "kde"] +[[bench]] +name = "beta" +harness = false +required-features = ["std"] + [features] default = ["std", "nalgebra", "rand"] std = ["approx/std", "num-traits/std", "nalgebra?/std", "rand?/std"] diff --git a/LICENSE-BOOST.md b/LICENSE-BOOST.md new file mode 100644 index 00000000..36b7cd93 --- /dev/null +++ b/LICENSE-BOOST.md @@ -0,0 +1,23 @@ +Boost Software License - Version 1.0 - August 17th, 2003 + +Permission is hereby granted, free of charge, to any person or organization +obtaining a copy of the software and accompanying documentation covered by +this license (the "Software") to use, reproduce, display, distribute, +execute, and transmit the Software, and to prepare derivative works of the +Software, and to permit third-parties to whom the Software is furnished to +do so, all subject to the following: + +The copyright notices in the Software and this entire statement, including +the above license grant, this restriction and the following disclaimer, +must be included in all copies of the Software, in whole or in part, and +all derivative works of the Software, unless such copies or derivative +works are solely in the form of machine-executable object code generated by +a source language processor. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/README.md b/README.md index 540f2642..701a4cbe 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # statrs ![tests][actions-test-badge] -[![MIT licensed][license-badge]](./LICENSE.md) +[![MIT and BSL-1.0 licensed][license-badge]](./THIRD_PARTY_NOTICES.md) [![Crate][crates-badge]][crates-url] [![docs.rs][docsrs-badge]][docs-url] [![codecov-statrs][codecov-badge]][codecov-url] @@ -10,7 +10,7 @@ [actions-test-badge]: https://github.com/statrs-dev/statrs/actions/workflows/test.yml/badge.svg [crates-badge]: https://img.shields.io/crates/v/statrs.svg [crates-url]: https://crates.io/crates/statrs -[license-badge]: https://img.shields.io/badge/license-MIT-blue.svg +[license-badge]: https://img.shields.io/badge/license-MIT%20AND%20BSL--1.0-blue.svg [docsrs-badge]: https://img.shields.io/docsrs/statrs [docs-url]: https://docs.rs/statrs/*/statrs [codecov-badge]: https://codecov.io/gh/statrs-dev/statrs/graph/badge.svg?token=XtMSMYXvIf diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md new file mode 100644 index 00000000..9441840b --- /dev/null +++ b/THIRD_PARTY_NOTICES.md @@ -0,0 +1,24 @@ +# Third-party notices + +## Boost.Math + +Portions of `src/function/beta/{bgrat,forward,fraction,log_forward,recurrence,series}.rs` are adapted from Boost.Math 1.90.0, `include/boost/math/special_functions/beta.hpp`. + +Portions of `src/function/beta/normal_tail.rs` are adapted from Boost.Math 1.90.0, `include/boost/math/special_functions/erf.hpp`. + +Copyright John Maddock 2006. +Copyright Matt Borland 2024. + +The Boost-derived portions are licensed under the Boost Software License 1.0; see `LICENSE-BOOST.md`. Statrs modifications are licensed under MIT, so these files are subject to both licenses. + +Source: https://github.com/boostorg/math/blob/e0fcd19f7227d81391770ea46015acc3c80af810/include/boost/math/special_functions/beta.hpp + +Source: https://github.com/boostorg/math/blob/e0fcd19f7227d81391770ea46015acc3c80af810/include/boost/math/special_functions/erf.hpp + +## special 0.8.1 + +The initial inverse-beta estimate in `src/function/beta/inverse/initial.rs` is adapted from `special` 0.8.1 under its MIT license option. + +Copyright 2014–2019 The special Developers. + +Source: https://github.com/stainless-steel/special/blob/c64902430bd50e8c8225c7c8b410334ffedf2f15/src/beta.rs diff --git a/benches/beta.rs b/benches/beta.rs new file mode 100644 index 00000000..928618b2 --- /dev/null +++ b/benches/beta.rs @@ -0,0 +1,70 @@ +use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main}; +use statrs::function::beta::{beta_reg, inv_beta_reg}; +use std::hint::black_box; + +fn bench_beta_reg(c: &mut Criterion) { + let mut group = c.benchmark_group("beta_reg"); + for (name, a, b, x) in [ + ("typical", 2.0, 5.0, 0.3), + ( + "large_symmetric_adjacent", + 1e8, + 1e8, + f64::from_bits(0.5_f64.to_bits() + 1), + ), + ( + "moderate_fraction", + 25.32628846940565, + 3.1028101710805442, + 0.9276950604606229, + ), + ] { + group.bench_with_input( + BenchmarkId::new("cdf", name), + &(a, b, x), + |bencher, input| { + bencher + .iter(|| beta_reg(black_box(input.0), black_box(input.1), black_box(input.2))); + }, + ); + } + group.finish(); +} + +fn bench_inv_beta_reg(c: &mut Criterion) { + let mut group = c.benchmark_group("inv_beta_reg"); + for (name, a, b, probability) in [ + ("typical", 2.0, 5.0, 0.3), + ("nontermination_regression", 200.0, 2.0, 1e-60), + ("panic_regression", 200.0, 2.0, 1e-165), + ("tiny_quantile", 0.1, 500.0, 1e-30), + ( + "subnormal_shape_two", + 0.5, + 2.0, + f64::from_bits(0x1e72_f942_2c23_c47c), + ), + ("subnormal_unit_shape", 1.0, 10.0, f64::from_bits(5)), + ("shape_two_large_median", 2.0, 1e308, 0.5), + ( + "shape_two_zero_cell", + 2.0, + 1e200, + f64::from_bits(0x0c8b_4ec7_f919_73ff), + ), + ] { + group.bench_with_input( + BenchmarkId::new("quantile", name), + &(a, b, probability), + |bencher, input| { + bencher.iter(|| { + inv_beta_reg(black_box(input.0), black_box(input.1), black_box(input.2)) + }); + }, + ); + } + group.finish(); +} + +criterion_group!(benches, bench_beta_reg, bench_inv_beta_reg); +criterion_main!(benches); diff --git a/src/distribution/beta.rs b/src/distribution/beta.rs index 38ae1517..1fc3ce9c 100644 --- a/src/distribution/beta.rs +++ b/src/distribution/beta.rs @@ -129,6 +129,10 @@ impl ContinuousCDF for Beta { /// Calculates the cumulative distribution function for the beta /// distribution at `x`. /// + /// # Panics + /// + /// If the numerical method does not converge. + /// /// # Formula /// /// ```text @@ -151,6 +155,10 @@ impl ContinuousCDF for Beta { /// Calculates the survival function for the beta distribution at `x`. /// + /// # Panics + /// + /// If the numerical method does not converge. + /// /// # Formula /// /// ```text @@ -167,12 +175,9 @@ impl ContinuousCDF for Beta { } else if self.shape_a == 1.0 && self.shape_b == 1.0 { 1. - x } else if x < (self.shape_a + 1.0) / (self.shape_a + self.shape_b + 2.0) { - // Below the continued fraction split point of `beta_reg`, - // `beta_reg(b, a, 1 - x)` reduces to `1 - beta_reg(a, b, x)`; - // computing the complement here instead avoids `1.0 - x` - // rounding to 1.0 for tiny x (< ~1.1e-16), which would lose - // the lower tail entirely. See #432 - 1.0 - beta::beta_reg(self.shape_a, self.shape_b, x) + beta::checked_ln_beta_reg_complement(self.shape_a, self.shape_b, x) + .unwrap() + .exp() } else { beta::beta_reg(self.shape_b, self.shape_a, 1.0 - x) } @@ -183,7 +188,7 @@ impl ContinuousCDF for Beta { /// /// # Panics /// - /// If x is not in `[0, 1]`. + /// If x is not in `[0, 1]` or the numerical method does not converge. /// /// # Formula /// @@ -208,6 +213,10 @@ impl ContinuousCDF for Beta { /// /// If x is not in `[0, 1]`. /// + /// # Panics + /// + /// If the numerical method does not converge. + /// /// # Formula /// /// ```text @@ -651,6 +660,18 @@ mod tests { } } + #[test] + fn test_cdf_large_symmetric_shapes() { + for shape in [1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8] { + let dist = Beta::new(shape, shape).unwrap(); + let cdf = dist.cdf(0.5); + let sf = dist.sf(0.5); + assert_eq!(cdf, 0.5); + assert_eq!(sf, 0.5); + assert_eq!(cdf + sf, 1.0); + } + } + #[test] fn test_sf() { let sf = |arg: f64| move |x: Beta| x.sf(arg); @@ -684,6 +705,17 @@ mod tests { } } + #[test] + fn test_sf_tiny_shape_preserves_representable_tail() { + let distribution = Beta::new( + f64::from_bits(0x00000000000007e8), + f64::from_bits(0x4040000000000000), + ) + .unwrap(); + let x = f64::from_bits(0x01556e1fc2f8f359); + assert_eq!(distribution.sf(x).to_bits(), 0x0000000000155101); + } + #[test] fn test_inverse_cdf() { // let inverse_cdf = |arg: f64| move |x: Beta| x.inverse_cdf(arg); @@ -705,6 +737,18 @@ mod tests { } } + #[test] + fn test_inverse_cdf_extreme_lower_tail() { + let dist = Beta::new(200.0, 2.0).unwrap(); + let actual = dist.inverse_cdf(1e-170); + let expected = 0.13765877485659653; + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-13, + "actual {actual}, expected {expected}" + ); + } + #[test] fn test_cdf_input_lt_0() { let cdf = |arg: f64| move |x: Beta| x.cdf(arg); diff --git a/src/distribution/binomial/mod.rs b/src/distribution/binomial/mod.rs index acf3d1f4..038dc17c 100644 --- a/src/distribution/binomial/mod.rs +++ b/src/distribution/binomial/mod.rs @@ -514,7 +514,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.9717524751, 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/function/beta.rs b/src/function/beta.rs deleted file mode 100644 index 1f879819..00000000 --- a/src/function/beta.rs +++ /dev/null @@ -1,654 +0,0 @@ -//! Provides the [beta](https://en.wikipedia.org/wiki/Beta_function) and related -//! function -//! -//! This module sets the default precision more tightly than crate defaults for `DEFAULT_EPS` - -use crate::function::gamma; -use crate::prec; -#[cfg(not(feature = "std"))] -use num_traits::Float as _; - -/// sample case of module level precision -const MODULE_EPS: f64 = 1e-15; - -/// Represents the errors that can occur when computing the natural logarithm -/// of the beta function or the regularized lower incomplete beta function. -#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)] -#[non_exhaustive] -pub enum BetaFuncError { - /// `a` is zero or less than zero. - ANotGreaterThanZero, - - /// `b` is zero or less than zero. - BNotGreaterThanZero, - - /// `x` is not in `[0, 1]`. - XOutOfRange, -} - -impl core::fmt::Display for BetaFuncError { - fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result { - match self { - BetaFuncError::ANotGreaterThanZero => write!(f, "a is zero or less than zero"), - BetaFuncError::BNotGreaterThanZero => write!(f, "b is zero or less than zero"), - BetaFuncError::XOutOfRange => write!(f, "x is not in [0, 1]"), - } - } -} - -impl core::error::Error for BetaFuncError {} - -/// Computes the natural logarithm -/// of the beta function -/// where `a` is the first beta parameter -/// and `b` is the second beta parameter -/// and `a > 0`, `b > 0`. -/// -/// # Panics -/// -/// if `a <= 0.0` or `b <= 0.0` -pub fn ln_beta(a: f64, b: f64) -> f64 { - checked_ln_beta(a, b).unwrap() -} - -/// Computes the natural logarithm -/// of the beta function -/// where `a` is the first beta parameter -/// and `b` is the second beta parameter -/// and `a > 0`, `b > 0`. -/// -/// # Errors -/// -/// if `a <= 0.0` or `b <= 0.0` -pub fn checked_ln_beta(a: f64, b: f64) -> Result { - if a <= 0.0 { - Err(BetaFuncError::ANotGreaterThanZero) - } else if b <= 0.0 { - Err(BetaFuncError::BNotGreaterThanZero) - } else { - Ok(gamma::ln_gamma(a) + gamma::ln_gamma(b) - gamma::ln_gamma(a + b)) - } -} - -/// Computes the beta function -/// where `a` is the first beta parameter -/// and `b` is the second beta parameter. -/// -/// -/// # Panics -/// -/// if `a <= 0.0` or `b <= 0.0` -pub fn beta(a: f64, b: f64) -> f64 { - checked_beta(a, b).unwrap() -} - -/// Computes the beta function -/// where `a` is the first beta parameter -/// and `b` is the second beta parameter. -/// -/// -/// # Errors -/// -/// if `a <= 0.0` or `b <= 0.0` -pub fn checked_beta(a: f64, b: f64) -> Result { - checked_ln_beta(a, b).map(|x| x.exp()) -} - -/// Computes the lower incomplete (unregularized) beta function -/// `B(a,b,x) = int(t^(a-1)*(1-t)^(b-1),t=0..x)` for `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 -/// -/// # Panics -/// -/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` -pub fn beta_inc(a: f64, b: f64, x: f64) -> f64 { - checked_beta_inc(a, b, x).unwrap() -} - -/// Computes the lower incomplete (unregularized) beta function -/// `B(a,b,x) = int(t^(a-1)*(1-t)^(b-1),t=0..x)` for `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 -/// -/// # Errors -/// -/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` -pub fn checked_beta_inc(a: f64, b: f64, x: f64) -> Result { - checked_beta_reg(a, b, x).and_then(|x| checked_beta(a, b).map(|y| x * y)) -} - -/// 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. -/// -/// # Panics -/// -/// if `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` -pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 { - checked_beta_reg(a, b, x).unwrap() -} - -/// 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. -/// -/// # Errors -/// -/// if `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` -pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { - if a <= 0.0 { - return Err(BetaFuncError::ANotGreaterThanZero); - } - - if b <= 0.0 { - return Err(BetaFuncError::BNotGreaterThanZero); - } - - if !(0.0..=1.0).contains(&x) { - return Err(BetaFuncError::XOutOfRange); - } - - 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() - }; - let symm_transform = x >= (a + 1.0) / (a + b + 2.0); - let eps = prec::F64_PREC; - let fpmin = f64::MIN_POSITIVE / eps; - - let mut a = a; - let mut b = b; - let mut x = x; - if symm_transform { - let swap = a; - x = 1.0 - x; - a = b; - b = swap; - } - - let qab = a + b; - let qap = a + 1.0; - let qam = a - 1.0; - let mut c = 1.0; - let mut d = 1.0 - qab * x / qap; - - if d.abs() < fpmin { - d = fpmin; - } - d = 1.0 / d; - let mut h = d; - - for m in 1..141 { - let m = f64::from(m); - let m2 = m * 2.0; - let mut aa = m * (b - m) * x / ((qam + m2) * (a + m2)); - d = 1.0 + aa * d; - - if d.abs() < fpmin { - d = fpmin; - } - - c = 1.0 + aa / c; - if c.abs() < fpmin { - c = fpmin; - } - - d = 1.0 / d; - h = h * d * c; - aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); - d = 1.0 + aa * d; - - if d.abs() < fpmin { - d = fpmin; - } - - c = 1.0 + aa / c; - - if c.abs() < fpmin { - c = fpmin; - } - - d = 1.0 / d; - let del = d * c; - h *= del; - - if (del - 1.0).abs() <= eps { - return if symm_transform { - Ok(1.0 - bt * h / a) - } else { - Ok(bt * h / a) - }; - } - } - - if symm_transform { - Ok(1.0 - bt * h / a) - } else { - Ok(bt * h / a) - } -} - -/// Computes the inverse of the regularized incomplete beta function -// This code is based on the implementation in the ["special"][1] crate, -// which in turn is based on a [C implementation][2] by John Burkardt. The -// original algorithm was published in Applied Statistics and is known as -// [Algorithm AS 64][3] and [Algorithm AS 109][4]. -// -// [1]: https://docs.rs/special/0.8.1/ -// [2]: http://people.sc.fsu.edu/~jburkardt/c_src/asa109/asa109.html -// [3]: http://www.jstor.org/stable/2346798 -// [4]: http://www.jstor.org/stable/2346887 -// -// > Copyright 2014–2019 The special Developers -// > -// > Permission is hereby granted, free of charge, to any person obtaining a copy of -// > this software and associated documentation files (the "Software"), to deal in -// > the Software without restriction, including without limitation the rights to -// > use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of -// > the Software, and to permit persons to whom the Software is furnished to do so, -// > subject to the following conditions: -// > -// > The above copyright notice and this permission notice shall be included in all -// > copies or substantial portions of the Software. -// > -// > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS -// > FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR -// > COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER -// > IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN -// > CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -pub fn inv_beta_reg(mut a: f64, mut b: f64, mut x: f64) -> f64 { - // Algorithm AS 64 - // http://www.jstor.org/stable/2346798 - // - // An approximation x₀ to x if found from (cf. Scheffé and Tukey, 1944) - // - // 1 + x₀ 4p + 2q - 2 - // ------ = ----------- - // 1 - x₀ χ²(α) - // - // where χ²(α) is the upper α point of the χ² distribution with 2q - // degrees of freedom and is obtained from Wilson and Hilferty's - // approximation (cf. Wilson and Hilferty, 1931) - // - // χ²(α) = 2q (1 - 1 / (9q) + y(α) sqrt(1 / (9q)))^3, - // - // y(α) being Hastings' approximation (cf. Hastings, 1955) for the upper - // α point of the standard normal distribution. If χ²(α) < 0, then - // - // x₀ = 1 - ((1 - α)q B(p, q))^(1 / q). - // - // Again if (4p + 2q - 2) / χ²(α) does not exceed 1, x₀ is obtained from - // - // x₀ = (αp B(p, q))^(1 / p). - // - // The final solution is obtained by the Newton–Raphson method from the - // relation - // - // f(x[i - 1]) - // x[i] = x[i - 1] - ------------ - // f'(x[i - 1]) - // - // where - // - // f(x) = I(x, p, q) - α. - let ln_beta = ln_beta(a, b); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - const SAE: i32 = -30; - const FPU: f64 = 1e-30; // 10^SAE - - debug_assert!((0.0..=1.0).contains(&x) && a > 0.0 && b > 0.0); - - if x == 0.0 { - return 0.0; - } - if x == 1.0 { - return 1.0; - } - - let mut p; - let mut q; - - let flip = 0.5 < x; - if flip { - p = a; - a = b; - b = p; - x = 1.0 - x; - } - - p = (-(x * x).ln()).sqrt(); - q = p - (2.30753 + 0.27061 * p) / (1.0 + (0.99229 + 0.04481 * p) * p); - - if 1.0 < a && 1.0 < b { - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - // - // For a and b > 1, the approximation given by Carter (1947), which - // improves the Fisher–Cochran formula, is generally better. For - // other values of a and b en empirical investigation has shown that - // the approximation given in AS 64 is adequate. - let r = (q * q - 3.0) / 6.0; - let s = 1.0 / (2.0 * a - 1.0); - let t = 1.0 / (2.0 * b - 1.0); - let h = 2.0 / (s + t); - let w = q * (h + r).sqrt() / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h)); - p = a / (a + b * (2.0 * w).exp()); - } else { - let mut t = 1.0 / (9.0 * b); - t = 2.0 * b * (1.0 - t + q * t.sqrt()).powf(3.0); - if t <= 0.0 { - p = 1.0 - ((((1.0 - x) * b).ln() + ln_beta) / b).exp(); - } else { - t = 2.0 * (2.0 * a + b - 1.0) / t; - if t <= 1.0 { - p = (((x * a).ln() + ln_beta) / a).exp(); - } else { - p = 1.0 - 2.0 / (t + 1.0); - } - } - } - - p = p.clamp(0.0001, 0.9999); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - let e = (-5.0 / a / a - 1.0 / x.powf(0.2) - 13.0) as i32; - let acu = if e > SAE { f64::powi(10.0, e) } else { FPU }; - - let mut pnext; - let mut qprev = 0.0; - let mut sq = 1.0; - let mut prev = 1.0; - - 'outer: loop { - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - q = beta_reg(a, b, p); - q = (q - x) * (ln_beta + (1.0 - a) * p.ln() + (1.0 - b) * (1.0 - p).ln()).exp(); - - // Remark AS R83 - // http://www.jstor.org/stable/2347779 - if q * qprev <= 0.0 { - prev = if sq > FPU { sq } else { FPU }; - } - - // Remark AS R19 and Algorithm AS 109 - // http://www.jstor.org/stable/2346887 - let mut g = 1.0; - loop { - loop { - let adj = g * q; - sq = adj * adj; - - if sq < prev { - pnext = p - adj; - if (0.0..=1.0).contains(&pnext) { - break; - } - } - g /= 3.0; - } - - if prev <= acu || q * q <= acu { - p = pnext; - break 'outer; - } - - if pnext != 0.0 && pnext != 1.0 { - break; - } - - g /= 3.0; - } - - if pnext == p { - break; - } - - p = pnext; - qprev = q; - } - - if flip { 1.0 - p } else { p } -} - -#[cfg(test)] -mod tests { - use super::*; - use crate::prec; - use core::f64::consts as f64_consts; - const MODULE_RELATIVE_ACC: f64 = 1e-14; - - fn beta_assert_relative_eq(a: f64, b: f64) { - prec::assert_relative_eq!( - a, - b, - epsilon = MODULE_EPS, - max_relative = MODULE_RELATIVE_ACC - ); - } - - fn beta_assert_abs_diff_eq(a: f64, b: f64) { - prec::assert_abs_diff_eq!(a, b, epsilon = MODULE_EPS); - } - - #[test] - fn test_ln_beta() { - beta_assert_relative_eq(ln_beta(0.5, 0.5), 1.144729885849400174144); - beta_assert_relative_eq(ln_beta(1.0, 0.5), f64_consts::LN_2); - beta_assert_relative_eq(ln_beta(2.5, 0.5), 0.163900632837673937284); - beta_assert_relative_eq(ln_beta(0.5, 1.0), f64_consts::LN_2); - beta_assert_relative_eq(ln_beta(1.0, 1.0), 0.0); - beta_assert_relative_eq(ln_beta(2.5, 1.0), -0.9162907318741550651835); - beta_assert_relative_eq(ln_beta(0.5, 2.5), 0.163900632837673937284); - beta_assert_relative_eq(ln_beta(1.0, 2.5), -0.9162907318741550651835); - beta_assert_relative_eq(ln_beta(2.5, 2.5), -2.608688089402107300388); - } - - #[test] - #[should_panic] - fn test_ln_beta_a_lte_0() { - ln_beta(0.0, 0.5); - } - - #[test] - #[should_panic] - fn test_ln_beta_b_lte_0() { - ln_beta(0.5, 0.0); - } - - #[test] - fn test_checked_ln_beta_a_lte_0() { - assert!(checked_ln_beta(0.0, 0.5).is_err()); - } - - #[test] - fn test_checked_ln_beta_b_lte_0() { - assert!(checked_ln_beta(0.5, 0.0).is_err()); - } - - #[test] - #[should_panic] - fn test_beta_a_lte_0() { - beta(0.0, 0.5); - } - - #[test] - #[should_panic] - fn test_beta_b_lte_0() { - beta(0.5, 0.0); - } - - #[test] - fn test_checked_beta_a_lte_0() { - assert!(checked_beta(0.0, 0.5).is_err()); - } - - #[test] - fn test_checked_beta_b_lte_0() { - assert!(checked_beta(0.5, 0.0).is_err()); - } - - #[test] - fn test_beta() { - beta_assert_relative_eq(beta(0.5, 0.5), f64_consts::PI); - beta_assert_relative_eq(beta(1.0, 0.5), 2.0); - beta_assert_relative_eq(beta(2.5, 0.5), 1.17809724509617246442); - beta_assert_relative_eq(beta(0.5, 1.0), 2.0); - beta_assert_relative_eq(beta(1.0, 1.0), 1.0); - beta_assert_relative_eq(beta(2.5, 1.0), 0.4); - beta_assert_relative_eq(beta(0.5, 2.5), 1.17809724509617246442); - beta_assert_relative_eq(beta(1.0, 2.5), 0.4); - beta_assert_relative_eq(beta(2.5, 2.5), 0.073631077818510779026); - } - - #[test] - fn test_beta_inc() { - beta_assert_relative_eq(beta_inc(0.5, 0.5, 0.5), f64_consts::FRAC_PI_2); - beta_assert_relative_eq(beta_inc(0.5, 0.5, 1.0), f64_consts::PI); - beta_assert_relative_eq(beta_inc(1.0, 0.5, 0.5), 0.5857864376269049511983); - beta_assert_relative_eq(beta_inc(1.0, 0.5, 1.0), 2.0); - beta_assert_relative_eq(beta_inc(2.5, 0.5, 0.5), 0.0890486225480862322117); - beta_assert_relative_eq(beta_inc(2.5, 0.5, 1.0), 1.17809724509617246442); - beta_assert_relative_eq(beta_inc(0.5, 1.0, 0.5), f64_consts::SQRT_2); - beta_assert_relative_eq(beta_inc(0.5, 1.0, 1.0), 2.0); - beta_assert_relative_eq(beta_inc(1.0, 1.0, 0.5), 0.5); - beta_assert_relative_eq(beta_inc(1.0, 1.0, 1.0), 1.0); - beta_assert_relative_eq(beta_inc(2.5, 1.0, 0.5), 0.0707106781186547524401); - beta_assert_relative_eq(beta_inc(2.5, 1.0, 1.0), 0.4); - beta_assert_relative_eq(beta_inc(0.5, 2.5, 0.5), 1.08904862254808623221); - beta_assert_relative_eq(beta_inc(0.5, 2.5, 1.0), 1.17809724509617246442); - beta_assert_relative_eq(beta_inc(1.0, 2.5, 0.5), 0.32928932188134524756); - beta_assert_relative_eq(beta_inc(1.0, 2.5, 1.0), 0.4); - beta_assert_relative_eq(beta_inc(2.5, 2.5, 0.5), 0.03681553890925538951323); - beta_assert_relative_eq(beta_inc(2.5, 2.5, 1.0), 0.073631077818510779026); - } - - #[test] - #[should_panic] - fn test_beta_inc_a_lte_0() { - beta_inc(0.0, 1.0, 1.0); - } - - #[test] - #[should_panic] - fn test_beta_inc_b_lte_0() { - beta_inc(1.0, 0.0, 1.0); - } - - #[test] - #[should_panic] - fn test_beta_inc_x_lt_0() { - beta_inc(1.0, 1.0, -1.0); - } - - #[test] - #[should_panic] - fn test_beta_inc_x_gt_1() { - beta_inc(1.0, 1.0, 2.0); - } - - #[test] - fn test_checked_beta_inc_a_lte_0() { - assert!(checked_beta_inc(0.0, 1.0, 1.0).is_err()); - } - - #[test] - fn test_checked_beta_inc_b_lte_0() { - assert!(checked_beta_inc(1.0, 0.0, 1.0).is_err()); - } - - #[test] - fn test_checked_beta_inc_x_lt_0() { - assert!(checked_beta_inc(1.0, 1.0, -1.0).is_err()); - } - - #[test] - fn test_checked_beta_inc_x_gt_1() { - assert!(checked_beta_inc(1.0, 1.0, 2.0).is_err()); - } - - #[test] - fn test_beta_reg() { - beta_assert_abs_diff_eq(beta_reg(0.5, 0.5, 0.5), 0.5); - assert_eq!(beta_reg(0.5, 0.5, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(1.0, 0.5, 0.5), 0.292893218813452475599); - assert_eq!(beta_reg(1.0, 0.5, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(2.5, 0.5, 0.5), 0.07558681842161243795); - assert_eq!(beta_reg(2.5, 0.5, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(0.5, 1.0, 0.5), f64_consts::FRAC_1_SQRT_2); - assert_eq!(beta_reg(0.5, 1.0, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(1.0, 1.0, 0.5), 0.5); - assert_eq!(beta_reg(1.0, 1.0, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(2.5, 1.0, 0.5), 0.1767766952966368811); - assert_eq!(beta_reg(2.5, 1.0, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(0.5, 2.5, 0.5), 0.92441318157838756205); - assert_eq!(beta_reg(0.5, 2.5, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(1.0, 2.5, 0.5), 0.8232233047033631189); - assert_eq!(beta_reg(1.0, 2.5, 1.0), 1.0); - beta_assert_abs_diff_eq(beta_reg(2.5, 2.5, 0.5), 0.5); - assert_eq!(beta_reg(2.5, 2.5, 1.0), 1.0); - } - - #[test] - #[should_panic] - fn test_beta_reg_a_lte_0() { - beta_reg(0.0, 1.0, 1.0); - } - - #[test] - #[should_panic] - fn test_beta_reg_b_lte_0() { - beta_reg(1.0, 0.0, 1.0); - } - - #[test] - #[should_panic] - fn test_beta_reg_x_lt_0() { - beta_reg(1.0, 1.0, -1.0); - } - - #[test] - #[should_panic] - fn test_beta_reg_x_gt_1() { - beta_reg(1.0, 1.0, 2.0); - } - - #[test] - fn test_checked_beta_reg_a_lte_0() { - assert!(checked_beta_reg(0.0, 1.0, 1.0).is_err()); - } - - #[test] - fn test_checked_beta_reg_b_lte_0() { - assert!(checked_beta_reg(1.0, 0.0, 1.0).is_err()); - } - - #[test] - fn test_checked_beta_reg_x_lt_0() { - assert!(checked_beta_reg(1.0, 1.0, -1.0).is_err()); - } - - #[test] - fn test_checked_beta_reg_x_gt_1() { - assert!(checked_beta_reg(1.0, 1.0, 2.0).is_err()); - } - - #[test] - fn test_error_is_sync_send() { - fn assert_sync_send() {} - assert_sync_send::(); - } -} diff --git a/src/function/beta/api.rs b/src/function/beta/api.rs new file mode 100644 index 00000000..4bf6c433 --- /dev/null +++ b/src/function/beta/api.rs @@ -0,0 +1,63 @@ +use super::*; + +/// Computes the beta function +/// where `a` is the first beta parameter +/// and `b` is the second beta parameter. +/// +/// +/// # Panics +/// +/// if `a <= 0.0` or `b <= 0.0` +pub fn beta(a: f64, b: f64) -> f64 { + checked_beta(a, b).unwrap() +} + +/// Computes the beta function +/// where `a` is the first beta parameter +/// and `b` is the second beta parameter. +/// +/// +/// # Errors +/// +/// if `a <= 0.0` or `b <= 0.0` +pub fn checked_beta(a: f64, b: f64) -> Result { + checked_ln_beta(a, b).map(|x| x.exp()) +} + +/// Computes the lower incomplete (unregularized) beta function +/// `B(a,b,x) = int(t^(a-1)*(1-t)^(b-1),t=0..x)` for `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 +/// +/// # Panics +/// +/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` +pub fn beta_inc(a: f64, b: f64, x: f64) -> f64 { + checked_beta_inc(a, b, x).unwrap() +} + +/// Computes the lower incomplete (unregularized) beta function +/// `B(a,b,x) = int(t^(a-1)*(1-t)^(b-1),t=0..x)` for `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 +/// +/// # Errors +/// +/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, or `x > 1.0` +pub fn checked_beta_inc(a: f64, b: f64, x: f64) -> Result { + checked_beta_reg(a, b, x).and_then(|x| checked_beta(a, b).map(|y| x * y)) +} + +/// 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. +/// +/// # Panics +/// +/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, `x > 1.0`, or the numerical method +/// does not converge. +pub fn beta_reg(a: f64, b: f64, x: f64) -> f64 { + checked_beta_reg(a, b, x).unwrap() +} diff --git a/src/function/beta/asymptotic.rs b/src/function/beta/asymptotic.rs new file mode 100644 index 00000000..afd43c3b --- /dev/null +++ b/src/function/beta/asymptotic.rs @@ -0,0 +1,150 @@ +use super::*; + +pub(super) fn beta_shape_statistics(a: f64, b: f64) -> (f64, f64, f64, f64) { + let scale = a.max(b); + let scaled_a = a / scale; + let scaled_b = b / scale; + let scaled_sum = scaled_a + scaled_b; + let mean = scaled_a / scaled_sum; + let complement = scaled_b / scaled_sum; + let log_sum = scale.ln() + scaled_sum.ln(); + let root_sum = scale.sqrt() * scaled_sum.sqrt(); + (mean, complement, log_sum, root_sum) +} + +pub(super) fn beta_log_ratio(a: f64, b: f64, x: f64) -> (f64, f64) { + let residual = x.mul_add(b, -((1.0 - x) * a)); + let log_ratio = a * log1pmx(residual / a) + b * log1pmx(-residual / b); + (residual, log_ratio) +} + +fn beta_asymptotic_log_ratio(a: f64, b: f64, x: f64) -> (f64, (f64, f64)) { + let complement = two_sum(1.0, -x); + let left = dd_mul((x, 0.0), (b, 0.0)); + let right = dd_mul(complement, (a, 0.0)); + let residual_parts = dd_add(left, (-right.0, -right.1)); + let residual = residual_parts.0 + residual_parts.1; + let left_ratio = dd_div((residual_parts.0, residual_parts.1), (a, 0.0)); + let right_ratio = dd_div((-residual_parts.0, -residual_parts.1), (b, 0.0)); + let left_log = dd_add( + accurate_ln_one_plus_dd(left_ratio), + (-left_ratio.0, -left_ratio.1), + ); + let right_log = dd_add( + accurate_ln_one_plus_dd(right_ratio), + (-right_ratio.0, -right_ratio.1), + ); + let log_ratio = dd_add(dd_mul((a, 0.0), left_log), dd_mul((b, 0.0), right_log)); + (residual, log_ratio) +} + +fn beta_reg_symmetric_central(a: f64, b: f64, x: f64) -> Option { + if a != b || a < 100.0 { + return None; + } + + let delta = x - 0.5; + let delta_squared = delta * delta; + if 4.0 * a * delta_squared > 0.5 { + return None; + } + + // DLMF 8.17.1 and 5.5.5 give the symmetric-center integral; DLMF + // 5.11.13 supplies the asymptotic gamma ratio in its normalization. + let inverse_a = 1.0 / a; + let mut gamma_ratio: f64 = 869.0 / 4_194_304.0; + for coefficient in [ + -399.0 / 262_144.0, + -21.0 / 32_768.0, + 5.0 / 1_024.0, + 1.0 / 128.0, + -1.0 / 8.0, + 1.0, + ] { + gamma_ratio = gamma_ratio.mul_add(inverse_a, coefficient); + } + let central_density = 2.0 * (a / core::f64::consts::PI).sqrt() * gamma_ratio; + + let mut term = delta; + let mut integral = term; + for index in 1..=32 { + let n = f64::from(index); + term *= -4.0 * (a - n) * delta_squared * (2.0 * n - 1.0) / (n * (2.0 * n + 1.0)); + let previous = integral; + integral += term; + if integral == previous { + break; + } + } + Some(central_density.mul_add(integral, 0.5)) +} + +pub(super) fn beta_reg_asymptotic(a: f64, b: f64, x: f64) -> Option { + if let Some(result) = beta_reg_symmetric_central(a, b, x) { + return Some(result); + } + + let (mean, complement, _, root_sum) = beta_shape_statistics(a, b); + if a < ASYMPTOTIC_MIN_SUM && b < ASYMPTOTIC_MIN_SUM - a { + return None; + } + + if mean.min(complement) < 0.01 + || (mean.min(complement) < 0.1 && a.min(b) < ASYMPTOTIC_MIN_SHAPE) + { + return None; + } + + let (residual, log_ratio) = beta_asymptotic_log_ratio(a, b, x); + let scaled_deviance_parts = (-log_ratio.0, -log_ratio.1); + let scaled_deviance = scaled_deviance_parts.0 + scaled_deviance_parts.1; + if scaled_deviance > ASYMPTOTIC_MAX_DEVIANCE { + if scaled_deviance > -f64::from_bits(1).ln() { + return Some(if residual < 0.0 { 0.0 } else { 1.0 }); + } + return None; + } + + let use_centered = mean.min(complement) >= 0.01 && residual.abs() < 0.05 * a.min(b); + let (c0, c1) = if use_centered { + let (series_mean, delta_parts) = temme_delta(a, b, x); + temme_coefficients(series_mean, delta_parts.0 + delta_parts.1) + } else { + let eta = ((2.0 * scaled_deviance).sqrt() / root_sum).copysign(residual); + let c0 = 1.0 / eta - a.sqrt() * b.sqrt() / residual; + (c0, 0.0) + }; + let normal_argument = if scaled_deviance == 0.0 { + (0.0, 0.0) + } else { + let normal_root = scaled_deviance_parts.0.sqrt(); + let normal_root_error = (scaled_deviance_parts.1 + + (-normal_root).mul_add(normal_root, scaled_deviance_parts.0)) + / (2.0 * normal_root); + if residual < 0.0 { + (normal_root, normal_root_error) + } else { + (-normal_root, -normal_root_error) + } + }; + let absolute_normal_argument = if normal_argument.0 >= 0.0 { + normal_argument + } else { + (-normal_argument.0, -normal_argument.1) + }; + let tail = normal_tail(absolute_normal_argument); + let coefficient = (-c1 / (root_sum * root_sum)).mul_add(1.0, c0); + let correction = dd_exp(log_ratio) * coefficient / (consts::SQRT_2PI * root_sum); + let result_parts = if normal_argument.0 >= 0.0 { + dd_add((tail, 0.0), (correction, 0.0)) + } else { + let complement = dd_add((tail, 0.0), (-correction, 0.0)); + dd_add((1.0, 0.0), (-complement.0, -complement.1)) + }; + let result = result_parts.0 + result_parts.1; + if (0.0..=1.0).contains(&result) { + Some(result) + } else { + None + } +} diff --git a/src/function/beta/bgrat.rs b/src/function/beta/bgrat.rs new file mode 100644 index 00000000..1386d7b6 --- /dev/null +++ b/src/function/beta/bgrat.rs @@ -0,0 +1,219 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +const ACCURATE_SMALL_B_MAX_RATIO: f64 = 1e-4; + +pub(super) fn use_beta_small_b_shifted_accurate(a: f64, b: f64, y: f64) -> bool { + (0.0..1.0).contains(&a) && b <= ACCURATE_SMALL_B_MAX_RATIO * a && y < 0.3 +} + +pub(super) fn beta_small_b_large_a_factor( + a: f64, + b: f64, + x: f64, + y: f64, +) -> Result<(f64, f64), BetaFuncError> { + let bm1 = b - 1.0; + let t = a + 0.5 * bm1; + let lx = if y < 0.35 { (-y).ln_1p() } else { x.ln() }; + let u = -t * lx; + let log_h = b * u.ln() - u - ln_gamma_stable(b); + let log_prefix = log_h + ln_gamma_delta(a, b) - b * t.ln(); + + let mut odd_factorials = [1.0; 30]; + let mut factorial = 1.0; + for k in 1..=59 { + factorial *= k as f64; + if k >= 3 && k % 2 == 1 { + odd_factorials[(k - 3) as usize / 2] = factorial; + } + } + + let mut coefficients = [0.0; 30]; + coefficients[0] = 1.0; + let mut j = if u >= SCALED_GAMMA_MIN_X { + upper_gamma_scaled_asymptotic(b, u)? + } else if u > 1.0 { + upper_gamma_scaled_continued_fraction(b, u)? + } else if b <= 1e-4 && u <= 1.0 { + upper_gamma_scaled_small_shape(b, u)? + } else { + gamma::gamma_ur(b, u) / log_h.exp() + }; + let mut sum = j; + let mut compensation = 0.0_f64; + let lx2 = (0.5 * lx) * (0.5 * lx); + let mut lx_power = 1.0; + let t4 = 4.0 * t * t; + let mut b_plus_2n = b; + let mut converged = false; + + for n in 1..30 { + let n_f64 = n as f64; + let mut coefficient = 0.0; + for m in 1..n { + coefficient += (m as f64 * b - n_f64) * coefficients[n - m] / odd_factorials[m - 1]; + } + coefficient /= n_f64; + coefficient += bm1 / odd_factorials[n - 1]; + coefficients[n] = coefficient; + + j = (b_plus_2n * (b_plus_2n + 1.0) * j + (u + b_plus_2n + 1.0) * lx_power) / t4; + lx_power *= lx2; + b_plus_2n += 2.0; + let term = coefficient * j; + let corrected = term - compensation; + let next = sum + corrected; + compensation = (next - sum) - corrected; + sum = next; + if term.abs() <= prec::F64_PREC * sum.abs() { + converged = true; + break; + } + } + + if converged && sum > 0.0 { + Ok((log_prefix, sum)) + } else { + Err(BetaFuncError::ConvergenceFailed) + } +} + +pub(super) fn beta_small_b_large_a_series( + a: f64, + b: f64, + x: f64, + y: f64, + initial: f64, +) -> Result { + let (log_prefix, factor) = beta_small_b_large_a_factor(a, b, x, y)?; + let sum = initial + log_prefix.exp() * factor; + if (0.0..=1.0).contains(&sum) { + Ok(sum) + } else { + Err(BetaFuncError::ConvergenceFailed) + } +} + +pub(super) fn beta_small_b_large_a_series_log( + a: f64, + b: f64, + x: f64, + y: f64, + initial: f64, +) -> Result { + let (log_prefix, factor) = beta_small_b_large_a_factor(a, b, x, y)?; + let tail = log_prefix + factor.ln(); + if initial == 0.0 { + Ok(tail) + } else { + let initial = initial.ln(); + let maximum = initial.max(tail); + Ok(maximum + (initial.min(tail) - maximum).exp().ln_1p()) + } +} + +pub(super) fn beta_reg_small_b_shifted_log( + a: f64, + b: f64, + x: f64, + y: f64, + log_beta: (f64, f64), +) -> Result { + let steps = (10.0 - a).ceil() as usize; + let shifted = a + steps as f64; + let shifted_log = beta_small_b_large_a_series_log(shifted, b, x, y, 0.0)?; + let recurrence_log = beta_a_step_log(a, b, x, steps, log_beta); + let maximum = shifted_log.max(recurrence_log); + Ok(maximum + (shifted_log.min(recurrence_log) - maximum).exp().ln_1p()) +} + +pub(super) fn beta_reg_small_b_shifted_accurate( + a: f64, + b: f64, + x: f64, + y: f64, +) -> Result<(f64, (f64, f64)), BetaFuncError> { + let steps = (10.0 - a).ceil() as usize; + let shifted = a + steps as f64; + let (_, factor) = beta_small_b_large_a_factor(shifted, b, x, y)?; + let t = shifted + 0.5 * (b - 1.0); + let lx = if y < 0.35 { (-y).ln_1p() } else { x.ln() }; + let u = -t * lx; + let log_b = accurate_ln(b); + let log_gamma_b = dd_add((ln_gamma_one_plus_series(b), 0.0), (-log_b.0, -log_b.1)); + let log_h = dd_add( + dd_add(dd_mul((b, 0.0), accurate_ln(u)), (-u, 0.0)), + (-log_gamma_b.0, -log_gamma_b.1), + ); + let gamma_delta = ln_gamma_delta_accurate_parts(shifted, b); + let b_log_t = dd_mul((b, 0.0), accurate_ln(t)); + let log_prefix = dd_add(dd_add(log_h, gamma_delta), (-b_log_t.0, -b_log_t.1)); + let shifted_value = dd_exp(log_prefix) * factor; + let shifted_log = dd_add(log_prefix, accurate_ln(factor)); + let gamma = ln_gamma_accurate_parts(b); + let gamma_delta = ln_gamma_delta_accurate_parts(a, b); + let log_beta = dd_add(gamma, (-gamma_delta.0, -gamma_delta.1)); + let recurrence_log = beta_a_step_log_parts(a, b, x, steps, log_beta); + let recurrence_value = dd_exp(recurrence_log); + let value = dd_add((shifted_value, 0.0), (recurrence_value, 0.0)); + let (maximum, minimum) = if shifted_log.0 > recurrence_log.0 { + (shifted_log, recurrence_log) + } else { + (recurrence_log, shifted_log) + }; + let difference = dd_add(minimum, (-maximum.0, -maximum.1)); + let correction = accurate_ln_one_plus_dd((dd_exp(difference), 0.0)); + let logarithm = dd_add(maximum, correction); + Ok((value.0 + value.1, logarithm)) +} + +pub(super) fn beta_reg_small_b_large_a( + a: f64, + b: f64, + x: f64, + y: f64, +) -> Result, BetaFuncError> { + if a < 10.0 || b >= 40.0 || y >= 0.3 { + return Ok(None); + } + let mut steps = b.floor() as usize; + if b == steps as f64 { + steps -= 1; + } + let reduced_b = b - steps as f64; + let initial = if steps == 0 { + 0.0 + } else { + beta_a_step(reduced_b, a, y, steps) + }; + beta_small_b_large_a_series(a, reduced_b, x, y, initial).map(Some) +} + +pub(super) fn beta_reg_small_b_large_a_log( + a: f64, + b: f64, + x: f64, + y: f64, +) -> Result, BetaFuncError> { + if a < 10.0 || b >= 40.0 || y >= 0.3 { + return Ok(None); + } + let mut steps = b.floor() as usize; + if b == steps as f64 { + steps -= 1; + } + let reduced_b = b - steps as f64; + let initial = if steps == 0 { + 0.0 + } else { + beta_a_step(reduced_b, a, y, steps) + }; + beta_small_b_large_a_series_log(a, reduced_b, x, y, initial).map(Some) +} diff --git a/src/function/beta/dd.rs b/src/function/beta/dd.rs new file mode 100644 index 00000000..257f4368 --- /dev/null +++ b/src/function/beta/dd.rs @@ -0,0 +1,185 @@ +#[cfg(all(not(feature = "std"), not(test)))] +use super::Float; + +pub(super) fn accurate_ln_dd(value: (f64, f64)) -> (f64, f64) { + let logarithm = accurate_ln(value.0); + dd_add(logarithm, ((value.1 / value.0).ln_1p(), 0.0)) +} + +pub(super) fn accurate_ln_one_plus_dd(value: (f64, f64)) -> (f64, f64) { + if value.0 == 0.0 && value.1 == 0.0 { + return (0.0, 0.0); + } + if value.0.abs() > 0.5 { + return accurate_ln_dd(dd_add((1.0, 0.0), value)); + } + let ratio = dd_div(value, dd_add((2.0, 0.0), value)); + let ratio_squared = dd_mul(ratio, ratio); + let mut term = ratio; + let mut sum = ratio; + for index in 1..=24 { + term = dd_mul(term, ratio_squared); + if term.0 == 0.0 && term.1 == 0.0 { + break; + } + sum = dd_add(sum, dd_div_f64(term, f64::from(2 * index + 1))); + if term.0.abs() <= f64::EPSILON * f64::EPSILON * sum.0.abs() { + break; + } + } + dd_mul((2.0, 0.0), sum) +} + +pub(super) fn accurate_ln_one_minus_dd(value: f64) -> (f64, f64) { + if value <= 0.5 { + accurate_ln_one_plus_dd((-value, 0.0)) + } else { + let complement = two_sum(1.0, -value); + accurate_ln_dd(complement) + } +} + +pub(super) fn log1pmx(x: f64) -> f64 { + if x.abs() > 0.01 { + return x.ln_1p() - x; + } + + let mut term = -0.5 * x * x; + let mut sum = term; + for n in 3..=64 { + term *= -x * f64::from(n - 1) / f64::from(n); + sum += term; + } + sum +} + +pub(super) fn two_sum(left: f64, right: f64) -> (f64, f64) { + let sum = left + right; + let virtual_right = sum - left; + let error = (left - (sum - virtual_right)) + (right - virtual_right); + (sum, error) +} + +pub(super) fn dd_add( + (left, left_error): (f64, f64), + (right, right_error): (f64, f64), +) -> (f64, f64) { + let (sum, error) = two_sum(left, right); + two_sum(sum, error + left_error + right_error) +} + +pub(super) fn dd_mul( + (left, left_error): (f64, f64), + (right, right_error): (f64, f64), +) -> (f64, f64) { + let product = left * right; + let error = left.mul_add(right, -product) + + left * right_error + + left_error * right + + left_error * right_error; + two_sum(product, error) +} + +pub(super) fn dd_div_f64((numerator, numerator_error): (f64, f64), denominator: f64) -> (f64, f64) { + let quotient = numerator / denominator; + let remainder = (-quotient).mul_add(denominator, numerator) + numerator_error; + two_sum(quotient, remainder / denominator) +} + +pub(super) fn dd_div(numerator: (f64, f64), denominator: (f64, f64)) -> (f64, f64) { + let quotient = numerator.0 / denominator.0; + let product = dd_mul((quotient, 0.0), denominator); + let remainder = dd_add(numerator, (-product.0, -product.1)); + dd_add( + (quotient, 0.0), + ((remainder.0 + remainder.1) / denominator.0, 0.0), + ) +} + +pub(super) fn dd_exp((value, error): (f64, f64)) -> f64 { + let combined = value + error; + if combined < f64::from_bits(1).ln() - core::f64::consts::LN_2 { + return 0.0; + } + let exponential = value.exp(); + let error_expm1 = error.exp_m1(); + if exponential == 0.0 || !error_expm1.is_finite() { + return combined.exp(); + } + exponential.mul_add(error_expm1, exponential) +} + +pub(super) fn dd_negative_expm1((value, error): (f64, f64)) -> f64 { + let combined = value + error; + if combined < f64::from_bits(1).ln() - core::f64::consts::LN_2 { + return 1.0; + } + let exponential = value.exp(); + let error_expm1 = error.exp_m1(); + if exponential == 0.0 || !error_expm1.is_finite() { + return -combined.exp_m1(); + } + -value.exp_m1() - exponential * error_expm1 +} + +pub(super) fn accurate_ln(value: f64) -> (f64, f64) { + if value == 0.0 { + return (f64::NEG_INFINITY, 0.0); + } + if value == 1.0 { + return (0.0, 0.0); + } + let mut scaled = value; + let mut exponent_adjustment = 0_i32; + if scaled < f64::MIN_POSITIVE { + scaled *= 18_014_398_509_481_984.0; + exponent_adjustment = -54; + } + let value_bits = scaled.to_bits(); + let mut exponent = ((value_bits >> 52) & 0x7ff) as i32 - 1023 + exponent_adjustment; + let mut mantissa = f64::from_bits((value_bits & 0x000f_ffff_ffff_ffff) | (1023_u64 << 52)); + if mantissa > core::f64::consts::SQRT_2 { + mantissa *= 0.5; + exponent += 1; + } + let numerator = dd_add((mantissa, 0.0), (-1.0, 0.0)); + let denominator = dd_add((mantissa, 0.0), (1.0, 0.0)); + let ratio = dd_div(numerator, denominator); + let ratio_squared = dd_mul(ratio, ratio); + let mut term = ratio; + let mut sum = ratio; + for index in 1..=24 { + term = dd_mul(term, ratio_squared); + sum = dd_add(sum, dd_div_f64(term, f64::from(2 * index + 1))); + if term.0.abs() <= f64::EPSILON * f64::EPSILON * sum.0.abs() { + break; + } + } + let log_mantissa = dd_mul((2.0, 0.0), sum); + let log_two = (core::f64::consts::LN_2, 2.3190468138462996e-17); + dd_add(dd_mul((f64::from(exponent), 0.0), log_two), log_mantissa) +} + +pub(super) fn accurate_ln_one_minus(value: f64) -> (f64, f64) { + accurate_ln_one_minus_dd(value) +} + +pub(super) fn compensated_ln(value: f64) -> (f64, f64) { + let high = value.ln(); + let low = if value >= f64::MIN_POSITIVE && !(0.5..=2.0).contains(&value) { + value.mul_add((-high).exp(), -1.0).ln_1p() + } else { + 0.0 + }; + (high, low) +} + +pub(super) fn compensated_ln_one_minus(value: f64) -> (f64, f64) { + if value <= 0.5 { + ((-value).ln_1p(), 0.0) + } else { + let (complement, complement_error) = two_sum(1.0, -value); + let (high, low) = compensated_ln(complement); + (high, low + (complement_error / complement).ln_1p()) + } +} diff --git a/src/function/beta/forward.rs b/src/function/beta/forward.rs new file mode 100644 index 00000000..328eff94 --- /dev/null +++ b/src/function/beta/forward.rs @@ -0,0 +1,164 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +/// 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. +/// +/// # Errors +/// +/// If `a <= 0.0`, `b <= 0.0`, `x < 0.0`, `x > 1.0`, or the numerical method +/// does not converge. +pub fn checked_beta_reg(a: f64, b: f64, x: f64) -> Result { + if a <= 0.0 { + return Err(BetaFuncError::ANotGreaterThanZero); + } + + if b <= 0.0 { + return Err(BetaFuncError::BNotGreaterThanZero); + } + + if !(0.0..=1.0).contains(&x) { + return Err(BetaFuncError::XOutOfRange); + } + + if x == 0.0 { + return Ok(0.0); + } + if x == 1.0 { + return Ok(1.0); + } + if a == b && x == 0.5 { + return Ok(0.5); + } + if b == 1.0 { + return Ok(x.powf(a)); + } + if a == 1.0 { + return Ok(-(b * (-x).ln_1p()).exp_m1()); + } + let y = 1.0 - x; + if let Some((log_result, invert)) = beta_small_shapes_series_log(a, b, x, y)? { + let result = if invert { + -log_result.exp_m1() + } else { + log_result.exp() + }; + return if (0.0..=1.0).contains(&result) { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + }; + } + if let Some(result) = beta_reg_asymptotic(a, b, x) { + return Ok(result); + } + if a.mul_add(y, -(b * x)) >= 0.0 + && let Some(result) = beta_reg_small_b_large_a(a, b, x, y)? + { + return Ok(result); + } + if use_beta_small_b_shifted_accurate(a, b, y) { + return beta_reg_small_b_shifted_accurate(a, b, x, y).map(|result| result.0); + } + if (1.0..10.0).contains(&a) && b < 1.0 && y < 0.3 { + let result = beta_reg_small_b_shifted_log(a, b, x, y, ln_beta_accurate_parts(a, b))?.exp(); + return if (0.0..=1.0).contains(&result) { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + }; + } + let symm_transform = + !use_beta_power_series_before_symmetry(a, b, x) && use_beta_symmetry(a, b, x); + let (transformed_a, transformed_b, transformed_x, transformed_y) = if symm_transform { + (b, a, y, x) + } else { + (a, b, x, y) + }; + if !use_exact_complement_continued_fraction(a, b, symm_transform) + && let Some(tail) = + beta_reg_small_b_large_a(transformed_a, transformed_b, transformed_x, transformed_y)? + { + return Ok(if symm_transform { 1.0 - tail } else { tail }); + } + if use_beta_power_series(transformed_a, transformed_b, transformed_x) { + let log_result = beta_power_series_log_parts(transformed_a, transformed_b, transformed_x)?; + let result = if symm_transform { + dd_negative_expm1(log_result) + } else { + (log_result.0 + log_result.1).exp() + }; + return if (0.0..=1.0).contains(&result) { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + }; + } + + let lanczos_power = use_beta_reg_lanczos_power(transformed_a, transformed_b).then(|| { + let coordinates = if symm_transform { + (two_sum(1.0, -x), (x, 0.0)) + } else { + ((x, 0.0), two_sum(1.0, -x)) + }; + beta_reg_lanczos_power(transformed_a, transformed_b, coordinates.0, coordinates.1) + }); + let log_power = lanczos_power + .is_none() + .then(|| beta_reg_log_power_parts(a, b, x)); + let power = lanczos_power + .map(|value| value.0 + value.1) + .unwrap_or_else(|| { + let value = log_power.unwrap(); + (value.0 + value.1).exp() + }); + if power == 0.0 { + return Ok(if symm_transform { 1.0 } else { 0.0 }); + } + let fraction = beta_fraction_for_transformed_tail( + a, + b, + x, + transformed_a, + transformed_b, + transformed_x, + symm_transform, + )?; + let accurate_fraction = + 1.0 - transformed_x == 1.0 || use_exact_complement_continued_fraction(a, b, symm_transform); + let result = if let Some(lanczos_power) = lanczos_power { + let tail = dd_div(lanczos_power, fraction); + if symm_transform { + dd_add((1.0, 0.0), (-tail.0, -tail.1)).0 + } else { + tail.0 + tail.1 + } + } else if accurate_fraction { + let log_power = log_power.unwrap(); + let log_fraction = accurate_ln_dd(fraction); + let log_result = dd_add(log_power, (-log_fraction.0, -log_fraction.1)); + if symm_transform { + dd_negative_expm1(log_result) + } else { + dd_exp(log_result) + } + } else if symm_transform { + 1.0 - power / (fraction.0 + fraction.1) + } else { + power / (fraction.0 + fraction.1) + }; + if (0.0..=1.0).contains(&result) { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + } +} diff --git a/src/function/beta/fraction.rs b/src/function/beta/fraction.rs new file mode 100644 index 00000000..9956852f --- /dev/null +++ b/src/function/beta/fraction.rs @@ -0,0 +1,140 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +pub(super) fn beta_continued_fraction(a: f64, b: f64, x: f64) -> Result { + let y = 1.0 - x; + let tiny = 16.0 * f64::MIN_POSITIVE; + let mut fraction = a * (a * y - b * x + 1.0) / (a + 1.0); + if fraction == 0.0 { + fraction = tiny; + } + let mut c = fraction; + let mut d = 0.0; + + for m in 1..=MAX_BETA_REG_ITERATIONS { + let m = f64::from(m); + let denominator = a + 2.0 * m - 1.0; + let numerator = + (m * (a + m - 1.0) / denominator) * ((a + b + m - 1.0) / denominator) * (b - m) * x * x; + let denominator_term = m + + m * (b - m) * x / denominator + + (a + m) * (a * y - b * x + 1.0 + m * (2.0 - x)) / (a + 2.0 * m + 1.0); + + d = denominator_term + numerator * d; + if d == 0.0 { + d = tiny; + } + c = denominator_term + numerator / c; + if c == 0.0 { + c = tiny; + } + d = 1.0 / d; + let delta = c * d; + fraction *= delta; + + if (delta - 1.0).abs() <= prec::F64_PREC { + return Ok(fraction); + } + } + + Err(BetaFuncError::ConvergenceFailed) +} + +pub(super) fn beta_continued_fraction_dd( + a: f64, + b: f64, + x: (f64, f64), +) -> Result<(f64, f64), BetaFuncError> { + let y = dd_add((1.0, 0.0), (-x.0, -x.1)); + let mut residual = dd_mul((a, 0.0), y); + residual = dd_add(residual, dd_mul((-b, 0.0), x)); + residual = dd_add(residual, (1.0, 0.0)); + let tiny = 16.0 * f64::MIN_POSITIVE; + let mut fraction = dd_div_f64(dd_mul((a, 0.0), residual), a + 1.0); + if fraction.0 == 0.0 { + fraction = (tiny, 0.0); + } + let mut c = fraction; + let mut d = (0.0, 0.0); + + for integer in 1..=MAX_BETA_REG_ITERATIONS { + let m = f64::from(integer); + let denominator = a + 2.0 * m - 1.0; + let mut numerator = dd_div_f64(dd_mul((m, 0.0), (a + m - 1.0, 0.0)), denominator); + let a_plus_b_plus_m_minus_one = dd_add((b, 0.0), dd_add((a, 0.0), (m - 1.0, 0.0))); + numerator = dd_mul( + numerator, + dd_div_f64(dd_mul(a_plus_b_plus_m_minus_one, x), denominator), + ); + let b_minus_m = dd_add((b, 0.0), (-m, 0.0)); + numerator = dd_mul(numerator, dd_mul(b_minus_m, x)); + + let first = dd_div_f64(dd_mul((m, 0.0), dd_mul(b_minus_m, x)), denominator); + let inner = dd_add(residual, dd_mul((m, 0.0), dd_add((2.0, 0.0), (-x.0, -x.1)))); + let second = dd_div_f64(dd_mul((a + m, 0.0), inner), a + 2.0 * m + 1.0); + let denominator_term = dd_add((m, 0.0), dd_add(first, second)); + + let mut scaled_denominator = dd_add(denominator_term, dd_mul(numerator, d)); + if scaled_denominator.0 == 0.0 { + scaled_denominator = (tiny, 0.0); + } + d = dd_div((1.0, 0.0), scaled_denominator); + c = dd_add(denominator_term, dd_div(numerator, c)); + if c.0 == 0.0 { + c = (tiny, 0.0); + } + let delta = dd_mul(c, d); + fraction = dd_mul(fraction, delta); + let convergence = dd_add(delta, (-1.0, 0.0)); + if (convergence.0 + convergence.1).abs() <= f64::EPSILON { + return Ok(fraction); + } + } + + Err(BetaFuncError::ConvergenceFailed) +} + +pub(super) fn selected_beta_continued_fraction( + a: f64, + b: f64, + x: f64, +) -> Result<(f64, f64), BetaFuncError> { + if x <= f64::EPSILON { + beta_continued_fraction_dd(a, b, (x, 0.0)) + } else { + beta_continued_fraction(a, b, x).map(|fraction| (fraction, 0.0)) + } +} + +pub(super) fn use_exact_complement_continued_fraction( + a: f64, + b: f64, + symm_transform: bool, +) -> bool { + symm_transform && a >= 1.0 && b >= 2.0 * (a + 1.0) +} + +pub(super) fn beta_fraction_for_transformed_tail( + a: f64, + b: f64, + x: f64, + transformed_a: f64, + transformed_b: f64, + transformed_x: f64, + symm_transform: bool, +) -> Result<(f64, f64), BetaFuncError> { + if use_exact_complement_continued_fraction(a, b, symm_transform) { + beta_continued_fraction_dd(transformed_a, transformed_b, two_sum(1.0, -x)) + } else { + selected_beta_continued_fraction(transformed_a, transformed_b, transformed_x) + } +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/fraction/tests.rs b/src/function/beta/fraction/tests.rs new file mode 100644 index 00000000..9eaee9f9 --- /dev/null +++ b/src/function/beta/fraction/tests.rs @@ -0,0 +1,10 @@ +use super::*; + +#[test] +fn double_double_fraction_handles_zero_initial_denominator() { + let fraction = beta_continued_fraction_dd(1.0, 3.0, (0.5, 0.0)).unwrap(); + + assert!(fraction.0.is_finite()); + assert!(fraction.1.is_finite()); + assert!(fraction.0 + fraction.1 > 0.0); +} diff --git a/src/function/beta/inverse/initial.rs b/src/function/beta/inverse/initial.rs new file mode 100644 index 00000000..44f2d19b --- /dev/null +++ b/src/function/beta/inverse/initial.rs @@ -0,0 +1,95 @@ +use super::super::*; + +// The initial estimate in `inverse_beta_initial` is derived from the +// implementation in the `special` crate, which in turn follows John Burkardt's +// implementation of Applied Statistics Algorithms AS 64 and AS 109: +// +// - https://docs.rs/special/0.8.1/ +// - https://people.sc.fsu.edu/~jburkardt/c_src/asa109/asa109.html +// - https://www.jstor.org/stable/2346798 +// - https://www.jstor.org/stable/2346887 +// +// Copyright 2014–2019 The special Developers +// +// Permission is hereby granted, free of charge, to any person obtaining a copy of +// this software and associated documentation files (the "Software"), to deal in +// the Software without restriction, including without limitation the rights to +// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +// the Software, and to permit persons to whom the Software is furnished to do so, +// subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + +pub(super) fn lower_tail_initial(a: f64, b: f64, probability: f64, ln_beta: f64) -> (f64, f64) { + let log_initial = (probability.ln() + a.ln() + ln_beta) / a; + let initial = log_initial.exp(); + let initial = if initial == 0.0 { + 0.0 + } else if initial < 1.0 { + initial + } else { + let (mean, _, _, _) = beta_shape_statistics(a, b); + if mean < 1.0 { + mean + } else { + f64::from_bits(1.0_f64.to_bits() - 1) + } + }; + (initial, log_initial) +} + +pub(super) fn lower_tail_initial_accurate( + a: f64, + probability: f64, + log_beta: (f64, f64), +) -> (f64, (f64, f64)) { + let mut logarithm = accurate_ln(probability); + logarithm = dd_add(logarithm, accurate_ln(a)); + logarithm = dd_add(logarithm, log_beta); + logarithm = dd_div_f64(logarithm, a); + (dd_exp(logarithm), logarithm) +} + +pub(super) fn lower_tail_initial_from_log_normalizer( + a: f64, + probability: f64, + log_normalizer: (f64, f64), +) -> (f64, (f64, f64)) { + let logarithm = dd_div_f64(dd_add(accurate_ln(probability), log_normalizer), a); + (dd_exp(logarithm), logarithm) +} + +pub(super) fn inverse_beta_initial(a: f64, b: f64, probability: f64, ln_beta: f64) -> (f64, f64) { + if a > 1.0 && b > 1.0 && (probability >= 1e-4 || a.min(b) >= STIRLING_MIN) { + let normal_tail = (-2.0 * probability.ln()).sqrt(); + let normal_quantile = normal_tail + - (2.30753 + 0.27061 * normal_tail) + / (1.0 + (0.99229 + 0.04481 * normal_tail) * normal_tail); + let correction = (normal_quantile * normal_quantile - 3.0) / 6.0; + let reciprocal_a = 1.0 / (2.0 * a - 1.0); + let reciprocal_b = 1.0 / (2.0 * b - 1.0); + let scale = 2.0 / (reciprocal_a + reciprocal_b); + let w = normal_quantile * (scale + correction).sqrt() / scale + - (reciprocal_b - reciprocal_a) * (correction + 5.0 / 6.0 - 2.0 / (3.0 * scale)); + let log_ratio = b.ln() - a.ln() + 2.0 * w; + let initial = if log_ratio > 0.0 { + let reciprocal = (-log_ratio).exp(); + reciprocal / (1.0 + reciprocal) + } else { + 1.0 / (1.0 + log_ratio.exp()) + }; + if initial > 0.0 && initial < 1.0 { + return (initial, f64::NAN); + } + } + + lower_tail_initial(a, b, probability, ln_beta) +} diff --git a/src/function/beta/inverse/mod.rs b/src/function/beta/inverse/mod.rs new file mode 100644 index 00000000..235e4849 --- /dev/null +++ b/src/function/beta/inverse/mod.rs @@ -0,0 +1,122 @@ +mod initial; +mod shape_two; +mod solve; + +use super::*; +use initial::*; +use shape_two::*; +use solve::*; + +// Near one, the reflected solver preserves upper-tail information needed to round to 1.0. +const SHAPE_TWO_SPECIALIZATION_MAX: f64 = 0.999_999_999; + +fn small_first_shape_lower_tail(a: f64, b: f64, probability: f64) -> Option { + if probability > 0.5 || a > 0.125 || b < STIRLING_MIN { + return None; + } + let log_normalizer = ln_a_beta_small_first_shape_parts(a, b); + let current = lower_tail_initial_from_log_normalizer(a, probability, log_normalizer).0; + if current <= 0.0 || current >= 1.0 { + return None; + } + let first_correction = ((b - 1.0).abs() / (a + 1.0)) * current; + let remainder_ratio = (b - 2.0).abs() * current; + (first_correction <= f64::EPSILON / 32.0 && remainder_ratio <= 0.5).then_some(current) +} + +/// Computes the inverse of the regularized incomplete beta function. +/// +/// # Panics +/// +/// Panics for arguments outside `a > 0`, `b > 0`, and `0 <= probability <= 1` +/// in debug builds, or if the numerical method does not converge. +pub fn inv_beta_reg(a: f64, b: f64, probability: f64) -> f64 { + debug_assert!((0.0..=1.0).contains(&probability) && a > 0.0 && b > 0.0); + + if probability == 0.0 { + return 0.0; + } + if probability == 1.0 { + return 1.0; + } + if a == b && probability == 0.5 { + return 0.5; + } + if b == 2.0 && probability > SHAPE_TWO_SPECIALIZATION_MAX { + let initial = ((probability.ln() - (a + 1.0).ln()) / a).exp(); + if let Some(quantile) = shape_two_lower_endpoint_result(a, b, probability, initial) { + return quantile; + } + } + if let Some(quantile) = beta_concentrated_quantile(a, b, probability) { + return quantile; + } + if b == 1.0 { + return probability.powf(1.0 / a); + } + if a == 1.0 + && probability <= f64::MIN_POSITIVE + && b == b.trunc() + && b <= 9_007_199_254_740_992.0 + { + let probability_bits = probability.to_bits(); + let divisor = b as u64; + let quotient = probability_bits / divisor; + let remainder = probability_bits % divisor; + return f64::from_bits(quotient + u64::from(2 * remainder >= divisor)); + } + if a == 1.0 { + return -((-probability).ln_1p() / b).exp_m1(); + } + if (a == 2.0 || b == 2.0) && probability <= SHAPE_TWO_SPECIALIZATION_MAX { + return inverse_beta_shape_two(a, b, probability); + } + if let Some(quantile) = small_first_shape_lower_tail(a, b, probability) { + return quantile; + } + + let log_beta = ln_beta_inverse_parts(a, b); + let flip = inverse_beta_reflect(a, b, probability, log_beta); + let (a, b, target) = if flip { + (b, a, 1.0 - probability) + } else { + (a, b, probability) + }; + let ln_beta = log_beta.0 + log_beta.1; + let (mut current, mut log_initial) = inverse_beta_initial(a, b, target, ln_beta); + let smaller = a.min(b); + let larger = a.max(b); + if log_initial.is_finite() + && larger >= STIRLING_MIN + && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) + { + let accurate_initial = lower_tail_initial_accurate(a, target, log_beta); + if accurate_initial.0 > 0.0 && accurate_initial.0 < 1.0 { + current = accurate_initial.0; + log_initial = accurate_initial.1.0 + accurate_initial.1.1; + let first_correction = ((b - 1.0).abs() / (a + 1.0)) * current; + let remainder_ratio = (b - 2.0).abs() * current; + if first_correction <= f64::EPSILON / 32.0 && remainder_ratio <= 0.5 { + return if flip { 1.0 - current } else { current }; + } + } + } + let min_subnormal = f64::from_bits(1); + if current == 0.0 && log_initial < min_subnormal.ln() - core::f64::consts::LN_2 { + return if flip { 1.0 } else { 0.0 }; + } + let first_correction = ((b - 1.0).abs() / (a + 1.0)) * current; + let remainder_ratio = (b - 2.0).abs() * current; + if first_correction <= f64::EPSILON / 32.0 && remainder_ratio <= 0.5 { + return if flip { 1.0 - current } else { current }; + } + if current < f64::MIN_POSITIVE { + let relative_correction = b * current / (a + 1.0); + let relative_half_ulp = 0.5 * (min_subnormal / current); + if relative_correction < 0.25 * relative_half_ulp { + return if flip { 1.0 - current } else { current }; + } + } + let result = inverse_beta_log_tail(a, b, target, current, log_beta, ln_beta); + if flip { 1.0 - result } else { result } +} diff --git a/src/function/beta/inverse/shape_two.rs b/src/function/beta/inverse/shape_two.rs new file mode 100644 index 00000000..f35b9eed --- /dev/null +++ b/src/function/beta/inverse/shape_two.rs @@ -0,0 +1,114 @@ +//! Shape-two identities follow from DLMF 8.17.7--8.17.8 and symmetry 8.17.4. +//! The large-`b` endpoint path expands `I_x(2,b) = 1 - (1-x)^b(1+bx)` in +//! `y = bx` and uses outward fixed-point intervals to certify rounding cells. + +use super::super::*; +use adjacent::adjacent_result; +use endpoint::lower_endpoint_result; +use value::{fast_cdf_and_pdf, log_cdf}; + +mod adjacent; +mod endpoint; +mod endpoint_beta_two; +mod endpoint_certified; +mod endpoint_fixed; +mod value; + +#[cfg(test)] +use value::log_cdf_parts; + +pub(super) fn shape_two_lower_endpoint_result( + a: f64, + b: f64, + probability: f64, + initial: f64, +) -> Option { + lower_endpoint_result(a, b, probability, initial) +} + +pub(super) fn inverse_beta_shape_two(a: f64, b: f64, probability: f64) -> f64 { + let mut current = if b == 2.0 { + ((probability.ln() - (a + 1.0).ln()) / a).exp() + } else { + (0.5 * (probability.ln() + core::f64::consts::LN_2 - b.ln() - (b + 1.0).ln())).exp() + }; + if let Some(result) = shape_two_lower_endpoint_result(a, b, probability, current) { + return result; + } + current = current.clamp(f64::from_bits(1), f64::from_bits(1.0_f64.to_bits() - 1)); + let mut lower = 0.0; + let mut upper = 1.0; + for _ in 0..128 { + if let Some((cdf, pdf)) = fast_cdf_and_pdf(a, b, current) { + let error = cdf - probability; + if error == 0.0 { + return adjacent_result(a, b, probability, current); + } + if error < 0.0 { + lower = current; + } else { + upper = current; + } + if upper.to_bits().abs_diff(lower.to_bits()) == 1 { + return adjacent_result(a, b, probability, current); + } + let step = error / pdf; + let pdf_ratio = (a - 1.0) / current - (b - 1.0) / (1.0 - current); + let denominator = 1.0 - 0.5 * step * pdf_ratio; + let candidate = current - step / denominator; + if candidate == current { + return adjacent_result(a, b, probability, current); + } + let next = if denominator > 0.0 && candidate > lower && candidate < upper { + candidate + } else { + lower + 0.5 * (upper - lower) + }; + if next == current { + return adjacent_result(a, b, probability, current); + } + current = next; + continue; + } + let target = accurate_ln(probability); + let log_value = log_cdf(a, b, current); + let error = log_value - (target.0 + target.1); + if error == 0.0 { + return adjacent_result(a, b, probability, current); + } + if error < 0.0 { + lower = current; + } else { + upper = current; + } + if upper.to_bits().abs_diff(lower.to_bits()) == 1 { + return adjacent_result(a, b, probability, current); + } + let log_pdf = if b == 2.0 { + (a - 1.0).mul_add(current.ln(), a.ln() + (a + 1.0).ln() + (-current).ln_1p()) + } else { + (b - 1.0).mul_add((-current).ln_1p(), b.ln() + (b + 1.0).ln() + current.ln()) + }; + let inverse_derivative = (log_value - log_pdf).exp(); + let step = error * inverse_derivative; + let log_pdf_derivative = (a - 1.0) / current - (b - 1.0) / (1.0 - current); + let denominator = 1.0 - 0.5 * step * (log_pdf_derivative - 1.0 / inverse_derivative); + let candidate = current - step / denominator; + if candidate == current { + return adjacent_result(a, b, probability, current); + } + let next = if denominator > 0.0 && candidate > lower && candidate < upper { + candidate + } else { + lower + 0.5 * (upper - lower) + }; + if next == current { + return adjacent_result(a, b, probability, current); + } + current = next; + } + panic!("shape-two inverse did not converge for a={a}, b={b}, probability={probability}") +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/inverse/shape_two/adjacent.rs b/src/function/beta/inverse/shape_two/adjacent.rs new file mode 100644 index 00000000..76c408a2 --- /dev/null +++ b/src/function/beta/inverse/shape_two/adjacent.rs @@ -0,0 +1,182 @@ +use super::super::super::*; +use super::super::inverse_beta_adjacent_result; +use super::value::{direct_cdf_and_pdf, log_cdf_parts}; + +#[derive(Copy, Clone, PartialEq, Eq)] +enum ErrorScale { + Boundary, + Series, + Tail, + Power, + Log, +} + +#[derive(Copy, Clone)] +struct Endpoint { + value: f64, + error: f64, + scale: ErrorScale, +} + +#[derive(Copy, Clone)] +struct Evaluation { + endpoint: Endpoint, + pdf: Option, +} + +fn get_log_target(probability: f64, target: &mut Option<(f64, f64)>) -> (f64, f64) { + *target.get_or_insert_with(|| accurate_ln(probability)) +} + +fn evaluate( + a: f64, + b: f64, + value: f64, + probability: f64, + target: &mut Option<(f64, f64)>, +) -> Evaluation { + let (error, pdf, scale) = if let Some((cdf, pdf)) = direct_cdf_and_pdf(a, b, value) { + let scale = if b == 2.0 { + ErrorScale::Power + } else if value < 0.5 && b * value < 0.5 { + ErrorScale::Series + } else { + ErrorScale::Tail + }; + (dd_add(cdf, (-probability, 0.0)), Some(pdf), scale) + } else { + let target = get_log_target(probability, target); + ( + dd_add(log_cdf_parts(a, b, value), (-target.0, -target.1)), + None, + ErrorScale::Log, + ) + }; + Evaluation { + endpoint: Endpoint { + value, + error: error.0 + error.1, + scale, + }, + pdf, + } +} + +fn pair_result( + a: f64, + b: f64, + mut lower: Endpoint, + mut upper: Endpoint, + probability: f64, + target: &mut Option<(f64, f64)>, +) -> f64 { + if lower.scale != upper.scale + || lower.scale == ErrorScale::Boundary + || upper.scale == ErrorScale::Boundary + { + let target = get_log_target(probability, target); + for endpoint in [&mut lower, &mut upper] { + let error = dd_add(log_cdf_parts(a, b, endpoint.value), (-target.0, -target.1)); + endpoint.error = error.0 + error.1; + } + } + inverse_beta_adjacent_result(lower.value, upper.value, lower.error, upper.error) +} + +fn neighboring_result( + a: f64, + b: f64, + probability: f64, + current: Endpoint, + neighbor: f64, + target: &mut Option<(f64, f64)>, +) -> Option { + let neighbor = evaluate(a, b, neighbor, probability, target).endpoint; + if current.error * neighbor.error > 0.0 { + return None; + } + let (lower, upper) = if neighbor.value < current.value { + (neighbor, current) + } else { + (current, neighbor) + }; + Some(pair_result(a, b, lower, upper, probability, target)) +} + +pub(super) fn adjacent_result(a: f64, b: f64, probability: f64, mut current: f64) -> f64 { + let mut log_target = None; + let mut lower = Endpoint { + value: 0.0, + error: f64::NEG_INFINITY, + scale: ErrorScale::Boundary, + }; + let mut upper = Endpoint { + value: 1.0, + error: f64::INFINITY, + scale: ErrorScale::Boundary, + }; + for _ in 0..64 { + if current == 0.0 || current == 1.0 { + return current; + } + let evaluation = evaluate(a, b, current, probability, &mut log_target); + let endpoint = evaluation.endpoint; + if endpoint.error < 0.0 { + lower = endpoint; + } else { + upper = endpoint; + } + if upper.value.to_bits().abs_diff(lower.value.to_bits()) == 1 { + return pair_result(a, b, lower, upper, probability, &mut log_target); + } + let step = if let Some(pdf) = evaluation.pdf { + endpoint.error / pdf + } else { + let target = get_log_target(probability, &mut log_target); + let log_pdf = if b == 2.0 { + (a - 1.0).mul_add(current.ln(), a.ln() + (a + 1.0).ln() + (-current).ln_1p()) + } else { + (b - 1.0).mul_add((-current).ln_1p(), b.ln() + (b + 1.0).ln() + current.ln()) + }; + endpoint.error * ((target.0 + target.1) - log_pdf).exp() + }; + let candidate = current - step; + if candidate == current { + let neighbor = f64::from_bits(if endpoint.error > 0.0 { + current.to_bits() - 1 + } else { + current.to_bits() + 1 + }); + if let Some(result) = + neighboring_result(a, b, probability, endpoint, neighbor, &mut log_target) + { + return result; + } + current = neighbor; + continue; + } + let next = if candidate.is_finite() && candidate > lower.value && candidate < upper.value { + candidate + } else { + lower.value + 0.5 * (upper.value - lower.value) + }; + if next == current { + let neighbor = f64::from_bits(if endpoint.error > 0.0 { + current.to_bits() - 1 + } else { + current.to_bits() + 1 + }); + if let Some(result) = + neighboring_result(a, b, probability, endpoint, neighbor, &mut log_target) + { + return result; + } + current = lower.value + 0.5 * (upper.value - lower.value); + } else { + current = next; + } + } + panic!( + "shape-two inverse did not resolve adjacent values for a={a}, b={b}, probability={probability}" + ) +} diff --git a/src/function/beta/inverse/shape_two/endpoint.rs b/src/function/beta/inverse/shape_two/endpoint.rs new file mode 100644 index 00000000..64179439 --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint.rs @@ -0,0 +1,165 @@ +use super::super::super::*; +use super::endpoint_beta_two; +use super::endpoint_certified; +use super::endpoint_certified::Certificate; + +const MIN_SUBNORMAL: f64 = f64::from_bits(1); +const MIN_NORMAL_BITS: u64 = f64::MIN_POSITIVE.to_bits(); +const FIRST_NORMAL_BIN_END_BITS: u64 = f64::MIN_POSITIVE.to_bits() + (1_u64 << 52) - 1; +const BETA_TWO_ENDPOINT_SHAPE: f64 = f64::from_bits(1507_u64 << 52); +const LOG_TWO: (f64, f64) = (core::f64::consts::LN_2, 2.3190468138462996e-17); + +fn log_cdf_at_min_subnormal_multiple(a: f64, multiple: (f64, f64)) -> (f64, f64) { + let log_x = dd_add(accurate_ln_dd(multiple), dd_mul((-1074.0, 0.0), LOG_TWO)); + if a > f64::MAX / -log_x.0 { + return (f64::NEG_INFINITY, 0.0); + } + let ax = dd_mul(dd_mul((a, 0.0), (MIN_SUBNORMAL, 0.0)), multiple); + let factor = dd_add(dd_add((1.0, 0.0), (a, 0.0)), (-ax.0, -ax.1)); + dd_add(dd_mul((a, 0.0), log_x), accurate_ln_dd(factor)) +} + +fn compare_logs(left: (f64, f64), right: (f64, f64)) -> core::cmp::Ordering { + if right.0 == f64::NEG_INFINITY { + return core::cmp::Ordering::Greater; + } + let difference = dd_add(left, (-right.0, -right.1)); + let difference = difference.0 + difference.1; + if difference < 0.0 { + core::cmp::Ordering::Less + } else if difference > 0.0 { + core::cmp::Ordering::Greater + } else { + core::cmp::Ordering::Equal + } +} + +fn midpoint_log_cdf(a: f64, lower_bits: u64) -> (f64, f64) { + let multiple = dd_add((lower_bits as f64, 0.0), (0.5, 0.0)); + log_cdf_at_min_subnormal_multiple(a, multiple) +} + +pub(super) fn overlap_result(lower_bits: u64) -> f64 { + let result_bits = if lower_bits & 1 == 0 { + lower_bits + } else { + lower_bits + 1 + }; + f64::from_bits(result_bits) +} + +fn midpoint_order( + a: f64, + b: f64, + probability: f64, + target: (f64, f64), + lower_bits: u64, +) -> Option> { + if a == 2.0 { + if let Some(order) = endpoint_beta_two::certified_midpoint_order(b, probability, lower_bits) + { + return Some(Ok(order)); + } + return Some( + match endpoint_certified::midpoint_certificate(b, probability, lower_bits).ok()? { + Certificate::Ordered(order) => Ok(order), + Certificate::Overlap => Err(overlap_result(lower_bits)), + }, + ); + } + let order = compare_logs(target, midpoint_log_cdf(a, lower_bits)); + if order.is_eq() { + let even = if lower_bits & 1 == 0 { + lower_bits + } else { + lower_bits + 1 + }; + Some(Err(f64::from_bits(even))) + } else { + Some(Ok(order)) + } +} + +pub(super) fn lower_endpoint_result(a: f64, b: f64, probability: f64, initial: f64) -> Option { + let last_candidate = if a == 2.0 { + FIRST_NORMAL_BIN_END_BITS + } else { + MIN_NORMAL_BITS - 1 + }; + if !b.is_finite() + || a == 2.0 && b < BETA_TWO_ENDPOINT_SHAPE + || a != 2.0 && (b != 2.0 || initial >= f64::MIN_POSITIVE) + { + return None; + } + let target = accurate_ln(probability); + let candidate = if a == 2.0 { + endpoint_beta_two::initial_candidate(b, probability) + } else { + initial + } + .to_bits() + .min(last_candidate); + let (mut lower, mut upper) = match midpoint_order(a, b, probability, target, candidate)? { + Err(result) => return Some(result), + Ok(core::cmp::Ordering::Less) => { + let mut upper = candidate; + let mut step = 1_u64; + loop { + if upper == 0 { + if a == 2.0 { + break (0, 0); + } + return Some(0.0); + } + let probe = upper.saturating_sub(step); + match midpoint_order(a, b, probability, target, probe)? { + Err(result) => return Some(result), + Ok(core::cmp::Ordering::Less) => { + upper = probe; + step = step.saturating_mul(2); + } + Ok(core::cmp::Ordering::Greater) => break (probe + 1, upper), + Ok(core::cmp::Ordering::Equal) => unreachable!(), + } + } + } + Ok(core::cmp::Ordering::Greater) => { + let mut lower = candidate + 1; + let mut probe = candidate; + let mut step = 1_u64; + loop { + if probe == last_candidate { + if a == 2.0 { + break (last_candidate + 1, last_candidate + 1); + } + return None; + } + probe = probe.saturating_add(step).min(last_candidate); + match midpoint_order(a, b, probability, target, probe)? { + Err(result) => return Some(result), + Ok(core::cmp::Ordering::Less) => break (lower, probe), + Ok(core::cmp::Ordering::Greater) => { + lower = probe + 1; + step = step.saturating_mul(2); + } + Ok(core::cmp::Ordering::Equal) => unreachable!(), + } + } + } + Ok(core::cmp::Ordering::Equal) => unreachable!(), + }; + while lower < upper { + let midpoint = lower + (upper - lower) / 2; + match midpoint_order(a, b, probability, target, midpoint)? { + Err(result) => return Some(result), + Ok(core::cmp::Ordering::Less) => upper = midpoint, + Ok(core::cmp::Ordering::Greater) => lower = midpoint + 1, + Ok(core::cmp::Ordering::Equal) => unreachable!(), + } + } + (lower <= last_candidate).then(|| f64::from_bits(lower)) +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/inverse/shape_two/endpoint/tests.rs b/src/function/beta/inverse/shape_two/endpoint/tests.rs new file mode 100644 index 00000000..59187642 --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint/tests.rs @@ -0,0 +1,12 @@ +use super::overlap_result; + +#[test] +fn overlap_policy_is_monotone_and_ties_to_even() { + let cases = [(0_u64, 0_u64), (1, 2), (2, 2), (3, 4)]; + let results = cases.map(|(lower, expected)| { + let result = overlap_result(lower); + assert_eq!(result.to_bits(), expected); + result + }); + assert!(results.windows(2).all(|pair| pair[0] <= pair[1])); +} diff --git a/src/function/beta/inverse/shape_two/endpoint_beta_two.rs b/src/function/beta/inverse/shape_two/endpoint_beta_two.rs new file mode 100644 index 00000000..cdf2dced --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint_beta_two.rs @@ -0,0 +1,149 @@ +use super::super::super::*; + +const CERTIFIED_DD_RADIUS: f64 = f64::from_bits((1023_u64 - 79) << 52); +const CERTIFIED_MEDIAN_RADIUS: f64 = f64::from_bits((1023_u64 - 95) << 52); +const ENDPOINT_LIMIT_BITS: u64 = 0x0020_0000_0000_0000; +const MINIMUM_SHAPE_BITS: u64 = (1023_u64 + 484) << 52; +const MEDIAN_SCALE: f64 = f64::from_bits(0x3ffada825f9762b2); +const MEDIAN_SCALE_LOW: f64 = f64::from_bits(0x3c9f4a493534d79d); + +fn normal_power_of_two(exponent: i32) -> f64 { + f64::from_bits(((exponent + 1023) as u64) << 52) +} + +fn binary_exponent(value: f64) -> i32 { + ((value.to_bits() >> 52) & 0x7ff) as i32 - 1023 +} + +fn normalized_probability(value: f64) -> (f64, i32) { + let (scaled, adjustment) = if value < f64::MIN_POSITIVE { + (value * 18_014_398_509_481_984.0, -54) + } else { + (value, 0) + }; + let bits = scaled.to_bits(); + let exponent = ((bits >> 52) & 0x7ff) as i32 - 1023 + adjustment; + let mantissa = f64::from_bits((bits & 0x000f_ffff_ffff_ffff) | (1023_u64 << 52)); + (mantissa, exponent) +} + +fn scale_probability(value: f64, adjustment: i32) -> Option { + let (mantissa, exponent) = normalized_probability(value); + let exponent = exponent + adjustment; + if exponent > 1023 { + None + } else if exponent < -1022 { + Some(0.0) + } else { + Some(mantissa * normal_power_of_two(exponent)) + } +} + +pub(super) fn initial_candidate(b: f64, probability: f64) -> f64 { + if probability == 0.5 { + return MEDIAN_SCALE / b; + } + let log_tail = (-probability).ln_1p(); + let tail_scale = -log_tail; + let mut scaled = (2.0 * probability).sqrt().max(tail_scale + tail_scale.ln()); + for _ in 0..3 { + if scaled < 0.0001 { + break; + } + let log_survival = -scaled + scaled.ln_1p(); + let survival = log_survival.exp(); + let cdf = 1.0 - survival; + let derivative = scaled * survival / (1.0 + scaled); + let step = (cdf - probability) / derivative; + let derivative_ratio = 1.0 / scaled - 1.0; + let next = scaled - step / (1.0 - 0.5 * step * derivative_ratio); + if next == scaled || !next.is_finite() || next <= 0.0 { + break; + } + scaled = next; + } + scaled / b +} + +fn limiting_correction(scaled_x: (f64, f64)) -> (f64, f64) { + let mut term = (0.5, 0.0); + let mut tail = (0.0, 0.0); + for index in 1..=64 { + let index = f64::from(index); + let coefficient = dd_div_f64((-(index + 1.0), 0.0), index * (index + 2.0)); + term = dd_mul(dd_mul(term, scaled_x), coefficient); + tail = dd_add(tail, term); + if index >= 8.0 && term.0.abs() <= f64::EPSILON * f64::EPSILON * tail.0.abs() { + break; + } + } + tail +} + +fn scaled_midpoint(b: f64, lower_bits: u64) -> (f64, f64) { + let lower = f64::from_bits(lower_bits); + let upper = f64::from_bits(lower_bits + 1); + let scaled_lower = dd_mul((b, 0.0), (lower, 0.0)); + let scaled_step = dd_mul((b, 0.0), (upper - lower, 0.0)); + dd_add(scaled_lower, (0.5 * scaled_step.0, 0.5 * scaled_step.1)) +} + +fn midpoint_difference(b: f64, probability: f64, lower_bits: u64) -> f64 { + let scaled_x = scaled_midpoint(b, lower_bits); + let exponent = binary_exponent(scaled_x.0); + let scale = normal_power_of_two(-exponent); + let normalized_x = (scaled_x.0 * scale, scaled_x.1 * scale); + let normalized_square = dd_mul(normalized_x, normalized_x); + let Some(normalized_probability) = scale_probability(probability, -2 * exponent) else { + return f64::INFINITY; + }; + if normalized_probability == 0.0 { + return f64::NEG_INFINITY; + } + let high_square = dd_mul((normalized_x.0, 0.0), (normalized_x.0, 0.0)); + let mut numerator = dd_add( + (normalized_probability, 0.0), + (-0.5 * high_square.0, -0.5 * high_square.1), + ); + let cross = dd_mul((normalized_x.0, 0.0), (normalized_x.1, 0.0)); + numerator = dd_add(numerator, (-cross.0, -cross.1)); + let low_square = dd_mul((normalized_x.1, 0.0), (normalized_x.1, 0.0)); + numerator = dd_add(numerator, (-0.5 * low_square.0, -0.5 * low_square.1)); + let target = dd_div(numerator, normalized_square); + let limiting = limiting_correction(scaled_x); + let difference = dd_add(target, (-limiting.0, -limiting.1)); + difference.0 + difference.1 +} + +fn certified_median_midpoint_order(b: f64, lower_bits: u64) -> Option { + let scaled_x = scaled_midpoint(b, lower_bits); + let difference = dd_add(scaled_x, (-MEDIAN_SCALE, -MEDIAN_SCALE_LOW)); + let difference = difference.0 + difference.1; + if difference > CERTIFIED_MEDIAN_RADIUS { + Some(core::cmp::Ordering::Less) + } else if difference < -CERTIFIED_MEDIAN_RADIUS { + Some(core::cmp::Ordering::Greater) + } else { + None + } +} + +pub(super) fn certified_midpoint_order( + b: f64, + probability: f64, + lower_bits: u64, +) -> Option { + if !b.is_finite() + || b.is_sign_negative() + || b.to_bits() < MINIMUM_SHAPE_BITS + || !(0.0..1.0).contains(&probability) + || lower_bits >= ENDPOINT_LIMIT_BITS + { + return None; + } + if probability == 0.5 { + return certified_median_midpoint_order(b, lower_bits); + } + let difference = midpoint_difference(b, probability, lower_bits); + (difference.abs() > CERTIFIED_DD_RADIUS).then(|| difference.total_cmp(&0.0)) +} diff --git a/src/function/beta/inverse/shape_two/endpoint_certified.rs b/src/function/beta/inverse/shape_two/endpoint_certified.rs new file mode 100644 index 00000000..a4558c6f --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint_certified.rs @@ -0,0 +1,249 @@ +use super::endpoint_fixed::Fixed; +use core::cmp::Ordering; + +const ENDPOINT_LIMIT_BITS: u64 = 0x0020_0000_0000_0000; +const MINIMUM_SHAPE_BITS: u64 = (1023_u64 + 484) << 52; +const FRACTION_LIMBS: usize = 32; +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub(super) enum Certificate { + Ordered(Ordering), + Overlap, +} + +#[derive(Clone, Copy, Debug, Eq, PartialEq)] +pub(super) enum CertifierError { + InvalidInput, + ArithmeticInvariant, +} + +enum Stage { + Ordered(Ordering), + Overlap, + Refine, +} + +#[derive(Clone, Copy)] +struct Interval { + lower: Fixed, + upper: Fixed, + cutoff: usize, +} + +impl Interval { + fn exact(value: Fixed, cutoff: usize) -> Option { + Some(Self { + lower: value.quantize_floor(cutoff)?, + upper: value.quantize_ceil(cutoff)?, + cutoff, + }) + } + + fn zero(cutoff: usize) -> Self { + Self { + lower: Fixed::zero(), + upper: Fixed::zero(), + cutoff, + } + } + + fn same_precision(self, other: Self) -> Option { + (self.cutoff == other.cutoff).then_some(self.cutoff) + } + + fn checked_add(self, other: Self) -> Option { + let cutoff = self.same_precision(other)?; + Some(Self { + lower: self.lower.checked_add(other.lower, cutoff)?, + upper: self.upper.checked_add(other.upper, cutoff)?, + cutoff, + }) + } + + fn checked_sub(self, other: Self) -> Option { + let cutoff = self.same_precision(other)?; + Some(Self { + lower: self.lower.checked_sub(other.upper, cutoff)?, + upper: self.upper.checked_sub(other.lower, cutoff)?, + cutoff, + }) + } + + fn checked_mul(self, other: Self) -> Option { + let cutoff = self.same_precision(other)?; + Some(Self { + lower: self.lower.mul_floor(other.lower, cutoff)?, + upper: self.upper.mul_ceil(other.upper, cutoff)?, + cutoff, + }) + } + + fn checked_mul_small(self, factor: u64) -> Option { + Some(Self { + lower: self.lower.checked_mul_small(factor, self.cutoff)?, + upper: self.upper.checked_mul_small(factor, self.cutoff)?, + cutoff: self.cutoff, + }) + } + + fn checked_div_small(self, divisor: u64) -> Option { + Some(Self { + lower: self.lower.div_small_floor(divisor, self.cutoff), + upper: self.upper.div_small_ceil(divisor, self.cutoff)?, + cutoff: self.cutoff, + }) + } +} + +fn next_term(term: Interval, factor: Interval, index: u64) -> Option { + term.checked_mul(factor)? + .checked_mul_small(index + 1)? + .checked_div_small(index.checked_mul(index + 2)?) +} + +fn remainder_interval( + positive: Interval, + negative: Interval, + term: Interval, + index: u64, +) -> Option { + let cutoff = positive.same_precision(negative)?; + (cutoff == term.cutoff).then_some(())?; + let lower = positive.lower.checked_sub(negative.upper, cutoff)?; + let upper = positive.upper.checked_sub(negative.lower, cutoff)?; + if index & 1 == 0 { + Some(Interval { + lower, + upper: upper.checked_add(term.upper, cutoff)?, + cutoff, + }) + } else { + Some(Interval { + lower: lower.checked_sub(term.upper, cutoff)?, + upper, + cutoff, + }) + } +} + +fn series_interval(y: Interval, x: Interval) -> Option { + let cutoff = y.same_precision(x)?; + let half = Interval::exact(Fixed::half()?, cutoff)?; + let quantum = Fixed::quantum(cutoff)?; + let mut positive = half; + let mut negative = Interval::zero(cutoff); + let mut term = half; + for index in 1_u64..=513 { + let factor = y.checked_sub(x.checked_mul_small(index)?)?; + term = next_term(term, factor, index)?; + if index >= 8 && term.upper <= quantum { + return remainder_interval(positive, negative, term, index); + } + if index & 1 == 0 { + positive = positive.checked_add(term)?; + } else { + negative = negative.checked_add(term)?; + } + } + None +} + +fn order_at_precision(b: f64, probability: f64, lower_bits: u64, active_limbs: usize) -> Stage { + let Some(cutoff) = FRACTION_LIMBS.checked_sub(active_limbs) else { + return Stage::Refine; + }; + let Some(x) = Fixed::midpoint(lower_bits).and_then(|value| Interval::exact(value, cutoff)) + else { + return Stage::Refine; + }; + let Some(y) = + Fixed::scaled_midpoint(b, lower_bits).and_then(|value| Interval::exact(value, cutoff)) + else { + return Stage::Refine; + }; + let Some(eight) = Fixed::integer(8) else { + return Stage::Refine; + }; + if y.upper >= eight { + return Stage::Refine; + } + let Some(series) = series_interval(y, x) else { + return Stage::Refine; + }; + let Some(cdf) = y + .checked_add(x) + .and_then(|sum| y.checked_mul(sum)) + .and_then(|prefactor| prefactor.checked_mul(series)) + else { + return Stage::Refine; + }; + let Some(probability) = + Fixed::from_f64(probability).and_then(|value| Interval::exact(value, cutoff)) + else { + return Stage::Refine; + }; + if probability.upper < cdf.lower { + Stage::Ordered(Ordering::Less) + } else if probability.lower > cdf.upper { + Stage::Ordered(Ordering::Greater) + } else { + Stage::Overlap + } +} + +fn exponent(value: f64) -> Option { + let bits = value.to_bits(); + let encoded = ((bits >> 52) & 0x7ff) as i32; + if encoded != 0 { + return (encoded != 0x7ff).then_some(encoded - 1023); + } + let mantissa = bits & 0x000f_ffff_ffff_ffff; + if mantissa == 0 { + return None; + } + let leading = i32::try_from(63 - mantissa.leading_zeros()).ok()?; + Some(leading - 1074) +} + +fn initial_active_limbs(probability: f64) -> Option { + let bits = exponent(probability)?.checked_neg()?.checked_add(128)?; + let rounded = bits.checked_add(63)?.checked_div(64)?; + usize::try_from(rounded) + .ok() + .map(|limbs| limbs.clamp(4, FRACTION_LIMBS)) +} + +pub(super) fn midpoint_certificate( + b: f64, + probability: f64, + lower_bits: u64, +) -> Result { + if !b.is_finite() + || b.is_sign_negative() + || b.to_bits() < MINIMUM_SHAPE_BITS + || !(probability > 0.0 && probability < 1.0) + || lower_bits >= ENDPOINT_LIMIT_BITS + { + return Err(CertifierError::InvalidInput); + } + let initial = initial_active_limbs(probability).ok_or(CertifierError::InvalidInput)?; + let stages = [initial, 8, 16, 24]; + let mut previous = 0; + for active in stages { + let active = active.max(initial).min(FRACTION_LIMBS); + if active == previous { + continue; + } + if let Stage::Ordered(order) = order_at_precision(b, probability, lower_bits, active) { + return Ok(Certificate::Ordered(order)); + } + previous = active; + } + match order_at_precision(b, probability, lower_bits, FRACTION_LIMBS) { + Stage::Ordered(order) => Ok(Certificate::Ordered(order)), + Stage::Overlap => Ok(Certificate::Overlap), + Stage::Refine => Err(CertifierError::ArithmeticInvariant), + } +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/inverse/shape_two/endpoint_certified/tests.rs b/src/function/beta/inverse/shape_two/endpoint_certified/tests.rs new file mode 100644 index 00000000..5a5ca55a --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint_certified/tests.rs @@ -0,0 +1,60 @@ +use super::*; + +fn exact_integer(value: u64) -> Interval { + Interval::exact(Fixed::integer(value).unwrap(), 0).unwrap() +} + +#[test] +fn alternating_remainder_uses_the_correct_parity() { + let positive = exact_integer(10); + let negative = exact_integer(3); + let term = exact_integer(1); + let even = remainder_interval(positive, negative, term, 8).unwrap(); + let odd = remainder_interval(positive, negative, term, 9).unwrap(); + + assert!(even.lower == Fixed::integer(7).unwrap()); + assert!(even.upper == Fixed::integer(8).unwrap()); + assert!(odd.lower == Fixed::integer(6).unwrap()); + assert!(odd.upper == Fixed::integer(7).unwrap()); +} + +#[test] +fn full_precision_reaches_an_enclosure_at_domain_extremes() { + let shapes = [ + f64::from_bits(MINIMUM_SHAPE_BITS), + f64::from_bits(0x7fef_ffff_ffff_ffff), + ]; + let probabilities = [ + f64::from_bits(1), + 0.5, + f64::from_bits(1.0_f64.to_bits() - 1), + ]; + let lower_bits = [0, 1, ENDPOINT_LIMIT_BITS - 1]; + + for shape in shapes { + for probability in probabilities { + for lower in lower_bits { + assert!(!matches!( + order_at_precision(shape, probability, lower, FRACTION_LIMBS), + Stage::Refine + )); + } + } + } +} + +#[test] +fn invalid_inputs_are_rejected() { + for shape in [f64::NAN, f64::INFINITY, -1.0, 1.0] { + assert_eq!( + midpoint_certificate(shape, 0.5, 0), + Err(CertifierError::InvalidInput) + ); + } + for probability in [f64::NAN, 0.0, 1.0] { + assert_eq!( + midpoint_certificate(f64::from_bits(MINIMUM_SHAPE_BITS), probability, 0), + Err(CertifierError::InvalidInput) + ); + } +} diff --git a/src/function/beta/inverse/shape_two/endpoint_fixed.rs b/src/function/beta/inverse/shape_two/endpoint_fixed.rs new file mode 100644 index 00000000..646d3535 --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint_fixed.rs @@ -0,0 +1,246 @@ +use core::cmp::Ordering; + +const LIMBS: usize = 34; +const WIDE_LIMBS: usize = 68; +const FRACTION_BITS: i32 = 2048; +const FRACTION_LIMBS: usize = 32; + +#[derive(Clone, Copy, Eq, PartialEq)] +pub(super) struct Fixed([u64; LIMBS]); + +impl Ord for Fixed { + fn cmp(&self, other: &Self) -> Ordering { + for index in (0..LIMBS).rev() { + match self.0[index].cmp(&other.0[index]) { + Ordering::Equal => {} + order => return order, + } + } + Ordering::Equal + } +} + +impl PartialOrd for Fixed { + fn partial_cmp(&self, other: &Self) -> Option { + Some(self.cmp(other)) + } +} + +impl Fixed { + pub(super) const fn zero() -> Self { + Self([0; LIMBS]) + } + + fn from_shifted_u128(value: u128, shift: usize) -> Option { + let mut result = Self::zero(); + let word = shift / 64; + let offset = shift % 64; + for (part_index, part) in [value as u64, (value >> 64) as u64].into_iter().enumerate() { + if part == 0 { + continue; + } + let index = word + part_index; + if index >= LIMBS { + return None; + } + result.0[index] |= part << offset; + if offset != 0 { + if index + 1 >= LIMBS { + return None; + } + result.0[index + 1] |= part >> (64 - offset); + } + } + Some(result) + } + + pub(super) fn integer(value: u64) -> Option { + Self::from_shifted_u128(value as u128, FRACTION_BITS as usize) + } + + pub(super) fn half() -> Option { + Self::from_shifted_u128(1, FRACTION_BITS as usize - 1) + } + + pub(super) fn quantum(cutoff: usize) -> Option { + (cutoff <= FRACTION_LIMBS).then(|| { + let mut result = Self::zero(); + result.0[cutoff] = 1; + result + }) + } + + pub(super) fn quantize_floor(mut self, cutoff: usize) -> Option { + (cutoff <= FRACTION_LIMBS).then(|| { + self.0[..cutoff].fill(0); + self + }) + } + + pub(super) fn quantize_ceil(mut self, cutoff: usize) -> Option { + if cutoff > FRACTION_LIMBS { + return None; + } + let discarded = self.0[..cutoff].iter().any(|&limb| limb != 0); + self.0[..cutoff].fill(0); + if discarded { + self = self.checked_add(Self::quantum(cutoff)?, cutoff)?; + } + Some(self) + } + + pub(super) fn from_f64(value: f64) -> Option { + let bits = value.to_bits(); + let encoded_exponent = ((bits >> 52) & 0x7ff) as i32; + let (mantissa, power) = if encoded_exponent == 0 { + (bits & 0x000f_ffff_ffff_ffff, -1074) + } else { + ( + (bits & 0x000f_ffff_ffff_ffff) | (1_u64 << 52), + encoded_exponent - 1023 - 52, + ) + }; + let shift = FRACTION_BITS.checked_add(power)?; + Self::from_shifted_u128(mantissa as u128, usize::try_from(shift).ok()?) + } + + pub(super) fn midpoint(lower_bits: u64) -> Option { + let numerator = (2_u128).checked_mul(lower_bits as u128)?.checked_add(1)?; + Self::from_shifted_u128(numerator, (FRACTION_BITS - 1075) as usize) + } + + pub(super) fn scaled_midpoint(b: f64, lower_bits: u64) -> Option { + let bits = b.to_bits(); + let encoded_exponent = ((bits >> 52) & 0x7ff) as i32; + if encoded_exponent == 0 || encoded_exponent == 0x7ff { + return None; + } + let mantissa = (bits & 0x000f_ffff_ffff_ffff) | (1_u64 << 52); + let numerator = (2_u128) + .checked_mul(lower_bits as u128)? + .checked_add(1)? + .checked_mul(mantissa as u128)?; + let exponent = encoded_exponent - 1023; + let shift = FRACTION_BITS.checked_add(exponent)?.checked_sub(1127)?; + Self::from_shifted_u128(numerator, usize::try_from(shift).ok()?) + } + + pub(super) fn checked_add(self, other: Self, cutoff: usize) -> Option { + let mut result = Self::zero(); + let mut carry = 0_u128; + for index in cutoff..LIMBS { + let sum = self.0[index] as u128 + other.0[index] as u128 + carry; + result.0[index] = sum as u64; + carry = sum >> 64; + } + (carry == 0).then_some(result) + } + + pub(super) fn checked_sub(self, other: Self, cutoff: usize) -> Option { + if self < other { + return None; + } + let mut result = Self::zero(); + let mut borrow = 0_u128; + for index in cutoff..LIMBS { + let subtrahend = other.0[index] as u128 + borrow; + let value = self.0[index] as u128; + result.0[index] = value.wrapping_sub(subtrahend) as u64; + borrow = u128::from(value < subtrahend); + } + (borrow == 0).then_some(result) + } + + pub(super) fn checked_mul_small(self, factor: u64, cutoff: usize) -> Option { + let mut result = Self::zero(); + let mut carry = 0_u128; + for index in cutoff..LIMBS { + let product = self.0[index] as u128 * factor as u128 + carry; + result.0[index] = product as u64; + carry = product >> 64; + } + (carry == 0).then_some(result) + } + + fn div_small(self, divisor: u64, cutoff: usize) -> (Self, u64) { + let mut result = Self::zero(); + let mut remainder = 0_u128; + for index in (cutoff..LIMBS).rev() { + let numerator = (remainder << 64) | self.0[index] as u128; + result.0[index] = (numerator / divisor as u128) as u64; + remainder = numerator % divisor as u128; + } + (result, remainder as u64) + } + + pub(super) fn div_small_floor(self, divisor: u64, cutoff: usize) -> Self { + self.div_small(divisor, cutoff).0 + } + + pub(super) fn div_small_ceil(self, divisor: u64, cutoff: usize) -> Option { + let (mut result, remainder) = self.div_small(divisor, cutoff); + if remainder != 0 { + result = result.checked_add(Self::quantum(cutoff)?, cutoff)?; + } + Some(result) + } + + fn product(self, other: Self, cutoff: usize) -> Option<[u64; WIDE_LIMBS]> { + let mut wide = [0_u64; WIDE_LIMBS]; + for left in cutoff..LIMBS { + if self.0[left] == 0 { + continue; + } + let mut carry = 0_u128; + for right in cutoff..LIMBS { + let index = left + right; + let product = + self.0[left] as u128 * other.0[right] as u128 + wide[index] as u128 + carry; + wide[index] = product as u64; + carry = product >> 64; + } + let mut index = left + LIMBS; + while carry != 0 { + if index == WIDE_LIMBS { + return None; + } + let sum = wide[index] as u128 + carry; + wide[index] = sum as u64; + carry = sum >> 64; + index += 1; + } + } + Some(wide) + } + + fn scaled_product(self, other: Self, cutoff: usize) -> Option<(Self, bool)> { + let wide = self.product(other, cutoff)?; + if wide[FRACTION_LIMBS + LIMBS..].iter().any(|&limb| limb != 0) { + return None; + } + let mut result = Self::zero(); + result + .0 + .copy_from_slice(&wide[FRACTION_LIMBS..FRACTION_LIMBS + LIMBS]); + result.0[..cutoff].fill(0); + let discarded = wide[..FRACTION_LIMBS + cutoff] + .iter() + .any(|&limb| limb != 0); + Some((result, discarded)) + } + + pub(super) fn mul_floor(self, other: Self, cutoff: usize) -> Option { + Some(self.scaled_product(other, cutoff)?.0) + } + + pub(super) fn mul_ceil(self, other: Self, cutoff: usize) -> Option { + let (mut result, discarded) = self.scaled_product(other, cutoff)?; + if discarded { + result = result.checked_add(Self::quantum(cutoff)?, cutoff)?; + } + Some(result) + } +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/inverse/shape_two/endpoint_fixed/tests.rs b/src/function/beta/inverse/shape_two/endpoint_fixed/tests.rs new file mode 100644 index 00000000..b8d9343b --- /dev/null +++ b/src/function/beta/inverse/shape_two/endpoint_fixed/tests.rs @@ -0,0 +1,32 @@ +use super::Fixed; + +#[test] +fn multiplication_bounds_contain_the_full_precision_product() { + let cases = [(0.1, 0.3), (1.5, 2.25), (7.75, 7.875)]; + for cutoff in [0, 8, 16, 24, 31] { + for (left, right) in cases { + let left = Fixed::from_f64(left).unwrap(); + let right = Fixed::from_f64(right).unwrap(); + let exact = left.mul_floor(right, 0).unwrap(); + let lower = left.mul_floor(right, cutoff).unwrap(); + let upper = left.mul_ceil(right, cutoff).unwrap(); + assert!(lower <= exact && exact <= upper); + } + } +} + +#[test] +fn division_bounds_are_outward_and_detect_discarded_bits() { + let value = Fixed::from_f64(1.0).unwrap(); + for cutoff in [0, 8, 16, 24, 31] { + let fine_lower = value.div_small_floor(3, 0); + let fine_upper = value.div_small_ceil(3, 0).unwrap(); + let lower = value.div_small_floor(3, cutoff); + let upper = value.div_small_ceil(3, cutoff).unwrap(); + assert!(lower <= fine_lower && fine_upper <= upper); + assert!(lower < upper); + } + + let exact = Fixed::from_f64(1.5).unwrap(); + assert!(exact.div_small_floor(2, 0) == exact.div_small_ceil(2, 0).unwrap()); +} diff --git a/src/function/beta/inverse/shape_two/tests.rs b/src/function/beta/inverse/shape_two/tests.rs new file mode 100644 index 00000000..f8ac70ab --- /dev/null +++ b/src/function/beta/inverse/shape_two/tests.rs @@ -0,0 +1,259 @@ +use super::*; + +#[test] +fn shape_two_scalar_phase_preserves_accurate_rounding() { + let cases = [ + (200.0, 2.0, 1e-165, 0x3fc2aa4f7f316421_u64), + (79.0, 2.0, 8.048559608467247e-58, 0x3fc6fca7645f9501_u64), + (2.0, 732.0, 6.925147302269184e-72, 0x37fba96fc46c76d6_u64), + ]; + for (a, b, probability, expected) in cases { + assert_eq!( + inverse_beta_shape_two(a, b, probability).to_bits(), + expected + ); + } +} + +#[test] +fn shape_two_log_values_match_references() { + let cases: [(f64, f64, f64, f64); 3] = [ + (200.0_f64, 2.0, 0.48970503636005447, -138.15510557964274), + (200.0, 2.0, 0.14582246504394994, -379.9265403440175), + (2.0, 5.0, 0.18180347131894917, -1.203972804325936), + ]; + for (a, b, x, expected) in cases { + let actual = log_cdf_parts(a, b, x); + assert!( + (actual.0 + actual.1).to_bits().abs_diff(expected.to_bits()) <= 2, + "a={a} b={b} actual={:?} expected={expected:?}", + actual.0 + actual.1 + ); + } +} + +#[test] +fn shape_two_subnormal_probabilities_are_monotone() { + for (a, b) in [(2.0, 200.0), (200.0, 2.0)] { + let values = + [1_u64, 2, 3, 4].map(|bits| inverse_beta_shape_two(a, b, f64::from_bits(bits))); + assert!(values.windows(2).all(|pair| pair[0] <= pair[1])); + assert!(values[0] > 0.0); + } +} + +#[test] +fn shape_two_adjacent_selection_uses_one_error_scale() { + let probability = f64::from_bits(0x3fb7_8ac0_9e9f_630f); + let expected = 0x3f7f_81f8_1f81_f820_u64; + let actual = inverse_beta_shape_two(2.0, 65.0, probability).to_bits(); + assert!( + actual.abs_diff(expected) <= 1, + "actual={actual:#018x} expected={expected:#018x}" + ); +} + +#[test] +fn shape_two_large_shape_rounds_at_unit_boundary() { + let shape = 18_446_744_073_709_551_616.0; + assert_eq!(inverse_beta_shape_two(shape, 2.0, 0.5), 1.0); + assert_eq!( + inverse_beta_shape_two(2.0, shape, 0.5).to_bits(), + 0x3bfa_da82_5f97_62b2 + ); +} + +#[test] +fn shape_two_large_second_shape_does_not_overflow() { + assert_eq!( + crate::function::beta::inv_beta_reg(2.0, 1e308, 0.5).to_bits(), + 0x000c_1190_8513_0dd9 + ); +} + +#[test] +fn real_shape_two_inputs_match_500_digit_references() { + let cases = [ + (2.0, 5.5, 0.3, 0x3fc5_7cb4_a4fa_6544_u64), + (5.5, 2.0, 0.3, 0x3fe5_4170_6a98_75f7), + (2.0, 65.5, 0.3, 0x3f90_e2ac_2408_d260), + (65.5, 2.0, 0.3, 0x3fee_d6c0_af00_418b), + (2.0, 200.5, 0.3, 0x3f76_4d12_9f5f_2879), + (200.5, 2.0, 0.3, 0x3fef_9d2f_de46_a3e6), + (2.0, 200.5, 0.999_999_999, 0x3fbc_bed1_fc37_f6e0), + (200.5, 2.0, 0.999_999_999, 0x3fef_ffff_888c_10e7), + (2.0, 1_000_000.5, 1e-60, 0x387e_13b4_0c5d_6bc5), + (1_000_000.5, 2.0, 1e-60, 0x3fef_fed3_dd80_3e82), + (2.0, 1e308, 0.5, 0x000c_1190_8513_0dd9), + (1e308, 2.0, 0.5, 0x3ff0_0000_0000_0000), + ]; + for (a, b, probability, expected) in cases { + let actual = crate::function::beta::inv_beta_reg(a, b, probability).to_bits(); + assert!( + actual.abs_diff(expected) <= 1, + "a={a} b={b} probability={probability} actual={actual:#018x} expected={expected:#018x}" + ); + } +} + +#[test] +fn real_shape_two_subnormal_quantiles_match_500_digit_references() { + let cases = [ + (0x2007_9905_9deb_7818, 0x000f_77d8_b988_842e), + (0x2007_9905_9deb_7819, 0x000f_77d8_b988_842f), + (0x2007_9905_9deb_781a, 0x000f_77d8_b988_8431), + ]; + let values = cases.map(|(probability, expected)| { + let actual = crate::function::beta::inv_beta_reg(0.5, 2.0, f64::from_bits(probability)); + assert_eq!(actual.to_bits(), expected); + actual + }); + assert!(values.windows(2).all(|pair| pair[0] <= pair[1])); +} + +#[test] +fn real_shape_two_zero_cell_is_rounded_and_monotone() { + let cases = [(0x1668_7e92_154e_f7ac, 0_u64), (0x1e6c_cb05_3660_8d61, 1)]; + let values = cases.map(|(probability, expected)| { + let actual = crate::function::beta::inv_beta_reg(0.5, 2.0, f64::from_bits(probability)); + assert_eq!(actual.to_bits(), expected); + actual + }); + assert!(values[0] <= values[1]); + + let tiny_shape = f64::from_bits(1); + for probability in [0.1, 0.5, 0.9, 0.999_999_999] { + assert_eq!( + crate::function::beta::inv_beta_reg(tiny_shape, 2.0, probability), + 0.0 + ); + } +} + +#[test] +fn real_shape_two_zero_cell_matches_512_bit_reference() { + let cases = [ + (1_u64, 0_u64), + (0x0c8b_4ec7_f919_73fe, 0_u64), + (0x0c8b_4ec7_f919_73ff, 1_u64), + (0x0cbe_b8a0_f83c_a27e, 1_u64), + (0x0cbe_b8a0_f83c_a27f, 2_u64), + (0x0cd5_558c_3a9b_e29f, 2_u64), + (0x0cd5_558c_3a9b_e2a0, 3_u64), + ]; + let values = cases.map(|(probability, expected)| { + let actual = crate::function::beta::inv_beta_reg(2.0, 1e200, f64::from_bits(probability)); + assert_eq!(actual.to_bits(), expected); + actual + }); + assert!(values[0] <= values[1]); +} + +#[test] +fn real_shape_two_near_midpoint_matches_512_bit_reference() { + let shape = f64::from_bits(0x6570_0000_0000_0000); + let cases = [ + (0x0480_0000_0000_0000, 1_u64), + (0x0adf_ffff_ffff_fff8, 0x0003_ffff_ffff_ffff), + ]; + for (probability, expected) in cases { + let actual = + crate::function::beta::inv_beta_reg(2.0, shape, f64::from_bits(probability)).to_bits(); + assert_eq!(actual, expected, "probability={probability:#018x}"); + } +} + +#[test] +fn real_shape_two_first_normal_boundary_matches_500_bit_reference() { + let actual = crate::function::beta::inv_beta_reg( + 2.0, + f64::from_bits(0x6040_0000_0000_0000), + f64::from_bits(0x00bf_ffff_ffff_fffe), + ); + assert_eq!(actual.to_bits(), 0x000f_ffff_ffff_ffff); +} + +#[test] +fn real_shape_two_certified_midpoint_cases_match_1100_digit_references() { + let cases = [ + (0x61a0_0000_0000_0000, 2, 1), + (0x6570_0000_0000_0000, 0x04c9_0000_0000_0000, 3), + (0x6570_0000_0000_0000, 0x04d8_8000_0000_0000, 4), + (0x6180_0000_0000_0000, 1, 1), + (0x6180_0000_0000_0000, 2, 2), + (0x5e36_a09e_667f_3bcd, 1, 0x001f_ffff_ffff_ffff), + ( + 0x7fe0_0000_0000_0000, + 0x3fe3_0200_0530_5ea7, + 0x0010_0000_0000_0000, + ), + ( + 0x7fef_ffff_ffff_ffff, + 0x3fed_11ca_9b3a_ce79, + 0x000f_ffff_ffff_ffff, + ), + ]; + for (shape, probability, expected) in cases { + let actual = crate::function::beta::inv_beta_reg( + 2.0, + f64::from_bits(shape), + f64::from_bits(probability), + ); + assert_eq!( + actual.to_bits(), + expected, + "shape={shape:#018x} probability={probability:#018x}" + ); + } +} + +#[test] +fn real_shape_two_subnormal_cells_match_550_digit_references() { + let cases = [ + (0x3fdf_ffff_ffff_ffff, 0x1e6d_64d5_1e0d_b31c, 0x2_u64), + (0x3fdf_ffff_ffff_ffff, 0x1e87_9f3c_f0a9_fc5f, 0x10), + (0x3fdf_ffff_ffff_ffff, 0x1e91_1a46_0cc5_8ef7, 0x21), + (0x3fe0_0000_0000_0000, 0x1e72_f942_2c23_c47c, 0x2), + (0x3fe0_0000_0000_0000, 0x1e88_5f42_f152_1410, 0x10), + (0x3fe0_0000_0000_0000, 0x1e90_d663_a9ca_a39a, 0x20), + (0x3fe0_0000_0000_0001, 0x1e6d_64d5_1e0d_af1b, 0x2), + (0x3fe0_0000_0000_0001, 0x1e89_198c_8b83_0581, 0x11), + (0x3fe0_0000_0000_0001, 0x1e90_916b_2b5f_fe37, 0x1f), + ]; + for (a, probability, expected) in cases { + assert_eq!( + crate::function::beta::inv_beta_reg( + f64::from_bits(a), + 2.0, + f64::from_bits(probability), + ) + .to_bits(), + expected, + "a={a:#018x} probability={probability:#018x}" + ); + } +} + +#[test] +fn subnormal_first_shape_rounds_to_zero_above_specialization_gate() { + let probability = f64::from_bits(1.0_f64.to_bits() - 1); + for shape_bits in [2_u64, 3, 0x10, 0x100, 0x0010_0000_0000_0000] { + assert_eq!( + crate::function::beta::inv_beta_reg(f64::from_bits(shape_bits), 2.0, probability), + 0.0, + "shape_bits={shape_bits:#018x}" + ); + } +} + +#[test] +fn unit_first_shape_rounds_subnormal_quantiles_upward() { + for probability_bits in 1_u64..=1024 { + assert_eq!( + crate::function::beta::inv_beta_reg(1.0, 2.0, f64::from_bits(probability_bits),) + .to_bits(), + probability_bits.div_ceil(2), + "probability_bits={probability_bits:#018x}" + ); + } +} diff --git a/src/function/beta/inverse/shape_two/value.rs b/src/function/beta/inverse/shape_two/value.rs new file mode 100644 index 00000000..55f26e06 --- /dev/null +++ b/src/function/beta/inverse/shape_two/value.rs @@ -0,0 +1,246 @@ +use super::super::super::*; + +pub(super) fn log_cdf(a: f64, b: f64, x: f64) -> f64 { + if b == 2.0 { + a.mul_add(x.ln(), a.mul_add(1.0 - x, 1.0).ln()) + } else if x < 0.5 && b * x < 0.5 { + let mut term = 0.5; + let mut sum = term; + for k in 1..64 { + let k = f64::from(k); + term *= -(b - k) / k * x * (k + 1.0) / (k + 2.0); + sum += term; + if term.abs() <= f64::EPSILON * sum.abs() { + return b.ln() + (b + 1.0).ln() + 2.0 * x.ln() + sum.ln(); + } + } + panic!("shape-two beta series did not converge for b={b}, x={x}") + } else { + let log_tail = b.mul_add((-x).ln_1p(), b.mul_add(x, 1.0).ln()); + log1mexp(log_tail) + } +} + +pub(super) fn log_cdf_parts(a: f64, b: f64, x: f64) -> (f64, f64) { + if b == 2.0 { + let complement = two_sum(1.0, -x); + let factor = dd_add((1.0, 0.0), dd_mul((a, 0.0), complement)); + dd_add( + dd_mul((a, 0.0), shape_two_ln((x, 0.0))), + shape_two_ln(factor), + ) + } else if x < 0.5 && b * x < 0.5 { + let sum = series_sum_dd(b, x) + .unwrap_or_else(|| panic!("shape-two beta series did not converge for b={b}, x={x}")); + dd_add( + dd_add( + dd_add(accurate_ln_dd((b, 0.0)), accurate_ln_dd((b + 1.0, 0.0))), + dd_mul((2.0, 0.0), accurate_ln_dd((x, 0.0))), + ), + accurate_ln_dd(sum), + ) + } else { + accurate_ln_dd(tail_cdf_parts(b, x)) + } +} + +fn tail_cdf_parts(b: f64, x: f64) -> (f64, f64) { + tail_cdf_parts_dd(b, (x, 0.0)) +} + +fn tail_cdf_parts_dd(b: f64, x: (f64, f64)) -> (f64, f64) { + let log_tail = dd_add( + dd_mul((b, 0.0), shape_two_ln_one_plus((-x.0, -x.1))), + shape_two_ln_one_plus(dd_mul((b, 0.0), x)), + ); + let exponential = log_tail.0.exp(); + let exponential_error = exponential * log_tail.1.exp_m1(); + if log_tail.0 < -core::f64::consts::LN_2 { + dd_add((1.0, 0.0), (-exponential, -exponential_error)) + } else { + two_sum(-log_tail.0.exp_m1(), -exponential_error) + } +} + +fn shape_two_ln_one_plus(value: (f64, f64)) -> (f64, f64) { + if value.0.abs() > 0.5 { + return shape_two_ln(dd_add((1.0, 0.0), value)); + } + let ratio = dd_div(value, dd_add((2.0, 0.0), value)); + let ratio_squared = dd_mul(ratio, ratio); + let mut term = ratio; + let mut sum = ratio; + for index in 1..=24 { + term = dd_mul(term, ratio_squared); + sum = dd_add(sum, dd_div_f64(term, f64::from(2 * index + 1))); + if term.0.abs() <= f64::EPSILON * f64::EPSILON * sum.0.abs() { + break; + } + } + dd_mul((2.0, 0.0), sum) +} + +fn shape_two_ln(value: (f64, f64)) -> (f64, f64) { + let mut scaled = value; + let mut exponent_adjustment = 0_i32; + if scaled.0 < f64::MIN_POSITIVE { + scaled.0 *= 18_014_398_509_481_984.0; + scaled.1 *= 18_014_398_509_481_984.0; + exponent_adjustment = -54; + } + let mut exponent = ((scaled.0.to_bits() >> 52) & 0x7ff) as i32 - 1023; + let mut mantissa = + f64::from_bits((scaled.0.to_bits() & 0x000f_ffff_ffff_ffff) | (1023_u64 << 52)); + if mantissa > core::f64::consts::SQRT_2 { + mantissa *= 0.5; + exponent += 1; + } + let scale = 2.0_f64.powi(exponent); + let mantissa = (mantissa, scaled.1 / scale); + dd_add( + dd_mul( + (f64::from(exponent + exponent_adjustment), 0.0), + (core::f64::consts::LN_2, 2.3190468138462996e-17), + ), + shape_two_ln_one_plus(dd_add(mantissa, (-1.0, 0.0))), + ) +} + +fn series_sum_dd(b: f64, x: f64) -> Option<(f64, f64)> { + let mut term = (0.5, 0.0); + let mut sum = term; + for k in 1..64 { + let k = f64::from(k); + let coefficient = dd_div_f64(dd_mul(dd_add((b, 0.0), (-k, 0.0)), (x, 0.0)), k); + let coefficient = dd_mul(coefficient, (-(k + 1.0) / (k + 2.0), 0.0)); + term = dd_mul(term, coefficient); + sum = dd_add(sum, term); + if term.0.abs() <= f64::EPSILON * sum.0.abs() { + return Some(sum); + } + } + None +} + +fn integer_power(value: (f64, f64), exponent: f64) -> Option<(f64, f64)> { + if exponent != exponent.trunc() || !(1.0..18_446_744_073_709_551_616.0).contains(&exponent) { + return None; + } + let mut exponent = exponent as u64; + let mut factor = value; + let mut result = (1.0, 0.0); + while exponent != 0 { + if exponent & 1 != 0 { + result = dd_mul(result, factor); + } + exponent >>= 1; + if exponent != 0 { + factor = dd_mul(factor, factor); + } + } + (result.0 >= f64::MIN_POSITIVE && result.0.is_finite()).then_some(result) +} + +pub(super) fn direct_cdf_and_pdf(a: f64, b: f64, x: f64) -> Option<((f64, f64), f64)> { + if a == 2.0 && (x < 0.5 && b * x < 0.5 || (b == b.trunc() && b <= 64.0 && b * x <= 1.0)) { + let prefactor = dd_mul(dd_mul((b, 0.0), (x, 0.0)), dd_mul((b + 1.0, 0.0), (x, 0.0))); + let sum = series_sum_dd(b, x) + .unwrap_or_else(|| panic!("shape-two beta series did not converge for b={b}, x={x}")); + let cdf = dd_mul(prefactor, sum); + let pdf = ((b * x) * (1.0 - x).powf(b - 1.0)) * (b + 1.0); + if cdf.0 >= f64::MIN_POSITIVE && cdf.0 < 1.0 && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + if a == 2.0 && !(x < 0.5 && b * x < 0.5) { + let cdf = if let Some(power) = integer_power(two_sum(1.0, -x), b) { + let factor = dd_add((1.0, 0.0), dd_mul((b, 0.0), (x, 0.0))); + let tail = dd_mul(power, factor); + dd_add((1.0, 0.0), (-tail.0, -tail.1)) + } else { + tail_cdf_parts(b, x) + }; + let tail = 1.0 - (cdf.0 + cdf.1); + let pdf = tail * (b * x) * ((b + 1.0) / ((1.0 - x) * (1.0 + b * x))); + if cdf.0 >= f64::MIN_POSITIVE && cdf.0 < 1.0 && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + if b == 2.0 { + let power = integer_power((x, 0.0), a)?; + let complement = two_sum(1.0, -x); + let factor = dd_add((1.0, 0.0), dd_mul((a, 0.0), complement)); + let cdf = dd_mul(power, factor); + let pdf = dd_mul( + dd_mul((a * (a + 1.0), 0.0), power), + dd_div(complement, (x, 0.0)), + ); + let pdf = pdf.0 + pdf.1; + if cdf.0 >= f64::MIN_POSITIVE && cdf.0 < 1.0 && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + None +} + +fn integer_power_scalar(value: f64, exponent: f64) -> Option { + if exponent != exponent.trunc() || !(1.0..18_446_744_073_709_551_616.0).contains(&exponent) { + return None; + } + let mut exponent = exponent as u64; + let mut factor = value; + let mut result = 1.0; + while exponent != 0 { + if exponent & 1 != 0 { + result *= factor; + } + exponent >>= 1; + if exponent != 0 { + factor *= factor; + } + } + (result >= f64::MIN_POSITIVE && result.is_finite()).then_some(result) +} + +fn series_sum(b: f64, x: f64) -> Option { + let mut term = 0.5; + let mut sum = term; + for k in 1..64 { + let k = f64::from(k); + term *= -(b - k) / k * x * (k + 1.0) / (k + 2.0); + sum += term; + if term.abs() <= f64::EPSILON * sum.abs() { + return Some(sum); + } + } + None +} + +pub(super) fn fast_cdf_and_pdf(a: f64, b: f64, x: f64) -> Option<(f64, f64)> { + if a == 2.0 && x < 0.5 && b * x < 0.5 { + let cdf = (b * x) * ((b + 1.0) * x) * series_sum(b, x)?; + let pdf = ((b * x) * (1.0 - x).powf(b - 1.0)) * (b + 1.0); + if (f64::MIN_POSITIVE..1.0).contains(&cdf) && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + if a == 2.0 && x >= f64::MIN_POSITIVE && !(x < 0.5 && b * x < 0.5) && b <= 1.0 / f64::EPSILON { + let log_tail = b.mul_add((-x).ln_1p(), (b * x).ln_1p()); + let cdf = -log_tail.exp_m1(); + let tail = 1.0 - cdf; + let pdf = tail * (b * x) * ((b + 1.0) / ((1.0 - x) * (1.0 + b * x))); + if (f64::MIN_POSITIVE..1.0).contains(&cdf) && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + if b == 2.0 { + let power = integer_power_scalar(x, a).unwrap_or_else(|| x.powf(a)); + let complement = 1.0 - x; + let cdf = power * (1.0 + a * complement); + let pdf = (a * power) * ((a + 1.0) * complement / x); + if (f64::MIN_POSITIVE..1.0).contains(&cdf) && pdf > 0.0 && pdf.is_finite() { + return Some((cdf, pdf)); + } + } + None +} diff --git a/src/function/beta/inverse/solve.rs b/src/function/beta/inverse/solve.rs new file mode 100644 index 00000000..7a65c12e --- /dev/null +++ b/src/function/beta/inverse/solve.rs @@ -0,0 +1,182 @@ +use super::super::*; + +pub(super) fn inverse_beta_midpoint(lower: f64, upper: f64) -> f64 { + let arithmetic = lower + 0.5 * (upper - lower); + let candidate = if upper < 0.5 { + let positive_lower = if lower == 0.0 { + f64::from_bits(1) + } else { + lower + }; + (0.5 * (positive_lower.ln() + upper.ln())).exp() + } else if lower > 0.5 { + let lower_complement = 1.0 - lower; + let upper_complement = if upper == 1.0 { + f64::from_bits(1) + } else { + 1.0 - upper + }; + 1.0 - (0.5 * (lower_complement.ln() + upper_complement.ln())).exp() + } else { + arithmetic + }; + if candidate > lower && candidate < upper { + candidate + } else { + arithmetic + } +} + +pub(super) fn inverse_beta_adjacent_result( + lower: f64, + upper: f64, + lower_error: f64, + upper_error: f64, +) -> f64 { + if !lower_error.is_finite() { + return upper; + } + let fraction = -lower_error / (upper_error - lower_error); + if fraction < 0.5 { + lower + } else if fraction > 0.5 || upper.to_bits() & 1 == 0 { + upper + } else { + lower + } +} + +pub(super) fn inverse_beta_log_value_parts( + a: f64, + b: f64, + x: f64, + log_beta: (f64, f64), + accurate_log_beta: Option<(f64, f64)>, +) -> Result<(f64, f64), BetaFuncError> { + if (0.01..10.0).contains(&a) && b < 1.0 && 1.0 - x < 0.3 { + let accurate_log_beta = + accurate_log_beta.unwrap_or_else(|| ln_beta_inverse_accurate_parts(a, b)); + return beta_reg_small_b_shifted_log(a, b, x, 1.0 - x, accurate_log_beta) + .map(|value| (value, 0.0)); + } + if (10.0..1e15).contains(&a) + && b < 1.0 + && 1.0 - x < 0.3 + && let Some(value) = beta_reg_small_b_large_a_log(a, b, x, 1.0 - x)? + { + return Ok((value, 0.0)); + } + if use_beta_power_series(a, b, x) + && (!use_beta_symmetry(a, b, x) || use_beta_power_series_before_symmetry(a, b, x)) + { + beta_power_series_log_parts_with_log_beta(a, b, x, Some(log_beta)) + } else { + checked_ln_beta_reg_with_log_beta(a, b, x, Some(log_beta)).map(|value| (value, 0.0)) + } +} + +pub(super) fn inverse_beta_log_tail( + a: f64, + b: f64, + target: f64, + mut current: f64, + log_beta: (f64, f64), + ln_beta: f64, +) -> f64 { + const FAST_ITERATIONS: usize = 64; + const MAX_ITERATIONS: usize = 256; + + let (log_target, log_target_correction) = accurate_ln(target); + let mut lower = 0.0; + let mut upper = 1.0; + let mut lower_error = f64::NEG_INFINITY; + let mut upper_error = -log_target - log_target_correction; + let accurate_log_beta = if (0.01..10.0).contains(&a) && b < 1.0 { + Some(ln_beta_inverse_accurate_parts(a, b)) + } else { + None + }; + + for iteration in 0..MAX_ITERATIONS { + let log_value = inverse_beta_log_value_parts(a, b, current, log_beta, accurate_log_beta) + .unwrap_or_else(|error| { + panic!("inv_beta_reg evaluation failed at x={current:?}: {error}") + }); + let error_parts = dd_add(log_value, (-log_target, -log_target_correction)); + let error = error_parts.0 + error_parts.1; + if error_parts.0 == 0.0 && error_parts.1 == 0.0 { + return current; + } + + if error < 0.0 { + lower = current; + lower_error = error; + } else { + upper = current; + upper_error = error; + } + + let midpoint = inverse_beta_midpoint(lower, upper); + if midpoint == lower || midpoint == upper { + return inverse_beta_adjacent_result(lower, upper, lower_error, upper_error); + } + + let log_pdf = (a - 1.0) * current.ln() + (b - 1.0) * (-current).ln_1p() - ln_beta; + let step = error * (log_value.0 + log_value.1 - log_pdf).exp(); + let newton = current - step; + let next = if iteration < FAST_ITERATIONS + && newton.is_finite() + && ((newton > lower && newton < upper) || newton == current) + { + newton + } else { + midpoint + }; + + if next == current { + let neighbor = if error > 0.0 { + f64::from_bits(current.to_bits() - 1) + } else { + f64::from_bits(current.to_bits() + 1) + }; + let neighbor_value = + inverse_beta_log_value_parts(a, b, neighbor, log_beta, accurate_log_beta) + .unwrap_or_else(|evaluation_error| { + panic!("inv_beta_reg evaluation failed: {evaluation_error}") + }); + let neighbor_error = dd_add(neighbor_value, (-log_target, -log_target_correction)); + let neighbor_error = neighbor_error.0 + neighbor_error.1; + if error * neighbor_error <= 0.0 { + return if error > 0.0 { + inverse_beta_adjacent_result(neighbor, current, neighbor_error, error) + } else { + inverse_beta_adjacent_result(current, neighbor, error, neighbor_error) + }; + } + current = if neighbor_error.abs() <= error.abs() { + neighbor + } else { + midpoint + }; + } else { + current = next; + } + } + + panic!("inv_beta_reg did not converge for a={a}, b={b}, probability={target}") +} + +pub(super) fn inverse_beta_reflect(a: f64, b: f64, probability: f64, log_beta: (f64, f64)) -> bool { + if probability <= 0.5 { + false + } else if a >= b { + true + } else { + let midpoint_log_probability = checked_ln_beta_reg_with_log_beta(a, b, 0.5, Some(log_beta)) + .unwrap_or_else(|error| panic!("inv_beta_reg evaluation failed: {error}")); + midpoint_log_probability < probability.ln() + } +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/inverse/solve/tests.rs b/src/function/beta/inverse/solve/tests.rs new file mode 100644 index 00000000..795efadd --- /dev/null +++ b/src/function/beta/inverse/solve/tests.rs @@ -0,0 +1,20 @@ +use super::*; + +#[test] +fn log_value_computes_missing_accurate_beta_parts() { + let a = 2.5; + let b = 0.5; + let x = 0.8; + let log_beta = ln_beta_inverse_parts(a, b); + let expected = inverse_beta_log_value_parts( + a, + b, + x, + log_beta, + Some(ln_beta_inverse_accurate_parts(a, b)), + ) + .unwrap(); + let actual = inverse_beta_log_value_parts(a, b, x, log_beta, None).unwrap(); + + assert_eq!(actual, expected); +} diff --git a/src/function/beta/lanczos.rs b/src/function/beta/lanczos.rs new file mode 100644 index 00000000..5ad574b5 --- /dev/null +++ b/src/function/beta/lanczos.rs @@ -0,0 +1,76 @@ +use super::*; + +const SHIFT: f64 = 10.400511; +const BETA_SCALE: (f64, f64) = (3.06725258552748459, -1.3709967328337378e-16); +const LOW_TOTAL_MAX: f64 = 8.0; +const IMBALANCE_RATIO: f64 = 8.0; + +// Pugh, "An Analysis of the Lanczos Gamma Approximation", Table 8.5 and Eq. 6.14. +// These polynomials are the Horner form of statrs' existing partial-fraction sum. +fn sum(x: f64) -> f64 { + const NUMERATOR: [f64; 11] = [ + 2.48574089138753550e-5, + 2.59434050880906703e-3, + 1.21848070364446573e-1, + 3.39136624401530806, + 6.19452889142209600e1, + 7.75877940545563547e2, + 6.74876752593456695e3, + 4.02538353814263901e4, + 1.57567999493601179e5, + 3.65505352696257003e5, + 3.81540663397352677e5, + ]; + const DENOMINATOR: [f64; 10] = [ + 1.0, 45.0, 870.0, 9450.0, 63273.0, 269325.0, 723680.0, 1172700.0, 1026576.0, 362880.0, + ]; + let reciprocal = 1.0 / x; + let numerator = NUMERATOR[..10] + .iter() + .rev() + .fold(NUMERATOR[10], |value, coefficient| { + value.mul_add(reciprocal, *coefficient) + }); + let denominator = DENOMINATOR[..9] + .iter() + .rev() + .fold(DENOMINATOR[9], |value, coefficient| { + value.mul_add(reciprocal, *coefficient) + }); + numerator / denominator +} + +pub(super) fn use_beta_reg_lanczos_power(a: f64, b: f64) -> bool { + let smaller = a.min(b); + let larger = a.max(b); + let total = a + b; + smaller >= 2.0 + && total <= STIRLING_MIN + && (total <= LOW_TOTAL_MAX || larger >= IMBALANCE_RATIO * smaller) +} + +pub(super) fn beta_reg_lanczos_power(a: f64, b: f64, x: (f64, f64), y: (f64, f64)) -> (f64, f64) { + debug_assert!(use_beta_reg_lanczos_power(a, b)); + + let total = two_sum(a, b); + let shifted_a = dd_add((a, 0.0), (SHIFT, 0.0)); + let shifted_b = dd_add((b, 0.0), (SHIFT, 0.0)); + let shifted_total = dd_add(total, (SHIFT, 0.0)); + let delta_a = dd_div( + dd_add(dd_mul(shifted_total, x), (-shifted_a.0, -shifted_a.1)), + shifted_a, + ); + let delta_b = dd_div( + dd_add(dd_mul(shifted_total, y), (-shifted_b.0, -shifted_b.1)), + shifted_b, + ); + let log_a = (delta_a.0.ln_1p(), delta_a.1 / (1.0 + delta_a.0)); + let log_b = (delta_b.0.ln_1p(), delta_b.1 / (1.0 + delta_b.0)); + let exponent = dd_add(dd_mul((a, 0.0), log_a), dd_mul((b, 0.0), log_b)); + let shifted_scale = dd_div(dd_mul(shifted_a, shifted_b), shifted_total); + let lanczos_scale = sum(total.0 + total.1) / (sum(a) * sum(b)); + let sqrt_scale = shifted_scale.0.sqrt(); + let sqrt_scale = (sqrt_scale, shifted_scale.1 / (2.0 * sqrt_scale)); + let scale = dd_div(dd_mul((lanczos_scale, 0.0), sqrt_scale), BETA_SCALE); + dd_mul(scale, (dd_exp(exponent), 0.0)) +} diff --git a/src/function/beta/log_beta.rs b/src/function/beta/log_beta.rs new file mode 100644 index 00000000..ba6c2c89 --- /dev/null +++ b/src/function/beta/log_beta.rs @@ -0,0 +1,240 @@ +use super::*; + +/// Computes the natural logarithm of the beta function +/// where `a` is the first beta parameter +/// and `b` is the second beta parameter +/// and `a > 0`, `b > 0`. +/// +/// # Panics +/// +/// if `a <= 0.0` or `b <= 0.0` +pub fn ln_beta(a: f64, b: f64) -> f64 { + checked_ln_beta(a, b).unwrap() +} + +/// Computes the natural logarithm +/// of the beta function +/// where `a` is the first beta parameter +/// and `b` is the second beta parameter +/// and `a > 0`, `b > 0`. +/// +/// # Errors +/// +/// if `a <= 0.0` or `b <= 0.0` +pub fn checked_ln_beta(a: f64, b: f64) -> Result { + if a <= 0.0 { + Err(BetaFuncError::ANotGreaterThanZero) + } else if b <= 0.0 { + Err(BetaFuncError::BNotGreaterThanZero) + } else { + Ok(ln_beta_stable(a, b)) + } +} + +pub(super) fn stirling_correction(x: f64) -> f64 { + let reciprocal = 1.0 / x; + let x2 = reciprocal * reciprocal; + reciprocal + * (1.0 / 12.0 + + x2 * (-1.0 / 360.0 + + x2 * (1.0 / 1260.0 + + x2 * (-1.0 / 1680.0 + x2 * (1.0 / 1188.0 - x2 * 691.0 / 360360.0))))) +} + +pub(super) fn stirling_correction_log(log_x: f64) -> f64 { + let reciprocal = (-log_x).exp(); + let x2 = reciprocal * reciprocal; + reciprocal + * (1.0 / 12.0 + + x2 * (-1.0 / 360.0 + + x2 * (1.0 / 1260.0 + + x2 * (-1.0 / 1680.0 + x2 * (1.0 / 1188.0 - x2 * 691.0 / 360360.0))))) +} + +pub(super) fn ln_gamma_delta(base: f64, delta: f64) -> f64 { + let log_ratio = (delta / base).ln_1p(); + let log_sum = base.ln() + log_ratio; + delta * base.ln() + base.mul_add(log_ratio, (delta - 0.5) * log_ratio) - delta + + stirling_correction_log(log_sum) + - stirling_correction(base) +} + +pub(super) fn ln_gamma_stable(x: f64) -> f64 { + if x < 0.5 { + gamma::ln_gamma(1.0 + x) - x.ln() + } else { + gamma::ln_gamma(x) + } +} + +pub(super) fn ln_gamma_stirling_parts(value: (f64, f64)) -> (f64, f64) { + let shifted = dd_add(value, (-0.5, 0.0)); + let mut result = dd_mul(shifted, accurate_ln_dd(value)); + result = dd_add(result, (-value.0, -value.1)); + result = dd_add(result, (consts::LN_SQRT_2PI, -3.8782941580672414e-17)); + dd_add(result, (stirling_correction(value.0), 0.0)) +} + +pub(super) fn ln_gamma_accurate_parts(x: f64) -> (f64, f64) { + if x == 1.0 || x == 2.0 { + return (0.0, 0.0); + } + if x <= 0.125 { + let mut result = dd_add((x, 0.0), (1.0, 0.0)); + let mut recurrence = (0.0, 0.0); + while result.0 < STIRLING_MIN { + recurrence = dd_add(recurrence, accurate_ln_dd(result)); + result = dd_add(result, (1.0, 0.0)); + } + let gamma_one_plus = dd_add( + ln_gamma_stirling_parts(result), + (-recurrence.0, -recurrence.1), + ); + let logarithm = accurate_ln(x); + return dd_add(gamma_one_plus, (-logarithm.0, -logarithm.1)); + } + + let mut shifted = (x, 0.0); + let mut recurrence = (0.0, 0.0); + while shifted.0 < STIRLING_MIN { + recurrence = dd_add(recurrence, accurate_ln_dd(shifted)); + shifted = dd_add(shifted, (1.0, 0.0)); + } + let result = ln_gamma_stirling_parts(shifted); + dd_add(result, (-recurrence.0, -recurrence.1)) +} + +pub(super) fn ln_gamma_fast_accurate(x: f64) -> f64 { + if x <= 0.125 { + ln_gamma_one_plus_series(x) - x.ln() + } else { + ln_gamma_stable(x) + } +} + +pub(super) fn ln_gamma_delta_parts(base: f64, delta: f64) -> (f64, f64) { + let base_log = accurate_ln(base); + let ratio = dd_div_f64((delta, 0.0), base); + let log_ratio = accurate_ln_one_plus_dd(ratio); + let mut result = dd_mul((delta, 0.0), base_log); + result = dd_add(result, dd_mul((base, 0.0), log_ratio)); + result = dd_add(result, dd_mul((delta - 0.5, 0.0), log_ratio)); + result = dd_add(result, (-delta, 0.0)); + result = dd_add(result, (stirling_correction(base + delta), 0.0)); + dd_add(result, (-stirling_correction(base), 0.0)) +} + +pub(super) fn ln_gamma_delta_accurate_parts(base: f64, delta: f64) -> (f64, f64) { + if base >= STIRLING_MIN { + return ln_gamma_delta_parts(base, delta); + } + let steps = (STIRLING_MIN - base).ceil() as usize; + let mut result = ln_gamma_delta_parts(base + steps as f64, delta); + for step in 0..steps { + let ratio = dd_div_f64((delta, 0.0), base + step as f64); + let logarithm = accurate_ln_one_plus_dd(ratio); + result = dd_add(result, (-logarithm.0, -logarithm.1)); + } + result +} + +fn ln_beta_accurate_parts_impl(a: f64, b: f64, fast_small_gamma: bool) -> (f64, f64) { + let gamma_parts = |x| { + if fast_small_gamma && x <= 0.125 { + ln_gamma_small_accurate_parts(x) + } else { + ln_gamma_accurate_parts(x) + } + }; + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) { + let gamma = gamma_parts(smaller); + let delta = ln_gamma_delta_parts(larger, smaller); + return dd_add(gamma, (-delta.0, -delta.1)); + } + if a + b == f64::INFINITY { + return (ln_beta_stable(a, b), 0.0); + } + let gamma_a = gamma_parts(a); + let gamma_b = gamma_parts(b); + let gamma_sum = gamma_parts(a + b); + dd_add(dd_add(gamma_a, gamma_b), (-gamma_sum.0, -gamma_sum.1)) +} + +pub(super) fn ln_beta_accurate_parts(a: f64, b: f64) -> (f64, f64) { + ln_beta_accurate_parts_impl(a, b, false) +} + +pub(super) fn ln_beta_inverse_accurate_parts(a: f64, b: f64) -> (f64, f64) { + ln_beta_accurate_parts_impl(a, b, true) +} + +pub(super) fn ln_a_beta_small_first_shape_parts(a: f64, b: f64) -> (f64, f64) { + let gamma_one_plus = ln_gamma_one_plus_inverse_parts(a); + let gamma_delta = ln_gamma_delta_parts(b, a); + dd_add(gamma_one_plus, (-gamma_delta.0, -gamma_delta.1)) +} + +pub(super) fn ln_beta_stable_parts(a: f64, b: f64) -> (f64, f64) { + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) { + ln_beta_accurate_parts(a, b) + } else { + (ln_beta_stable(a, b), 0.0) + } +} + +pub(super) fn ln_beta_inverse_parts(a: f64, b: f64) -> (f64, f64) { + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) { + ln_beta_inverse_accurate_parts(a, b) + } else { + (ln_beta_stable(a, b), 0.0) + } +} + +pub(super) fn imbalanced_ln_beta(a: f64, b: f64) -> Option { + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && smaller < STIRLING_MIN { + Some(ln_gamma_stable(smaller) - ln_gamma_delta(larger, smaller)) + } else if smaller <= 1e-8 * larger { + Some(ln_gamma_stable(smaller) - smaller * gamma::digamma(larger)) + } else { + None + } +} + +pub(super) fn ln_beta_stable(a: f64, b: f64) -> f64 { + if a.min(b) <= 0.125 { + if a.max(b) >= STIRLING_MIN { + let result = ln_beta_accurate_parts(a, b); + return result.0 + result.1; + } + if a.max(b) <= 0.125 { + return (a + b).ln() - a.ln() - b.ln() + + ln_gamma_one_plus_series(a) + + ln_gamma_one_plus_series(b) + - ln_gamma_one_plus_series(a + b); + } + return ln_gamma_fast_accurate(a) + ln_gamma_fast_accurate(b) - ln_gamma_stable(a + b); + } + if let Some(ln_beta) = imbalanced_ln_beta(a, b) { + return ln_beta; + } + if a < STIRLING_MIN || b < STIRLING_MIN { + return ln_gamma_stable(a) + ln_gamma_stable(b) - ln_gamma_stable(a + b); + } + + let (mean, complement, log_sum, _) = beta_shape_statistics(a, b); + a * mean.ln() + + b * complement.ln() + + consts::LN_SQRT_2PI + + 0.5 * (log_sum - a.ln() - b.ln()) + + stirling_correction(a) + + stirling_correction(b) + - stirling_correction_log(log_sum) +} diff --git a/src/function/beta/log_forward.rs b/src/function/beta/log_forward.rs new file mode 100644 index 00000000..976bcae7 --- /dev/null +++ b/src/function/beta/log_forward.rs @@ -0,0 +1,176 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +pub(super) fn log1mexp(x: f64) -> f64 { + if x < -core::f64::consts::LN_2 { + (-x.exp()).ln_1p() + } else { + (-x.exp_m1()).ln() + } +} + +pub(crate) fn checked_ln_beta_reg_complement(a: f64, b: f64, x: f64) -> Result { + if a <= 0.0 { + return Err(BetaFuncError::ANotGreaterThanZero); + } + if b <= 0.0 { + return Err(BetaFuncError::BNotGreaterThanZero); + } + if !(0.0..=1.0).contains(&x) { + return Err(BetaFuncError::XOutOfRange); + } + if x == 1.0 { + return Ok(f64::NEG_INFINITY); + } + if x == 0.0 { + return Ok(0.0); + } + if a <= f64::EPSILON.sqrt() && b >= STIRLING_MIN && x.powf(a) > 0.5 { + let log_cdf = checked_ln_beta_reg(a, b, x)?; + return Ok(log1mexp(log_cdf)); + } + if use_beta_symmetry(a, b, x) { + let y = 1.0 - x; + if use_beta_power_series(b, a, y) { + return beta_power_series_log(b, a, y); + } + } + let log_cdf = checked_ln_beta_reg(a, b, x)?; + if log_cdf < -core::f64::consts::LN_2 { + Ok(log1mexp(log_cdf)) + } else { + checked_ln_beta_reg(b, a, 1.0 - x) + } +} + +pub(super) fn checked_ln_beta_reg_with_log_beta( + a: f64, + b: f64, + x: f64, + log_beta: Option<(f64, f64)>, +) -> Result { + if a <= 0.0 { + return Err(BetaFuncError::ANotGreaterThanZero); + } + if b <= 0.0 { + return Err(BetaFuncError::BNotGreaterThanZero); + } + if !(0.0..=1.0).contains(&x) { + return Err(BetaFuncError::XOutOfRange); + } + if x == 0.0 { + return Ok(f64::NEG_INFINITY); + } + if x == 1.0 { + return Ok(0.0); + } + if a == b && x == 0.5 { + return Ok(-core::f64::consts::LN_2); + } + if b == 1.0 { + return Ok(a * x.ln()); + } + if a == 1.0 { + return Ok((-(b * (-x).ln_1p()).exp_m1()).ln()); + } + let y = 1.0 - x; + if let Some((log_result, invert)) = + beta_small_shapes_series_log_with_log_beta(a, b, x, y, log_beta)? + { + return Ok(if invert { + log1mexp(log_result) + } else { + log_result + }); + } + if let Some(result) = beta_reg_asymptotic(a, b, x) { + return Ok(result.ln()); + } + if a.mul_add(y, -(b * x)) >= 0.0 + && let Some(result) = beta_reg_small_b_large_a_log(a, b, x, y)? + { + return Ok(result); + } + if use_beta_small_b_shifted_accurate(a, b, y) { + let result = beta_reg_small_b_shifted_accurate(a, b, x, y)?.1; + return Ok(result.0 + result.1); + } + if (1.0..10.0).contains(&a) && b < 1.0 && y < 0.3 { + return beta_reg_small_b_shifted_log(a, b, x, y, ln_beta_accurate_parts(a, b)); + } + let symm_transform = + !use_beta_power_series_before_symmetry(a, b, x) && use_beta_symmetry(a, b, x); + let (transformed_a, transformed_b, transformed_x, transformed_y) = if symm_transform { + (b, a, y, x) + } else { + (a, b, x, y) + }; + if !use_exact_complement_continued_fraction(a, b, symm_transform) + && let Some(log_tail) = beta_reg_small_b_large_a_log( + transformed_a, + transformed_b, + transformed_x, + transformed_y, + )? + { + return Ok(if symm_transform { + log1mexp(log_tail) + } else { + log_tail + }); + } + if use_beta_power_series(transformed_a, transformed_b, transformed_x) { + let log_result = beta_power_series_log_parts_with_log_beta( + transformed_a, + transformed_b, + transformed_x, + log_beta, + )?; + let log_result = log_result.0 + log_result.1; + return Ok(if symm_transform { + log1mexp(log_result) + } else { + log_result + }); + } + + let log_power = if let Some(log_beta) = log_beta { + beta_reg_log_power_parts_with_log_beta(a, b, x, log_beta) + } else { + beta_reg_log_power_parts(a, b, x) + }; + if symm_transform && (log_power.0 + log_power.1).exp() == 0.0 { + return Ok(0.0); + } + let fraction = beta_fraction_for_transformed_tail( + a, + b, + x, + transformed_a, + transformed_b, + transformed_x, + symm_transform, + )?; + let smaller = a.min(b); + let larger = a.max(b); + let log_fraction = if fraction.1 != 0.0 + || (larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger)) + { + accurate_ln_dd(fraction) + } else { + (fraction.0.ln(), 0.0) + }; + let log_result = dd_add(log_power, (-log_fraction.0, -log_fraction.1)); + let log_result = log_result.0 + log_result.1; + if symm_transform { + Ok(log1mexp(log_result)) + } else { + Ok(log_result) + } +} diff --git a/src/function/beta/mod.rs b/src/function/beta/mod.rs new file mode 100644 index 00000000..64382948 --- /dev/null +++ b/src/function/beta/mod.rs @@ -0,0 +1,106 @@ +//! Provides the [beta](https://en.wikipedia.org/wiki/Beta_function) and related +//! function +//! +//! This module sets the default precision more tightly than crate defaults for `DEFAULT_EPS` +//! +//! The implementation combines the incomplete-beta series, recurrences, and +//! continued fraction in [DLMF 8.17](https://dlmf.nist.gov/8.17), the BGRAT +//! expansion in [Algorithm 708](https://doi.org/10.1145/131766.131776), the +//! large-parameter expansion in [DLMF 8.18](https://dlmf.nist.gov/8.18), and +//! the small-argument log-gamma series in [DLMF 5.7.3](https://dlmf.nist.gov/5.7.E3). +//! Its direct beta prefactor uses [Pugh's Lanczos approximation, Table 8.5](https://web.viu.ca/pughg/phdThesis/phdThesis.pdf). + +mod api; +mod asymptotic; +mod bgrat; +mod dd; +mod forward; +mod fraction; +mod inverse; +mod lanczos; +mod log_beta; +mod log_forward; +mod normal_tail; +mod prefactor; +mod quantile; +mod recurrence; +mod scaled_gamma; +mod series; +mod small_gamma; +mod temme; + +pub use api::{beta, beta_inc, beta_reg, checked_beta, checked_beta_inc}; +pub use forward::checked_beta_reg; +pub use inverse::inv_beta_reg; +pub use log_beta::{checked_ln_beta, ln_beta}; +pub(crate) use log_forward::checked_ln_beta_reg_complement; + +use asymptotic::*; +use bgrat::*; +use dd::*; +use fraction::*; +use lanczos::*; +use log_beta::*; +use log_forward::*; +use normal_tail::*; +use prefactor::*; +use quantile::*; +use recurrence::*; +use scaled_gamma::*; +use series::*; +use small_gamma::*; +use temme::*; + +use crate::consts; +use crate::function::{erf, gamma}; +use crate::prec; +#[cfg(all(not(feature = "std"), not(test)))] +use num_traits::Float; + +/// sample case of module level precision +#[cfg(test)] +const MODULE_EPS: f64 = 1e-15; +const STIRLING_MIN: f64 = 32.0; +const SCALED_GAMMA_MIN_X: f64 = 64.0; +const MAX_BETA_REG_ITERATIONS: u32 = 100_000; +const ASYMPTOTIC_MIN_SUM: f64 = 1e7; +const ASYMPTOTIC_MIN_SHAPE: f64 = 1e6; +const ASYMPTOTIC_MAX_DEVIANCE: f64 = 9.0; + +/// Represents the errors that can occur when computing the natural logarithm +/// of the beta function or the regularized lower incomplete beta function. +#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)] +#[non_exhaustive] +pub enum BetaFuncError { + /// `a` is zero or less than zero. + ANotGreaterThanZero, + + /// `b` is zero or less than zero. + BNotGreaterThanZero, + + /// `x` is not in `[0, 1]`. + XOutOfRange, + + /// The numerical method did not converge. + ConvergenceFailed, +} + +impl core::fmt::Display for BetaFuncError { + fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result { + match self { + BetaFuncError::ANotGreaterThanZero => write!(f, "a is zero or less than zero"), + BetaFuncError::BNotGreaterThanZero => write!(f, "b is zero or less than zero"), + BetaFuncError::XOutOfRange => write!(f, "x is not in [0, 1]"), + BetaFuncError::ConvergenceFailed => write!(f, "computation did not converge"), + } + } +} + +impl core::error::Error for BetaFuncError {} + +pub(crate) fn checked_ln_beta_reg(a: f64, b: f64, x: f64) -> Result { + checked_ln_beta_reg_with_log_beta(a, b, x, None) +} + +#[cfg(test)] +mod tests; diff --git a/src/function/beta/normal_tail.rs b/src/function/beta/normal_tail.rs new file mode 100644 index 00000000..646432fd --- /dev/null +++ b/src/function/beta/normal_tail.rs @@ -0,0 +1,181 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) +// Adapted from the 53-bit erfc approximation in Boost.Math 1.90 erf.hpp. + +#[cfg(all(not(feature = "std"), not(test)))] +use super::Float; +use super::two_sum; + +fn evaluate(argument: f64, coefficients: &[f64]) -> f64 { + coefficients.iter().rev().fold(0.0, |value, coefficient| { + value.mul_add(argument, *coefficient) + }) +} + +pub(super) fn normal_tail(argument: (f64, f64)) -> f64 { + let (argument_value, argument_error) = two_sum(argument.0, argument.1); + let value = if argument_value < 0.5 { + let squared = argument_value * argument_value; + let numerator = evaluate( + squared, + &[ + 0.08343058921465318, + -0.33816513445936094, + -0.050999073514677746, + -0.007727583458021333, + -0.0003227801209646057, + ], + ); + let denominator = evaluate( + squared, + &[ + 1.0, + 0.455004033050794, + 0.08752226001422525, + 0.008585719250744063, + 0.000370900071787748, + ], + ); + let erf = argument_value * (1.0449485778808594 + numerator / denominator); + 0.5 - 0.5 * erf + } else { + normal_tail_rational(argument_value) + }; + let derivative = -(-argument_value * argument_value).exp() / core::f64::consts::PI.sqrt(); + let corrected = two_sum(value, derivative * argument_error); + corrected.0 + corrected.1 +} + +fn normal_tail_rational(argument: f64) -> f64 { + let (offset, numerator, denominator, constant) = if argument < 1.5 { + ( + argument - 0.5, + &[ + -0.09809059221628124, + 0.17811466584112034, + 0.19100369579677543, + 0.08889003689678845, + 0.01950490012512188, + 0.0018042453829701422, + ][..], + &[ + 1.0, + 1.8475907098300222, + 1.4262800484551132, + 0.5780528048899024, + 0.12385097467900864, + 0.011338523357700142, + 0.0000033751147248309468, + ][..], + 0.40593576431274414, + ) + } else if argument < 2.5 { + ( + argument - 1.5, + &[ + -0.024350047620769844, + 0.03865403750357072, + 0.04394818964209516, + 0.01756794363118021, + 0.0032396240629084213, + 0.00023583911559688072, + ][..], + &[ + 1.0, + 1.5399149494855245, + 0.9824037091579202, + 0.32573292478244445, + 0.056392183742047816, + 0.004103697239789046, + ][..], + 0.5067281723022461, + ) + } else { + ( + argument - 3.5, + &[ + 0.0029527671653097166, + 0.013738442589635533, + 0.008408076155555854, + 0.0021282562091461865, + 0.00025026996154479463, + 0.000011321240664884757, + ][..], + &[ + 1.0, + 1.0421781416693842, + 0.4425976594815631, + 0.09584927263010614, + 0.010598290648487653, + 0.0004794112695217145, + ][..], + 0.5405750274658203, + ) + }; + let rational = constant + evaluate(offset, numerator) / evaluate(offset, denominator); + let squared = argument * argument; + let squared_error = argument.mul_add(argument, -squared); + 0.5 * rational * (-squared).exp() * (-squared_error).exp() / argument +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn boundaries_match_500_digit_references() { + let cases = [ + (0x3fdfffffffffffff_u64, 0x3fceb02147ce245d_u64), + (0x3fe0000000000000, 0x3fceb02147ce245c), + (0x3fe0000000000001, 0x3fceb02147ce245a), + (0x3ff7ffffffffffff, 0x3f915aaa8ec85209), + (0x3ff8000000000000, 0x3f915aaa8ec85205), + (0x3ff8000000000001, 0x3f915aaa8ec85201), + (0x4003ffffffffffff, 0x3f2aab859b20acb0), + (0x4004000000000000, 0x3f2aab859b20ac9e), + (0x4004000000000001, 0x3f2aab859b20ac8d), + ]; + for (argument, expected) in cases { + let actual = normal_tail((f64::from_bits(argument), 0.0)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "argument={argument:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + } + + #[test] + fn double_double_arguments_match_500_digit_references() { + let cases = [ + ( + 0x3fd0000000000000_u64, + 2.0_f64.powi(-55), + 0x3fd728558ee694fb_u64, + ), + (0x3fd0000000000000, -2.0_f64.powi(-55), 0x3fd728558ee694fc), + (0x3ff8000000000000, 2.0_f64.powi(-53), 0x3f915aaa8ec85203), + ( + 0x400408614bd1f138, + f64::from_bits(0xbcb38087e4245eb0), + 0x3f2a178104a215c2, + ), + ( + 0x4006a2ae09a6ce40, + f64::from_bits(0xbc8cd2b297d889bc), + 0x3f0081598e5e54ed, + ), + (0x4006af973ca608ed, 0.0, 0x3effc9d9ca2b4546), + ]; + for (argument, error, expected) in cases { + let actual = normal_tail((f64::from_bits(argument), error)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "argument={argument:#018x}, error={error:?}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + } +} diff --git a/src/function/beta/prefactor.rs b/src/function/beta/prefactor.rs new file mode 100644 index 00000000..3d2ea07e --- /dev/null +++ b/src/function/beta/prefactor.rs @@ -0,0 +1,107 @@ +use super::*; + +pub(super) fn beta_reg_central_log_power_parts(a: f64, b: f64, x: f64) -> Option<(f64, f64)> { + if a >= STIRLING_MIN && b >= STIRLING_MIN && 1.0 - x < 1.0 { + let (residual, log_ratio) = beta_log_ratio(a, b, x); + if residual.abs() <= 0.1 * a.min(b) { + let (_, _, log_sum, _) = beta_shape_statistics(a, b); + let log_scale = consts::LN_SQRT_2PI + + 0.5 * (log_sum - a.ln() - b.ln()) + + stirling_correction(a) + + stirling_correction(b) + - stirling_correction_log(log_sum); + return Some(two_sum(log_ratio, -log_scale)); + } + } + None +} + +pub(super) fn beta_reg_log_power_parts_with_log_x( + a: f64, + b: f64, + (log_x, log_x_error): (f64, f64), + (log_y, log_y_error): (f64, f64), + (log_beta, log_beta_error): (f64, f64), +) -> (f64, f64) { + let a_log_x = a * log_x; + let a_log_x_error = a.mul_add(log_x, -a_log_x) + a * log_x_error; + let b_log_y = b * log_y; + let b_log_y_error = b.mul_add(log_y, -b_log_y) + b * log_y_error; + let (variable, variable_error) = two_sum(a_log_x, b_log_y); + let variable_error = variable_error + a_log_x_error + b_log_y_error; + let (result, result_error) = two_sum(variable, -log_beta); + (result, result_error + variable_error - log_beta_error) +} + +pub(super) fn beta_reg_log_power_parts(a: f64, b: f64, x: f64) -> (f64, f64) { + beta_reg_central_log_power_parts(a, b, x).unwrap_or_else(|| { + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) { + return beta_reg_log_power_parts_with_log_x( + a, + b, + accurate_ln(x), + accurate_ln_one_minus(x), + ln_beta_accurate_parts(a, b), + ); + } + beta_reg_log_power_parts_with_log_x( + a, + b, + compensated_ln(x), + compensated_ln_one_minus(x), + ln_beta_stable_parts(a, b), + ) + }) +} + +pub(super) fn beta_reg_log_power_parts_with_log_beta( + a: f64, + b: f64, + x: f64, + log_beta: (f64, f64), +) -> (f64, f64) { + beta_reg_central_log_power_parts(a, b, x).unwrap_or_else(|| { + let smaller = a.min(b); + let larger = a.max(b); + if larger >= STIRLING_MIN && (smaller < STIRLING_MIN || smaller <= 0.25 * larger) { + beta_reg_log_power_parts_with_log_x( + a, + b, + accurate_ln(x), + accurate_ln_one_minus(x), + log_beta, + ) + } else { + beta_reg_log_power_parts_with_log_x( + a, + b, + compensated_ln(x), + compensated_ln_one_minus(x), + log_beta, + ) + } + }) +} + +pub(super) fn beta_reg_log_power_parts_accurate(a: f64, b: f64, x: f64) -> (f64, f64) { + beta_reg_log_power_parts_accurate_with_log_beta(a, b, x, ln_beta_accurate_parts(a, b)) +} + +pub(super) fn beta_reg_log_power_parts_accurate_with_log_beta( + a: f64, + b: f64, + x: f64, + log_beta: (f64, f64), +) -> (f64, f64) { + beta_reg_central_log_power_parts(a, b, x).unwrap_or_else(|| { + beta_reg_log_power_parts_with_log_x( + a, + b, + accurate_ln(x), + accurate_ln_one_minus(x), + log_beta, + ) + }) +} diff --git a/src/function/beta/quantile.rs b/src/function/beta/quantile.rs new file mode 100644 index 00000000..417b30f1 --- /dev/null +++ b/src/function/beta/quantile.rs @@ -0,0 +1,60 @@ +use super::*; + +mod bounds; + +use bounds::{central_cell_is_certified, extreme_ratio_cell_is_certified}; + +fn beta_mean(a: f64, b: f64, mean: f64) -> (f64, f64) { + let scale = a.max(b); + let scaled_a = a / scale; + let scaled_b = b / scale; + let scaled_sum = scaled_a + scaled_b; + let scaled_a_error = (-scaled_a).mul_add(scale, a) / scale; + let scaled_b_error = (-scaled_b).mul_add(scale, b) / scale; + let virtual_scaled_b = scaled_sum - scaled_a; + let scaled_sum_error = (scaled_a - (scaled_sum - virtual_scaled_b)) + + (scaled_b - virtual_scaled_b) + + scaled_a_error + + scaled_b_error; + let product = mean * scaled_sum; + let product_error = mean.mul_add(scaled_sum, -product); + let difference = scaled_a - product; + let virtual_product = difference - scaled_a; + let difference_error = + (scaled_a - (difference - virtual_product)) + (-product - virtual_product); + let residual = + difference + (difference_error + scaled_a_error - product_error - mean * scaled_sum_error); + (mean, residual / scaled_sum) +} + +pub(super) fn beta_concentrated_quantile(a: f64, b: f64, probability: f64) -> Option { + let (mean, complement, _, root_sum) = beta_shape_statistics(a, b); + let spacing = if mean == 0.0 { + f64::from_bits(1) + } else if mean == 1.0 { + 1.0 - f64::from_bits(1.0_f64.to_bits() - 1) + } else { + let lower_spacing = mean - f64::from_bits(mean.to_bits() - 1); + let upper_spacing = f64::from_bits(mean.to_bits() + 1) - mean; + lower_spacing.min(upper_spacing) + }; + let standard_deviation = (mean * complement).sqrt() / root_sum; + if 64.0 * standard_deviation >= 0.5 * spacing { + return None; + } + + let mean = beta_mean(a, b, mean); + let normal_quantile = -core::f64::consts::SQRT_2 * erf::erfc_inv(2.0 * probability); + let reciprocal_sum = (1.0 / root_sum) / root_sum; + let skew_correction = + (complement - mean.0) * normal_quantile.mul_add(normal_quantile, -1.0) * reciprocal_sum + / 3.0; + let offset = normal_quantile.mul_add(standard_deviation, mean.1 + skew_correction); + let candidate = (mean.0 + offset).clamp(0.0, 1.0); + + if a.min(b) >= ASYMPTOTIC_MIN_SHAPE && mean.0.min(complement) >= 0.1 { + central_cell_is_certified(a, b, probability, mean, candidate).then_some(candidate) + } else { + extreme_ratio_cell_is_certified(a, b, probability, mean, candidate).then_some(candidate) + } +} diff --git a/src/function/beta/quantile/bounds.rs b/src/function/beta/quantile/bounds.rs new file mode 100644 index 00000000..ddec464a --- /dev/null +++ b/src/function/beta/quantile/bounds.rs @@ -0,0 +1,217 @@ +use super::super::*; + +const LOG_EXTREME_MARGIN: f64 = 16.0 * core::f64::consts::LN_2; +const LOG_CENTRAL_MARGIN: f64 = 4.0 * core::f64::consts::LN_2; + +type Double = (f64, f64); + +struct ScaledParts { + scale: f64, + a: Double, + b: Double, + sum: Double, +} + +fn scaled_parts(a: f64, b: f64) -> ScaledParts { + let scale = a.max(b); + let scaled_a = dd_div_f64((a, 0.0), scale); + let scaled_b = dd_div_f64((b, 0.0), scale); + ScaledParts { + scale, + a: scaled_a, + b: scaled_b, + sum: dd_add(scaled_a, scaled_b), + } +} + +fn log_sum(a: f64, b: f64) -> (f64, f64) { + let parts = scaled_parts(a, b); + dd_add(accurate_ln(parts.scale), accurate_ln_dd(parts.sum)) +} + +fn log_variance(a: f64, b: f64) -> (f64, f64) { + let sum = log_sum(a, b); + let mut result = dd_add(accurate_ln(a), accurate_ln(b)); + result = dd_add(result, dd_mul((-3.0, 0.0), sum)); + let reciprocal_sum = dd_exp((-sum.0, -sum.1)); + let correction = accurate_ln_one_plus_dd((reciprocal_sum, 0.0)); + dd_add(result, (-correction.0, -correction.1)) +} + +fn log_cantelli_bound(variance: (f64, f64), distance: (f64, f64)) -> f64 { + if distance.0 + distance.1 <= 0.0 { + return 0.0; + } + let ratio = dd_add( + dd_mul((2.0, 0.0), accurate_ln_dd(distance)), + (-variance.0, -variance.1), + ); + let ratio = ratio.0 + ratio.1; + if ratio > 0.0 { + -ratio - (-ratio).exp().ln_1p() + } else { + -ratio.exp().ln_1p() + } +} + +fn relative_logs(a: f64, b: f64, point: (f64, f64)) -> ((f64, f64), (f64, f64), f64) { + let parts = scaled_parts(a, b); + let residual = dd_add(dd_mul(point, parts.sum), (-parts.a.0, -parts.a.1)); + let relative_a = dd_div(residual, parts.a); + let relative_b = dd_div((-residual.0, -residual.1), parts.b); + let log_a = accurate_ln_one_plus_dd(relative_a); + let log_b = accurate_ln_one_plus_dd(relative_b); + let centered_a = dd_add(log_a, (-relative_a.0, -relative_a.1)); + let centered_b = dd_add(log_b, (-relative_b.0, -relative_b.1)); + let exponent = dd_add(dd_mul(parts.a, centered_a), dd_mul(parts.b, centered_b)); + let exponent = exponent.0 + exponent.1; + let exponent = if exponent < -f64::MAX / parts.scale { + f64::NEG_INFINITY + } else { + parts.scale * exponent + }; + (log_a, log_b, exponent) +} + +// Beta = Ga / (Ga + Gb); Markov applied to exp(-t((1-x)Ga-xGb)) gives the KL/Chernoff tail bound. +fn log_chernoff_bound(a: f64, b: f64, point: (f64, f64)) -> f64 { + relative_logs(a, b, point).2 +} + +fn upper_bound_is_below(bound: f64, target: (f64, f64), margin: f64) -> bool { + let proof = dd_add(target, (-bound - margin, 0.0)); + proof.0 + proof.1 > 0.0 +} + +fn log_density(a: f64, b: f64, point: (f64, f64)) -> f64 { + let sum = log_sum(a, b); + let mut center = dd_mul((1.5, 0.0), sum); + center = dd_add(center, dd_mul((-0.5, 0.0), accurate_ln(a))); + center = dd_add(center, dd_mul((-0.5, 0.0), accurate_ln(b))); + center = dd_add(center, (-consts::LN_SQRT_2PI, 3.8782941580672414e-17)); + center = dd_add( + center, + ( + -stirling_correction(a) - stirling_correction(b) + + stirling_correction_log(sum.0 + sum.1), + 0.0, + ), + ); + let (log_a, log_b, exponent) = relative_logs(a, b, point); + let result = dd_add(center, (exponent, 0.0)); + let result = dd_add(result, (-log_a.0, -log_a.1)); + let result = dd_add(result, (-log_b.0, -log_b.1)); + result.0 + result.1 +} + +fn mode(a: f64, b: f64) -> (f64, f64) { + let parts = scaled_parts(a, b); + let reciprocal_scale = 1.0 / parts.scale; + dd_div( + dd_add(parts.a, (-reciprocal_scale, 0.0)), + dd_add(parts.sum, (-2.0 * reciprocal_scale, 0.0)), + ) +} + +fn lower_bound_is_above( + a: f64, + b: f64, + point: (f64, f64), + target: (f64, f64), + variance: (f64, f64), +) -> bool { + let mode_distance = dd_add(mode(a, b), (-point.0, -point.1)); + if mode_distance.0 + mode_distance.1 < 0.0 { + return false; + } + let width = dd_exp(dd_mul((0.5, 0.0), variance)).min(0.5 * (point.0 + point.1)); + if width <= 0.0 { + return false; + } + let left = dd_add(point, (-width, 0.0)); + let bound = width.ln() + log_density(a, b, left) - LOG_CENTRAL_MARGIN; + let proof = dd_add((bound, 0.0), (-target.0, -target.1)); + proof.0 + proof.1 > 0.0 +} + +pub(super) fn extreme_ratio_cell_is_certified( + a: f64, + b: f64, + probability: f64, + mean: (f64, f64), + candidate: f64, +) -> bool { + let variance = log_variance(a, b); + if candidate > 0.0 { + let previous = f64::from_bits(candidate.to_bits() - 1); + let midpoint = dd_mul(dd_add((candidate, 0.0), (previous, 0.0)), (0.5, 0.0)); + let distance = dd_add(mean, (-midpoint.0, -midpoint.1)); + let bound = log_cantelli_bound(variance, distance); + if !upper_bound_is_below(bound, accurate_ln(probability), LOG_EXTREME_MARGIN) { + return false; + } + } + if candidate < 1.0 { + let next = f64::from_bits(candidate.to_bits() + 1); + let midpoint = dd_mul(dd_add((candidate, 0.0), (next, 0.0)), (0.5, 0.0)); + let distance = dd_add(midpoint, (-mean.0, -mean.1)); + let bound = log_cantelli_bound(variance, distance); + if !upper_bound_is_below( + bound, + accurate_ln_one_minus(probability), + LOG_EXTREME_MARGIN, + ) { + return false; + } + } + true +} + +pub(super) fn central_cell_is_certified( + a: f64, + b: f64, + probability: f64, + mean: (f64, f64), + candidate: f64, +) -> bool { + let variance = log_variance(a, b); + let log_probability = accurate_ln(probability); + let log_complement = accurate_ln_one_minus(probability); + if candidate > 0.0 { + let previous = f64::from_bits(candidate.to_bits() - 1); + let midpoint = dd_mul(dd_add((candidate, 0.0), (previous, 0.0)), (0.5, 0.0)); + let distance = dd_add(mean, (-midpoint.0, -midpoint.1)); + let proven = if distance.0 + distance.1 > 0.0 { + upper_bound_is_below( + log_chernoff_bound(a, b, midpoint), + log_probability, + LOG_CENTRAL_MARGIN, + ) + } else { + let reflected = dd_add((1.0, 0.0), (-midpoint.0, -midpoint.1)); + lower_bound_is_above(b, a, reflected, log_complement, variance) + }; + if !proven { + return false; + } + } + if candidate < 1.0 { + let next = f64::from_bits(candidate.to_bits() + 1); + let midpoint = dd_mul(dd_add((candidate, 0.0), (next, 0.0)), (0.5, 0.0)); + let distance = dd_add(midpoint, (-mean.0, -mean.1)); + let proven = if distance.0 + distance.1 > 0.0 { + let reflected = dd_add((1.0, 0.0), (-midpoint.0, -midpoint.1)); + upper_bound_is_below( + log_chernoff_bound(b, a, reflected), + log_complement, + LOG_CENTRAL_MARGIN, + ) + } else { + lower_bound_is_above(a, b, midpoint, log_probability, variance) + }; + if !proven { + return false; + } + } + true +} diff --git a/src/function/beta/recurrence.rs b/src/function/beta/recurrence.rs new file mode 100644 index 00000000..f2d29739 --- /dev/null +++ b/src/function/beta/recurrence.rs @@ -0,0 +1,46 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +pub(super) fn beta_a_step(a: f64, b: f64, x: f64, steps: usize) -> f64 { + let power = beta_reg_log_power_parts(a, b, x); + (power.0 + power.1 + beta_a_step_log_sum(a, b, x, steps) - a.ln()).exp() +} + +pub(super) fn beta_a_step_log_sum(a: f64, b: f64, x: f64, steps: usize) -> f64 { + let mut log_sum = 0.0_f64; + let mut log_term = 0.0_f64; + let log_x = x.ln(); + for i in 0..steps.saturating_sub(1) { + let i = i as f64; + log_term += (a + b + i).ln() + log_x - (a + i + 1.0).ln(); + let maximum = log_sum.max(log_term); + log_sum = maximum + (log_sum.min(log_term) - maximum).exp().ln_1p(); + } + log_sum +} + +pub(super) fn beta_a_step_log_parts( + a: f64, + b: f64, + x: f64, + steps: usize, + log_beta: (f64, f64), +) -> (f64, f64) { + let power = beta_reg_log_power_parts_accurate_with_log_beta(a, b, x, log_beta); + let log_a = accurate_ln(a); + dd_add( + dd_add(power, (beta_a_step_log_sum(a, b, x, steps), 0.0)), + (-log_a.0, -log_a.1), + ) +} + +pub(super) fn beta_a_step_log(a: f64, b: f64, x: f64, steps: usize, log_beta: (f64, f64)) -> f64 { + let result = beta_a_step_log_parts(a, b, x, steps, log_beta); + result.0 + result.1 +} diff --git a/src/function/beta/scaled_gamma.rs b/src/function/beta/scaled_gamma.rs new file mode 100644 index 00000000..1375026f --- /dev/null +++ b/src/function/beta/scaled_gamma.rs @@ -0,0 +1,118 @@ +//! Scaled upper-gamma evaluations follow DLMF 8.9.2 and 8.11.2. The bounded +//! continued-fraction recurrence is adapted from the existing `statrs` gamma +//! implementation; the small-shape form is the cancellation-safe DLMF 8.7.3 series. + +use super::*; + +pub(super) fn upper_gamma_scaled_asymptotic(shape: f64, x: f64) -> Result { + let mut term = 1.0_f64; + let mut sum = 1.0_f64; + for n in 1..=64 { + term *= (shape - f64::from(n)) / x; + sum += term; + if term.abs() <= prec::F64_PREC * sum.abs() { + return Ok(sum / x); + } + } + Err(BetaFuncError::ConvergenceFailed) +} + +pub(super) fn upper_gamma_scaled_continued_fraction( + shape: f64, + x: f64, +) -> Result { + const BIG: f64 = 4_503_599_627_370_496.0; + const BIG_INVERSE: f64 = 2.220446049250313e-16; + + let mut y = 1.0 - shape; + let mut z = x + y + 1.0; + let mut c = 0.0; + let mut pkm2 = 1.0; + let mut qkm2 = x; + let mut pkm1 = x + 1.0; + let mut qkm1 = z * x; + let mut result = pkm1 / qkm1; + for _ in 0..256 { + y += 1.0; + z += 2.0; + c += 1.0; + let yc = y * c; + let pk = pkm1 * z - pkm2 * yc; + let qk = qkm1 * z - qkm2 * yc; + + pkm2 = pkm1; + pkm1 = pk; + qkm2 = qkm1; + qkm1 = qk; + + if pk.abs() > BIG { + pkm2 *= BIG_INVERSE; + pkm1 *= BIG_INVERSE; + qkm2 *= BIG_INVERSE; + qkm1 *= BIG_INVERSE; + } + + if qk != 0.0 { + let next = pk / qk; + let relative_change = ((result - next) / next).abs(); + result = next; + if relative_change <= 4.0 * prec::F64_PREC { + return if result > 0.0 && result.is_finite() { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + }; + } + } + } + Err(BetaFuncError::ConvergenceFailed) +} + +pub(super) fn expm1c(x: f64) -> f64 { + if x.abs() < 1e-5 { + 1.0 + x * (0.5 + x * (1.0 / 6.0 + x * (1.0 / 24.0 + x / 120.0))) + } else { + x.exp_m1() / x + } +} + +pub(super) fn ln_gamma_one_plus_over_x(x: f64) -> f64 { + if x <= 1e-4 { + -consts::EULER_MASCHERONI + + x * (0.8224670334241132 + + x * (-0.40068563438653143 + + x * (0.27058080842778455 + + x * (-0.20738555102867398 + x * 0.1695571769974082)))) + } else { + gamma::ln_gamma(1.0 + x) / x + } +} + +pub(super) fn upper_gamma_scaled_small_shape(shape: f64, x: f64) -> Result { + let log_x = x.ln(); + let log_gamma_ratio = ln_gamma_one_plus_over_x(shape); + let difference = log_x - log_gamma_ratio; + let scaled_difference = shape * difference; + let mut term = -x / (shape + 1.0); + let mut sum = term; + let mut compensation = 0.0_f64; + for n in 2..=128 { + let n = f64::from(n); + term *= (-x / n) * (shape + n - 1.0) / (shape + n); + let corrected = term - compensation; + let next = sum + corrected; + compensation = (next - sum) - corrected; + sum = next; + if term.abs() <= prec::F64_PREC * sum.abs() { + let upper_gamma = + -difference * expm1c(scaled_difference) - scaled_difference.exp() * sum; + let result = upper_gamma * (x - scaled_difference).exp(); + return if result > 0.0 && result.is_finite() { + Ok(result) + } else { + Err(BetaFuncError::ConvergenceFailed) + }; + } + } + Err(BetaFuncError::ConvergenceFailed) +} diff --git a/src/function/beta/series.rs b/src/function/beta/series.rs new file mode 100644 index 00000000..7d5e4901 --- /dev/null +++ b/src/function/beta/series.rs @@ -0,0 +1,138 @@ +// (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. +// SPDX-License-Identifier: MIT AND BSL-1.0 +// Use, modification and distribution are subject to the Boost Software +// License, Version 1.0. (See accompanying file LICENSE-BOOST.md or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +use super::*; + +pub(super) fn beta_power_series_log_parts_with_log_beta( + a: f64, + b: f64, + x: f64, + log_beta: Option<(f64, f64)>, +) -> Result<(f64, f64), BetaFuncError> { + let scaled_b = b * x; + let scaled_b = (scaled_b, b.mul_add(x, -scaled_b)); + let a_minus_one = dd_add((a, 0.0), (-1.0, 0.0)); + let mut term = (1.0_f64, 0.0_f64); + let mut sum = (1.0_f64, 0.0_f64); + for n in 1..=MAX_BETA_REG_ITERATIONS { + let n = f64::from(n); + let shape_numerator = dd_add(a_minus_one, (n, 0.0)); + let scaled_numerator = dd_mul(shape_numerator, (x, 0.0)); + let factor = dd_div_f64(dd_add(scaled_numerator, scaled_b), a + n); + term = dd_mul(term, factor); + sum = dd_add(sum, term); + if term.0.abs() <= f64::EPSILON * f64::EPSILON * sum.0.abs() { + if sum.0 <= 0.0 { + return Err(BetaFuncError::ConvergenceFailed); + } + let (log_sum, log_sum_error) = accurate_ln(sum.0); + let log_sum_error = log_sum_error + (sum.1 / sum.0).ln_1p(); + if use_beta_gamma_limit(a, b, scaled_b.0) { + let (log_scaled_b, log_scaled_b_error) = accurate_ln(scaled_b.0); + let log_scaled_b_error = log_scaled_b_error + (scaled_b.1 / scaled_b.0).ln_1p(); + let mut result = dd_mul((a, 0.0), (log_scaled_b, log_scaled_b_error)); + result = dd_add(result, (-scaled_b.0, -scaled_b.1)); + let log_gamma = if a <= 1e-4 { + a * ln_gamma_one_plus_over_x(a) + } else { + gamma::ln_gamma(1.0 + a) + }; + result = dd_add(result, (-log_gamma, 0.0)); + return Ok(dd_add(result, (log_sum, log_sum_error))); + } + let (log_power, log_power_error) = if let Some(log_beta) = log_beta { + beta_reg_log_power_parts_accurate_with_log_beta(a, b, x, log_beta) + } else { + beta_reg_log_power_parts_accurate(a, b, x) + }; + let (variable, variable_error) = two_sum(log_power, log_sum); + let log_a = accurate_ln(a); + return Ok(dd_add( + (variable, variable_error + log_power_error + log_sum_error), + (-log_a.0, -log_a.1), + )); + } + } + Err(BetaFuncError::ConvergenceFailed) +} + +pub(super) fn beta_power_series_log_parts( + a: f64, + b: f64, + x: f64, +) -> Result<(f64, f64), BetaFuncError> { + beta_power_series_log_parts_with_log_beta(a, b, x, None) +} + +pub(super) fn beta_power_series_log(a: f64, b: f64, x: f64) -> Result { + beta_power_series_log_parts(a, b, x).map(|(result, error)| result + error) +} + +pub(super) fn beta_small_shapes_series_log( + a: f64, + b: f64, + x: f64, + y: f64, +) -> Result, BetaFuncError> { + beta_small_shapes_series_log_with_log_beta(a, b, x, y, None) +} + +pub(super) fn beta_small_shapes_series_log_with_log_beta( + a: f64, + b: f64, + x: f64, + y: f64, + log_beta: Option<(f64, f64)>, +) -> Result, BetaFuncError> { + if a.max(b) > 1.0 { + return Ok(None); + } + let invert = !(a >= 0.2_f64.min(b) || x.powf(a) <= 0.9); + let (transformed_a, transformed_b, transformed_x) = if invert { (b, a, y) } else { (a, b, x) }; + if transformed_x > 0.9 { + return Ok(None); + } + beta_power_series_log_parts_with_log_beta(transformed_a, transformed_b, transformed_x, log_beta) + .map(|result| Some((result.0 + result.1, invert))) +} + +pub(super) fn use_beta_gamma_limit(a: f64, b: f64, scaled_x: f64) -> bool { + let correction_scale = a + scaled_x + 1.0; + correction_scale.is_finite() && correction_scale / b.sqrt() <= 0.25 * f64::EPSILON.sqrt() +} + +pub(super) fn use_beta_power_series(a: f64, b: f64, x: f64) -> bool { + let scaled_x = b * x; + x < 1.0 + && ((scaled_x <= 0.7 && x <= 0.95) + || (a <= f64::EPSILON.sqrt() && scaled_x <= 2.0 && x < beta_symmetry_split(a, b)) + || (a <= 0.3 && b >= 32.0 && scaled_x <= 2.0) + || (a <= 40.0 && b >= 32.0 && x < beta_symmetry_split(a, b)) + || (use_beta_gamma_limit(a, b, scaled_x) && scaled_x <= 64.0)) +} + +pub(super) fn use_beta_power_series_before_symmetry(a: f64, b: f64, x: f64) -> bool { + let scaled_x = b * x; + x < 1.0 + && !(a <= f64::EPSILON.sqrt() && b >= STIRLING_MIN && x.powf(a) > 0.5) + && ((a <= f64::EPSILON.sqrt() && scaled_x <= 2.0 && x < beta_symmetry_split(a, b)) + || (a <= 0.3 && b >= 32.0 && scaled_x <= 2.0) + || (a <= 40.0 && b >= 32.0 && x < beta_symmetry_split(a, b)) + || (use_beta_gamma_limit(a, b, scaled_x) && scaled_x <= 64.0)) +} + +pub(super) fn beta_symmetry_split(a: f64, b: f64) -> f64 { + let a1 = a + 1.0; + let b1 = b + 1.0; + let scale = a1.max(b1); + (a1 / scale) / (a1 / scale + b1 / scale) +} + +pub(super) fn use_beta_symmetry(a: f64, b: f64, x: f64) -> bool { + a < 1.0 && a <= f64::EPSILON.sqrt() && b >= STIRLING_MIN && x.powf(a) > 0.5 + || (a < 1.0 || x > f64::EPSILON) && 1.0 - x < 1.0 && x >= beta_symmetry_split(a, b) +} diff --git a/src/function/beta/small_gamma.rs b/src/function/beta/small_gamma.rs new file mode 100644 index 00000000..e0307859 --- /dev/null +++ b/src/function/beta/small_gamma.rs @@ -0,0 +1,89 @@ +use super::*; + +// DLMF 5.7.3: https://dlmf.nist.gov/5.7.E3 +const LN_GAMMA_ONE_PLUS_COEFFICIENTS: [(f64, f64); 31] = [ + (0.8224670334241132, 1.520336175199238e-17), + (-0.40068563438653143, 2.250747042487504e-18), + (0.27058080842778455, 1.1871280107138412e-17), + (-0.20738555102867398, -4.099767328621813e-18), + (0.1695571769974082, 2.2393851330167238e-18), + (-0.14404989676884612, -9.623140085232555e-18), + (0.12550966952474304, -2.5214685384672305e-18), + (-0.11133426586956469, -4.643990572582924e-18), + (0.10009945751278181, 2.6102404859583283e-18), + (-0.09095401714582904, -8.306705457691885e-19), + (0.083353840546109, 2.963832603652642e-19), + (-0.0769325164113522, 3.2900356019181198e-18), + (0.07143294629536133, 6.278806024191499e-18), + (-0.06666870588242047, -3.2295860759966306e-18), + (0.06250095514121304, 2.551099464019315e-18), + (-0.058823978658684585, 2.6912901341966357e-18), + (0.055555767627403614, -3.0261864849830964e-18), + (-0.05263167937961666, -2.523843702471215e-18), + (0.05000004769810169, 2.7894418264458796e-19), + (-0.047619070330142226, -2.4796342684293355e-18), + (0.04545455629320467, 4.382931774550076e-19), + (-0.04347826605304026, 1.8462229880395943e-18), + (0.04166666915034121, 2.308174687248266e-18), + (-0.04000000119214014, -3.145690613937729e-18), + (0.03846153903467518, 3.3927204223959168e-18), + (-0.037037037312989324, -1.7709932414949877e-18), + (0.035714285847333355, 3.3772026865595416e-18), + (-0.034482758684919304, 3.2599869270595477e-18), + (0.03333333336437758, -2.2936827368961794e-18), + (-0.03225806453115042, 1.9360221160020273e-18), + (0.03125000000727597, 2.9882678459447273e-18), +]; + +fn ln_gamma_one_plus_series_prefix_parts(x: f64, count: usize) -> (f64, f64) { + let coefficients = &LN_GAMMA_ONE_PLUS_COEFFICIENTS[..count]; + let (mut high, mut low) = *coefficients.last().unwrap(); + for coefficient in coefficients[..count - 1].iter().rev() { + let product = high * x; + let product_error = high.mul_add(x, -product) + low * x; + let (sum, sum_error) = two_sum(product, coefficient.0); + (high, low) = two_sum(sum, product_error + sum_error + coefficient.1); + } + dd_mul( + (x, 0.0), + dd_add( + (-consts::EULER_MASCHERONI, 4.942915152430645e-18), + dd_mul((x, 0.0), (high, low)), + ), + ) +} + +pub(super) fn ln_gamma_one_plus_series_parts(x: f64) -> (f64, f64) { + ln_gamma_one_plus_series_prefix_parts(x, LN_GAMMA_ONE_PLUS_COEFFICIENTS.len()) +} + +pub(super) fn ln_gamma_one_plus_inverse_parts(x: f64) -> (f64, f64) { + // DLMF 5.7.3. For 0 < x <= 1/8, the k=22 onward remainder is below + // zeta(22) * x^22 / (22 * (1-x)); after division by x it is below 5.7e-21. + ln_gamma_one_plus_series_prefix_parts(x, 20) +} + +pub(super) fn ln_gamma_small_accurate_parts(x: f64) -> (f64, f64) { + let gamma_one_plus = ln_gamma_one_plus_series_parts(x); + let logarithm = accurate_ln(x); + dd_add(gamma_one_plus, (-logarithm.0, -logarithm.1)) +} + +pub(super) fn ln_gamma_one_plus_series(x: f64) -> f64 { + let mut polynomial = LN_GAMMA_ONE_PLUS_COEFFICIENTS.last().unwrap().0; + for coefficient in LN_GAMMA_ONE_PLUS_COEFFICIENTS[..30].iter().rev() { + polynomial = polynomial.mul_add(x, coefficient.0); + } + x * (-consts::EULER_MASCHERONI + x * polynomial) +} + +#[test] +fn inverse_prefix_tracks_the_full_series_on_its_domain() { + for index in 1..=4096 { + let x = 0.125 * f64::from(index) / 4096.0; + let full = ln_gamma_one_plus_series_parts(x); + let prefix = ln_gamma_one_plus_inverse_parts(x); + let difference = dd_add(full, (-prefix.0, -prefix.1)); + assert!((difference.0 + difference.1).abs() < 7.1e-22); + } +} diff --git a/src/function/beta/temme.rs b/src/function/beta/temme.rs new file mode 100644 index 00000000..5c650cb5 --- /dev/null +++ b/src/function/beta/temme.rs @@ -0,0 +1,136 @@ +#[cfg(all(not(feature = "std"), not(test)))] +use super::Float; +use super::{dd_add, dd_div, dd_mul}; + +// DLMF 8.18.9--12 gives the uniform expansion and c0. The centered series and +// c1 recurrence follow Temme, Special Functions (1996), section 11.3.3.2. + +const SERIES_DEGREE: usize = 6; + +fn series_multiply( + left: &[f64; SERIES_DEGREE + 1], + right: &[f64; SERIES_DEGREE + 1], +) -> [f64; SERIES_DEGREE + 1] { + let mut result = [0.0; SERIES_DEGREE + 1]; + for index in 0..=SERIES_DEGREE { + for offset in 0..=index { + result[index] = left[offset].mul_add(right[index - offset], result[index]); + } + } + result +} + +fn series_reciprocal(value: &[f64; SERIES_DEGREE + 1]) -> [f64; SERIES_DEGREE + 1] { + let mut result = [0.0; SERIES_DEGREE + 1]; + result[0] = 1.0 / value[0]; + for index in 1..=SERIES_DEGREE { + let mut sum = 0.0; + for offset in 1..=index { + sum = value[offset].mul_add(result[index - offset], sum); + } + result[index] = -sum / value[0]; + } + result +} + +fn series_sqrt(value: &[f64; SERIES_DEGREE + 1]) -> [f64; SERIES_DEGREE + 1] { + let mut result = [0.0; SERIES_DEGREE + 1]; + result[0] = value[0].sqrt(); + for index in 1..=SERIES_DEGREE { + let mut sum = 0.0; + for offset in 1..index { + sum = result[offset].mul_add(result[index - offset], sum); + } + result[index] = (value[index] - sum) / (2.0 * result[0]); + } + result +} + +fn evaluate(coefficients: &[f64], argument: f64) -> f64 { + coefficients.iter().rev().fold(0.0, |value, coefficient| { + value.mul_add(argument, *coefficient) + }) +} + +pub(super) fn temme_coefficients(mean: f64, delta: f64) -> (f64, f64) { + let complement = 1.0 - mean; + let variance_root = (mean * complement).sqrt(); + let mut deviance = [0.0; SERIES_DEGREE + 1]; + let mut mean_power = mean; + let mut complement_power = complement; + for (index, coefficient) in deviance.iter_mut().enumerate() { + let order = index + 2; + let sign = if order & 1 == 0 { 1.0 } else { -1.0 }; + *coefficient = + 2.0 * mean * complement * (sign * mean_power.recip() + complement_power.recip()) + / order as f64; + mean_power *= mean; + complement_power *= complement; + } + + let eta_over_delta = series_sqrt(&deviance); + let delta_over_eta = series_reciprocal(&eta_over_delta); + let mut c0_coefficients = [0.0; SERIES_DEGREE + 1]; + for index in 0..SERIES_DEGREE { + c0_coefficients[index] = variance_root * delta_over_eta[index + 1]; + } + + let mut eta_derivative = [0.0; SERIES_DEGREE + 1]; + for index in 0..=SERIES_DEGREE { + eta_derivative[index] = (index + 1) as f64 * eta_over_delta[index]; + } + let mut eta_derivative_reciprocal = series_reciprocal(&eta_derivative); + for coefficient in &mut eta_derivative_reciprocal { + *coefficient *= variance_root; + } + let mut c0_delta_derivative = [0.0; SERIES_DEGREE + 1]; + for index in 0..SERIES_DEGREE { + c0_delta_derivative[index] = (index + 1) as f64 * c0_coefficients[index + 1]; + } + let c0_eta_derivative = series_multiply(&c0_delta_derivative, &eta_derivative_reciprocal); + let mut divided_numerator = [0.0; SERIES_DEGREE + 1]; + for index in 0..SERIES_DEGREE { + divided_numerator[index] = -c0_eta_derivative[index + 1]; + } + let mut c1_coefficients = series_multiply(&delta_over_eta, ÷d_numerator); + let stirling = (1.0 - mean.recip() - complement.recip()) / 12.0; + for index in 0..=SERIES_DEGREE { + c1_coefficients[index] = + variance_root.mul_add(c1_coefficients[index], -stirling * c0_coefficients[index]); + } + + ( + evaluate(&c0_coefficients, delta), + evaluate(&c1_coefficients, delta), + ) +} + +pub(super) fn temme_delta(a: f64, b: f64, x: f64) -> (f64, (f64, f64)) { + let scale = a.max(b); + let scaled_a = (a / scale, (-a / scale).mul_add(scale, a) / scale); + let scaled_b = (b / scale, (-b / scale).mul_add(scale, b) / scale); + let scaled_sum = dd_add(scaled_a, scaled_b); + let mean = dd_div(scaled_a, scaled_sum); + let numerator = dd_add(dd_mul((x, 0.0), scaled_sum), (-scaled_a.0, -scaled_a.1)); + let delta = dd_div(numerator, scaled_sum); + (mean.0 + mean.1, delta) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn coefficients_match_symbolic_reference() { + let (c0, c1) = temme_coefficients(1.0 / 3.0, -2.7216552652253867e-5); + assert!((c0 - 0.23571910018211537).abs() < 2e-15, "c0={c0:?}"); + assert!((c1 - -0.0360076116745462).abs() < 2e-14, "c1={c1:?}"); + } + + #[test] + fn compensated_delta_matches_exact_input_reference() { + let (mean, delta) = temme_delta(40_000_000.0, 80_000_000.0, 1.0 / 3.0); + assert_eq!(mean, 1.0 / 3.0); + assert!((delta.0 + delta.1 - -1.850371707708594e-17).abs() < 1e-32); + } +} diff --git a/src/function/beta/tests.rs b/src/function/beta/tests.rs new file mode 100644 index 00000000..21c7d8a7 --- /dev/null +++ b/src/function/beta/tests.rs @@ -0,0 +1,2046 @@ +use super::*; +use crate::prec; +use core::f64::consts as f64_consts; +const MODULE_RELATIVE_ACC: f64 = 1e-14; +const INVERSE_REFERENCE_MAX_ULPS: u64 = 4; + +// Tests named `*_500_digit_*` use exact decimal encodings of the input `f64`s +// with `cpp_dec_float<500>` and one final round to binary64. Inverse references +// were independently checked with mpmath at 550 decimal digits. + +fn beta_assert_relative_eq(a: f64, b: f64) { + prec::assert_relative_eq!( + a, + b, + epsilon = MODULE_EPS, + max_relative = MODULE_RELATIVE_ACC + ); +} + +fn beta_assert_abs_diff_eq(a: f64, b: f64) { + prec::assert_abs_diff_eq!(a, b, epsilon = MODULE_EPS); +} + +#[test] +fn test_ln_beta() { + beta_assert_relative_eq(ln_beta(0.5, 0.5), 1.144729885849400174144); + beta_assert_relative_eq(ln_beta(1.0, 0.5), f64_consts::LN_2); + beta_assert_relative_eq(ln_beta(2.5, 0.5), 0.163900632837673937284); + beta_assert_relative_eq(ln_beta(0.5, 1.0), f64_consts::LN_2); + beta_assert_relative_eq(ln_beta(1.0, 1.0), 0.0); + beta_assert_relative_eq(ln_beta(2.5, 1.0), -0.9162907318741550651835); + beta_assert_relative_eq(ln_beta(0.5, 2.5), 0.163900632837673937284); + beta_assert_relative_eq(ln_beta(1.0, 2.5), -0.9162907318741550651835); + beta_assert_relative_eq(ln_beta(2.5, 2.5), -2.608688089402107300388); +} + +#[test] +#[should_panic] +fn test_ln_beta_a_lte_0() { + ln_beta(0.0, 0.5); +} + +#[test] +#[should_panic] +fn test_ln_beta_b_lte_0() { + ln_beta(0.5, 0.0); +} + +#[test] +fn test_checked_ln_beta_a_lte_0() { + assert!(checked_ln_beta(0.0, 0.5).is_err()); +} + +#[test] +fn test_checked_ln_beta_b_lte_0() { + assert!(checked_ln_beta(0.5, 0.0).is_err()); +} + +#[test] +fn test_checked_ln_beta_reg_validates_inputs_and_boundaries() { + assert_eq!( + checked_ln_beta_reg(0.0, 1.0, 0.5), + Err(BetaFuncError::ANotGreaterThanZero) + ); + assert_eq!( + checked_ln_beta_reg(1.0, 0.0, 0.5), + Err(BetaFuncError::BNotGreaterThanZero) + ); + assert_eq!( + checked_ln_beta_reg(1.0, 1.0, f64::NAN), + Err(BetaFuncError::XOutOfRange) + ); + assert_eq!(checked_ln_beta_reg(2.0, 3.0, 0.0), Ok(f64::NEG_INFINITY)); + assert_eq!(checked_ln_beta_reg(2.0, 3.0, 1.0), Ok(0.0)); + assert_eq!(checked_ln_beta_reg(2.0, 2.0, 0.5), Ok(-f64_consts::LN_2)); + assert_eq!(checked_ln_beta_reg(2.0, 1.0, 0.25), Ok(-2.772588722239781)); + let actual = checked_ln_beta_reg(0.5, 1e-8, 0.95).unwrap(); + let expected = f64::from_bits(0xc030f2f1c4aab2eb); + assert!( + actual.to_bits().abs_diff(expected.to_bits()) <= 1, + "actual={actual:e}, expected={expected:e}" + ); +} + +#[test] +fn test_checked_ln_beta_reg_complement_validates_inputs_and_boundaries() { + assert_eq!( + checked_ln_beta_reg_complement(0.0, 1.0, 0.5), + Err(BetaFuncError::ANotGreaterThanZero) + ); + assert_eq!( + checked_ln_beta_reg_complement(1.0, 0.0, 0.5), + Err(BetaFuncError::BNotGreaterThanZero) + ); + assert_eq!( + checked_ln_beta_reg_complement(1.0, 1.0, f64::NAN), + Err(BetaFuncError::XOutOfRange) + ); + assert_eq!(checked_ln_beta_reg_complement(2.0, 3.0, 0.0), Ok(0.0)); + assert_eq!( + checked_ln_beta_reg_complement(2.0, 3.0, 1.0), + Ok(f64::NEG_INFINITY) + ); + let actual = checked_ln_beta_reg_complement(2.0, 3.0, 0.75).unwrap(); + let expected = f64::from_bits(0xc007d781d2c2e78d); + assert!( + actual.to_bits().abs_diff(expected.to_bits()) <= 1, + "actual={actual:e}, expected={expected:e}" + ); +} + +#[test] +#[should_panic] +fn test_beta_a_lte_0() { + beta(0.0, 0.5); +} + +#[test] +#[should_panic] +fn test_beta_b_lte_0() { + beta(0.5, 0.0); +} + +#[test] +fn test_checked_beta_a_lte_0() { + assert!(checked_beta(0.0, 0.5).is_err()); +} + +#[test] +fn test_checked_beta_b_lte_0() { + assert!(checked_beta(0.5, 0.0).is_err()); +} + +#[test] +fn test_beta() { + beta_assert_relative_eq(beta(0.5, 0.5), f64_consts::PI); + beta_assert_relative_eq(beta(1.0, 0.5), 2.0); + beta_assert_relative_eq(beta(2.5, 0.5), 1.17809724509617246442); + beta_assert_relative_eq(beta(0.5, 1.0), 2.0); + beta_assert_relative_eq(beta(1.0, 1.0), 1.0); + beta_assert_relative_eq(beta(2.5, 1.0), 0.4); + beta_assert_relative_eq(beta(0.5, 2.5), 1.17809724509617246442); + beta_assert_relative_eq(beta(1.0, 2.5), 0.4); + beta_assert_relative_eq(beta(2.5, 2.5), 0.073631077818510779026); +} + +#[test] +fn test_beta_inc() { + beta_assert_relative_eq(beta_inc(0.5, 0.5, 0.5), f64_consts::FRAC_PI_2); + beta_assert_relative_eq(beta_inc(0.5, 0.5, 1.0), f64_consts::PI); + beta_assert_relative_eq(beta_inc(1.0, 0.5, 0.5), 0.5857864376269049511983); + beta_assert_relative_eq(beta_inc(1.0, 0.5, 1.0), 2.0); + beta_assert_relative_eq(beta_inc(2.5, 0.5, 0.5), 0.0890486225480862322117); + beta_assert_relative_eq(beta_inc(2.5, 0.5, 1.0), 1.17809724509617246442); + beta_assert_relative_eq(beta_inc(0.5, 1.0, 0.5), f64_consts::SQRT_2); + beta_assert_relative_eq(beta_inc(0.5, 1.0, 1.0), 2.0); + beta_assert_relative_eq(beta_inc(1.0, 1.0, 0.5), 0.5); + beta_assert_relative_eq(beta_inc(1.0, 1.0, 1.0), 1.0); + beta_assert_relative_eq(beta_inc(2.5, 1.0, 0.5), 0.0707106781186547524401); + beta_assert_relative_eq(beta_inc(2.5, 1.0, 1.0), 0.4); + beta_assert_relative_eq(beta_inc(0.5, 2.5, 0.5), 1.08904862254808623221); + beta_assert_relative_eq(beta_inc(0.5, 2.5, 1.0), 1.17809724509617246442); + beta_assert_relative_eq(beta_inc(1.0, 2.5, 0.5), 0.32928932188134524756); + beta_assert_relative_eq(beta_inc(1.0, 2.5, 1.0), 0.4); + beta_assert_relative_eq(beta_inc(2.5, 2.5, 0.5), 0.03681553890925538951323); + beta_assert_relative_eq(beta_inc(2.5, 2.5, 1.0), 0.073631077818510779026); +} + +#[test] +#[should_panic] +fn test_beta_inc_a_lte_0() { + beta_inc(0.0, 1.0, 1.0); +} + +#[test] +#[should_panic] +fn test_beta_inc_b_lte_0() { + beta_inc(1.0, 0.0, 1.0); +} + +#[test] +#[should_panic] +fn test_beta_inc_x_lt_0() { + beta_inc(1.0, 1.0, -1.0); +} + +#[test] +#[should_panic] +fn test_beta_inc_x_gt_1() { + beta_inc(1.0, 1.0, 2.0); +} + +#[test] +fn test_checked_beta_inc_a_lte_0() { + assert!(checked_beta_inc(0.0, 1.0, 1.0).is_err()); +} + +#[test] +fn test_checked_beta_inc_b_lte_0() { + assert!(checked_beta_inc(1.0, 0.0, 1.0).is_err()); +} + +#[test] +fn test_checked_beta_inc_x_lt_0() { + assert!(checked_beta_inc(1.0, 1.0, -1.0).is_err()); +} + +#[test] +fn test_checked_beta_inc_x_gt_1() { + assert!(checked_beta_inc(1.0, 1.0, 2.0).is_err()); +} + +#[test] +fn test_beta_reg() { + beta_assert_abs_diff_eq(beta_reg(0.5, 0.5, 0.5), 0.5); + assert_eq!(beta_reg(0.5, 0.5, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(1.0, 0.5, 0.5), 0.292893218813452475599); + assert_eq!(beta_reg(1.0, 0.5, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(2.5, 0.5, 0.5), 0.07558681842161243795); + assert_eq!(beta_reg(2.5, 0.5, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(0.5, 1.0, 0.5), f64_consts::FRAC_1_SQRT_2); + assert_eq!(beta_reg(0.5, 1.0, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(1.0, 1.0, 0.5), 0.5); + assert_eq!(beta_reg(1.0, 1.0, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(2.5, 1.0, 0.5), 0.1767766952966368811); + assert_eq!(beta_reg(2.5, 1.0, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(0.5, 2.5, 0.5), 0.92441318157838756205); + assert_eq!(beta_reg(0.5, 2.5, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(1.0, 2.5, 0.5), 0.8232233047033631189); + assert_eq!(beta_reg(1.0, 2.5, 1.0), 1.0); + beta_assert_abs_diff_eq(beta_reg(2.5, 2.5, 0.5), 0.5); + assert_eq!(beta_reg(2.5, 2.5, 1.0), 1.0); +} + +#[test] +fn test_beta_reg_large_parameters_against_reference() { + let cases = [ + (1e6, 2e6, 0.333, 0.11032283951664962), + (1e6, 2e6, 1.0 / 3.0, 0.5000542891707268), + (1e6, 2e6, 0.334, 0.9928335645421132), + (1e8, 2e8, 0.3333, 0.11033439854811466), + (1e8, 2e8, 1.0 / 3.0, 0.5000054289165304), + (1e8, 2e8, 0.3334, 0.992845709515461), + (1e5, 1e5, 0.49, 1.8571347290404196e-19), + (1e5, 1e5, 0.499, 0.18554674455755675), + (1e5, 1e5, 0.501, 0.8144532554424433), + (40.0, 32.0, 1e-8, 1.2676414050441584e-300), + (32.0, 40.0, 1e-8, 1.5845516362868252e-236), + (0.1, 1e8, 1e-8, 0.9758726562930068), + (0.1, 1e8, 1e-9, 0.8275517592836537), + (2.0, 1e8, 1e-8, 0.2642411213359098), + (10.0, 1e8, 1e-7, 0.5420704043826821), + (1e13, 9.9e14, 0.01, 0.5000000414451727), + ( + 1.098252731340299, + 1.780042655540735e17, + 5.235783704840033e-17, + 0.999881646675342, + ), + ( + 7_627_209.761, + 11.3319, + 0.9999105965110135, + 1.6790000011611638e-274, + ), + (99_999.0, 11.3319, 0.9998, 0.013667998876668642), + (100_001.0, 11.3319, 0.9998, 0.013665136770782414), + (100_000.0, 10.0, 0.992653308338289, 1.0000000000029653e-300), + ]; + + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + let error = (actual - expected).abs(); + let tolerance = 5e-12 * expected.max(1e-300); + assert!( + error <= tolerance, + "beta_reg({a}, {b}, {x}) = {actual}, expected {expected}, error {error}" + ); + } +} + +#[test] +fn test_beta_reg_large_asymmetric_temme_focus_against_500_digit_references() { + let cases = [ + ( + 100_000_000.0, + 200_000_000.0, + 0x3fd554e32dc84e59_u64, + 0x3fc44ed0bc353f04_u64, + ), + ( + 12_000_000.0, + 108_000_000.0, + 0x3fb99926bbf81f29, + 0x3fd9af470acc030a, + ), + ( + 40_000_000.0, + 80_000_000.0, + 0x3fd5555555555555, + 0x3fe00012006e3f11, + ), + ( + 108_000_000.0, + 12_000_000.0, + 0x3fecccbe71189d7f, + 0x3fd9ae504ae8a2e1, + ), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "a={a:?}, b={b:?}, x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + } +} + +#[test] +fn test_beta_reg_large_asymmetric_temme_grid_against_500_digit_references() { + let cases = [ + ( + 100_000_000.0, + 900_000_000.0, + 0x3fb998fa6ff66eb1_u64, + 0x3fc44ed0bb21a1af_u64, + ), + ( + 100_000_000.0, + 900_000_000.0, + 0x3fb999999999999a, + 0x3fe00017846dc7c6, + ), + ( + 100_000_000.0, + 900_000_000.0, + 0x3fb99a38c33cc483, + 0x3feaec4bd137a086, + ), + ( + 333_333_333.333_333_3, + 666_666_666.666_666_7, + 0x3fd55516ceef6eb5, + 0x3fc44ed0bbb46ae2, + ), + ( + 333_333_333.333_333_3, + 666_666_666.666_666_7, + 0x3fd5555555555555, + 0x3fe000063c68549a, + ), + ( + 333_333_333.333_333_3, + 666_666_666.666_666_7, + 0x3fd55593dbbb3bf5, + 0x3feaec4bd112fd8f, + ), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "a={a:?}, b={b:?}, x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + } +} + +#[test] +fn test_beta_reg_asymptotic_deviance_gate_against_500_digit_references() { + let groups = [ + ( + 1_000_000.0, + 9_000_000.0, + [ + (0x3fb97fe3e2209f30_u64, 0x3ef23a83b513c998_u64), + (0x3fb97fe3e2209f31, 0x3ef23a83b513d668), + (0x3fb97fe3e2209f32, 0x3ef23a83b513e337), + (0x3fb97fe3e2209f33, 0x3ef23a83b513f006), + (0x3fb97fe3e2209f34, 0x3ef23a83b513fcd6), + ], + ), + ( + 1_000_100.0, + 98_999_900.0, + [ + (0x3f846555255c20b9_u64, 0x3ee7d0fb4a5ec8c2_u64), + (0x3f846555255c20ba, 0x3ee7d0fb4a5edd23), + (0x3f846555255c20bb, 0x3ee7d0fb4a5ef184), + (0x3f846555255c20bc, 0x3ee7d0fb4a5f05e5), + (0x3f846555255c20bd, 0x3ee7d0fb4a5f1a46), + ], + ), + ]; + for (a, b, cases) in groups { + let mut previous = 0_u64; + for (x, expected) in cases { + let actual = beta_reg(a, b, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "a={a:?}, b={b:?}, x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + assert!(actual > previous, "a={a:?}, b={b:?}, x={x:#018x}"); + previous = actual; + } + } +} + +#[test] +fn test_beta_reg_asymptotic_shape_gates_against_500_digit_references() { + let groups = [ + ( + "sum", + [ + ( + 3_333_333.0, + 6_666_666.0, + 0x3fd5555555555555_u64, + 0x3fe0003e5c12c87a_u64, + 1_024_u64, + ), + ( + 3_333_333.0, + 6_666_667.0, + 0x3fd55555318abc87, + 0x3fe0003e5c138137, + 4, + ), + ( + 3_333_334.0, + 6_666_667.0, + 0x3fd55555791fede8, + 0x3fe0003e5c117848, + 4, + ), + ], + ), + ( + "minimum shape", + [ + ( + 999_999.0, + 19_000_001.0, + 0x3fa99997ec1a6fee, + 0x3fe0010183687a50, + 256, + ), + ( + 1_000_000.0, + 19_000_000.0, + 0x3fa999999999999a, + 0x3fe00101835e9c34, + 4, + ), + ( + 1_000_001.0, + 18_999_999.0, + 0x3fa9999b4718c345, + 0x3fe001018354bc19, + 4, + ), + ], + ), + ( + "mean", + [ + ( + 999_999.0, + 99_000_001.0, + 0x3f847adff0152658, + 0x3fe00112ae2480be, + 128, + ), + ( + 1_000_000.0, + 99_000_000.0, + 0x3f847ae147ae147b, + 0x3fe00112ae1b39b3, + 4, + ), + ( + 1_000_001.0, + 98_999_999.0, + 0x3f847ae29f47029e, + 0x3fe00112ae11f2a9, + 4, + ), + ], + ), + ]; + for (gate, cases) in groups { + for (a, b, x, expected, max_ulp) in cases { + let actual = beta_reg(a, b, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= max_ulp, + "gate={gate}, a={a:?}, b={b:?}, x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + } +} + +#[test] +fn test_beta_reg_large_symmetric_adjacent_against_500_digit_references() { + // cpp_dec_float<500>, with each f64 input converted from its exact binary ratio. + let lower = f64::from_bits(0.5_f64.to_bits() - 1); + let upper = f64::from_bits(0.5_f64.to_bits() + 1); + let cases = [ + (1e2, 0x3fdffffffffffff5_u64, 0x3fe000000000000b_u64), + (1e3, 0x3fdfffffffffffdc, 0x3fe0000000000024), + (1e4, 0x3fdfffffffffff8f, 0x3fe0000000000071), + (1e5, 0x3fdffffffffffe9b, 0x3fe0000000000165), + (1e6, 0x3fdffffffffffb98, 0x3fe0000000000468), + (1e7, 0x3fdffffffffff210, 0x3fe0000000000df0), + (1e8, 0x3fdfffffffffd3ec, 0x3fe0000000002c14), + ]; + for (shape, lower_expected, upper_expected) in cases { + for (x, expected) in [(lower, lower_expected), (upper, upper_expected)] { + let actual = beta_reg(shape, shape, x).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "shape={shape:?}, x={x:?}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + } +} + +#[test] +fn test_beta_reg_symmetric_central_boundary_against_500_digit_references() { + // mpmath at 550 digits using exact f64 ratios and the equivalent gamma/hyp2f1 identity. + let cases = [ + ( + 100.0, + 0x3fddbcbcf5c0139f_u64, + 0x3fc44eca5b83b728_u64, + 0x3fe121a1851ff630_u64, + 0x3feaec4d691f1233_u64, + ), + ( + 100.125, + 0x3fddbd198e18a036, + 0x3fc44eca5f98d22f, + 0x3fe1217338f3afe5, + 0x3feaec4d6819cb74, + ), + ( + 1e6, + 0x3fdffa3516f00033, + 0x3fc44ed0bb7cb51c, + 0x3fe002e57487ffe6, + 0x3feaec4bd120d163, + ), + ( + 1e12, + 0x3fdffffe845ffbe8, + 0x3fc44ed0bb87ad45, + 0x3fe00000bdd0020c, + 0x3feaec4bd11e14af, + ), + ]; + for (shape, lower_x, lower_expected, upper_x, upper_expected) in cases { + assert_eq!(beta_reg(shape, shape, 0.5), 0.5); + for (x, expected) in [(lower_x, lower_expected), (upper_x, upper_expected)] { + let actual = beta_reg(shape, shape, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "shape={shape:?}, x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + } +} + +#[test] +fn test_beta_reg_extreme_ratio_central_value_against_reference() { + let cases: [(f64, f64, f64, f64); 2] = [ + ( + 1.2e7, + 1.2000000000000001e307, + 9.999999999999999e-301, + 0.50003838823874907, + ), + (1.2e7, 1e308, 1.2e-301, 0.50003838823881181), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + assert!( + actual.to_bits().abs_diff(expected.to_bits()) <= 1024, + "beta_reg({a}, {b}, {x}) = {actual}, expected {expected}" + ); + } +} + +#[test] +fn test_beta_reg_overflowing_shape_sum() { + let lower = f64::from_bits(0.5_f64.to_bits() - 1); + let upper = f64::from_bits(0.5_f64.to_bits() + 1); + assert_eq!(beta_reg(1e308, 1e308, lower), 0.0); + assert_eq!(beta_reg(1e308, 1e308, 0.5), 0.5); + assert_eq!(beta_reg(1e308, 1e308, upper), 1.0); + let actual = checked_ln_beta(1e308, 1e308).unwrap(); + assert!(actual.is_finite()); + assert!((actual / 1e308 + 2.0 * core::f64::consts::LN_2).abs() <= 2e-15); + let expected = -2.0007184997951635e301; + let actual = checked_ln_beta(f64::MAX, 1e300).unwrap(); + assert!(((actual - expected) / expected).abs() <= 3e-10); + + let mean = f64::from_bits(0x3fe5555555555555); + assert_eq!(beta_reg(1e308, 5e307, mean), 0.0); + assert_eq!( + beta_reg(1e308, 5e307, f64::from_bits(mean.to_bits() + 1)), + 1.0 + ); +} + +#[test] +fn test_beta_reg_algorithm_boundaries_against_reference() { + let cases = [ + (39_999_999.0, 79_999_999.0, 0.33335, 0.6507629787874431), + (40_000_001.0, 80_000_001.0, 0.33335, 0.6507151999304125), + (29_999_999.0, 270_000_001.0, 0.10001, 0.7182251069092127), + (30_000_001.0, 269_999_999.0, 0.10001, 0.7180951316317142), + (1e8, 2e8, 0.33328635138267637, 0.042150859881784875), + (1e8, 2e8, 0.33328603712606697, 0.04112293252416181), + ]; + + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 2e-12, + "beta_reg({a}, {b}, {x}) = {actual}, expected {expected}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_beta_reg_large_a_small_b_subnormal() { + let cases = [ + ( + 1e18, + 39.9, + f64::from_bits(0x3feffffffffffff8), + f64::from_bits(0x1520b9), + ), + (1e8, 0.9, 0.99999284, f64::from_bits(0xfce148c723)), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + assert!( + actual.to_bits().abs_diff(expected.to_bits()) <= 4, + "beta_reg({a}, {b}, {x}) = {actual:e} ({:#x}), expected {expected:e} ({:#x})", + actual.to_bits(), + expected.to_bits() + ); + } +} + +#[test] +fn test_beta_reg_large_a_tiny_b_rounded_complement() { + let x = f64::from_bits(1.0_f64.to_bits() - 1); + let cases = [ + ( + 1.7492718718060828e16, + 1.7529350052864036e-11, + f64::from_bits(0x3d7057be8b9ff83b), + ), + ( + 2.6496319847741348e16, + 3.8997923472821135e-12, + f64::from_bits(0x3d2edb1e5cecbc3f), + ), + ( + 1.3443603650606364e16, + 3.8682302848162155e-10, + f64::from_bits(0x3dc581e85bf535df), + ), + ( + 1.4398454548018444e16, + 1.1381500822684144e-10, + f64::from_bits(0x3da5a5b386b28adf), + ), + ( + 9_288_475_808_954_264.0, + 5.299156768316511e-9, + f64::from_bits(0x3e12f55b03b79471), + ), + ( + 1.6977806187270128e16, + 5.491909396055591e-12, + f64::from_bits(0x3d562f56c473937b), + ), + ]; + for (a, b, expected) in cases { + let actual = beta_reg(a, b, x); + assert!( + actual.to_bits().abs_diff(expected.to_bits()) <= 64, + "beta_reg({a}, {b}, {x}) = {actual:e}, expected {expected:e}" + ); + } +} + +#[test] +fn test_beta_reg_small_shape_upper_gamma_against_reference() { + let cases = [ + ( + 112_176_097_488.593_9, + 1.3959752253898728e-12, + f64::from_bits(0x3fefffffffff851d), + [ + 9.999569151432288e-13, + 9.999869047052781e-13, + 1.000016895594139e-12, + ], + ), + ( + 238_641_107_383.443_27, + 1.799146819367202e-12, + f64::from_bits(0x3fefffffffffb5cc), + [ + 9.999141915967447e-13, + 9.999714463464230e-13, + 1.000028705627258e-12, + ], + ), + ( + 246.932962952654, + 1.1953991131275682e-12, + f64::from_bits(0x3feffffee33a9e66), + [ + 9.999999999706979e-12, + 9.999999999957152e-12, + 1.000000000020733e-11, + ], + ), + ]; + for (a, b, x, expected) in cases { + for (offset, expected) in [-1_i64, 0, 1].into_iter().zip(expected) { + let x = f64::from_bits(x.to_bits().wrapping_add_signed(offset)); + let actual = beta_reg(a, b, x); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-13, + "beta_reg({a}, {b}, {x}) = {actual:e}, expected {expected:e}, relative error {relative_error}" + ); + } + } +} + +#[test] +fn test_ln_beta_reg_tiny_shape_scaled_gamma_against_reference() { + let cases = [ + ( + f64::from_bits(0x3feffffbce423b02), + [ + f64::from_bits(0xc085ae5914154dec), + f64::from_bits(0xc085ae59141548a5), + f64::from_bits(0xc085ae591415435d), + ], + ), + ( + f64::from_bits(0x3fefffeb0750a667), + [ + f64::from_bits(0xc085f9547c11ffd4), + f64::from_bits(0xc085f9547c11fbaa), + f64::from_bits(0xc085f9547c11f77f), + ], + ), + ( + f64::from_bits(0x3fefff7be22e5816), + [ + f64::from_bits(0xc087af793037cd6d), + f64::from_bits(0xc087af793037c98d), + f64::from_bits(0xc087af793037c5ad), + ], + ), + ]; + for (x, expected) in cases { + for (offset, expected) in [-1_i64, 0, 1].into_iter().zip(expected) { + let x = f64::from_bits(x.to_bits().wrapping_add_signed(offset)); + let actual = checked_ln_beta_reg(1e6, 1e-300, x).unwrap(); + assert!((actual - expected).abs() <= 2e-13); + } + } +} + +#[test] +fn test_ln_beta_reg_power_series_is_locally_monotone() { + let cases = [ + ( + 9.11327743985456, + 133_525_174_076_797.34, + f64::from_bits(0x3cf3d7e149ac36dd), + ), + ( + 6.078046923216118, + 31_131_628_187_944.344, + f64::from_bits(0x3cf592225b607c93), + ), + ]; + for (a, b, root) in cases { + let mut previous = f64::NEG_INFINITY; + for offset in -100_i64..=100 { + let x = f64::from_bits(root.to_bits().wrapping_add_signed(offset)); + let value = checked_ln_beta_reg(a, b, x).unwrap(); + assert!(value >= previous, "a={a}, b={b}, x={x}"); + previous = value; + } + } +} + +#[test] +fn test_beta_reg_power_series_is_locally_monotone() { + let cases: [(f64, f64, f64); 7] = [ + ( + 0.47937889777569664, + 390_713_368_494_940.25, + 5.842150555453333e-16, + ), + ( + 0.20713927131052443, + 1_264_447_072_006_281.8, + 7.355559632987759e-17, + ), + ( + 0.5883286844875396, + 53_930_034_336_347.77, + 2.7619798816617607e-15, + ), + ( + 0.3047929367901273, + 258_195_370_359_324.8, + 1.5384576649827977e-15, + ), + ( + 0.21280081734067854, + 54_626_561.16286868, + 4.363878090733803e-9, + ), + ( + 42.51394493556042, + 2_256_890_178_438.929, + 1.0526514336858459e-13, + ), + ( + 77.54913939933753, + 14_481_621_713.827797, + 2.8605493321691776e-11, + ), + ]; + for (a, b, center) in cases { + let mut previous = 0.0; + for offset in -64_i64..=64 { + let x = f64::from_bits(center.to_bits().wrapping_add_signed(offset)); + let value = beta_reg(a, b, x); + assert!(value >= previous, "a={a}, b={b}, x={x}"); + previous = value; + } + } +} + +#[test] +fn test_beta_reg_power_series_subnormal_result_against_reference() { + let actual = beta_reg( + 147.13149557601173, + 1.6465152935404156e16, + f64::from_bits(0x3c78ef1d912aaa46), + ); + assert_eq!(actual.to_bits(), 4); +} + +#[test] +fn test_beta_reg_power_series_boundary_against_reference() { + let (log_beta, log_beta_error) = ln_beta_accurate_parts(10.0, 32.0); + assert_eq!(log_beta.to_bits(), 0xc03723e193251f2a); + assert!((log_beta_error - f64::from_bits(0xbcd496eeab49e82c)).abs() <= 2e-19); + let cases = [ + (0x3fcfffffffffff7f, 0x3fe30d694d7fb0f1), + (0x3fcfffffffffff80, 0x3fe30d694d7fb0f2), + (0x3fcfffffffffff81, 0x3fe30d694d7fb0f4), + (0x3fcfffffffffff82, 0x3fe30d694d7fb0f5), + ]; + let mut previous = 0; + for (x, expected) in cases { + let actual = beta_reg(10.0, 32.0, f64::from_bits(x)).to_bits(); + assert!( + actual.abs_diff(expected) <= 2, + "x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + assert!( + actual > previous, + "x={x:#018x}, actual={actual:#018x}, previous={previous:#018x}" + ); + previous = actual; + } +} + +#[test] +fn test_beta_reg_near_one_moderate_shapes_converges() { + let x = f64::from_bits(1.0_f64.to_bits() - 1); + for (a, b) in [(39.9, 40.0), (40.0, 40.0), (40.0, 41.0)] { + let actual = checked_beta_reg(a, b, x).unwrap(); + assert!( + (0.0..=1.0).contains(&actual), + "a={a}, b={b}, actual={actual:?}" + ); + } +} + +#[test] +fn test_beta_reg_near_one_uses_convergent_power_series() { + let x = f64::from_bits(1.0_f64.to_bits() - 1); + let cases = [ + (217348.9453342118, 7.083729216298346e17), + (74.50754210941346, 4.6813710928374765e17), + (13.940004463756644, 5.294575065065153e17), + ]; + for (a, b) in cases { + let actual = checked_beta_reg(a, b, x).unwrap(); + assert!( + (0.0..=1.0).contains(&actual), + "a={a}, b={b}, actual={actual:?}" + ); + } +} + +#[test] +fn test_beta_reg_tiny_first_shape_remains_monotone_below_split() { + let a = 2.1856409177373306e-11; + let b = 18.619031676940928; + let references = [ + (0x3ea669742f6d91e9_u64, 0x3fefffffffdfb936_u64), + (0x3fa7d0724ba189c0_u64, 0x3fefffffffff2a9e_u64), + ]; + let mut previous = 0_u64; + for (x, expected) in references { + let actual = checked_beta_reg(a, b, f64::from_bits(x)).unwrap().to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + assert!(actual >= previous); + previous = actual; + } +} + +#[test] +fn test_beta_reg_exact_complement_fraction_against_reference() { + let center = 0x3ee7118258b21dd3_u64; + let references = [ + (-128_i64, 0x3fe51a846b074d53_u64), + (-64, 0x3fe51a846b074dbd), + (-1, 0x3fe51a846b074e25), + (0, 0x3fe51a846b074e27), + (1, 0x3fe51a846b074e29), + (64, 0x3fe51a846b074e91), + (128, 0x3fe51a846b074efb), + ]; + for (offset, expected) in references { + let x = f64::from_bits(center.wrapping_add_signed(offset)); + let actual = checked_beta_reg(10.0, 1e6, x).unwrap().to_bits(); + assert!( + actual.abs_diff(expected) <= 3, + "offset={offset}, actual={actual:#018x}, expected={expected:#018x}" + ); + } + let mut previous = 0.0; + for bits in center - 128..=center + 128 { + let actual = checked_beta_reg(10.0, 1e6, f64::from_bits(bits)).unwrap(); + assert!( + actual >= previous, + "bits={bits:#018x}, previous={previous:?}, actual={actual:?}" + ); + previous = actual; + } +} + +#[test] +fn test_beta_reg_continued_fraction_adjacent_reference() { + let a = 1833.469197457969; + let b = 648975.2550258434; + let cases = [ + (0x3f63feb8f2cd8c97, 0x3e112e0be826bc4b, 0xc034b927f32c0140), + (0x3f63feb8f2cd8c98, 0x3e112e0be826bd23, 0xc034b927f32c0133), + ]; + let mut previous = 0; + for (x, expected, expected_log) in cases { + let x = f64::from_bits(x); + assert_eq!( + checked_ln_beta_reg(a, b, x).unwrap().to_bits(), + expected_log + ); + let actual = beta_reg(a, b, x).to_bits(); + let log_power = beta_reg_log_power_parts(a, b, x); + let fraction = beta_continued_fraction(a, b, x).unwrap(); + let direct = ((log_power.0 + log_power.1).exp() / fraction).to_bits(); + assert!( + actual.abs_diff(expected) <= 2, + "actual={actual:#018x}, direct={direct:#018x}, expected={expected:#018x}" + ); + assert!(actual > previous); + previous = actual; + } +} + +#[test] +fn test_beta_reg_moderate_fraction_against_500_digit_reference() { + let a = 25.32628846940565; + let b = 3.1028101710805442; + let x = 0.9276950604606229; + let expected = 0x3fe69562e02877e6_u64; + let actual = beta_reg(a, b, x).to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "actual={actual:#018x}, expected={expected:#018x}" + ); +} + +#[test] +fn test_beta_reg_typical_fraction_against_500_digit_reference() { + let actual = beta_reg(2.5, 3.5, 0.4).to_bits(); + let expected = 0x3fdf297032b8f5ac_u64; + assert!( + actual.abs_diff(expected) <= 4, + "actual={actual:#018x}, expected={expected:#018x}" + ); +} + +#[test] +fn test_inv_beta_reg_typical_against_500_digit_reference() { + let actual = inv_beta_reg(2.0, 5.0, 0.3).to_bits(); + let expected = 0x3fc745560dce9cd1_u64; + assert!( + actual.abs_diff(expected) <= 2, + "actual={actual:#018x}, expected={expected:#018x}" + ); +} + +#[test] +fn test_inv_beta_reg_shape_two_specialization_boundary() { + let probability = 0.999_999_999_f64; + let probabilities = [ + f64::from_bits(probability.to_bits() - 1), + probability, + f64::from_bits(probability.to_bits() + 1), + ]; + let cases = [ + ( + 2.0, + 200.0, + [0x3fbcd0193e77fee6, 0x3fbcd01940aae67a, 0x3fbcd01942ddce12], + ), + ( + 200.0, + 2.0, + [0x3fefffff883fcdff, 0x3fefffff883fce6e, 0x3fefffff883fcede], + ), + ]; + for (a, b, expected) in cases { + let actual = probabilities.map(|p| inv_beta_reg(a, b, p)); + assert!(actual.windows(2).all(|pair| pair[0] <= pair[1])); + for (value, reference) in actual.into_iter().zip(expected) { + assert!(value.to_bits().abs_diff(reference) <= 1); + } + } +} + +#[test] +fn test_beta_reg_tiny_b_against_500_digit_references() { + let cases = [ + ( + 0.8144818117006096, + 1.250857626649459e-12, + 0.9669920517519052, + 0x3d94af09e6a6b751_u64, + ), + ( + 0.2623971057030866, + 5.23256841817563e-12, + 0.9924817752047999, + 0x3dc7f760fcea90cd, + ), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x).to_bits(); + assert!( + actual.abs_diff(expected) <= 1, + "a={a:?}, b={b:?}, x={x:?}, actual={actual:#018x}, expected={expected:#018x}" + ); + } +} + +#[test] +fn test_beta_reg_tiny_b_boundary_against_500_digit_references() { + let cases = [ + ( + 0.015778004354037867, + 3.91414134306449e-9, + 0.9138081692744422, + 0x3e91430c6cd6e778_u64, + ), + (0.5, 5e-5, 0.95, 0x3f2c8c230a2377e9), + (0.1, 1e-5, 0.99999, 0x3f2bfe1c7f2d26f5), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x).to_bits(); + assert!( + actual.abs_diff(expected) <= 2, + "a={a:?}, b={b:?}, x={x:?}, actual={actual:#018x}, expected={expected:#018x}" + ); + } +} + +#[test] +fn test_beta_reg_tiny_x_large_b_against_reference() { + let cases: [(f64, f64, f64, u64); 2] = [ + (100.0, 1e308, 1.01e-306, 0x3fe1b153914c2fe1_u64), + (1e6, 1e308, 1.000001e-302, 0x3fe0045b85d90000_u64), + ]; + for (a, b, center, expected) in cases { + let center_bits = center.to_bits(); + let mut previous = 0.0; + for bits in center_bits - 128..=center_bits + 128 { + let actual = checked_beta_reg(a, b, f64::from_bits(bits)).unwrap(); + assert!( + actual >= previous, + "a={a}, b={b}, bits={bits:#018x}, previous={previous:?}, actual={actual:?}" + ); + previous = actual; + } + let actual = checked_beta_reg(a, b, center).unwrap().to_bits(); + assert!( + actual.abs_diff(expected) <= 4, + "a={a}, b={b}, actual={actual:#018x}, expected={expected:#018x}" + ); + } +} + +#[test] +fn test_beta_reg_tiny_x_continued_fraction_singularity() { + let references = [ + (0x3c9d1c7c0f1fd2c9_u64, 0x3fe1b153914c2fde_u64), + (0x3c9d1c7c0f1fd2ca_u64, 0x3fe1b153914c2fe2_u64), + (0x3c9d1c7c0f1fd2cb_u64, 0x3fe1b153914c2fe7_u64), + ]; + let mut previous = 0_u64; + for (x, expected) in references { + let actual = checked_beta_reg(100.0, 1e18, f64::from_bits(x)) + .unwrap() + .to_bits(); + assert!( + actual.abs_diff(expected) <= 8, + "x={x:#018x}, actual={actual:#018x}, expected={expected:#018x}" + ); + assert!(actual >= previous); + previous = actual; + } +} + +#[test] +fn test_beta_reg_tiny_x_does_not_lose_complement() { + let a = 40.0; + let b = 1e18; + let center = 0x3c87a28834d566b4_u64; + let mut previous = 0.0; + for bits in center - 128..=center + 128 { + let actual = checked_beta_reg(a, b, f64::from_bits(bits)).unwrap(); + assert!( + actual >= previous, + "bits={bits:#018x}, previous={previous:?}, actual={actual:?}" + ); + previous = actual; + } + let actual = checked_beta_reg(a, b, f64::from_bits(center)).unwrap(); + assert_eq!(actual.to_bits(), 0x3fe2a783c7380c04); +} + +#[test] +fn test_beta_reg_power_series_tiny_shape_boundary() { + let a = f64::from_bits(0x00000000000007e8); + let b = f64::from_bits(0x4040000000000000); + let x = f64::from_bits(0x01556e1fc2f8f359); + assert!(beta_power_series_log_parts(a, b, x).is_ok()); + for offset in -3_i64..=3 { + let x = f64::from_bits(x.to_bits().wrapping_add_signed(offset)); + assert_eq!(checked_beta_reg(a, b, x).unwrap(), 1.0); + assert_eq!( + checked_ln_beta_reg(a, b, x).unwrap().to_bits(), + 0x8000000000155101 + ); + } +} + +#[test] +fn test_beta_reg_power_series_tiny_shape_is_locally_monotone() { + let a = f64::from_bits(0x3d719799812dea11); + let b = f64::from_bits(0x43abc16d674ec800); + let x = f64::from_bits(0x3c32725dd1d243ac); + for offset in -2_i64..=3 { + let x = f64::from_bits(x.to_bits().wrapping_add_signed(offset)); + assert_eq!(beta_reg(a, b, x).to_bits(), 0x3feffffffffff848); + } +} + +#[test] +fn test_accurate_ln_against_multiprecision_reference() { + assert_eq!(accurate_ln(0.0), (f64::NEG_INFINITY, 0.0)); + let cases = [ + (0x0000000000000001, 0xc0874385446d71c3, 0xbd28e569fa8ee781), + (0x0010000000000000, 0xc086232bdd7abcd2, 0xbd1eef3fec1be37f), + (0x39b0000000000000, 0xc051542457337d43, 0x3cde3948c376279d), + (0x3fe8000000000000, 0xbfd269621134db92, 0xbc7e0efadd9db02b), + (0x3ff6a09e667f3bcc, 0x3fd62e42fefa39ee, 0xbc78d6e518e495a3), + (0x3ff6a09e667f3bcd, 0x3fd62e42fefa39f0, 0x3c7c2e0e1b1548c2), + (0x3ff6a09e667f3bce, 0x3fd62e42fefa39f3, 0x3c7133014f0f271f), + (0x3ff8000000000000, 0x3fd9f323ecbf984c, 0xbc4a92e513217f5c), + (0x4000000000000000, 0x3fe62e42fefa39ef, 0x3c7abc9e3b39803f), + (0x4630000000000000, 0x4051542457337d43, 0xbcde3948c376279d), + (0x7fefffffffffffff, 0x40862e42fefa39ef, 0x3d1a9c9e3b39803f), + ]; + for (input, expected_high, expected_low) in cases { + let (high, low) = accurate_ln(f64::from_bits(input)); + let expected_low = f64::from_bits(expected_low); + let magnitude = expected_low.abs(); + let spacing = f64::from_bits(magnitude.to_bits() + 1) - magnitude; + assert_eq!(high.to_bits(), expected_high); + assert!( + (low - expected_low).abs() <= 8.0 * spacing, + "input={input:#018x}, low={low:?}, expected={expected_low:?}" + ); + } +} + +#[test] +fn test_beta_reg_bgrat_lower_shape_boundary() { + let cases = [ + (31.999, 0.5, 0.9, f64::from_bits(0x3f83d8d11db5fecb)), + (32.001, 0.5, 0.9, f64::from_bits(0x3f83d79daec1916d)), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 1e-12, + "beta_reg({a}, {b}, {x}) = {actual}, expected {expected}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_beta_reg_scaled_gamma_boundary_against_reference() { + let cases = [ + (100_000.0, 10.1, 0.9996800497549934, 1.904358612390508e-6), + (100_000.0, 10.1, 0.9996200704814975, 2.132915725768903e-8), + (100_000.0, 10.1, 0.9996000781900461, 4.537230484132134e-9), + (1e8, 0.1, 0.9999993610002013, 4.358202373741317e-31), + (1e8, 0.1, 0.999999360000202, 3.9380016482795125e-31), + (1e8, 0.9, 0.9999928600254862, 3.9763309919351194e-311), + (1e8, 0.9, 0.9999928400256292, 5.37987584721e-312), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-12, + "beta_reg({a}, {b}, {x}) = {actual:e}, expected {expected:e}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_beta_reg_small_shapes_stays_in_range() { + let cases = [ + ( + 0.1350095402068847, + 2.522023373459552e-11, + 0.858047569045879, + 2.2760966295231215e-10, + ), + ( + 1.6182184909371272e-12, + 0.8611154417262772, + 0.2090095742796264, + 0.9999999999971043, + ), + ]; + for (a, b, x, expected) in cases { + let actual = beta_reg(a, b, x); + assert!((0.0..=1.0).contains(&actual)); + assert!( + (actual - expected).abs() <= 5e-15 * expected.max(1e-10), + "beta_reg({a}, {b}, {x}) = {actual}, expected {expected}" + ); + } +} + +#[test] +#[should_panic] +fn test_beta_reg_a_lte_0() { + beta_reg(0.0, 1.0, 1.0); +} + +#[test] +#[should_panic] +fn test_beta_reg_b_lte_0() { + beta_reg(1.0, 0.0, 1.0); +} + +#[test] +#[should_panic] +fn test_beta_reg_x_lt_0() { + beta_reg(1.0, 1.0, -1.0); +} + +#[test] +#[should_panic] +fn test_beta_reg_x_gt_1() { + beta_reg(1.0, 1.0, 2.0); +} + +#[test] +fn test_checked_beta_reg_a_lte_0() { + assert!(checked_beta_reg(0.0, 1.0, 1.0).is_err()); +} + +#[test] +fn test_checked_beta_reg_b_lte_0() { + assert!(checked_beta_reg(1.0, 0.0, 1.0).is_err()); +} + +#[test] +fn test_checked_beta_reg_x_lt_0() { + assert!(checked_beta_reg(1.0, 1.0, -1.0).is_err()); +} + +#[test] +fn test_checked_beta_reg_x_gt_1() { + assert!(checked_beta_reg(1.0, 1.0, 2.0).is_err()); +} + +#[test] +fn test_inv_beta_reg_extreme_probability_does_not_panic() { + assert_eq!( + inv_beta_reg(200.0, 2.0, 1e-165).to_bits(), + 0x3fc2_aa4f_7f31_6421 + ); +} + +#[test] +fn test_inv_beta_reg_extreme_probability_terminates() { + assert_eq!( + inv_beta_reg(200.0, 2.0, 1e-60).to_bits(), + 0x3fdf_5753_caf6_9652 + ); +} + +#[test] +fn test_inv_beta_reg_unit_first_shape_subnormal_integer_second_shape() { + for second_shape in [2_u64, 3, 10, 1000] { + for probability_bits in 1_u64..=1024 { + let quotient = probability_bits / second_shape; + let remainder = probability_bits % second_shape; + let expected = quotient + u64::from(2 * remainder >= second_shape); + assert_eq!( + inv_beta_reg(1.0, second_shape as f64, f64::from_bits(probability_bits),).to_bits(), + expected, + "second_shape={second_shape} probability_bits={probability_bits:#018x}" + ); + } + } +} + +#[test] +fn test_inv_beta_reg_small_shape_lower_tail() { + let cases = [ + (1e-33, 0.0), + (1e-32, f64::from_bits(2)), + (1e-31, 1.215703604971242e-313), + (1e-30, 1.2157036049544172e-303), + (1e-20, 1.2157036049544e-203), + (1e-10, 1.2157036049543856e-103), + (1e-4, 1.2157036049543764e-43), + (1e-2, 1.215703604954373e-23), + ]; + let mut previous = 0.0; + + for (probability, expected) in cases { + let actual = inv_beta_reg(0.1, 500.0, probability); + if expected == 0.0 { + assert_eq!(actual, expected); + continue; + } + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-14, + "inv_beta_reg(0.1, 500, {probability}) = {actual}, expected {expected}, relative error {relative_error}" + ); + assert!(actual >= previous); + previous = actual; + } +} + +#[test] +fn test_inv_beta_reg_small_shape_rounds_extreme_tail() { + let cases = [ + (1e-30, 0x0010aad919ea62cfa), + (1e-31, 0x00000005baa38454), + (1e-32, 0x0000000000000002), + ]; + for (probability, expected) in cases { + assert_eq!(inv_beta_reg(0.1, 500.0, probability).to_bits(), expected); + } +} + +#[test] +fn test_inv_beta_reg_early_tail_correction_against_reference() { + assert_eq!( + inv_beta_reg(10.0, 1e18, f64::from_bits(0x206b45a31ae6c90e),).to_bits(), + 0x392f275e33972f0c + ); +} + +#[test] +fn test_inv_beta_reg_large_a_tiny_b_lower_tail() { + let cases = [ + ( + 27.229198855436444, + 3.192251825919222e-12, + 1e-12, + 0x3fef0fdff94fb881, + ), + ( + 10.741694769633645, + 2.057645959850482e-10, + 5e-9, + 0x3fefffffffffca0f, + ), + ( + 3.791228906881053, + 3.2160853621997853e-9, + 5e-9, + 0x3feeb7bc46a5108f, + ), + ( + 0.07111267420172858, + 2.459402818189203e-11, + 1e-9, + 0x3fefffffffffa790, + ), + ( + 0.0715388852036888, + 3.187243980970482e-9, + 1e-7, + 0x3feffffff2a24e82, + ), + ( + 0.14934701587929067, + 4.564473364066682e-9, + 1e-7, + 0x3fefffff95a64b04, + ), + ]; + for (a, b, probability, expected) in cases { + let actual = inv_beta_reg(a, b, probability).to_bits(); + assert!( + actual.abs_diff(expected) <= 2, + "a={a}, b={b}, actual={actual:#x}, expected={expected:#x}" + ); + } +} + +#[test] +fn test_beta_reg_moderate_a_tiny_b_against_reference() { + let actual = beta_reg( + 6.333131463399467, + 1.3323977213610329e-11, + 0.9137396220685055, + ) + .to_bits(); + let expected = 0x3d9ef22640629504_u64; + assert!( + actual.abs_diff(expected) <= 4, + "actual={actual:#018x}, expected={expected:#018x}" + ); +} + +#[test] +fn test_beta_reg_small_shapes_near_one_against_reference() { + let actual = checked_beta_reg(0.8593272045160161, 0.9835139781033098, 0.9999999999999999) + .unwrap() + .to_bits(); + assert!(actual.abs_diff(0x3feffffffffffffe) <= 1); +} + +#[test] +fn test_ln_beta_accurate_parts_reference() { + let cases = [ + (0.1, 32.0, 0x3ffe85545aa95cd9, 0xbc8fef9442e0fba4), + (0.3, 1000.0, 0xbfef3edcaae7008a, 0xbc8237c135557682), + (10.0, 32.0, 0xc03723e193251f2a, 0xbcd496eeab49e82c), + ]; + for (a, b, high, low) in cases { + let actual = ln_beta_accurate_parts(a, b); + assert_eq!(actual.0.to_bits(), high); + let expected = f64::from_bits(low); + let high_value = f64::from_bits(high).abs(); + let spacing = f64::from_bits(high_value.to_bits() + 1) - high_value; + assert!( + (actual.1 - expected).abs() <= 0.01 * spacing, + "a={a}, b={b}, actual={:?}, expected={expected:?}", + actual.1 + ); + } + let gamma = ln_gamma_accurate_parts(0.1); + assert_eq!(gamma.0.to_bits(), 0x4002058e35f3deee); + assert!( + (gamma.1 - f64::from_bits(0xbc97ad885b23066b)).abs() <= 5e-19, + "gamma={gamma:?}" + ); + let gamma_cases = [ + (0.125, 0x400027c4cfd515b0, 0x3c91baac8949b315), + (0.03125, 0x400b968177c407c6, 0x3ca22ff84657c0bd), + (0.01, 0x401265de0d9b33c4, 0xbca5ae6a9ccd75b9), + (1e-4, 0x40226baa2b2b7f63, 0xbcc49a20e6676b4a), + (1e-8, 0x40326bb1bb9c8a88, 0xbcc0a09c9d154d84), + (1e-12, 0x403ba18a998ffefe, 0xbcc9835734ab358b), + (1e-16, 0x40426bb1bbb55516, 0xbcef9d9398a70762), + (f64::MIN_POSITIVE, 0x4086232bdd7abcd2, 0x3d1eef3fec1be37f), + (f64::from_bits(1), 0x40874385446d71c3, 0x3d28e569fa8ee781), + ]; + for (x, expected_high, expected_low) in gamma_cases { + let expected_low_value = f64::from_bits(expected_low); + let recurrence = ln_gamma_accurate_parts(x); + assert_eq!(recurrence.0.to_bits(), expected_high, "x={x:?}"); + assert!( + (recurrence.1 - expected_low_value).abs() <= 5e-19, + "x={x:?}, recurrence={recurrence:?}, expected_low={expected_low_value:?}" + ); + let series = ln_gamma_small_accurate_parts(x); + assert_eq!(series.0.to_bits(), expected_high, "x={x:?}"); + assert!( + series.1.to_bits().abs_diff(expected_low) <= 8, + "x={x:?}, series={series:?}, expected_low={expected_low_value:?}" + ); + } + let delta = ln_gamma_delta_parts(32.0, 0.1); + assert_eq!(delta.0.to_bits(), 0x3fd6172044f9840c); + assert!((delta.1 - f64::from_bits(0xbc7ed6f8e6ca2265)).abs() <= 5e-19); +} + +#[test] +fn test_inv_beta_reg_regular_shape_lower_tail() { + let cases = [ + (1e-300, 7.053456158585983e-153), + (1e-100, 7.053456158585983e-53), + (1e-40, 7.053456158585983e-23), + (1e-30, 7.053456158585999e-18), + (1e-20, 7.053456158916007e-13), + ]; + let mut previous = 0.0; + + for (probability, expected) in cases { + let actual = inv_beta_reg(2.0, 200.0, probability); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-12, + "inv_beta_reg(2, 200, {probability}) = {actual}, expected {expected}, relative error {relative_error}" + ); + assert!(actual > previous); + previous = actual; + } +} + +#[test] +fn test_inv_beta_reg_large_parameters() { + let cases = [(0.1, 0.3332984541555588), (0.9, 0.3333682129869408)]; + + for (probability, expected) in cases { + let actual = inv_beta_reg(1e8, 2e8, probability); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-12, + "inv_beta_reg(1e8, 2e8, {probability}) = {actual}, expected {expected}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_inv_beta_reg_overflowing_shape_sum() { + for shape in [1e307, 1e308] { + assert_eq!(inv_beta_reg(shape, shape, 0.1), 0.5); + assert_eq!(inv_beta_reg(shape, shape, 0.9), 0.5); + } + let expected = f64::from_bits(0x3fe5555555555555); + for probability in [0.1, 0.5, 0.9] { + assert_eq!(inv_beta_reg(1e308, 5e307, probability), expected); + } +} + +#[test] +fn test_inv_beta_reg_min_subnormal_large_a_tiny_b() { + let cases = [ + ( + 1.418970410722184e16, + 0.0001029663852090984, + f64::from_bits(0x3feffffffffffe31), + ), + ( + 4.674866848491979e16, + 1.8053488701439817e-11, + f64::from_bits(0x3fefffffffffff77), + ), + ( + 3.2111418342313892e16, + 0.004499324538510611, + f64::from_bits(0x3fefffffffffff33), + ), + ( + 3.117388966777583e17, + 0.00105319319351692, + f64::from_bits(0x3fefffffffffffeb), + ), + ( + 9.629243664883278e17, + 3.208469262232818e-5, + f64::from_bits(0x3feffffffffffff9), + ), + ( + 7.351984375091425e17, + 2.6812348495943197e-11, + f64::from_bits(0x3feffffffffffff7), + ), + ( + 1.7012222411445178e17, + 7.129120396546662e-6, + f64::from_bits(0x3fefffffffffffda), + ), + ( + 1.9543788953358486e17, + 1.1304448170316649e-12, + f64::from_bits(0x3fefffffffffffdf), + ), + ( + 9.996829742803416e17, + 1.410942501012109e-8, + f64::from_bits(0x3feffffffffffffa), + ), + ]; + for (a, b, expected) in cases { + let actual = inv_beta_reg(a, b, f64::from_bits(1)); + assert_eq!(actual, expected, "a={a}, b={b}"); + } +} + +#[test] +fn test_inv_beta_reg_small_shape_upper_gamma() { + let cases = [ + ( + 112_176_097_488.593_9, + 1.3959752253898728e-12, + 1e-12, + f64::from_bits(0x3fefffffffff851d), + ), + ( + 238_641_107_383.443_27, + 1.799146819367202e-12, + 1e-12, + f64::from_bits(0x3fefffffffffb5cc), + ), + ( + 246.932962952654, + 1.1953991131275682e-12, + 1e-11, + f64::from_bits(0x3feffffee33a9e66), + ), + ]; + for (a, b, probability, expected) in cases { + assert_eq!(inv_beta_reg(a, b, probability), expected); + } +} + +#[test] +fn test_inv_beta_reg_large_a_tiny_b_is_monotone() { + let cases = [ + (5.034263241208714e17, 1.8917307295846354e-5), + (7.663354755004902e17, 0.06629881964843289), + (9.703110430017175e17, 1.3592520602121614e-6), + (7.633216846220836e17, 0.04203941489807821), + (9.846275348488209e17, 7.919461066109182e-7), + (8.324653375999025e17, 5.050727538603147e-11), + (6.519274800253329e17, 1.3080952792915084e-9), + (9.600975622510844e17, 3.1549066745793863e-7), + (5.0359005294126995e17, 4.282989132250602e-6), + (8.523009112110578e17, 2.1697803811832315e-7), + ]; + for (a, b) in cases { + let lower = inv_beta_reg(a, b, 1e-310); + let upper = inv_beta_reg(a, b, 1e-300); + assert!(lower <= upper, "a={a}, b={b}, lower={lower}, upper={upper}"); + } +} + +#[test] +fn test_inv_beta_reg_log_solver_boundary_is_monotone() { + let probability = 1e-8_f64; + let probabilities = [ + f64::from_bits(probability.to_bits() - 1), + probability, + f64::from_bits(probability.to_bits() + 1), + ]; + let cases = [ + ( + 2.0, + 200.0, + [ + f64::from_bits(0x3ea7ab27fd13660a), + f64::from_bits(0x3ea7ab27fd13660b), + f64::from_bits(0x3ea7ab27fd13660b), + ], + ), + ( + 0.1, + 500.0, + [ + f64::from_bits(0x2eb79df9fcc6b8b8), + f64::from_bits(0x2eb79df9fcc6b8c3), + f64::from_bits(0x2eb79df9fcc6b8ce), + ], + ), + (3.508179849994976e17, 0.8360747930277879, [1.0; 3]), + ]; + for (a, b, expected) in cases { + let actual = probabilities.map(|p| inv_beta_reg(a, b, p)); + assert!( + actual[0] <= actual[1] && actual[1] <= actual[2], + "a={a}, b={b}, actual={actual:?}" + ); + for ((value, reference), probability) in actual.into_iter().zip(expected).zip(probabilities) + { + let ulp_error = value.to_bits().abs_diff(reference.to_bits()); + assert!( + ulp_error <= INVERSE_REFERENCE_MAX_ULPS, + "a={a}, b={b}, probability={probability}, value={value}, reference={reference}, ulp_error={ulp_error}" + ); + let quantile_relative_error = ((value - reference) / reference).abs(); + assert!( + quantile_relative_error <= 4e-14, + "a={a}, b={b}, probability={probability}, value={value}, reference={reference}, quantile_relative_error={quantile_relative_error}" + ); + if value > 0.0 && value < 1.0 { + let relative_error = ((beta_reg(a, b, value) - probability) / probability).abs(); + assert!( + relative_error <= 1e-14, + "a={a}, b={b}, probability={probability}, value={value}, relative_error={relative_error}" + ); + } + } + } +} + +#[test] +fn test_inv_beta_reg_adjacent_probability_is_monotone() { + let probability = 1e-8_f64; + let probabilities = [ + f64::from_bits(probability.to_bits() - 1), + probability, + f64::from_bits(probability.to_bits() + 1), + ]; + let cases = [ + ( + 9.11327743985456, + 133_525_174_076_797.34, + f64::from_bits(0x3cf3d7e149ac36dd), + ), + ( + 6.078046923216118, + 31_131_628_187_944.344, + f64::from_bits(0x3cf592225b607c93), + ), + ]; + for (a, b, expected) in cases { + let actual = probabilities.map(|p| inv_beta_reg(a, b, p)); + assert!( + actual[0] <= actual[1] && actual[1] <= actual[2], + "a={a}, b={b}, actual={actual:?}" + ); + for value in actual { + assert!(value.to_bits().abs_diff(expected.to_bits()) <= INVERSE_REFERENCE_MAX_ULPS); + } + } +} + +#[test] +fn test_inv_beta_reg_upper_adjacent_probability_is_monotone() { + let cases = [ + ( + 100.0, + 1e6, + [0x3feffffffffffff9, 0x3feffffffffffffa], + [0x3f2a6e8528d3e729, 0x3f2a78942066b3b0], + ), + ( + 1000.0, + 1e6, + [0x3feffffffffffff7, 0x3feffffffffffff8], + [0x3f54d1ec0e95e0f5, 0x3f54d42ffc3c17aa], + ), + ( + 1000.0, + 1e6, + [0x3feffffffffffffb, 0x3feffffffffffffc], + [0x3f54dd318598d8ed, 0x3f54e1735a4b5c03], + ), + ( + 1000.0, + 1e6, + [0x3feffffffffffffd, 0x3feffffffffffffe], + [0x3f54e6ebec74e0ca, 0x3f54ee997db90e85], + ), + ( + 1000.0, + 1e8, + [0x3feffffffffffff3, 0x3feffffffffffff4], + [0x3eeaa4df95604c33, 0x3eeaa6db7106f8eb], + ), + ]; + for (a, b, probability_bits, expected_bits) in cases { + let actual = probability_bits.map(|bits| inv_beta_reg(a, b, f64::from_bits(bits))); + assert!(actual[0] <= actual[1]); + for (value, expected) in actual.into_iter().zip(expected_bits.map(f64::from_bits)) { + let ulp_error = value.to_bits().abs_diff(expected.to_bits()); + assert!( + ulp_error <= INVERSE_REFERENCE_MAX_ULPS, + "a={a}, b={b}, value={value}, expected={expected}, ulp_error={ulp_error}" + ); + } + } +} + +#[test] +fn test_inv_beta_reg_orientation_preserves_tiny_quantiles() { + let cases = [ + (0.49, f64::from_bits(0x083429b7deb4de35)), + (0.5, f64::from_bits(0x0a0650cbd0bac729)), + (0.51, f64::from_bits(0x0bd08de62d4b3d17)), + (0.9, f64::from_bits(0x3f064452047719b0)), + (0.99, 1.0), + ]; + let mut previous = 0.0; + for (probability, expected) in cases { + let actual = inv_beta_reg(0.001, 0.01, probability); + assert!(actual >= previous); + if expected == 1.0 { + assert_eq!(actual, expected); + } else { + assert!(((actual - expected) / expected).abs() <= 1e-12); + } + previous = actual; + } + let actual = inv_beta_reg(0.01, 1e8, 0.51); + let expected = f64::from_bits(0x38260460ad60f7d3); + assert!(((actual - expected) / expected).abs() <= 1e-12); +} + +#[test] +fn test_inv_beta_reg_concentrated_quantiles_round_correctly() { + let cases = [ + ( + 5.6337457945398355e35, + 3.4148653071385907e36, + 0.1, + f64::from_bits(0x3fc2206894075924), + ), + ( + 5.6337457945398355e35, + 3.4148653071385907e36, + 0.9, + f64::from_bits(0x3fc2206894075924), + ), + ( + 7.778370008599511e35, + 3.99094171205976e36, + f64::from_bits(1), + f64::from_bits(0x3fc4e0cc7f8ea39f), + ), + ( + 7.778370008599511e35, + 3.99094171205976e36, + 0.1, + f64::from_bits(0x3fc4e0cc7f8ea3a0), + ), + ( + 7.778370008599511e35, + 3.99094171205976e36, + 0.9, + f64::from_bits(0x3fc4e0cc7f8ea3a0), + ), + ]; + for (a, b, probability, expected) in cases { + assert_eq!(inv_beta_reg(a, b, probability), expected); + } +} + +#[test] +fn test_inv_beta_reg_concentrated_extreme_ratios_round_correctly() { + let cases = [ + ( + 0x6260d489122d216a, + 0x40791cf3e3defb3a, + 0x0a87f70771690bbe, + 0x3ff0000000000000, + ), + ( + 0x6d8891e7fd8b2193, + 0x6a40bb8ae6715d29, + 0x245f1e03bdf9d74e, + 0x3fefffffffffffff, + ), + ]; + for (a, b, probability, expected) in cases { + assert_eq!( + inv_beta_reg( + f64::from_bits(a), + f64::from_bits(b), + f64::from_bits(probability), + ) + .to_bits(), + expected, + ); + } +} + +#[test] +fn test_inv_beta_reg_concentrated_gate_preserves_extreme_tail_rounding() { + let probability = f64::from_bits(1); + let cases = [ + (0x43bbc16d674ec800, 0x3feffffffffffffd), + (0x43e158e460913d00, 0x3fefffffffffffff), + ]; + for (a, expected) in cases { + assert_eq!( + inv_beta_reg(f64::from_bits(a), 2.0, probability).to_bits(), + expected, + ); + } + assert_eq!(inv_beta_reg(1e308, 2.0, 0.5), 1.0); +} + +#[test] +fn test_inv_beta_reg_concentrated_unit_endpoint_is_monotone() { + for shape in [1e100, 1e200, 1e308] { + let mut previous = 0.0; + for probability in [ + 0.5, + 0.9, + 0.99, + 1.0 - 1e-12, + f64::from_bits(1.0_f64.to_bits() - 1), + ] { + let actual = inv_beta_reg(shape, 2.0, probability); + assert!(actual >= previous); + assert_eq!(actual, 1.0); + previous = actual; + } + let lower = inv_beta_reg(shape, 2.0, f64::from_bits(1)); + assert!(lower <= inv_beta_reg(shape, 2.0, 0.5)); + let mut swapped_previous = 0.0; + for probability in [f64::from_bits(1), 0.1, 0.5, 0.9, 1.0 - 1e-12] { + let actual = inv_beta_reg(2.0, shape, probability); + assert!(actual >= swapped_previous && actual < 1.0); + swapped_previous = actual; + } + } +} + +#[test] +fn test_inv_beta_reg_concentrated_endpoint_uses_exact_mean() { + let actual = inv_beta_reg( + f64::from_bits(0x44ea_7843_79d9_9db4), + f64::from_bits(0x4193_da32_9b63_3647), + 0.5, + ); + assert_eq!(actual.to_bits(), 0x3fef_ffff_ffff_ffff); +} + +#[test] +fn test_inv_beta_reg_extreme_tail_balanced_shapes() { + let cases = [ + (f64::from_bits(1), 0.1384383837250825), + (1e-300, 0.14764444133469024), + ]; + for (probability, expected) in cases { + let actual = inv_beta_reg(1000.0, 1000.0, probability); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-13, + "probability {probability}, actual {actual}, expected {expected}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_inv_beta_reg_extreme_tail_imbalanced_shapes() { + let cases = [ + (200.0, 2.0, 1e-192, 0.10683857283574616), + (1000.0, 2.0, f64::from_bits(1), 0.47203081850113066), + (1000.0, 2.0, 1e-303, 0.49464719057284383), + (1000.0, 2.0, 1e-200, 0.627230829476228), + (1000.0, 2.0, 1e-100, 0.7900887907081466), + (1000.0, 10.0, f64::from_bits(1), 0.454569346824437), + (1000.0, 10.0, 1e-303, 0.47650393899531424), + (1000.0, 10.0, 1e-200, 0.6055787273511661), + (1000.0, 10.0, 1e-100, 0.7659557362087095), + (1000.0, 100.0, f64::from_bits(1), 0.356892489498544), + (1000.0, 100.0, 1e-303, 0.3750351205470552), + (1000.0, 100.0, 1e-200, 0.48455098775995836), + (1000.0, 100.0, 1e-100, 0.6303764215497716), + (7_627_209.761, 11.3319, 1.679e-274, 0.9999105965110135), + ]; + for (a, b, probability, expected) in cases { + let actual = inv_beta_reg(a, b, probability); + let relative_error = ((actual - expected) / expected).abs(); + assert!( + relative_error <= 5e-13, + "inv_beta_reg({a}, {b}, {probability}) = {actual}, expected {expected}, relative error {relative_error}" + ); + } +} + +#[test] +fn test_inv_beta_reg_subnormal_power_series_boundary() { + let a = f64::from_bits(0x4024000000000000); + let b = f64::from_bits(0x7e37e43c8800759c); + let probability = f64::from_bits(0x2df5ed8667733d64); + for offset in -2_i64..=2 { + let probability = f64::from_bits(probability.to_bits().wrapping_add_signed(offset)); + assert_eq!( + inv_beta_reg(a, b, probability).to_bits(), + 0x000730d67819e860, + "offset={offset}" + ); + } +} + +#[test] +fn test_error_is_sync_send() { + fn assert_sync_send() {} + assert_sync_send::(); +}