From 165c812470001ed1edd6f83042d2c44d3f0238b8 Mon Sep 17 00:00:00 2001 From: Greg von Nessi Date: Wed, 1 Jul 2026 00:15:15 +0100 Subject: [PATCH] Detect floating-point roundoff floor in adaptive integration Adaptive Gauss-Kronrod integration could spin to max_evals when asked for a tolerance tighter than f64 can deliver. Add QUADPACK-style roundoff detection so it stops at the floating-point floor instead: - Apply the QUADPACK 50*eps_mach*integral|f| floor to each Gauss-Kronrod panel's error estimate (dqk21.f), so a panel integrated to rounding level reports an honest error rather than an unachievably small one. - Detect the ier=2 roundoff stall in the adaptive loop: a run of consecutive worst-panel bisections that neither reduce the error nor move the integral (and are not merely under-resolved) signals the floating-point floor. A bespoke cap-saturation guard keeps merely under-resolved panels from being mistaken for floor-limited ones. - Surface the outcome on QuadratureResult via a new roundoff_limited field (with is_roundoff_limited()), mutually exclusive with converged: converged / roundoff-limited / budget-exhausted are the three distinct stopping states. converged now means strictly "tolerance met". QuadratureResult is now #[non_exhaustive] so future fields are non-breaking; the field addition itself is a breaking change, hence the 0.3.0 version bump. --- CHANGELOG.md | 27 +++++++ Cargo.toml | 2 +- src/adaptive.rs | 145 +++++++++++++++++++++++++++++++++--- src/cauchy_pv.rs | 36 +++++++++ src/cubature/adaptive.rs | 1 + src/cubature/monte_carlo.rs | 4 + src/gauss_kronrod.rs | 29 +++++++- src/oscillatory.rs | 2 + src/result.rs | 32 +++++++- src/tanh_sinh.rs | 3 + 10 files changed, 265 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 010a173..f64ec12 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,33 @@ 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.3.0] - Unreleased + +### Added + +- Adaptive integration now detects the floating-point roundoff floor (analogous + to QUADPACK's `ier = 2`): when the requested tolerance is tighter than the + accuracy achievable in `f64`, it stops subdividing instead of spinning to + `max_evals`. The per-panel Gauss-Kronrod error estimate gains the QUADPACK + `max(50·εmach·∫|f|, error)` floor so a saturated panel reports an honest error. + +### Changed + +- **Breaking:** `QuadratureResult` gains a `roundoff_limited: bool` field (with an + `is_roundoff_limited()` accessor) distinguishing "stopped at the floating-point + floor" from "exhausted the evaluation budget"; `converged` now means strictly + "requested tolerance met". The struct is now `#[non_exhaustive]`, so subsequent + field additions will not be breaking. + +### Fixed + +- This release also rolls up a number of correctness and robustness fixes: the + corrupted Sobol direction-number table, the `resasc`-based Gauss-Kronrod error + estimate, oscillatory frequency bounding, tanh-sinh handling of reversed + limits, the Jacobi recurrence coefficients, adaptive and cubature resource + bounds, transform-endpoint guards, stricter input validation, and + Golub-Welsch eigensolver guards. + ## [0.2.0] - 2026-03-13 ### Changed diff --git a/Cargo.toml b/Cargo.toml index 4cfd220..1f4f2e3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bilby" -version = "0.2.0" +version = "0.3.0" edition = "2021" rust-version = "1.93" description = "A high-performance numerical quadrature (integration) library for Rust" diff --git a/src/adaptive.rs b/src/adaptive.rs index 99e6581..b3cb66e 100644 --- a/src/adaptive.rs +++ b/src/adaptive.rs @@ -127,8 +127,11 @@ impl AdaptiveIntegrator { /// Adaptively integrate `f` over \[a, b\]. /// /// Returns a [`QuadratureResult`] with the best estimate and error bound. - /// If the tolerance could not be achieved within the evaluation budget, - /// `converged` will be `false` but the best estimate is still returned. + /// If the tolerance could not be achieved, `converged` will be `false` but + /// the best estimate is still returned; in that case + /// [`roundoff_limited`](QuadratureResult::roundoff_limited) distinguishes + /// reaching the floating-point floor (the tolerance was unachievable in + /// `f64`) from exhausting the evaluation budget. /// /// # Errors /// @@ -150,7 +153,10 @@ impl AdaptiveIntegrator { /// /// Break points split the interval at known discontinuities or singularities. /// The adaptive algorithm then refines each sub-interval independently in a - /// global priority queue. + /// global priority queue. As with [`integrate`](Self::integrate), the result + /// reports [`converged`](QuadratureResult::converged) and + /// [`roundoff_limited`](QuadratureResult::roundoff_limited) so a tolerance + /// that is unachievable in `f64` is distinguishable from budget exhaustion. /// /// # Errors /// @@ -247,7 +253,19 @@ impl AdaptiveIntegrator { }); } - // Adaptive refinement loop + // Adaptive refinement loop. `roundoff_count` tracks consecutive panels + // whose bisection neither reduced the error nor moved the integral (and + // were not merely under-resolved) — the QUADPACK ier=2 roundoff signal; + // after ROUNDOFF_STALL_LIMIT of them the floating-point floor is reached. + // QUADPACK's dqagse accumulates iroff1/iroff2/iroff3 across the whole + // run (flagging at 10/20); the consecutive-with-reset policy here is a + // deliberately more conservative simplification — any genuinely + // improving bisection resets the count, so a flag requires an unbroken + // run of floor-level steps. `no_error_reduction` also folds in dqagse's + // iroff3 "error failed to decrease" case. + const ROUNDOFF_STALL_LIMIT: usize = 6; + let mut roundoff_count = 0usize; + let mut roundoff_limited = false; while !self.tolerance_met(total_estimate, total_error) && num_evals + 2 * evals_per_call <= self.max_evals { @@ -263,6 +281,10 @@ impl AdaptiveIntegrator { // re-evaluating the same panel. Exact equality is orientation- // independent, so this also holds for reversed bounds (a > b). Its // error stays in total_error, so convergence is not falsely declared. + // A drop is itself an at-resolution-limit event: it leaves + // `roundoff_count` untouched (neither incremented nor reset), so the + // consecutive-stall run is transparent to interleaved drops — which + // is benign, since a drop cannot manufacture the increment conditions. #[allow(clippy::float_cmp)] if mid == worst.a || mid == worst.b { continue; @@ -272,9 +294,29 @@ impl AdaptiveIntegrator { let right = gk.integrate_detail(mid, worst.b, f); num_evals += 2 * evals_per_call; + // QUADPACK dqagse roundoff-stall detection (ier=2): a bisection that + // neither reduced the error (children error >= 0.99 * parent) nor + // moved the integral (area stable) signals the floating-point floor + // — UNLESS a child's error merely saturated its resasc cap + // (`error == resasc`), which means the panel is under-resolved, not + // floor-limited, and refinement must continue. Counting only + // consecutive such events avoids a false stall on a hard integrand. + let children_estimate = left.estimate + right.estimate; + let children_error = left.error + right.error; + #[allow(clippy::float_cmp)] + let child_saturated = left.error == left.resasc || right.error == right.resasc; + let no_error_reduction = children_error >= worst.error * 0.99; + let area_stable = + (children_estimate - worst.estimate).abs() <= 1.0e-5 * children_estimate.abs(); + if no_error_reduction && area_stable && !child_saturated { + roundoff_count += 1; + } else { + roundoff_count = 0; + } + // Update totals: remove old contribution, add new - total_estimate = total_estimate - worst.estimate + left.estimate + right.estimate; - total_error = total_error - worst.error + left.error + right.error; + total_estimate = total_estimate - worst.estimate + children_estimate; + total_error = total_error - worst.error + children_error; heap.push(Subinterval { a: worst.a, @@ -288,15 +330,27 @@ impl AdaptiveIntegrator { estimate: right.estimate, error: right.error, }); + + if roundoff_count >= ROUNDOFF_STALL_LIMIT { + roundoff_limited = true; + break; + } } let converged = self.tolerance_met(total_estimate, total_error); + // A final bisection can both bring the total error within tolerance and + // trip the roundoff counter on the same iteration (the tolerance check + // is at the loop head, the roundoff break at its tail). `converged` + // takes precedence, keeping the three outcomes mutually exclusive: + // converged, roundoff-limited, or budget-exhausted. + let roundoff_limited = roundoff_limited && !converged; Ok(QuadratureResult { value: total_estimate, error_estimate: total_error, num_evals, converged, + roundoff_limited, }) } @@ -356,7 +410,10 @@ fn build_intervals(a: f64, b: f64, breaks: &[f64]) -> Vec<(f64, f64)> { /// Adaptively integrate `f` over \[a, b\] with default settings. /// -/// Uses G7-K15 with absolute and relative tolerance both set to `tol`. +/// Uses G7-K15 with absolute and relative tolerance both set to `tol`. Inspect +/// [`converged`](QuadratureResult::converged) and +/// [`roundoff_limited`](QuadratureResult::roundoff_limited) on the result to +/// tell a met tolerance from the floating-point floor or budget exhaustion. /// /// # Errors /// @@ -570,7 +627,12 @@ mod tests { /// Builder with higher-order pair achieves tight tolerance. #[test] - fn builder_high_order() { + fn builder_high_order_roundoff_limited() { + // 1e-14 is tighter than the f64 roundoff floor (≈ 50·εmach·∫|sin| ≈ + // 2.2e-14 over [0,π]), so the requested tolerance cannot be certified. + // The integrator must stop roundoff-limited (QUADPACK ier=2) — value + // accurate to machine precision, converged=false, roundoff_limited=true + // — rather than spinning to max_evals. let r = AdaptiveIntegrator::default() .with_pair(GKPair::G15K31) .with_abs_tol(1e-14) @@ -578,8 +640,71 @@ mod tests { .with_max_evals(100_000) .integrate(0.0, core::f64::consts::PI, f64::sin) .unwrap(); - assert!(r.converged, "err={}", r.error_estimate); - assert!((r.value - 2.0).abs() < 1e-14, "value={}", r.value); + assert!((r.value - 2.0).abs() < 1e-13, "value={}", r.value); + assert!(r.roundoff_limited, "expected roundoff_limited"); + assert!(!r.converged, "expected not converged at unachievable 1e-14"); + assert!( + r.num_evals < 5_000, + "did not stop early: num_evals={}", + r.num_evals + ); + } + + #[test] + fn smooth_integrand_roundoff_limited_at_tight_tol() { + // exp over [0,1] at rel_tol 1e-15 (below the floor ≈ 50·εmach·∫|exp|): + // the value is accurate but the tolerance is unachievable. + let r = AdaptiveIntegrator::default() + .with_abs_tol(0.0) + .with_rel_tol(1e-15) + .integrate(0.0, 1.0, f64::exp) + .unwrap(); + let exact = core::f64::consts::E - 1.0; + assert!((r.value - exact).abs() < 1e-13, "value={}", r.value); + assert!(r.roundoff_limited && !r.converged); + assert!(r.num_evals < 5_000, "num_evals={}", r.num_evals); + } + + #[test] + fn under_resolved_oscillatory_not_falsely_converged() { + // sin(1000x) with a budget far too small to resolve ~159 oscillations: + // the panels saturate their resasc cap, so the roundoff guard must NOT + // declare success on the (inaccurate) value — converged/roundoff_limited + // may only be true if the value is in fact accurate. + let exact = (1.0 - 1000.0_f64.cos()) / 1000.0; + let r = AdaptiveIntegrator::default() + .with_abs_tol(1e-12) + .with_rel_tol(1e-12) + .with_max_evals(2_000) + .integrate(0.0, 1.0, |x: f64| (1000.0 * x).sin()) + .unwrap(); + assert!( + !(r.roundoff_limited || r.converged) || (r.value - exact).abs() < 1e-6, + "falsely reported success: value={}, exact={exact}, conv={}, ro={}", + r.value, + r.converged, + r.roundoff_limited + ); + } + + #[test] + fn discontinuity_without_break_not_falsely_converged() { + // A jump at x=2 with no break point. The straddling panels are + // under-resolved (cap-saturated); success may only be reported with an + // accurate value. + let step = |x: f64| if x < 2.0 { 1.0 } else { 3.0 }; + let r = AdaptiveIntegrator::default() + .with_abs_tol(1e-14) + .with_rel_tol(1e-14) + .integrate(0.0, 4.0, step) + .unwrap(); + assert!( + !(r.roundoff_limited || r.converged) || (r.value - 8.0).abs() < 1e-6, + "value={}, conv={}, ro={}", + r.value, + r.converged, + r.roundoff_limited + ); } #[test] diff --git a/src/cauchy_pv.rs b/src/cauchy_pv.rs index 10c4213..5721558 100644 --- a/src/cauchy_pv.rs +++ b/src/cauchy_pv.rs @@ -72,6 +72,11 @@ impl CauchyPV { /// /// The singularity c must be strictly inside (a, b). /// + /// The result's [`converged`](QuadratureResult::converged) and + /// [`roundoff_limited`](QuadratureResult::roundoff_limited) flags are + /// propagated from the underlying adaptive integration of the subtracted + /// integrand. + /// /// # Errors /// /// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is non-finite. @@ -133,6 +138,7 @@ impl CauchyPV { error_estimate: adaptive_result.error_estimate, num_evals: adaptive_result.num_evals + 1, // +1 for f(c) converged: adaptive_result.converged, + roundoff_limited: adaptive_result.roundoff_limited, }) } } @@ -254,6 +260,36 @@ mod tests { assert!(pv_integrate(|x| x, 0.0, 1.0, f64::NAN, 1e-10).is_err()); } + #[test] + fn roundoff_limited_propagates_through_wrapper() { + // The subtracted integrand [f(x) - f(c)]/(x - c) is smooth here (f is a + // polynomial), so requesting a purely relative tolerance below the f64 + // floor drives the inner adaptive integrator to the roundoff floor. The + // wrapper must surface that outcome rather than reporting a spurious + // `converged`, while still returning the accurate principal value. + let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln(); + let result = CauchyPV::default() + .with_abs_tol(0.0) + .with_rel_tol(1e-15) + .with_max_evals(100_000) + .integrate(0.0, 1.0, 0.3, |x| x * x) + .unwrap(); + assert!( + result.roundoff_limited, + "expected roundoff_limited to propagate; converged={}, error={}", + result.converged, result.error_estimate + ); + assert!( + !result.converged, + "tolerance below the f64 floor is unachievable" + ); + assert!( + (result.value - exact).abs() < 1e-8, + "value should still be accurate: value={}, expected={exact}", + result.value + ); + } + #[test] fn inf_inputs_rejected() { assert!(pv_integrate(|x| x, f64::INFINITY, 1.0, 0.5, 1e-10).is_err()); diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs index 245073d..4d447e5 100644 --- a/src/cubature/adaptive.rs +++ b/src/cubature/adaptive.rs @@ -200,6 +200,7 @@ impl AdaptiveCubature { error_estimate: global_error, num_evals: total_evals, converged, + roundoff_limited: false, }) } } diff --git a/src/cubature/monte_carlo.rs b/src/cubature/monte_carlo.rs index 3dc6022..38d7d05 100644 --- a/src/cubature/monte_carlo.rs +++ b/src/cubature/monte_carlo.rs @@ -194,6 +194,7 @@ impl MonteCarloIntegrator { error_estimate: std_error, num_evals: n, converged: true, + roundoff_limited: false, }) } @@ -253,6 +254,7 @@ impl MonteCarloIntegrator { error_estimate: error, num_evals: n + n / 2, converged: true, + roundoff_limited: false, }) } @@ -401,6 +403,7 @@ impl MonteCarloIntegrator { error_estimate: std_error, num_evals: total_n, converged: true, + roundoff_limited: false, }) } @@ -459,6 +462,7 @@ impl MonteCarloIntegrator { error_estimate: error, num_evals: n + half, converged: true, + roundoff_limited: false, }) } } diff --git a/src/gauss_kronrod.rs b/src/gauss_kronrod.rs index b821f9a..f5c0005 100644 --- a/src/gauss_kronrod.rs +++ b/src/gauss_kronrod.rs @@ -24,8 +24,13 @@ use num_traits::Float as _; pub(crate) struct GKDetail { /// Kronrod estimate of the integral. pub estimate: f64, - /// Error estimate (Kronrod - Gauss difference, scaled). + /// Error estimate (Kronrod - Gauss difference, scaled), with the QUADPACK + /// roundoff floor applied. pub error: f64, + /// `resasc` = integral of `|f - mean|` over the panel. The adaptive + /// integrator uses `error == resasc` to detect an under-resolved panel + /// (its error saturated the resasc cap) versus a roundoff-limited one. + pub resasc: f64, } /// Identifies which embedded Gauss-Kronrod pair to use. @@ -152,6 +157,7 @@ impl GaussKronrod { return GKDetail { estimate: 0.0, error: 0.0, + resasc: 0.0, }; } @@ -218,17 +224,23 @@ impl GaussKronrod { // The error is scaled against resasc — the variation of f about its // mean — NOT resabs (= abs_sum, ∫|f|): a constant offset is integrated // exactly by both rules and must not deflate the error estimate. + // resabs = ∫|f| (Σ wgk·|f|); resasc = ∫|f − f̄| (variation about the + // mean f̄ = reskh). Both come from the retained panel values in one pass. let reskh = kronrod_sum * 0.5; let mut resasc = 0.0_f64; + let mut resabs = 0.0_f64; for i in 0..n { let wk = self.wgk[i]; if self.xgk[i] == 0.0 { resasc += wk * (fv1[i] - reskh).abs(); + resabs += wk * fv1[i].abs(); } else { resasc += wk * ((fv1[i] - reskh).abs() + (fv2[i] - reskh).abs()); + resabs += wk * (fv1[i].abs() + fv2[i].abs()); } } resasc *= half_width.abs(); + resabs *= half_width.abs(); // QUADPACK error estimate (Piessens et al. 1983, dqk15.f): scale the // raw Gauss-Kronrod difference δ by resasc via the (200·δ/resasc)^1.5 @@ -239,8 +251,21 @@ impl GaussKronrod { if resasc != 0.0 && error > 0.0 { error = resasc * (200.0 * error / resasc).powf(1.5).min(1.0); } + // QUADPACK roundoff floor: the error cannot be tighter than the f64 + // rounding level of the panel, 50·εmach·resabs. Gated as in dqk15.f to + // avoid underflow in the product. This is a max(), so it can only raise + // the estimate — it cannot reintroduce the resabs-deflation bug that + // scaling against resasc (not resabs) fixed. + let epmach = f64::EPSILON; + if resabs > f64::MIN_POSITIVE / (50.0 * epmach) { + error = error.max(50.0 * epmach * resabs); + } - GKDetail { estimate, error } + GKDetail { + estimate, + error, + resasc, + } } /// Returns the Kronrod order. diff --git a/src/oscillatory.rs b/src/oscillatory.rs index 0f24156..1fdee68 100644 --- a/src/oscillatory.rs +++ b/src/oscillatory.rs @@ -123,6 +123,7 @@ impl OscillatoryIntegrator { error_estimate: 0.0, num_evals: 0, converged: true, + roundoff_limited: false, }); } @@ -233,6 +234,7 @@ impl OscillatoryIntegrator { error_estimate: error, num_evals, converged, + roundoff_limited: false, }) } diff --git a/src/result.rs b/src/result.rs index db8547c..1dfe100 100644 --- a/src/result.rs +++ b/src/result.rs @@ -5,9 +5,15 @@ use num_traits::Float; /// Outcome of a quadrature computation. /// /// Always returned by adaptive integration methods. Contains the best estimate -/// even when convergence was not achieved — check [`converged`](Self::converged) -/// to determine whether the requested tolerance was met. +/// even when convergence was not achieved. Check [`converged`](Self::converged) +/// to determine whether the requested tolerance was met; if it is `false`, +/// [`roundoff_limited`](Self::roundoff_limited) distinguishes hitting the +/// floating-point floor from exhausting the evaluation budget. +/// +/// Marked `#[non_exhaustive]` so future fields can be added without a breaking +/// change; construct one only through the library's integration methods. #[derive(Debug, Clone)] +#[non_exhaustive] pub struct QuadratureResult { /// The computed integral value. pub value: F, @@ -15,8 +21,21 @@ pub struct QuadratureResult { pub error_estimate: F, /// Number of integrand evaluations. pub num_evals: usize, - /// Whether the requested tolerance was achieved. + /// Whether the requested tolerance was achieved (i.e. `error_estimate` is + /// within the requested absolute/relative tolerance). pub converged: bool, + /// Whether integration stopped at the floating-point roundoff floor + /// (analogous to QUADPACK's `ier = 2`): the requested tolerance was tighter + /// than the accuracy achievable in `f64`, so further subdivision could not + /// reduce `error_estimate`. When this is `true` the `value` is typically as + /// accurate as floating-point allows. Only the one-dimensional Gauss-Kronrod + /// adaptive integrator and its delegators (Cauchy principal value, + /// one-dimensional cubature) ever set it; it is always `false` for the + /// Monte-Carlo, tanh-sinh, oscillatory, and Genz-Malik cubature methods. + /// Mutually exclusive with `converged`: if the tolerance was met this is + /// `false`, so `converged` / `roundoff_limited` / neither are the three + /// distinct stopping outcomes (the last being budget exhaustion). + pub roundoff_limited: bool, } impl QuadratureResult { @@ -25,4 +44,11 @@ impl QuadratureResult { pub fn is_converged(&self) -> bool { self.converged } + + /// Returns `true` if integration stopped at the floating-point roundoff + /// floor (the requested tolerance was unachievable in `f64`). + #[inline] + pub fn is_roundoff_limited(&self) -> bool { + self.roundoff_limited + } } diff --git a/src/tanh_sinh.rs b/src/tanh_sinh.rs index dc392cc..828f8c8 100644 --- a/src/tanh_sinh.rs +++ b/src/tanh_sinh.rs @@ -102,6 +102,7 @@ impl TanhSinh { error_estimate: 0.0, num_evals: 0, converged: true, + roundoff_limited: false, }); } @@ -235,6 +236,7 @@ impl TanhSinh { error_estimate: error, num_evals: total_evals, converged: true, + roundoff_limited: false, }); } } @@ -255,6 +257,7 @@ impl TanhSinh { error_estimate: error, num_evals: total_evals, converged: false, + roundoff_limited: false, }) } }