diff --git a/src/adaptive.rs b/src/adaptive.rs index 4a6b41b..99e6581 100644 --- a/src/adaptive.rs +++ b/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() { diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs index 6a17d73..245073d 100644 --- a/src/cubature/adaptive.rs +++ b/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( &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(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]; diff --git a/src/cubature/monte_carlo.rs b/src/cubature/monte_carlo.rs index c460b82..3dc6022 100644 --- a/src/cubature/monte_carlo.rs +++ b/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( @@ -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 = (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 = (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() { diff --git a/src/error.rs b/src/error.rs index 69fe2c3..9fe3937 100644 --- a/src/error.rs +++ b/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 {} diff --git a/src/gauss_hermite.rs b/src/gauss_hermite.rs index ae68ef3..cfd0533 100644 --- a/src/gauss_hermite.rs +++ b/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 { if n == 0 { return Err(QuadratureError::ZeroOrder); diff --git a/src/gauss_kronrod.rs b/src/gauss_kronrod.rs index e57df13..b821f9a 100644 --- a/src/gauss_kronrod.rs +++ b/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; diff --git a/src/gauss_legendre.rs b/src/gauss_legendre.rs index 6a5fcdf..bbaa3ca 100644 --- a/src/gauss_legendre.rs +++ b/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, Vec) { } 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(); diff --git a/src/gauss_radau.rs b/src/gauss_radau.rs index b16100a..e08b06e 100644 --- a/src/gauss_radau.rs +++ b/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 { 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 { if n == 0 { return Err(QuadratureError::ZeroOrder); @@ -100,7 +104,7 @@ fn compute_radau_left(n: usize) -> Result<(Vec, Vec), 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) diff --git a/src/golub_welsch.rs b/src/golub_welsch.rs index 5a4e7d4..90e2aca 100644 --- a/src/golub_welsch.rs +++ b/src/golub_welsch.rs @@ -46,6 +46,13 @@ pub(crate) fn golub_welsch( } let mut d = diag.to_vec(); + // A negative squared off-diagonal means a non-PSD (malformed) Jacobi matrix; + // its sqrt would be NaN and surface later as a misleading "did not converge". + if off_diag_sq.iter().any(|&b| b < 0.0) { + return Err(QuadratureError::InvalidInput( + "off-diagonal squared values must be non-negative (non-PSD Jacobi matrix)", + )); + } let mut e: Vec = off_diag_sq.iter().map(|&b| b.sqrt()).collect(); // z[k] tracks the first component of the k-th eigenvector. @@ -74,13 +81,20 @@ pub(crate) fn golub_welsch( /// - `diag`: diagonal elements α₀, ..., α_{n-1} (modified in place) /// - `off_diag_sq`: squared off-diagonal elements β₁, ..., β_{n-1} /// - `x0`: the prescribed node (endpoint) -pub(crate) fn radau_modify(diag: &mut [f64], off_diag_sq: &[f64], x0: f64) { +pub(crate) fn radau_modify( + diag: &mut [f64], + off_diag_sq: &[f64], + x0: f64, +) -> Result<(), QuadratureError> { let n = diag.len(); assert_eq!(off_diag_sq.len(), n.saturating_sub(1)); - if n <= 1 { + if n == 0 { + return Ok(()); + } + if n == 1 { diag[0] = x0; - return; + return Ok(()); } // Evaluate the continued fraction r = p_{n-1}(x0) / p_{n-2}(x0) @@ -92,15 +106,34 @@ pub(crate) fn radau_modify(diag: &mut [f64], off_diag_sq: &[f64], x0: f64) { // The ratio r_k = p_k / p_{k-1} satisfies: // r_1 = x0 - α_0 // r_k = x0 - α_{k-1} - off_diag_sq[k-2] / r_{k-1} for k >= 2 + // A vanishing ratio r means the characteristic polynomial p_{k-1}(x0) = 0, + // which would divide by zero and produce a NaN/Inf diagonal that corrupts + // the rule; surface it as a clear error instead. For the classical families + // with x0 = ±1 this never occurs in practice. + #[allow(clippy::float_cmp)] + fn vanished(r: f64) -> bool { + r == 0.0 + } let mut r = x0 - diag[0]; // r_1 for k in 2..n { + if vanished(r) { + return Err(QuadratureError::InvalidInput( + "Radau modification: characteristic polynomial vanished at the endpoint", + )); + } r = x0 - diag[k - 1] - off_diag_sq[k - 2] / r; } + if vanished(r) { + return Err(QuadratureError::InvalidInput( + "Radau modification: characteristic polynomial vanished at the endpoint", + )); + } // Now r = r_{n-1} = p_{n-1}(x0) / p_{n-2}(x0) // Choose α_{n-1} so that p_n(x0) = 0: // p_n = (x0 - α_{n-1}) p_{n-1} - β_{n-1} p_{n-2} = 0 // α_{n-1} = x0 - β_{n-1} / r diag[n - 1] = x0 - off_diag_sq[n - 2] / r; + Ok(()) } /// Symmetric tridiagonal eigenvalue decomposition via the implicit QL @@ -199,6 +232,23 @@ fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) -> bool { mod tests { use super::*; + #[test] + fn negative_off_diagonal_rejected() { + // A negative squared off-diagonal (non-PSD Jacobi matrix) must error, + // not produce NaN nodes via sqrt and a misleading "did not converge". + assert!(golub_welsch(&[0.0, 0.0], &[-1.0], 2.0).is_err()); + } + + #[test] + fn radau_modify_degenerate_orders_ok() { + // n == 0 must not index diag[0]; n == 1 sets the single diagonal. + let mut empty: [f64; 0] = []; + assert!(radau_modify(&mut empty, &[], -1.0).is_ok()); + let mut one = [5.0]; + assert!(radau_modify(&mut one, &[], -1.0).is_ok()); + assert_eq!(one[0], -1.0); + } + /// Golub-Welsch with Legendre recurrence should recover Gauss-Legendre. #[test] fn legendre_recovery() { @@ -242,7 +292,7 @@ mod tests { let off_diag_sq = vec![1.0 / 3.0]; // β_1 = 1/3 for Legendre let mu0 = 2.0; - radau_modify(&mut diag, &off_diag_sq, -1.0); + radau_modify(&mut diag, &off_diag_sq, -1.0).unwrap(); let (nodes, weights) = golub_welsch(&diag, &off_diag_sq, mu0).unwrap(); assert!((nodes[0] - (-1.0)).abs() < 1e-14); diff --git a/src/rule.rs b/src/rule.rs index b0aeaba..2badb88 100644 --- a/src/rule.rs +++ b/src/rule.rs @@ -67,6 +67,9 @@ impl QuadratureRule { where G: Fn(F) -> F, { + if n_panels == 0 { + return F::zero(); + } let n = F::from(n_panels).unwrap(); let panel_width = (b - a) / n; let mut total = F::zero(); @@ -81,8 +84,11 @@ impl QuadratureRule { /// Parallel composite integration over \[a, b\] with `n_panels` subintervals. /// - /// Identical to [`integrate_composite`](Self::integrate_composite) but evaluates - /// panels in parallel using rayon. Requires `F: Send + Sync` and `f: Fn(F) -> F + Sync`. + /// Equivalent to [`integrate_composite`](Self::integrate_composite) up to + /// floating-point summation order (rayon's tree reduction sums panels in a + /// different order, so results can differ in the last ULPs), but evaluates + /// panels in parallel using rayon. Requires `F: Send + Sync` and + /// `f: Fn(F) -> F + Sync`. /// /// # Panics /// @@ -97,6 +103,9 @@ impl QuadratureRule { { use rayon::prelude::*; + if n_panels == 0 { + return F::zero(); + } let n = F::from(n_panels).unwrap(); let panel_width = (b - a) / n; @@ -110,3 +119,16 @@ impl QuadratureRule { .reduce(F::zero, |a, b| a + b) } } + +#[cfg(test)] +mod tests { + use crate::GaussLegendre; + + #[test] + fn zero_panels_returns_zero() { + // n_panels == 0 must short-circuit instead of computing (b - a) / 0. + let gl = GaussLegendre::new(5).unwrap(); + let r = gl.rule().integrate_composite(0.0, 1.0, 0, |x: f64| x); + assert_eq!(r, 0.0); + } +} diff --git a/src/weighted.rs b/src/weighted.rs index 35c56c2..b8cf702 100644 --- a/src/weighted.rs +++ b/src/weighted.rs @@ -71,6 +71,10 @@ impl WeightedIntegrator { } /// Set the quadrature order. + /// + /// Note that `n == 0` is not rejected here (unlike [`Self::new`]); a zero + /// order will instead panic when [`Self::integrate`] builds the underlying + /// rule. #[must_use] pub fn with_order(mut self, n: usize) -> Self { self.order = n; @@ -87,10 +91,10 @@ impl WeightedIntegrator { /// /// # Panics /// - /// Panics if the [`WeightFunction::Jacobi`] parameters are invalid - /// (`alpha <= -1` or `beta <= -1`) or if the [`WeightFunction::Laguerre`] - /// parameter is invalid (`alpha <= -1`), since the underlying quadrature - /// rule construction will fail. + /// Panics if the order is zero (set via [`Self::with_order`]), or if the + /// [`WeightFunction::Jacobi`] parameters are invalid (`alpha <= -1` or + /// `beta <= -1`) or the [`WeightFunction::Laguerre`] parameter is invalid + /// (`alpha <= -1`), since the underlying quadrature rule construction fails. pub fn integrate(&self, f: G) -> f64 where G: Fn(f64) -> f64,