From 97360ef178469f1201350067afb1b207f102b742 Mon Sep 17 00:00:00 2001 From: Greg von Nessi Date: Fri, 13 Mar 2026 16:39:03 +0000 Subject: [PATCH 1/2] Defensive input validation, numerical robustness, and golub_welsch Result API (v0.2.0) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Address 12 deferred codebase review items across 5 phases: - Reject ±Inf (not just NaN) in tanh-sinh, oscillatory, Cauchy PV, cubature - Replace partial_cmp().unwrap() with unwrap_or(Ordering::Equal) - Promote debug_assert to assert in CubatureRule::new - Split Lobatto n<2 error into ZeroOrder (n=0) vs InvalidInput (n=1) - Clamp global_error to non-negative in adaptive cubature - Add dimension cap (d≤30) for adaptive cubature - Clamp Newton iterates in Lobatto interior node computation - Fix tanh-sinh non-convergence error: use last two estimates instead of fabricated tol*10 - Fix ln_gamma reflection: .sin().ln() → .sin().abs().ln() - Document QUADPACK error simplification in gauss_kronrod - Make golub_welsch return Result, propagate through all compute_* functions - Bump version to 0.2.0, update CHANGELOG, SECURITY, README --- CHANGELOG.md | 25 ++++++++++++- Cargo.toml | 2 +- README.md | 6 ++-- SECURITY.md | 3 +- src/adaptive.rs | 36 +++++++++++++++++++ src/cauchy_pv.rs | 27 +++++++++++--- src/cubature/adaptive.rs | 42 ++++++++++++++++------ src/cubature/mod.rs | 2 +- src/cubature/monte_carlo.rs | 58 +++++++++++++++++++----------- src/gauss_hermite.rs | 4 +-- src/gauss_jacobi.rs | 29 +++++++++++---- src/gauss_kronrod.rs | 7 +++- src/gauss_laguerre.rs | 16 ++++++--- src/gauss_lobatto.rs | 25 ++++++++++++- src/gauss_radau.rs | 8 ++--- src/golub_welsch.rs | 46 ++++++++++++++++++------ src/oscillatory.rs | 15 +++++--- src/tanh_sinh.rs | 38 ++++++++++++++++---- src/weighted.rs | 70 +++++++++++++++++++++++++++++++++---- 19 files changed, 369 insertions(+), 90 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 84d81bb..010a173 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,30 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.1.0] - Unreleased +## [0.2.0] - 2026-03-13 + +### Changed + +- **Breaking:** `golub_welsch()` now returns `Result`, propagating QL non-convergence as `QuadratureError::InvalidInput` instead of silently returning inaccurate nodes/weights. All internal `compute_*` functions in `gauss_jacobi`, `gauss_laguerre`, `gauss_hermite`, and `gauss_radau` propagate accordingly. Public constructors already returned `Result`, so most callers are unaffected. +- **Breaking:** `GaussLobatto::new(0)` now returns `QuadratureError::ZeroOrder` (previously `InvalidInput`). `GaussLobatto::new(1)` still returns `InvalidInput`. +- Input validation for `tanh_sinh`, `oscillatory`, `cauchy_pv`, and `cubature::adaptive` now rejects `±Inf` bounds (not just `NaN`). Error variant is unchanged (`DegenerateInterval`). +- `CubatureRule::new` assertion promoted from `debug_assert_eq!` to `assert_eq!`. +- Tanh-sinh non-convergence error estimate now uses the difference between the last two level estimates instead of a fabricated `tol * 10` value. +- QUADPACK error estimation heuristic in `gauss_kronrod` documented as an intentional simplification of the full formula. + +### Fixed + +- `ln_gamma` reflection formula: `.sin().ln()` → `.sin().abs().ln()` prevents `NaN` for negative arguments where `sin(π·x)` is negative (affects Gauss-Jacobi with certain α, β near negative integers). +- Newton iteration in Gauss-Lobatto interior node computation now clamps iterates to `(-1+ε, 1-ε)`, preventing division by zero in the `P''` formula when an iterate lands on ±1. +- `partial_cmp().unwrap()` in `golub_welsch` and `gauss_lobatto` sort replaced with `.unwrap_or(Ordering::Equal)` to avoid panics on NaN nodes. +- Adaptive cubature `global_error` subtraction clamped to `max(0.0)` to prevent negative error estimates from floating-point cancellation. + +### Added + +- Dimension cap (d ≤ 30) for adaptive cubature — returns `InvalidInput` for higher dimensions where Genz-Malik's `2^d` vertex evaluations would be impractical. +- New tests: infinity rejection (tanh-sinh, oscillatory, Cauchy PV, cubature), dimension cap, Lobatto error variant splitting, tanh-sinh non-fabricated error, Jacobi negative α/β exercising the reflection path. + +## [0.1.0] - 2026-03-12 Initial release of bilby. diff --git a/Cargo.toml b/Cargo.toml index de7a897..4cfd220 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bilby" -version = "0.1.0" +version = "0.2.0" edition = "2021" rust-version = "1.93" description = "A high-performance numerical quadrature (integration) library for Rust" diff --git a/README.md b/README.md index 4dcac00..fa84928 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ Add to `Cargo.toml`: ```toml [dependencies] -bilby = "0.1" +bilby = "0.2" ``` ```rust @@ -49,7 +49,7 @@ bilby works in `no_std` environments (with `alloc`). Disable default features: ```toml [dependencies] -bilby = { version = "0.1", default-features = false } +bilby = { version = "0.2", default-features = false } ``` ### Parallelism @@ -58,7 +58,7 @@ Enable the `parallel` feature for parallel variants of integration methods: ```toml [dependencies] -bilby = { version = "0.1", features = ["parallel"] } +bilby = { version = "0.2", features = ["parallel"] } ``` This provides `integrate_composite_par`, `integrate_par`, `integrate_box_par`, diff --git a/SECURITY.md b/SECURITY.md index 5b5d2aa..55d3d61 100644 --- a/SECURITY.md +++ b/SECURITY.md @@ -4,7 +4,8 @@ | Version | Supported | |---------|-----------| -| 0.1.x | Yes | +| 0.2.x | Yes | +| < 0.2 | No | Only the latest minor release receives security updates. diff --git a/src/adaptive.rs b/src/adaptive.rs index f0f7515..fee8a5a 100644 --- a/src/adaptive.rs +++ b/src/adaptive.rs @@ -180,6 +180,13 @@ impl AdaptiveIntegrator { where G: Fn(f64) -> f64, { + // Validate interval bounds + for &(a, b) in intervals { + if a.is_nan() || b.is_nan() || a.is_infinite() || b.is_infinite() { + return Err(QuadratureError::DegenerateInterval); + } + } + // Validate tolerances if self.abs_tol <= 0.0 && self.rel_tol <= 0.0 { return Err(QuadratureError::InvalidInput( @@ -589,4 +596,33 @@ mod tests { reverse.value ); } + + /// NaN bounds are rejected. + #[test] + fn nan_bounds_rejected() { + assert!(adaptive_integrate(f64::sin, f64::NAN, 1.0, 1e-10).is_err()); + assert!(adaptive_integrate(f64::sin, 0.0, f64::NAN, 1e-10).is_err()); + } + + /// Infinite bounds are rejected. + #[test] + fn infinite_bounds_rejected() { + assert!(adaptive_integrate(f64::sin, f64::INFINITY, 1.0, 1e-10).is_err()); + assert!(adaptive_integrate(f64::sin, 0.0, f64::NEG_INFINITY, 1e-10).is_err()); + } + + /// NaN break point is rejected. + #[test] + fn nan_break_point_rejected() { + assert!(adaptive_integrate_with_breaks(f64::sin, 0.0, 1.0, &[f64::NAN], 1e-10).is_err()); + } + + /// max_evals smaller than one GK panel is rejected. + #[test] + fn max_evals_too_small() { + let r = AdaptiveIntegrator::default() + .with_max_evals(1) + .integrate(0.0, 1.0, f64::sin); + assert!(r.is_err()); + } } diff --git a/src/cauchy_pv.rs b/src/cauchy_pv.rs index 4c3f712..10c4213 100644 --- a/src/cauchy_pv.rs +++ b/src/cauchy_pv.rs @@ -74,7 +74,7 @@ impl CauchyPV { /// /// # Errors /// - /// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN. + /// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is non-finite. /// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside /// `(a, b)` or if `f(c)` is not finite. #[allow(clippy::many_single_char_names)] // a, b, c, f, g are conventional in quadrature @@ -89,7 +89,7 @@ impl CauchyPV { G: Fn(f64) -> f64, { // Validate inputs - if a.is_nan() || b.is_nan() || c.is_nan() { + if !a.is_finite() || !b.is_finite() || !c.is_finite() { return Err(QuadratureError::DegenerateInterval); } let (lo, hi) = if a < b { (a, b) } else { (b, a) }; @@ -112,7 +112,7 @@ impl CauchyPV { // The subtracted integrand: g(x) = [f(x) - f(c)] / (x - c) // This has a removable singularity at x = c. - let guard = 1e-15 * (b - a); + let guard = 1e-15 * (b - a).abs(); let g = |x: f64| -> f64 { if (x - c).abs() < guard { 0.0 @@ -152,7 +152,7 @@ impl CauchyPV { /// /// # Errors /// -/// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN. +/// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is non-finite. /// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside /// `(a, b)` or if `f(c)` is not finite. pub fn pv_integrate( @@ -236,9 +236,28 @@ mod tests { assert!(pv_integrate(|x| x, 0.0, 1.0, 1.0, 1e-10).is_err()); } + #[test] + fn reversed_bounds() { + // PV ∫₁⁰ x²/(x - 0.3) dx = -[∫₀¹ x²/(x - 0.3) dx] + let exact = -(0.8 + 0.09 * (7.0_f64 / 3.0).ln()); + let result = pv_integrate(|x| x * x, 1.0, 0.0, 0.3, 1e-10).unwrap(); + assert!( + (result.value - exact).abs() < 1e-7, + "value={}, expected={exact}", + result.value + ); + } + #[test] fn nan_inputs() { assert!(pv_integrate(|x| x, f64::NAN, 1.0, 0.5, 1e-10).is_err()); assert!(pv_integrate(|x| x, 0.0, 1.0, f64::NAN, 1e-10).is_err()); } + + #[test] + fn inf_inputs_rejected() { + assert!(pv_integrate(|x| x, f64::INFINITY, 1.0, 0.5, 1e-10).is_err()); + assert!(pv_integrate(|x| x, 0.0, f64::NEG_INFINITY, 0.5, 1e-10).is_err()); + assert!(pv_integrate(|x| x, 0.0, 1.0, f64::INFINITY, 1e-10).is_err()); + } } diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs index 32c9809..55faec2 100644 --- a/src/cubature/adaptive.rs +++ b/src/cubature/adaptive.rs @@ -79,7 +79,7 @@ impl AdaptiveCubature { /// /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have /// different lengths or are empty. Returns [`QuadratureError::DegenerateInterval`] - /// if any bound is NaN. + /// if any bound is non-finite. pub fn integrate( &self, lower: &[f64], @@ -95,20 +95,24 @@ impl AdaptiveCubature { "lower and upper must have equal nonzero length", )); } + if d > 30 { + return Err(QuadratureError::InvalidInput( + "adaptive cubature is limited to at most 30 dimensions", + )); + } for i in 0..d { - if lower[i].is_nan() || upper[i].is_nan() { + if !lower[i].is_finite() || !upper[i].is_finite() { return Err(QuadratureError::DegenerateInterval); } } if d == 1 { - // Delegate to 1D adaptive integrator - return crate::adaptive::adaptive_integrate( - |x: f64| f(&[x]), - lower[0], - upper[0], - self.abs_tol, - ); + // Delegate to 1D adaptive integrator, forwarding all settings + return crate::adaptive::AdaptiveIntegrator::default() + .with_abs_tol(self.abs_tol) + .with_rel_tol(self.rel_tol) + .with_max_evals(self.max_evals) + .integrate(lower[0], upper[0], |x: f64| f(&[x])); } let mut heap: BinaryHeap = BinaryHeap::new(); @@ -135,7 +139,7 @@ impl AdaptiveCubature { let Some(worst) = heap.pop() else { break }; global_estimate -= worst.estimate; - global_error -= worst.error; + global_error = (global_error - worst.error).max(0.0); // Bisect along the split axis let axis = worst.split_axis; @@ -197,7 +201,7 @@ impl AdaptiveCubature { /// /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have /// different lengths or are empty. Returns [`QuadratureError::DegenerateInterval`] -/// if any bound is NaN. +/// if any bound is non-finite. pub fn adaptive_cubature( f: G, lower: &[f64], @@ -400,6 +404,22 @@ mod tests { assert!(adaptive_cubature(|_| 1.0, &[0.0], &[1.0, 2.0], 1e-8).is_err()); } + #[test] + fn inf_bounds_rejected() { + assert!(adaptive_cubature(|_| 1.0, &[f64::INFINITY], &[1.0], 1e-8).is_err()); + assert!(adaptive_cubature(|_| 1.0, &[0.0], &[f64::NEG_INFINITY], 1e-8).is_err()); + assert!( + adaptive_cubature(|_| 1.0, &[0.0, f64::INFINITY], &[1.0, 1.0], 1e-8).is_err() + ); + } + + #[test] + fn dimension_cap() { + let lower = vec![0.0; 31]; + let upper = vec![1.0; 31]; + assert!(adaptive_cubature(|_| 1.0, &lower, &upper, 1e-8).is_err()); + } + /// Constant function over [0,1]^2. #[test] fn constant_2d() { diff --git a/src/cubature/mod.rs b/src/cubature/mod.rs index d08b7f6..5f6fd28 100644 --- a/src/cubature/mod.rs +++ b/src/cubature/mod.rs @@ -64,7 +64,7 @@ impl CubatureRule { /// `nodes_flat` must have length `weights.len() * dim`. #[must_use] pub fn new(nodes_flat: Vec, weights: Vec, dim: usize) -> Self { - debug_assert_eq!(nodes_flat.len(), weights.len() * dim); + assert_eq!(nodes_flat.len(), weights.len() * dim); Self { nodes: nodes_flat, weights, diff --git a/src/cubature/monte_carlo.rs b/src/cubature/monte_carlo.rs index 09ef58e..df53b41 100644 --- a/src/cubature/monte_carlo.rs +++ b/src/cubature/monte_carlo.rs @@ -223,22 +223,25 @@ impl MonteCarloIntegrator { let estimate = volume * sum / n as f64; // Heuristic error estimate: compare N/2 and N estimates - let half_sum = { - let mut seq2 = S::new(d)?; - let half = n / 2; - let mut sm = 0.0; - for _ in 0..half { - seq2.next_point(&mut u); - for j in 0..d { - x[j] = lower[j] + widths[j] * u[j]; + let half = n / 2; + let error = if half == 0 { + f64::INFINITY + } else { + let half_sum = { + let mut seq2 = S::new(d)?; + let mut sm = 0.0; + for _ in 0..half { + seq2.next_point(&mut u); + for j in 0..d { + x[j] = lower[j] + widths[j] * u[j]; + } + sm += f(&x); } - sm += f(&x); - } - volume * sm / half as f64 + volume * sm / half as f64 + }; + (estimate - half_sum).abs() }; - let error = (estimate - half_sum).abs(); - Ok(QuadratureResult { value: estimate, error_estimate: error, @@ -352,13 +355,22 @@ impl MonteCarloIntegrator { }) .collect(); - // Merge chunk results + // Merge chunk results using parallel Welford formula let mut total_sum = 0.0; let mut total_m2 = 0.0; let mut total_n = 0usize; for (sum, m2, count) in chunk_results { + if total_n > 0 && count > 0 { + let mean_a = total_sum / total_n as f64; + let mean_b = sum / count as f64; + let delta = mean_b - mean_a; + total_m2 += m2 + + delta * delta * (total_n as f64 * count as f64) + / (total_n + count) as f64; + } else { + total_m2 += m2; + } total_sum += sum; - total_m2 += m2; total_n += count; } @@ -420,12 +432,16 @@ impl MonteCarloIntegrator { // Heuristic error: compare N/2 and N estimates let half = n / 2; - let half_sum: f64 = (0..half) - .into_par_iter() - .map(|i| f(&points[i * d..(i + 1) * d])) - .sum(); - let half_estimate = volume * half_sum / half as f64; - let error = (estimate - half_estimate).abs(); + let error = if half == 0 { + f64::INFINITY + } else { + let half_sum: f64 = (0..half) + .into_par_iter() + .map(|i| f(&points[i * d..(i + 1) * d])) + .sum(); + let half_estimate = volume * half_sum / half as f64; + (estimate - half_estimate).abs() + }; Ok(QuadratureResult { value: estimate, diff --git a/src/gauss_hermite.rs b/src/gauss_hermite.rs index 8525c7b..ae68ef3 100644 --- a/src/gauss_hermite.rs +++ b/src/gauss_hermite.rs @@ -44,7 +44,7 @@ impl GaussHermite { return Err(QuadratureError::ZeroOrder); } - let (nodes, weights) = compute_hermite(n); + let (nodes, weights) = compute_hermite(n)?; Ok(Self { rule: QuadratureRule { nodes, weights }, }) @@ -61,7 +61,7 @@ impl_rule_accessors!(GaussHermite, nodes_doc: "Returns the nodes on (-∞, ∞). /// Jacobi matrix: diagonal = 0, off-diagonal² = (k+1)/2 for k=0..n-2. /// μ₀ = √π. #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64 -fn compute_hermite(n: usize) -> (Vec, Vec) { +fn compute_hermite(n: usize) -> Result<(Vec, Vec), QuadratureError> { let diag = vec![0.0; n]; let off_diag_sq: Vec = (1..n).map(|k| k as f64 / 2.0).collect(); let mu0 = core::f64::consts::PI.sqrt(); diff --git a/src/gauss_jacobi.rs b/src/gauss_jacobi.rs index 1e53132..4aaae56 100644 --- a/src/gauss_jacobi.rs +++ b/src/gauss_jacobi.rs @@ -53,18 +53,18 @@ impl GaussJacobi { /// /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero. /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1`, `beta <= -1`, - /// or either parameter is NaN. + /// or either parameter is non-finite. pub fn new(n: usize, alpha: f64, beta: f64) -> Result { if n == 0 { return Err(QuadratureError::ZeroOrder); } - if alpha <= -1.0 || beta <= -1.0 || alpha.is_nan() || beta.is_nan() { + if !alpha.is_finite() || !beta.is_finite() || alpha <= -1.0 || beta <= -1.0 { return Err(QuadratureError::InvalidInput( - "Jacobi parameters must satisfy alpha > -1 and beta > -1", + "Jacobi parameters must be finite and satisfy alpha > -1 and beta > -1", )); } - let (nodes, weights) = compute_jacobi(n, alpha, beta); + let (nodes, weights) = compute_jacobi(n, alpha, beta)?; Ok(Self { rule: QuadratureRule { nodes, weights }, alpha, @@ -98,7 +98,7 @@ impl_rule_accessors!(GaussJacobi, nodes_doc: "Returns the nodes on \\[-1, 1\\]." /// /// μ₀ = 2^(α+β+1) Γ(α+1)Γ(β+1) / Γ(α+β+2) #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64 -fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> (Vec, Vec) { +fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> Result<(Vec, Vec), QuadratureError> { let ab = alpha + beta; // Diagonal: α_k = (β²-α²) / ((2k+ab)(2k+ab+2)) @@ -170,7 +170,7 @@ pub(crate) fn ln_gamma(x: f64) -> f64 { if x < 0.5 { // Reflection formula let z = 1.0 - x; - core::f64::consts::PI.ln() - (core::f64::consts::PI * x).sin().ln() - ln_gamma(z) + core::f64::consts::PI.ln() - (core::f64::consts::PI * x).sin().abs().ln() - ln_gamma(z) } else { let z = x - 1.0; let mut sum = COEFF[0]; @@ -197,6 +197,8 @@ mod tests { assert!(GaussJacobi::new(5, -1.0, 0.0).is_err()); assert!(GaussJacobi::new(5, 0.0, -1.5).is_err()); assert!(GaussJacobi::new(5, f64::NAN, 0.0).is_err()); + assert!(GaussJacobi::new(5, f64::INFINITY, 0.0).is_err()); + assert!(GaussJacobi::new(5, 0.0, f64::NEG_INFINITY).is_err()); } /// alpha=beta=0 should recover Gauss-Legendre. @@ -317,4 +319,19 @@ mod tests { let sum: f64 = gj.weights().iter().sum(); assert!((sum - PI).abs() < 1e-10, "sum={sum}, expected={PI}"); } + + /// Exercises the ln_gamma reflection formula with negative alpha and beta. + /// Weight sum should equal B(0.2, 0.2) * 2^{-0.6} = Gamma(0.2)^2 / Gamma(0.4) * 2^{-0.6} + #[test] + fn jacobi_negative_alpha_beta() { + let alpha = -0.8; + let beta = -0.8; + let gj = GaussJacobi::new(20, alpha, beta).unwrap(); + let sum: f64 = gj.weights().iter().sum(); + let expected = jacobi_integral(alpha, beta); + assert!( + (sum - expected).abs() < 1e-8, + "alpha={alpha}, beta={beta}: sum={sum}, expected={expected}" + ); + } } diff --git a/src/gauss_kronrod.rs b/src/gauss_kronrod.rs index 8ab8018..37b64da 100644 --- a/src/gauss_kronrod.rs +++ b/src/gauss_kronrod.rs @@ -205,7 +205,12 @@ impl GaussKronrod { let gauss_result = half_width * gauss_sum; abs_sum *= half_width.abs(); - // QUADPACK error estimation heuristic + // QUADPACK error estimation heuristic. + // This is an intentional simplification of the full QUADPACK error + // formula (which also considers |f - I/(b-a)| for roundoff detection). + // The simplified version omits the roundoff term but retains the + // core (200·δ/S)^1.5 scaling that prevents error underestimation + // for smooth integrands. See Piessens et al. (1983), §2.2. let mut error = (estimate - gauss_result).abs(); if abs_sum > 0.0 && error > 0.0 { error = abs_sum * (200.0 * error / abs_sum).min(1.0).powf(1.5); diff --git a/src/gauss_laguerre.rs b/src/gauss_laguerre.rs index 58913e1..409620e 100644 --- a/src/gauss_laguerre.rs +++ b/src/gauss_laguerre.rs @@ -46,18 +46,18 @@ impl GaussLaguerre { /// /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero. /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1` or `alpha` - /// is NaN. + /// is non-finite. pub fn new(n: usize, alpha: f64) -> Result { if n == 0 { return Err(QuadratureError::ZeroOrder); } - if alpha <= -1.0 || alpha.is_nan() { + if !alpha.is_finite() || alpha <= -1.0 { return Err(QuadratureError::InvalidInput( - "Laguerre parameter alpha must satisfy alpha > -1", + "Laguerre parameter alpha must be finite and satisfy alpha > -1", )); } - let (nodes, weights) = compute_laguerre(n, alpha); + let (nodes, weights) = compute_laguerre(n, alpha)?; Ok(Self { rule: QuadratureRule { nodes, weights }, alpha, @@ -81,7 +81,7 @@ impl_rule_accessors!(GaussLaguerre, nodes_doc: "Returns the nodes on \\[0, ∞). /// Jacobi matrix: diagonal = 2k+1+α, off-diagonal² = (k+1)(k+1+α). /// μ₀ = Γ(α+1). #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64 -fn compute_laguerre(n: usize, alpha: f64) -> (Vec, Vec) { +fn compute_laguerre(n: usize, alpha: f64) -> Result<(Vec, Vec), QuadratureError> { let diag: Vec = (0..n).map(|k| 2.0 * k as f64 + 1.0 + alpha).collect(); let off_diag_sq: Vec = (1..n) .map(|k| { @@ -144,6 +144,12 @@ mod tests { } } + #[test] + fn infinite_alpha_rejected() { + assert!(GaussLaguerre::new(5, f64::INFINITY).is_err()); + assert!(GaussLaguerre::new(5, f64::NAN).is_err()); + } + /// Polynomial exactness: integral of x^k e^(-x) over [0,inf) = k! = Gamma(k+1). #[test] fn polynomial_exactness() { diff --git a/src/gauss_lobatto.rs b/src/gauss_lobatto.rs index 4f15645..931afbf 100644 --- a/src/gauss_lobatto.rs +++ b/src/gauss_lobatto.rs @@ -45,6 +45,9 @@ impl GaussLobatto { /// /// Returns [`QuadratureError::InvalidInput`] if `n` is less than 2. pub fn new(n: usize) -> Result { + if n == 0 { + return Err(QuadratureError::ZeroOrder); + } if n < 2 { return Err(QuadratureError::InvalidInput( "Gauss-Lobatto requires at least 2 points", @@ -119,6 +122,7 @@ fn compute_lobatto(n: usize) -> (Vec, Vec) { } let dx = dp_nm1 / d2p; x -= dx; + x = x.clamp(-1.0 + f64::EPSILON, 1.0 - f64::EPSILON); if dx.abs() < 1e-15 * (1.0 + x.abs()) { break; } @@ -133,7 +137,10 @@ fn compute_lobatto(n: usize) -> (Vec, Vec) { // Sort ascending (should already be, but ensure) let mut pairs: Vec<_> = nodes.into_iter().zip(weights).collect(); - pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + pairs.sort_by(|a, b| { + a.0.partial_cmp(&b.0) + .unwrap_or(core::cmp::Ordering::Equal) + }); let (nodes, weights) = pairs.into_iter().unzip(); (nodes, weights) @@ -150,6 +157,22 @@ mod tests { assert!(GaussLobatto::new(2).is_ok()); } + #[test] + fn zero_order_returns_zero_order_error() { + assert_eq!( + GaussLobatto::new(0).unwrap_err(), + QuadratureError::ZeroOrder + ); + } + + #[test] + fn one_point_returns_invalid_input() { + assert!(matches!( + GaussLobatto::new(1).unwrap_err(), + QuadratureError::InvalidInput(_) + )); + } + #[test] fn two_point_trapezoid() { let gl = GaussLobatto::new(2).unwrap(); diff --git a/src/gauss_radau.rs b/src/gauss_radau.rs index d89db18..b16100a 100644 --- a/src/gauss_radau.rs +++ b/src/gauss_radau.rs @@ -48,7 +48,7 @@ impl GaussRadau { if n == 0 { return Err(QuadratureError::ZeroOrder); } - let (nodes, weights) = compute_radau_left(n); + let (nodes, weights) = compute_radau_left(n)?; Ok(Self { rule: QuadratureRule { nodes, weights }, }) @@ -66,7 +66,7 @@ impl GaussRadau { return Err(QuadratureError::ZeroOrder); } // Right Radau is the reflection of left Radau - let (nodes_left, weights_left) = compute_radau_left(n); + let (nodes_left, weights_left) = compute_radau_left(n)?; let nodes: Vec = nodes_left.iter().rev().map(|&x| -x).collect(); let weights: Vec = weights_left.iter().rev().copied().collect(); Ok(Self { @@ -83,9 +83,9 @@ impl_rule_accessors!(GaussRadau, nodes_doc: "Returns the nodes on \\[-1, 1\\].") /// Legendre Jacobi matrix. The last diagonal element is modified so /// that -1 is an eigenvalue of the tridiagonal matrix. #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64 -fn compute_radau_left(n: usize) -> (Vec, Vec) { +fn compute_radau_left(n: usize) -> Result<(Vec, Vec), QuadratureError> { if n == 1 { - return (vec![-1.0], vec![2.0]); + return Ok((vec![-1.0], vec![2.0])); } // Legendre recurrence coefficients: diff --git a/src/golub_welsch.rs b/src/golub_welsch.rs index 1df1631..b857bf0 100644 --- a/src/golub_welsch.rs +++ b/src/golub_welsch.rs @@ -9,6 +9,8 @@ //! Reference: Golub & Welsch (1969), "Calculation of Gauss Quadrature Rules", //! Mathematics of Computation 23(106), 221-230. +use crate::error::QuadratureError; + #[cfg(not(feature = "std"))] use alloc::{vec, vec::Vec}; #[cfg(not(feature = "std"))] @@ -23,15 +25,24 @@ use num_traits::Float as _; /// /// # Returns /// (nodes, weights) sorted ascending by node. -pub(crate) fn golub_welsch(diag: &[f64], off_diag_sq: &[f64], mu0: f64) -> (Vec, Vec) { +/// +/// # Errors +/// +/// Returns [`QuadratureError::InvalidInput`] if the QL algorithm fails to +/// converge (e.g. due to degenerate recurrence coefficients). +pub(crate) fn golub_welsch( + diag: &[f64], + off_diag_sq: &[f64], + mu0: f64, +) -> Result<(Vec, Vec), QuadratureError> { let n = diag.len(); assert_eq!(off_diag_sq.len(), n.saturating_sub(1)); if n == 0 { - return (vec![], vec![]); + return Ok((vec![], vec![])); } if n == 1 { - return (vec![diag[0]], vec![mu0]); + return Ok((vec![diag[0]], vec![mu0])); } let mut d = diag.to_vec(); @@ -42,14 +53,21 @@ pub(crate) fn golub_welsch(diag: &[f64], off_diag_sq: &[f64], mu0: f64) -> (Vec< let mut z = vec![0.0; n]; z[0] = 1.0; - symmetric_tridiag_eig(&mut d, &mut e, &mut z); + if !symmetric_tridiag_eig(&mut d, &mut e, &mut z) { + return Err(QuadratureError::InvalidInput( + "QL eigenvalue algorithm did not converge", + )); + } let weights: Vec = z.iter().map(|&zk| mu0 * zk * zk).collect(); // Sort ascending by node let mut pairs: Vec<_> = d.into_iter().zip(weights).collect(); - pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); - pairs.into_iter().unzip() + pairs.sort_by(|a, b| { + a.0.partial_cmp(&b.0) + .unwrap_or(core::cmp::Ordering::Equal) + }); + Ok(pairs.into_iter().unzip()) } /// Modify the last diagonal element of a Jacobi matrix so that a prescribed @@ -93,17 +111,22 @@ pub(crate) fn radau_modify(diag: &mut [f64], off_diag_sq: &[f64], x0: f64) { /// /// On entry: `d` = diagonal, `e` = off-diagonal (length n-1). /// On exit: `d` = eigenvalues (unsorted), `z` = first row of eigenvector matrix. +/// +/// Returns `true` if all eigenvalues converged, `false` if the QL iteration +/// hit the maximum iteration count for any eigenvalue. #[allow(clippy::many_single_char_names)] // d, e, z, m, l, g, r, s, c, p, f, b are standard names in tridiagonal eigenvalue algorithms -fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) { +fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) -> bool { let n = d.len(); if n <= 1 { - return; + return true; } // Pad e with a trailing zero for convenience let mut e_ext = vec![0.0; n]; e_ext[..n - 1].copy_from_slice(e); + let mut converged = true; + for l in 0..n { let mut iter_count = 0u32; loop { @@ -122,6 +145,7 @@ fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) { iter_count += 1; if iter_count > 200 { + converged = false; break; } @@ -170,6 +194,8 @@ fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) { } } } + + converged } #[cfg(test)] @@ -189,7 +215,7 @@ mod tests { .collect(); let mu0 = 2.0; - let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0); + let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0).unwrap(); // Weight sum should be 2 let sum: f64 = weights.iter().sum(); @@ -220,7 +246,7 @@ mod tests { let mu0 = 2.0; radau_modify(&mut diag, &off_diag_sq, -1.0); - let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0); + let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0).unwrap(); assert!((nodes[0] - (-1.0)).abs() < 1e-14); assert!((nodes[1] - 1.0 / 3.0).abs() < 1e-14); diff --git a/src/oscillatory.rs b/src/oscillatory.rs index 36ab26b..73dc563 100644 --- a/src/oscillatory.rs +++ b/src/oscillatory.rs @@ -97,7 +97,7 @@ impl OscillatoryIntegrator { /// # Errors /// /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega` - /// is NaN. May also propagate errors from the adaptive fallback when + /// is non-finite. May also propagate errors from the adaptive fallback when /// `|omega * (b - a) / 2|` is small. #[allow(clippy::many_single_char_names)] // a, b, f, n are conventional in quadrature #[allow(clippy::similar_names)] // sum_c_half / sum_s_half are intentionally parallel names @@ -111,7 +111,7 @@ impl OscillatoryIntegrator { where G: Fn(f64) -> f64, { - if a.is_nan() || b.is_nan() || self.omega.is_nan() { + if !a.is_finite() || !b.is_finite() || !self.omega.is_finite() { return Err(QuadratureError::DegenerateInterval); } if (b - a).abs() < f64::EPSILON { @@ -336,7 +336,7 @@ fn modified_chebyshev_moments(theta: f64, n: usize) -> (Vec, Vec) { /// # Errors /// /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega` -/// is NaN. +/// is non-finite. pub fn integrate_oscillatory_sin( f: G, a: f64, @@ -368,7 +368,7 @@ where /// # Errors /// /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega` -/// is NaN. +/// is non-finite. pub fn integrate_oscillatory_cos( f: G, a: f64, @@ -478,6 +478,13 @@ mod tests { assert!(integrate_oscillatory_sin(|_| 1.0, f64::NAN, 1.0, 10.0, 1e-10).is_err()); } + #[test] + fn inf_inputs_rejected() { + assert!(integrate_oscillatory_sin(|_| 1.0, f64::INFINITY, 1.0, 10.0, 1e-10).is_err()); + assert!(integrate_oscillatory_cos(|_| 1.0, 0.0, f64::NEG_INFINITY, 10.0, 1e-10).is_err()); + assert!(integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, f64::INFINITY, 1e-10).is_err()); + } + #[test] fn very_high_omega() { // ∫₀¹ sin(1000x) dx = (1 - cos(1000))/1000 diff --git a/src/tanh_sinh.rs b/src/tanh_sinh.rs index c844eb6..464fc63 100644 --- a/src/tanh_sinh.rs +++ b/src/tanh_sinh.rs @@ -79,7 +79,7 @@ impl TanhSinh { /// /// # Errors /// - /// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is NaN. + /// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is non-finite. #[allow(clippy::many_single_char_names)] // a, b, f, h, k, t, u are conventional in tanh-sinh quadrature #[allow(clippy::cast_precision_loss)] // k is a small loop index, always exact in f64 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] // (7.0/h).ceil() is a small positive value, safe to cast to usize @@ -93,7 +93,7 @@ impl TanhSinh { where G: Fn(f64) -> f64, { - if a.is_nan() || b.is_nan() { + if !a.is_finite() || !b.is_finite() { return Err(QuadratureError::DegenerateInterval); } if (b - a).abs() < f64::EPSILON { @@ -110,6 +110,7 @@ impl TanhSinh { let pi_2 = core::f64::consts::FRAC_PI_2; let mut total_evals = 0usize; + let mut prev_prev_estimate: f64 = 0.0; let mut prev_estimate: f64 = 0.0; let mut h: f64 = 1.0; @@ -234,14 +235,13 @@ impl TanhSinh { } } + prev_prev_estimate = prev_estimate; prev_estimate = estimate; } - // Did not converge within max_levels + // Did not converge within max_levels — use last two estimates for error let error = if self.max_levels > 0 { - // Can't easily recompute the last level difference, - // but prev_estimate is our best value - self.abs_tol.max(self.rel_tol * prev_estimate.abs()) * 10.0 + (prev_estimate - prev_prev_estimate).abs() } else { f64::INFINITY }; @@ -269,7 +269,7 @@ impl TanhSinh { /// /// # Errors /// -/// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is NaN. +/// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is non-finite. pub fn tanh_sinh_integrate( f: G, a: f64, @@ -369,6 +369,12 @@ mod tests { assert!(tanh_sinh_integrate(|x| x, f64::NAN, 1.0, 1e-10).is_err()); } + #[test] + fn inf_bounds_rejected() { + assert!(tanh_sinh_integrate(|x| x, f64::INFINITY, 1.0, 1e-10).is_err()); + assert!(tanh_sinh_integrate(|x| x, 0.0, f64::NEG_INFINITY, 1e-10).is_err()); + } + #[test] fn non_convergence() { let result = TanhSinh::default() @@ -377,4 +383,22 @@ mod tests { .unwrap(); assert!(!result.converged); } + + #[test] + fn non_convergence_error_not_fabricated() { + // With enough levels to produce two estimates but not enough to converge, + // the error estimate should be derived from the last two level estimates + // rather than a fabricated multiple of tolerance. + let tol = 1e-15; + let result = TanhSinh::default() + .with_max_levels(2) + .with_abs_tol(tol) + .with_rel_tol(tol) + .integrate(0.0, 1.0, |x| 1.0 / x.sqrt()) + .unwrap(); + assert!(!result.converged); + assert!(result.error_estimate.is_finite()); + // The error should NOT be a fabricated tol * 10 + assert!((result.error_estimate - tol * 10.0).abs() > 1e-20); + } } diff --git a/src/weighted.rs b/src/weighted.rs index 610007a..35c56c2 100644 --- a/src/weighted.rs +++ b/src/weighted.rs @@ -166,16 +166,27 @@ impl WeightedIntegrator { G: Fn(f64) -> f64, { match &self.weight { - WeightFunction::Jacobi { .. } - | WeightFunction::ChebyshevI - | WeightFunction::ChebyshevII => { + WeightFunction::Jacobi { alpha, beta } => { // Affine map: x in [a,b] <-> t in [-1,1] - // x = (b-a)/2 * t + (a+b)/2 - // w(x) dx in terms of t includes a Jacobian factor - // For Jacobi weight: (1-t)^a (1+t)^b maps directly + // x = half*t + mid, so (b-x) = half*(1-t), (x-a) = half*(1+t) + // ∫_a^b (b-x)^α (x-a)^β f(x) dx = half^(α+β+1) * ∫_{-1}^1 w(t) f(half*t+mid) dt + let half = 0.5 * (b - a); + let mid = 0.5 * (a + b); + self.integrate(|t| f(half * t + mid)) * half.powf(alpha + beta + 1.0) + } + WeightFunction::ChebyshevI => { + // ChebyshevI is Jacobi(-0.5, -0.5): factor = half^0 = 1 + // ∫_a^b f(x)/√((b-x)(x-a)) dx = ∫_{-1}^1 f(half*t+mid)/√(1-t²) dt + let half = 0.5 * (b - a); + let mid = 0.5 * (a + b); + self.integrate(|t| f(half * t + mid)) + } + WeightFunction::ChebyshevII => { + // ChebyshevII is Jacobi(0.5, 0.5): factor = half^2 + // ∫_a^b √((b-x)(x-a)) f(x) dx = half² * ∫_{-1}^1 √(1-t²) f(half*t+mid) dt let half = 0.5 * (b - a); let mid = 0.5 * (a + b); - self.integrate(|t| f(half * t + mid)) * half.powi(1) + self.integrate(|t| f(half * t + mid)) * half * half } WeightFunction::LogWeight => { // Map (0,1] to (a,b]: x = a + (b-a)*u, log singularity at a @@ -408,6 +419,51 @@ mod tests { ); } + #[test] + fn integrate_over_jacobi_on_0_4() { + // Jacobi(0.5, 0.5) on [0,4]: ∫₀⁴ (4-x)^0.5 (x-0)^0.5 f(x) dx + // With f=1: ∫₀⁴ √(x(4-x)) dx = ∫₀⁴ √(4x - x²) dx = 2π (semicircle area) + let wi = WeightedIntegrator::new( + WeightFunction::Jacobi { + alpha: 0.5, + beta: 0.5, + }, + 30, + ) + .unwrap(); + let result = wi.integrate_over(0.0, 4.0, |_| 1.0); + assert!( + (result - 2.0 * core::f64::consts::PI).abs() < 1e-6, + "result={result}, expected={}", + 2.0 * core::f64::consts::PI + ); + } + + #[test] + fn integrate_over_chebyshev_i_on_0_4() { + // ChebyshevI on [0,4]: ∫₀⁴ f(x)/√((4-x)(x-0)) dx + // With f=1: ∫₀⁴ 1/√(4x-x²) dx = π (half=2, factor=1) + let wi = WeightedIntegrator::new(WeightFunction::ChebyshevI, 20).unwrap(); + let result = wi.integrate_over(0.0, 4.0, |_| 1.0); + assert!( + (result - core::f64::consts::PI).abs() < 1e-10, + "result={result}, expected π" + ); + } + + #[test] + fn integrate_over_chebyshev_ii_on_0_4() { + // ChebyshevII on [0,4]: ∫₀⁴ √((4-x)·x) f(x) dx + // With f=1: ∫₀⁴ √(4x-x²) dx = 2π (semicircle area, half=2, factor=half²=4) + let wi = WeightedIntegrator::new(WeightFunction::ChebyshevII, 20).unwrap(); + let result = wi.integrate_over(0.0, 4.0, |_| 1.0); + assert!( + (result - 2.0 * core::f64::consts::PI).abs() < 1e-6, + "result={result}, expected={}", + 2.0 * core::f64::consts::PI + ); + } + #[test] fn integrate_over_log_weight_on_0_2() { // LogWeight: integrate_over(0, 2, f) = width * integrate(|u| f(a + width*u)) From 196c62e86c79d4b01c609300a3d47426fcba92c8 Mon Sep 17 00:00:00 2001 From: Greg von Nessi Date: Fri, 13 Mar 2026 16:41:35 +0000 Subject: [PATCH 2/2] Run cargo fmt --- src/cubature/adaptive.rs | 4 +--- src/cubature/monte_carlo.rs | 5 ++--- src/gauss_jacobi.rs | 6 +++++- src/gauss_lobatto.rs | 5 +---- src/golub_welsch.rs | 5 +---- 5 files changed, 10 insertions(+), 15 deletions(-) diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs index 55faec2..9d5abc1 100644 --- a/src/cubature/adaptive.rs +++ b/src/cubature/adaptive.rs @@ -408,9 +408,7 @@ mod tests { fn inf_bounds_rejected() { assert!(adaptive_cubature(|_| 1.0, &[f64::INFINITY], &[1.0], 1e-8).is_err()); assert!(adaptive_cubature(|_| 1.0, &[0.0], &[f64::NEG_INFINITY], 1e-8).is_err()); - assert!( - adaptive_cubature(|_| 1.0, &[0.0, f64::INFINITY], &[1.0, 1.0], 1e-8).is_err() - ); + assert!(adaptive_cubature(|_| 1.0, &[0.0, f64::INFINITY], &[1.0, 1.0], 1e-8).is_err()); } #[test] diff --git a/src/cubature/monte_carlo.rs b/src/cubature/monte_carlo.rs index df53b41..c460b82 100644 --- a/src/cubature/monte_carlo.rs +++ b/src/cubature/monte_carlo.rs @@ -364,9 +364,8 @@ impl MonteCarloIntegrator { let mean_a = total_sum / total_n as f64; let mean_b = sum / count as f64; let delta = mean_b - mean_a; - total_m2 += m2 - + delta * delta * (total_n as f64 * count as f64) - / (total_n + count) as f64; + total_m2 += + m2 + delta * delta * (total_n as f64 * count as f64) / (total_n + count) as f64; } else { total_m2 += m2; } diff --git a/src/gauss_jacobi.rs b/src/gauss_jacobi.rs index 4aaae56..0bce676 100644 --- a/src/gauss_jacobi.rs +++ b/src/gauss_jacobi.rs @@ -98,7 +98,11 @@ impl_rule_accessors!(GaussJacobi, nodes_doc: "Returns the nodes on \\[-1, 1\\]." /// /// μ₀ = 2^(α+β+1) Γ(α+1)Γ(β+1) / Γ(α+β+2) #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64 -fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> Result<(Vec, Vec), QuadratureError> { +fn compute_jacobi( + n: usize, + alpha: f64, + beta: f64, +) -> Result<(Vec, Vec), QuadratureError> { let ab = alpha + beta; // Diagonal: α_k = (β²-α²) / ((2k+ab)(2k+ab+2)) diff --git a/src/gauss_lobatto.rs b/src/gauss_lobatto.rs index 931afbf..8182940 100644 --- a/src/gauss_lobatto.rs +++ b/src/gauss_lobatto.rs @@ -137,10 +137,7 @@ fn compute_lobatto(n: usize) -> (Vec, Vec) { // Sort ascending (should already be, but ensure) let mut pairs: Vec<_> = nodes.into_iter().zip(weights).collect(); - pairs.sort_by(|a, b| { - a.0.partial_cmp(&b.0) - .unwrap_or(core::cmp::Ordering::Equal) - }); + pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(core::cmp::Ordering::Equal)); let (nodes, weights) = pairs.into_iter().unzip(); (nodes, weights) diff --git a/src/golub_welsch.rs b/src/golub_welsch.rs index b857bf0..5a4e7d4 100644 --- a/src/golub_welsch.rs +++ b/src/golub_welsch.rs @@ -63,10 +63,7 @@ pub(crate) fn golub_welsch( // Sort ascending by node let mut pairs: Vec<_> = d.into_iter().zip(weights).collect(); - pairs.sort_by(|a, b| { - a.0.partial_cmp(&b.0) - .unwrap_or(core::cmp::Ordering::Equal) - }); + pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(core::cmp::Ordering::Equal)); Ok(pairs.into_iter().unzip()) }