Defensive validation, numerical robustness, golub_welsch Result API (v0.2.0) (#13)

* Defensive input validation, numerical robustness, and golub_welsch Result API (v0.2.0)

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

* Run cargo fmt

Author: gvonness-apolitical

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.
 

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"

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`,

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.
 

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());
+    }
 }

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<G>(
@@ -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());
+    }
 }

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<G>(
         &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<SubRegion> = 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<G>(
     f: G,
     lower: &[f64],
@@ -400,6 +404,20 @@ 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() {

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<f64>, weights: Vec<f64>, 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,

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,21 @@ 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 +431,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,

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<f64>, Vec<f64>) {
+fn compute_hermite(n: usize) -> Result<(Vec<f64>, Vec<f64>), QuadratureError> {
     let diag = vec![0.0; n];
     let off_diag_sq: Vec<f64> = (1..n).map(|k| k as f64 / 2.0).collect();
     let mu0 = core::f64::consts::PI.sqrt();