Harden input validation and numerical robustness

src/adaptive.rs

@@ -53,10 +53,10 @@ impl PartialOrd for Subinterval {
 
 impl Ord for Subinterval {
     fn cmp(&self, other: &Self) -> Ordering {
-        // Max-heap by error. NaN sorts low.
-        self.error
-            .partial_cmp(&other.error)
-            .unwrap_or(Ordering::Equal)
+        // Max-heap by error, as a total order including NaN. The previous
+        // `partial_cmp(...).unwrap_or(Equal)` made a NaN error compare equal to
+        // every value, which is non-transitive; total_cmp is a proper total order.
+        self.error.total_cmp(&other.error)
     }
 }
 
@@ -187,7 +187,12 @@ impl AdaptiveIntegrator {
             }
         }
 
-        // Validate tolerances
+        // Validate tolerances. A NaN tolerance would pass the `<= 0.0` check
+        // (NaN compares false) and then make every `error <= tol` test false,
+        // spinning to max_evals instead of erroring.
+        if self.abs_tol.is_nan() || self.rel_tol.is_nan() {
+            return Err(QuadratureError::InvalidInput("tolerances must not be NaN"));
+        }
         if self.abs_tol <= 0.0 && self.rel_tol <= 0.0 {
             return Err(QuadratureError::InvalidInput(
                 "at least one tolerance must be positive",
@@ -577,6 +582,23 @@ mod tests {
         assert!((r.value - 2.0).abs() < 1e-14, "value={}", r.value);
     }
 
+    #[test]
+    fn nan_integrand_terminates_without_panic() {
+        // A NaN-producing integrand makes some panel errors NaN; the heap must
+        // stay well-ordered (total_cmp) and the call must return, not panic.
+        let r = AdaptiveIntegrator::default().integrate(0.0, 1.0, |_x: f64| f64::NAN);
+        assert!(r.is_ok());
+    }
+
+    #[test]
+    fn nan_tolerance_rejected() {
+        // A NaN tolerance must error, not slip past the `<= 0.0` check and spin.
+        let r = AdaptiveIntegrator::default()
+            .with_abs_tol(f64::NAN)
+            .integrate(0.0, 1.0, |x: f64| x);
+        assert!(r.is_err());
+    }
+
     /// Non-convergence within tight budget is signalled via converged flag.
     #[test]
     fn non_convergence() {

src/cubature/adaptive.rs

@@ -78,8 +78,8 @@ impl AdaptiveCubature {
     /// # Errors
     ///
     /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
-    /// different lengths or are empty. Returns [`QuadratureError::DegenerateInterval`]
-    /// if any bound is non-finite.
+    /// different lengths, are empty, or if `lower[i] > upper[i]` on any axis.
+    /// Returns [`QuadratureError::DegenerateInterval`] if any bound is non-finite.
     pub fn integrate<G>(
         &self,
         lower: &[f64],
@@ -104,6 +104,13 @@ impl AdaptiveCubature {
             if !lower[i].is_finite() || !upper[i].is_finite() {
                 return Err(QuadratureError::DegenerateInterval);
             }
+            // An inverted axis would give a negative half-width and a silently
+            // sign-flipped / garbage result; reject rather than guess intent.
+            if lower[i] > upper[i] {
+                return Err(QuadratureError::InvalidInput(
+                    "cubature requires lower[i] <= upper[i] on every axis",
+                ));
+            }
         }
 
         if d == 1 {
@@ -255,9 +262,10 @@ impl PartialOrd for SubRegion {
 
 impl Ord for SubRegion {
     fn cmp(&self, other: &Self) -> Ordering {
-        self.error
-            .partial_cmp(&other.error)
-            .unwrap_or(Ordering::Equal)
+        // Total order including NaN. `partial_cmp(...).unwrap_or(Equal)` made a
+        // NaN region error compare equal to every value, breaking transitivity;
+        // total_cmp sorts a NaN error consistently to an extreme.
+        self.error.total_cmp(&other.error)
     }
 }
 
@@ -279,8 +287,8 @@ struct GMDetail {
 /// 1. Center: 1 point
 /// 2. Along each axis at ±λ₂: 2d points
 /// 3. Along each axis at ±λ₃: 2d points
-/// 4. All vertices at (±λ₄, ..., ±λ₄): 2^d points
-/// 5. All pairs of axes at (±λ₅, ±λ₅, 0, ...): 2d(d-1) points
+/// 4. All vertices at (±λ₅, ..., ±λ₅): 2^d points
+/// 5. All pairs of axes at (±λ₄, ±λ₄, 0, ...): 2d(d-1) points
 ///
 /// Total: 1 + 4d + 2^d + 2d(d-1) = 2^d + 2d² + 2d + 1 evaluations.
 fn genz_malik_eval<G>(d: usize, lower: &[f64], upper: &[f64], f: &G) -> GMDetail
@@ -431,6 +439,12 @@ mod tests {
         assert!(adaptive_cubature(|_| 1.0, &[0.0, f64::INFINITY], &[1.0, 1.0], 1e-8).is_err());
     }
 
+    #[test]
+    fn inverted_box_rejected() {
+        // An axis with lower > upper would give a silently sign-flipped result.
+        assert!(adaptive_cubature(|_| 1.0, &[0.0, 1.0], &[1.0, 0.0], 1e-8).is_err());
+    }
+
     #[test]
     fn dimension_cap() {
         let lower = vec![0.0; 31];

src/cubature/monte_carlo.rs

@@ -107,7 +107,8 @@ impl MonteCarloIntegrator {
     /// # Errors
     ///
     /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
-    /// different lengths, are empty, or if the number of samples is zero.
+    /// different lengths, are empty, or if the number of samples is zero, or
+    /// [`QuadratureError::DegenerateInterval`] if any bound is non-finite.
     /// For Sobol and Halton methods, also returns an error if the dimension
     /// exceeds the supported limit (40 for Sobol, 100 for Halton).
     pub fn integrate<G>(
@@ -130,6 +131,9 @@ impl MonteCarloIntegrator {
                 "number of samples must be >= 1",
             ));
         }
+        if lower.iter().chain(upper).any(|x| !x.is_finite()) {
+            return Err(QuadratureError::DegenerateInterval);
+        }
 
         let widths: Vec<f64> = (0..d).map(|j| upper[j] - lower[j]).collect();
         let volume: f64 = widths.iter().product();
@@ -181,7 +185,9 @@ impl MonteCarloIntegrator {
         }
 
         let variance = if n > 1 { m2 / (n - 1) as f64 } else { 0.0 };
-        let std_error = (variance / n as f64).sqrt() * volume;
+        // volume may be negative for reversed bounds; the integral value carries
+        // that sign but a standard error must not.
+        let std_error = (variance / n as f64).sqrt() * volume.abs();
 
         Ok(QuadratureResult {
             value: volume * mean,
@@ -261,7 +267,8 @@ impl MonteCarloIntegrator {
     /// # Errors
     ///
     /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
-    /// different lengths, are empty, or if the number of samples is zero.
+    /// different lengths, are empty, or if the number of samples is zero, or
+    /// [`QuadratureError::DegenerateInterval`] if any bound is non-finite.
     /// For Sobol and Halton methods, also returns an error if the dimension
     /// exceeds the supported limit (40 for Sobol, 100 for Halton).
     #[cfg(feature = "parallel")]
@@ -285,6 +292,9 @@ impl MonteCarloIntegrator {
                 "number of samples must be >= 1",
             ));
         }
+        if lower.iter().chain(upper).any(|x| !x.is_finite()) {
+            return Err(QuadratureError::DegenerateInterval);
+        }
 
         let widths: Vec<f64> = (0..d).map(|j| upper[j] - lower[j]).collect();
         let volume: f64 = widths.iter().product();
@@ -382,7 +392,9 @@ impl MonteCarloIntegrator {
             0.0
         };
         #[allow(clippy::cast_precision_loss)]
-        let std_error = (variance / total_n as f64).sqrt() * volume;
+        // volume may be negative for reversed bounds; the integral value carries
+        // that sign but a standard error must not.
+        let std_error = (variance / total_n as f64).sqrt() * volume.abs();
 
         Ok(QuadratureResult {
             value: volume * mean,
@@ -540,6 +552,24 @@ mod tests {
         assert!((result.value - 1.0).abs() < 1e-10, "value={}", result.value);
     }
 
+    #[test]
+    fn non_finite_bounds_rejected() {
+        assert!(monte_carlo_integrate(|_| 1.0, &[f64::NAN], &[1.0], 100).is_err());
+        assert!(monte_carlo_integrate(|_| 1.0, &[0.0], &[f64::INFINITY], 100).is_err());
+    }
+
+    #[test]
+    fn reversed_bounds_std_error_nonnegative() {
+        // Reversed bounds make the volume negative; the value may carry the
+        // sign but the reported standard error must not.
+        let r = MonteCarloIntegrator::default()
+            .with_method(MCMethod::Plain)
+            .with_samples(1000)
+            .integrate(&[1.0], &[0.0], |x| x[0])
+            .unwrap();
+        assert!(r.error_estimate >= 0.0, "std_error={}", r.error_estimate);
+    }
+
     /// Integral of x over [0,1] = 0.5.
     #[test]
     fn linear_1d() {

src/error.rs

@@ -11,7 +11,8 @@ use core::fmt;
 pub enum QuadratureError {
     /// The requested order `n` is zero.
     ZeroOrder,
-    /// The integration interval is degenerate (NaN bounds).
+    /// The integration interval is degenerate (a bound is non-finite, i.e. NaN
+    /// or infinite).
     DegenerateInterval,
     /// An input parameter is invalid.
     InvalidInput(&'static str),
@@ -27,5 +28,6 @@ impl fmt::Display for QuadratureError {
     }
 }
 
-#[cfg(feature = "std")]
-impl std::error::Error for QuadratureError {}
+// `core::error::Error` (stable since Rust 1.81; MSRV is 1.93) is available
+// without `std`, so the impl is unconditional and `no_std` users get it too.
+impl core::error::Error for QuadratureError {}

src/gauss_hermite.rs

@@ -38,7 +38,8 @@ impl GaussHermite {
     ///
     /// # Errors
     ///
-    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero, or
+    /// [`QuadratureError::InvalidInput`] if the eigensolver fails to converge.
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);

src/gauss_kronrod.rs

@@ -156,7 +156,8 @@ impl GaussKronrod {
         }
 
         let half_width = 0.5 * (b - a);
-        let midpoint = 0.5 * (a + b);
+        // a + (b - a)/2 rather than (a + b)/2, which overflows for large |a|, |b|.
+        let midpoint = a + half_width;
         let n = self.xgk.len(); // number of (positive + center) abscissae
 
         let mut kronrod_sum = 0.0_f64;

src/gauss_legendre.rs

@@ -135,7 +135,9 @@ fn compute_newton(n: usize, m: usize, nodes: &mut [f64], weights: &mut [f64]) {
         }
 
         let (_, p_n_deriv) = legendre_eval(n, x);
-        let w = 2.0 / ((1.0 - x * x) * p_n_deriv * p_n_deriv);
+        // (1-x)(1+x) instead of 1-x*x: for x near ±1 the latter cancels (both
+        // terms ~1), while (1-x) is exact there by Sterbenz's lemma.
+        let w = 2.0 / ((1.0 - x) * (1.0 + x) * p_n_deriv * p_n_deriv);
 
         // Store: nodes are in increasing order (-1 to 1).
         // i=0 is the most negative node.
@@ -168,7 +170,8 @@ pub(crate) fn legendre_eval(n: usize, x: f64) -> (f64, f64) {
 
     // P_n'(x) = n * (x * P_n(x) - P_{n-1}(x)) / (x^2 - 1)
     let nf = n as f64;
-    let deriv = nf * (x * p_curr - p_prev) / (x * x - 1.0);
+    // -(1-x)(1+x) instead of x*x-1: avoids cancellation for x near ±1.
+    let deriv = nf * (x * p_curr - p_prev) / -((1.0 - x) * (1.0 + x));
 
     (p_curr, deriv)
 }
@@ -402,7 +405,9 @@ fn compute_gl_pair_par(n: usize) -> (Vec<f64>, Vec<f64>) {
                 }
 
                 let (_, p_n_deriv) = legendre_eval(n, x);
-                let w = 2.0 / ((1.0 - x * x) * p_n_deriv * p_n_deriv);
+                // (1-x)(1+x) instead of 1-x*x: for x near ±1 the latter cancels (both
+                // terms ~1), while (1-x) is exact there by Sterbenz's lemma.
+                let w = 2.0 / ((1.0 - x) * (1.0 + x) * p_n_deriv * p_n_deriv);
                 (i, x, w)
             })
             .collect();

src/gauss_radau.rs

@@ -43,7 +43,9 @@ impl GaussRadau {
     ///
     /// # Errors
     ///
-    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero, or
+    /// [`QuadratureError::InvalidInput`] if the endpoint-modification step or
+    /// the eigensolver fails.
     pub fn left(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -60,7 +62,9 @@ impl GaussRadau {
     ///
     /// # Errors
     ///
-    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero, or
+    /// [`QuadratureError::InvalidInput`] if the endpoint-modification step or
+    /// the eigensolver fails.
     pub fn right(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -100,7 +104,7 @@ fn compute_radau_left(n: usize) -> Result<(Vec<f64>, Vec<f64>), QuadratureError>
         .collect();
 
     // Modify the last diagonal element so that -1 is an eigenvalue
-    radau_modify(&mut diag, &off_diag_sq, -1.0);
+    radau_modify(&mut diag, &off_diag_sq, -1.0)?;
 
     let mu0 = 2.0; // integral of 1 over [-1, 1]
     golub_welsch(&diag, &off_diag_sq, mu0)