Fix quadrature correctness and bound resource use across the library

src/adaptive.rs

@@ -208,6 +208,18 @@ impl AdaptiveIntegrator {
         let mut total_error = 0.0;
         let mut num_evals = 0;
 
+        // Seeding every (non-degenerate) break-point subinterval must fit within
+        // the eval budget; otherwise the documented max_evals contract would be
+        // silently overrun on the integrate_with_breaks path. Compared via
+        // division to avoid overflow in n_seed * evals_per_call.
+        #[allow(clippy::float_cmp)]
+        let n_seed = intervals.iter().filter(|(ia, ib)| ia != ib).count();
+        if n_seed > self.max_evals / evals_per_call {
+            return Err(QuadratureError::InvalidInput(
+                "max_evals is too small to seed all break-point subintervals",
+            ));
+        }
+
         // Seed the priority queue with initial intervals
         for &(ia, ib) in intervals {
             // Exact comparison is intentional: degenerate zero-width
@@ -239,6 +251,18 @@ impl AdaptiveIntegrator {
             // Bisect the worst interval
             let mid = 0.5 * (worst.a + worst.b);
 
+            // If the midpoint rounds to an endpoint the interval has shrunk to
+            // the floating-point resolution limit (typically at an integrable
+            // singularity); bisecting it cannot reduce the error. Drop it and
+            // refine the next-worst interval rather than spinning to max_evals
+            // 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.
+            #[allow(clippy::float_cmp)]
+            if mid == worst.a || mid == worst.b {
+                continue;
+            }
+
             let left = gk.integrate_detail(worst.a, mid, f);
             let right = gk.integrate_detail(mid, worst.b, f);
             num_evals += 2 * evals_per_call;
@@ -479,6 +503,66 @@ mod tests {
         assert!(r.is_err());
     }
 
+    /// Too many break-point subintervals for the eval budget is rejected,
+    /// rather than silently overrunning max_evals during seeding.
+    #[test]
+    fn seed_budget_rejected_when_too_many_breaks() {
+        let breaks: Vec<f64> = (1..200).map(|i| f64::from(i) / 200.0).collect();
+        let r = AdaptiveIntegrator::default()
+            .with_max_evals(100)
+            .integrate_with_breaks(0.0, 1.0, &breaks, |x| x);
+        assert!(r.is_err());
+    }
+
+    /// An integrable endpoint singularity forces refinement toward the limit of
+    /// float resolution. The min-width guard drops unrefineable panels instead
+    /// of spinning to max_evals, so the call terminates with a finite value.
+    #[test]
+    fn singular_integrand_terminates() {
+        let r = AdaptiveIntegrator::default()
+            .with_abs_tol(1e-13)
+            .with_rel_tol(1e-13)
+            .integrate(
+                0.0,
+                1.0,
+                |x: f64| if x > 0.0 { 1.0 / x.sqrt() } else { 0.0 },
+            )
+            .unwrap();
+        assert!(r.value.is_finite());
+        assert!((r.value - 2.0).abs() < 0.05, "value={}", r.value);
+    }
+
+    /// Reversed bounds that require subdivision must still negate the forward
+    /// result. A peaked integrand forces the refinement loop (and the min-width
+    /// guard) to run, where an orientation-dependent guard would wrongly drop
+    /// every reversed (a > b) interval.
+    #[test]
+    fn reversed_bounds_needing_refinement_negate() {
+        let f = |x: f64| 1.0 / (1.0 + (x - 0.3).powi(2) * 1.0e4);
+        let fwd = AdaptiveIntegrator::default()
+            .with_abs_tol(1e-10)
+            .with_rel_tol(1e-10)
+            .integrate(0.0, 1.0, f)
+            .unwrap();
+        let rev = AdaptiveIntegrator::default()
+            .with_abs_tol(1e-10)
+            .with_rel_tol(1e-10)
+            .integrate(1.0, 0.0, f)
+            .unwrap();
+        assert!(
+            fwd.converged && rev.converged,
+            "fwd.conv={}, rev.conv={}",
+            fwd.converged,
+            rev.converged
+        );
+        assert!(
+            (fwd.value + rev.value).abs() < 1e-9,
+            "fwd={}, rev={}",
+            fwd.value,
+            rev.value
+        );
+    }
+
     /// Builder with higher-order pair achieves tight tolerance.
     #[test]
     fn builder_high_order() {

src/cubature/adaptive.rs

@@ -115,6 +115,18 @@ impl AdaptiveCubature {
                 .integrate(lower[0], upper[0], |x: f64| f(&[x]));
         }
 
+        // A single Genz-Malik region costs 2^d + 2d² + 2d + 1 evaluations, all
+        // spent before max_evals is next checked. Refuse if even one region
+        // cannot fit the budget, so a high dimension with a small budget fails
+        // cleanly instead of spending ~2^d evaluations (a de-facto hang at large
+        // d). 1 << d is safe: d <= 30 was enforced above.
+        let region_evals = (1usize << d) + 2 * d * d + 2 * d + 1;
+        if region_evals > self.max_evals {
+            return Err(QuadratureError::InvalidInput(
+                "max_evals is smaller than one Genz-Malik region; reduce the dimension or raise max_evals",
+            ));
+        }
+
         let mut heap: BinaryHeap<SubRegion> = BinaryHeap::new();
         let mut total_evals = 0usize;
 
@@ -315,6 +327,11 @@ where
     let mut sum_3 = 0.0;
     let mut fourth_diffs = vec![0.0; d]; // for split axis selection
 
+    // Genz-Malik fourth divided difference: the λ₂²/λ₄² ratio cancels the
+    // second-derivative term so the axis with the largest fourth-derivative
+    // (non-polynomial) variation is chosen. λ₂²/λ₄² = (9/70)/(9/10) = 1/7.
+    let l2_over_l4_sq = lambda2 * lambda2 / (lambda4 * lambda4);
+
     for i in 0..d {
         x[i] = midpoints[i] + lambda2 * half_widths[i];
         let f_plus2 = f(&x);
@@ -328,9 +345,12 @@ where
         let f_minus4 = f(&x);
         sum_3 += f_plus4 + f_minus4;
 
-        // Fourth difference for split axis selection
-        fourth_diffs[i] = (f_plus2 + f_minus2 - 2.0 * f_center).abs()
-            - (f_plus4 + f_minus4 - 2.0 * f_center).abs();
+        // Fourth divided difference for split-axis selection: a single absolute
+        // value of the ratio-weighted second differences (NOT a difference of
+        // two absolute values, which cancels no derivative).
+        let delta2 = f_plus2 + f_minus2 - 2.0 * f_center;
+        let delta4 = f_plus4 + f_minus4 - 2.0 * f_center;
+        fourth_diffs[i] = (delta2 - l2_over_l4_sq * delta4).abs();
 
         x[i] = midpoints[i]; // restore
         num_evals += 4;
@@ -418,6 +438,27 @@ mod tests {
         assert!(adaptive_cubature(|_| 1.0, &lower, &upper, 1e-8).is_err());
     }
 
+    #[test]
+    fn high_dim_small_budget_rejected() {
+        // A single Genz-Malik region at d=20 costs 2^20 ≈ 1.05e6 evals, above
+        // this budget. It must error cleanly rather than spending ~2^d
+        // evaluations (a de-facto hang) before the budget is checked.
+        let lower = vec![0.0; 20];
+        let upper = vec![1.0; 20];
+        let r = AdaptiveCubature::default()
+            .with_max_evals(100_000)
+            .integrate(&lower, &upper, |_| 1.0);
+        assert!(r.is_err());
+    }
+
+    #[test]
+    fn split_axis_targets_the_varying_dimension() {
+        // f varies only along axis 0; the fourth-difference criterion must pick
+        // axis 0. The former |δ₂|-|δ₄| form picked a constant axis for convex f.
+        let detail = genz_malik_eval(3, &[0.0; 3], &[1.0; 3], &|x: &[f64]| (3.0 * x[0]).exp());
+        assert_eq!(detail.split_axis, 0, "split_axis={}", detail.split_axis);
+    }
+
     /// Constant function over [0,1]^2.
     #[test]
     fn constant_2d() {