diff --git a/src/adaptive.rs b/src/adaptive.rs
index fee8a5a..4a6b41b 100644
--- a/src/adaptive.rs
+++ b/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() {
diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs
index 9d5abc1..6a17d73 100644
--- a/src/cubature/adaptive.rs
+++ b/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() {
diff --git a/src/cubature/sobol.rs b/src/cubature/sobol.rs
index a2b2d7e..0798b36 100644
--- a/src/cubature/sobol.rs
+++ b/src/cubature/sobol.rs
@@ -1,6 +1,6 @@
 //! Sobol low-discrepancy sequence.
 //!
-//! Quasi-random sequence using direction numbers from Joe & Kuo (2010).
+//! Quasi-random sequence using direction numbers from Joe & Kuo (2008).
 //! Better uniformity than Halton sequences, especially in high dimensions.
 //!
 //! Uses gray-code enumeration for O(1) per-point generation.
@@ -16,7 +16,7 @@ const BITS: u32 = 32;
 /// Sobol sequence generator.
 ///
 /// Generates quasi-random points in \[0, 1)^d using direction numbers
-/// from Joe & Kuo (2010). The first dimension uses the Van der Corput
+/// from Joe & Kuo (2008). The first dimension uses the Van der Corput
 /// sequence (base 2), higher dimensions use primitive polynomials.
 ///
 /// # Example
@@ -108,10 +108,19 @@ impl SobolSequence {
     ///
     /// # Panics
     ///
-    /// Panics if `point.len()` is less than the sequence dimension.
+    /// Panics if `point.len()` is less than the sequence dimension, or if the
+    /// sequence is exhausted: the 32-bit construction has period 2³², so at
+    /// most 2³² − 1 points can be generated.
     pub fn next_point(&mut self, point: &mut [f64]) {
         assert!(point.len() >= self.dim);
         self.index += 1;
+        // Beyond the period, (index - 1) has BITS trailing ones, so the
+        // gray-code bit index would reach BITS and read past the direction
+        // numbers. Fail with a clear message rather than an out-of-bounds panic.
+        assert!(
+            self.index < (1u64 << BITS),
+            "Sobol sequence exhausted: at most 2^32 - 1 points are available"
+        );
 
         // Find the rightmost zero bit of (index - 1) (gray-code position)
         let c = (self.index - 1).trailing_ones() as usize;
@@ -130,7 +139,16 @@ impl SobolSequence {
     ///
     /// After calling `skip(n)`, the next call to `next_point` produces
     /// the (n+1)-th point.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `n >= 2³² − 1`: the gray-code reconstruction only consumes the
+    /// low 32 bits, and the subsequent `next_point` must stay within the period.
     pub fn skip(&mut self, n: u64) {
+        assert!(
+            n < (1u64 << BITS) - 1,
+            "Sobol sequence exhausted: at most 2^32 - 1 points are available"
+        );
         let b = BITS as usize;
         // Reset state
         self.state.fill(0);
@@ -174,206 +192,216 @@ struct SobolEntry {
     m: [u32; 8],
 }
 
-/// Direction numbers from Joe & Kuo (2010) for dimensions 2..40.
-/// Each entry: (degree s, polynomial coefficients a, initial m-values).
+/// Direction numbers from Joe & Kuo (2008) for dimensions 2..=40.
 ///
-/// Source: <https://web.maths.unsw.edu.au/~fkuo/sobol/joe-kuo-old.1111>
+/// `SOBOL_TABLE[k]` holds the parameters for spatial dimension `k + 2` (the
+/// first dimension is the parameter-free Van der Corput sequence). Each entry
+/// is `(degree s, polynomial coefficient encoding a, initial m-values m_1..m_s)`
+/// for a distinct primitive polynomial over GF(2), ordered by increasing degree;
+/// every `m_i` is odd with `m_i < 2^i`. There is exactly one primitive polynomial
+/// of degree 1, so the degree-1 block has a single entry — two degree-1 rows would
+/// produce identical direction numbers and collapse two coordinates onto each other.
+///
+/// Source: `new-joe-kuo-6.21201` from the authors' repository
+/// <https://github.com/joe-kuo/sobol_data> (Joe & Kuo, "Constructing Sobol
+/// sequences with better two-dimensional projections", SIAM J. Sci. Comput.
+/// 30 (2008) 2635-2654).
 static SOBOL_TABLE: [SobolEntry; 39] = [
     SobolEntry {
         degree: 1,
         coeffs: 0,
         m: [1, 0, 0, 0, 0, 0, 0, 0],
-    },
-    SobolEntry {
-        degree: 1,
-        coeffs: 1,
-        m: [1, 0, 0, 0, 0, 0, 0, 0],
-    },
+    }, // d=2
     SobolEntry {
         degree: 2,
         coeffs: 1,
-        m: [1, 1, 0, 0, 0, 0, 0, 0],
-    },
+        m: [1, 3, 0, 0, 0, 0, 0, 0],
+    }, // d=3
     SobolEntry {
         degree: 3,
         coeffs: 1,
         m: [1, 3, 1, 0, 0, 0, 0, 0],
-    },
+    }, // d=4
     SobolEntry {
         degree: 3,
         coeffs: 2,
         m: [1, 1, 1, 0, 0, 0, 0, 0],
-    },
+    }, // d=5
     SobolEntry {
         degree: 4,
         coeffs: 1,
         m: [1, 1, 3, 3, 0, 0, 0, 0],
-    },
+    }, // d=6
     SobolEntry {
         degree: 4,
         coeffs: 4,
         m: [1, 3, 5, 13, 0, 0, 0, 0],
-    },
+    }, // d=7
     SobolEntry {
         degree: 5,
         coeffs: 2,
         m: [1, 1, 5, 5, 17, 0, 0, 0],
-    },
+    }, // d=8
     SobolEntry {
         degree: 5,
         coeffs: 4,
         m: [1, 1, 5, 5, 5, 0, 0, 0],
-    },
+    }, // d=9
     SobolEntry {
         degree: 5,
         coeffs: 7,
         m: [1, 1, 7, 11, 19, 0, 0, 0],
-    },
+    }, // d=10
     SobolEntry {
         degree: 5,
         coeffs: 11,
         m: [1, 1, 5, 1, 1, 0, 0, 0],
-    },
+    }, // d=11
     SobolEntry {
         degree: 5,
         coeffs: 13,
         m: [1, 1, 1, 3, 11, 0, 0, 0],
-    },
+    }, // d=12
     SobolEntry {
         degree: 5,
         coeffs: 14,
         m: [1, 3, 5, 5, 31, 0, 0, 0],
-    },
+    }, // d=13
     SobolEntry {
         degree: 6,
         coeffs: 1,
         m: [1, 3, 3, 9, 7, 49, 0, 0],
-    },
+    }, // d=14
     SobolEntry {
         degree: 6,
         coeffs: 13,
         m: [1, 1, 1, 15, 21, 21, 0, 0],
-    },
+    }, // d=15
     SobolEntry {
         degree: 6,
         coeffs: 16,
         m: [1, 3, 1, 13, 27, 49, 0, 0],
-    },
+    }, // d=16
     SobolEntry {
-        degree: 7,
+        degree: 6,
         coeffs: 19,
-        m: [1, 1, 1, 15, 7, 5, 127, 0],
-    },
+        m: [1, 1, 1, 15, 7, 5, 0, 0],
+    }, // d=17
     SobolEntry {
-        degree: 7,
+        degree: 6,
         coeffs: 22,
-        m: [1, 3, 3, 5, 19, 33, 65, 0],
-    },
+        m: [1, 3, 1, 15, 13, 25, 0, 0],
+    }, // d=18
     SobolEntry {
-        degree: 7,
+        degree: 6,
         coeffs: 25,
-        m: [1, 3, 7, 11, 29, 17, 85, 0],
-    },
+        m: [1, 1, 5, 5, 19, 61, 0, 0],
+    }, // d=19
     SobolEntry {
         degree: 7,
-        coeffs: 37,
-        m: [1, 1, 3, 7, 23, 55, 41, 0],
-    },
+        coeffs: 1,
+        m: [1, 3, 7, 11, 23, 15, 103, 0],
+    }, // d=20
     SobolEntry {
         degree: 7,
-        coeffs: 41,
-        m: [1, 3, 5, 1, 15, 17, 63, 0],
-    },
+        coeffs: 4,
+        m: [1, 3, 7, 13, 13, 15, 69, 0],
+    }, // d=21
     SobolEntry {
         degree: 7,
-        coeffs: 50,
-        m: [1, 1, 7, 9, 31, 29, 17, 0],
-    },
+        coeffs: 7,
+        m: [1, 1, 3, 13, 7, 35, 63, 0],
+    }, // d=22
     SobolEntry {
         degree: 7,
-        coeffs: 55,
-        m: [1, 3, 7, 7, 21, 61, 119, 0],
-    },
+        coeffs: 8,
+        m: [1, 3, 5, 9, 1, 25, 53, 0],
+    }, // d=23
     SobolEntry {
         degree: 7,
-        coeffs: 59,
-        m: [1, 1, 5, 3, 5, 49, 89, 0],
-    },
+        coeffs: 14,
+        m: [1, 3, 1, 13, 9, 35, 107, 0],
+    }, // d=24
     SobolEntry {
         degree: 7,
-        coeffs: 62,
-        m: [1, 3, 1, 1, 11, 3, 117, 0],
-    },
+        coeffs: 19,
+        m: [1, 3, 1, 5, 27, 61, 31, 0],
+    }, // d=25
     SobolEntry {
-        degree: 8,
-        coeffs: 14,
-        m: [1, 3, 3, 3, 15, 17, 17, 193],
-    },
+        degree: 7,
+        coeffs: 21,
+        m: [1, 1, 5, 11, 19, 41, 61, 0],
+    }, // d=26
     SobolEntry {
-        degree: 8,
-        coeffs: 52,
-        m: [1, 1, 3, 7, 31, 29, 67, 45],
-    },
+        degree: 7,
+        coeffs: 28,
+        m: [1, 3, 5, 3, 3, 13, 69, 0],
+    }, // d=27
     SobolEntry {
-        degree: 8,
-        coeffs: 56,
-        m: [1, 3, 1, 7, 23, 59, 55, 165],
-    },
+        degree: 7,
+        coeffs: 31,
+        m: [1, 1, 7, 13, 1, 19, 1, 0],
+    }, // d=28
     SobolEntry {
-        degree: 8,
-        coeffs: 67,
-        m: [1, 3, 5, 9, 21, 37, 7, 51],
-    },
+        degree: 7,
+        coeffs: 32,
+        m: [1, 3, 7, 5, 13, 19, 59, 0],
+    }, // d=29
     SobolEntry {
-        degree: 8,
-        coeffs: 69,
-        m: [1, 1, 7, 5, 17, 13, 81, 251],
-    },
+        degree: 7,
+        coeffs: 37,
+        m: [1, 1, 3, 9, 25, 29, 41, 0],
+    }, // d=30
     SobolEntry {
-        degree: 8,
-        coeffs: 70,
-        m: [1, 1, 3, 15, 29, 47, 49, 143],
-    },
+        degree: 7,
+        coeffs: 41,
+        m: [1, 3, 5, 13, 23, 1, 55, 0],
+    }, // d=31
     SobolEntry {
-        degree: 8,
-        coeffs: 79,
-        m: [1, 3, 7, 3, 21, 39, 29, 45],
-    },
+        degree: 7,
+        coeffs: 42,
+        m: [1, 3, 7, 3, 13, 59, 17, 0],
+    }, // d=32
     SobolEntry {
-        degree: 8,
-        coeffs: 81,
-        m: [1, 1, 5, 7, 29, 7, 37, 67],
-    },
+        degree: 7,
+        coeffs: 50,
+        m: [1, 3, 1, 3, 5, 53, 69, 0],
+    }, // d=33
     SobolEntry {
-        degree: 8,
-        coeffs: 84,
-        m: [1, 1, 5, 5, 11, 57, 97, 175],
-    },
+        degree: 7,
+        coeffs: 55,
+        m: [1, 1, 5, 5, 23, 33, 13, 0],
+    }, // d=34
     SobolEntry {
-        degree: 8,
-        coeffs: 87,
-        m: [1, 3, 3, 13, 29, 25, 29, 15],
-    },
+        degree: 7,
+        coeffs: 56,
+        m: [1, 1, 7, 7, 1, 61, 123, 0],
+    }, // d=35
     SobolEntry {
-        degree: 8,
-        coeffs: 90,
-        m: [1, 1, 7, 7, 31, 37, 105, 189],
-    },
+        degree: 7,
+        coeffs: 59,
+        m: [1, 1, 7, 9, 13, 61, 49, 0],
+    }, // d=36
+    SobolEntry {
+        degree: 7,
+        coeffs: 62,
+        m: [1, 3, 3, 5, 3, 55, 33, 0],
+    }, // d=37
     SobolEntry {
         degree: 8,
-        coeffs: 97,
-        m: [1, 3, 5, 3, 25, 37, 5, 187],
-    },
+        coeffs: 14,
+        m: [1, 3, 1, 15, 31, 13, 49, 245],
+    }, // d=38
     SobolEntry {
         degree: 8,
-        coeffs: 103,
-        m: [1, 3, 1, 7, 21, 43, 81, 37],
-    },
+        coeffs: 21,
+        m: [1, 3, 5, 15, 31, 59, 63, 97],
+    }, // d=39
     SobolEntry {
         degree: 8,
-        coeffs: 115,
-        m: [1, 1, 3, 5, 1, 45, 47, 77],
-    },
+        coeffs: 22,
+        m: [1, 3, 1, 11, 11, 11, 77, 249],
+    }, // d=40
 ];
 
 #[cfg(test)]
@@ -431,6 +459,15 @@ mod tests {
         assert!(SobolSequence::new(41).is_err());
     }
 
+    #[test]
+    #[should_panic(expected = "exhausted")]
+    fn skip_past_period_panics() {
+        // The 32-bit construction has period 2^32; skipping to the boundary must
+        // panic with a clear message rather than silently aliasing high bits.
+        let mut sob = SobolSequence::new(2).unwrap();
+        sob.skip((1u64 << BITS) - 1);
+    }
+
     #[test]
     fn high_dim_points_in_unit_cube() {
         let dim = 20;
@@ -463,23 +500,126 @@ mod tests {
         }
     }
 
+    /// The first emitted points must match the unscrambled Joe-Kuo sequence
+    /// (`new-joe-kuo-6.21201`). bilby advances `index` before emitting, so it
+    /// drops the index-0 zero point: emitted point `n` equals the reference
+    /// sequence's index-`n` point. Reference oracle: SciPy `stats.qmc.Sobol`
+    /// with `scramble=false`, which uses the same 6.21201 direction numbers.
+    ///
+    /// This anchors the dimension-2 direction numbers to an external oracle; it
+    /// does NOT guard against the coordinate-collapse regression (which lived in
+    /// a later table entry) — that is covered by `all_dimensions_pairwise_distinct`
+    /// and `direction_table_structurally_valid`.
     #[test]
-    fn high_dim_constant_integral() {
-        // Quasi-Monte Carlo estimate of ∫_{[0,1]^d} 1 dx = 1
-        // Using Sobol with N points, the average of f=1 should be 1.
-        let dim = 20;
-        let n = 1024;
+    fn dim2_matches_reference_sequence() {
+        let expected = [
+            [0.5, 0.5],
+            [0.75, 0.25],
+            [0.25, 0.75],
+            [0.375, 0.375],
+            [0.875, 0.875],
+        ];
+        let mut sob = SobolSequence::new(2).unwrap();
+        let mut pt = [0.0; 2];
+        for (n, exp) in expected.iter().enumerate() {
+            sob.next_point(&mut pt);
+            for j in 0..2 {
+                assert!(
+                    (pt[j] - exp[j]).abs() < 1e-12,
+                    "point {}, dim {j}: {} vs {}",
+                    n + 1,
+                    pt[j],
+                    exp[j]
+                );
+            }
+        }
+    }
+
+    /// No two coordinates may follow the same sequence. Identical direction
+    /// numbers for two dimensions collapse their joint projection onto the
+    /// diagonal and silently bias any integrand that couples them.
+    #[test]
+    fn all_dimensions_pairwise_distinct() {
+        let dim = MAX_DIM;
+        let n = 256;
+        let mut sob = SobolSequence::new(dim).unwrap();
+        let mut cols: Vec<Vec<f64>> = (0..dim).map(|_| Vec::with_capacity(n)).collect();
+        let mut pt = vec![0.0; dim];
+        for _ in 0..n {
+            sob.next_point(&mut pt);
+            for (j, &x) in pt.iter().enumerate() {
+                cols[j].push(x);
+            }
+        }
+        for a in 0..dim {
+            for b in (a + 1)..dim {
+                assert!(cols[a] != cols[b], "dimensions {a} and {b} are identical");
+            }
+        }
+    }
+
+    /// The table must encode a valid set of distinct primitive polynomials over
+    /// GF(2): degrees non-decreasing, the per-degree counts matching the number
+    /// of primitive polynomials of each degree (1,1,2,2,6,6,18 for degrees 1-7,
+    /// plus the first 3 of the 16 degree-8 polynomials for dims 38-40), distinct
+    /// (degree, coeffs) pairs, and every initial direction number `m_i` odd with
+    /// `m_i < 2^i`. A duplicate or out-of-range entry silently collapses or
+    /// degrades the sequence.
+    #[test]
+    fn direction_table_structurally_valid() {
+        let mut counts = [0usize; 9]; // indexed by degree 1..=8
+        let mut prev_degree = 0u32;
+        for (k, e) in SOBOL_TABLE.iter().enumerate() {
+            assert!(e.degree >= prev_degree, "degree decreased at entry {k}");
+            prev_degree = e.degree;
+            let s = e.degree as usize;
+            assert!((1..=8).contains(&s), "degree {s} out of range at entry {k}");
+            counts[s] += 1;
+            for i in 0..s {
+                let m = e.m[i];
+                assert!(m % 2 == 1, "m[{i}]={m} not odd at entry {k}");
+                assert!(
+                    m < (1u32 << (i + 1)),
+                    "m[{i}]={m} >= 2^{} at entry {k}",
+                    i + 1
+                );
+            }
+            for i in s..8 {
+                assert_eq!(e.m[i], 0, "unused m[{i}] nonzero at entry {k}");
+            }
+        }
+        assert_eq!(&counts[1..=8], &[1, 1, 2, 2, 6, 6, 18, 3]);
+
+        // (degree, coeffs) pairs must be pairwise distinct.
+        for (k, e) in SOBOL_TABLE.iter().enumerate() {
+            for (l, f) in SOBOL_TABLE.iter().enumerate().skip(k + 1) {
+                assert!(
+                    !(e.degree == f.degree && e.coeffs == f.coeffs),
+                    "duplicate polynomial at entries {k} and {l}"
+                );
+            }
+        }
+    }
+
+    /// QMC estimate of `∫_{[0,1]^d} x_1 * x_2 dx` (0-indexed dims 1, 2) = 1/4.
+    /// If two dimensions collapse to the same sequence, the estimate drifts
+    /// toward `E[x^2] = 1/3` instead. (Replaces a former f≡1 test that never
+    /// inspected the generated points.)
+    #[test]
+    fn sobol_couples_dimensions() {
+        let dim = 5;
+        let n = 4096;
         let mut sob = SobolSequence::new(dim).unwrap();
         let mut pt = vec![0.0; dim];
         let mut sum = 0.0;
         for _ in 0..n {
             sob.next_point(&mut pt);
-            sum += 1.0; // f(x) = 1 for all x
+            sum += pt[1] * pt[2];
         }
         let estimate = sum / n as f64;
         assert!(
-            (estimate - 1.0).abs() < 1e-14,
-            "estimate={estimate}, expected 1.0"
+            (estimate - 0.25).abs() < 1e-3,
+            "estimate={estimate}, expected 0.25"
         );
     }
 
diff --git a/src/cubature/sparse_grid.rs b/src/cubature/sparse_grid.rs
index 801cc50..f943e34 100644
--- a/src/cubature/sparse_grid.rs
+++ b/src/cubature/sparse_grid.rs
@@ -55,13 +55,29 @@ impl SparseGrid {
     ///
     /// # Errors
     ///
-    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero.
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero, if `dim` or
+    /// `level` exceeds the supported maximum (32 and 12 respectively), or if the
+    /// resulting grid would exceed roughly one million points.
     pub fn new(dim: usize, level: usize, _basis: SparseGridBasis) -> Result<Self, QuadratureError> {
         if dim == 0 {
             return Err(QuadratureError::InvalidInput("dimension must be >= 1"));
         }
+        // Bound the inputs so the integer arithmetic stays safe (the `1 << level`
+        // in cc_order and the binomial coefficients cannot overflow) and a single
+        // 1-D Clenshaw-Curtis rule stays small. The combined point count is
+        // bounded separately inside build_smolyak.
+        if dim > MAX_DIM {
+            return Err(QuadratureError::InvalidInput(
+                "sparse grid dimension exceeds the supported maximum (32)",
+            ));
+        }
+        if level > MAX_LEVEL {
+            return Err(QuadratureError::InvalidInput(
+                "sparse grid level exceeds the supported maximum (12)",
+            ));
+        }
 
-        let rule = build_smolyak(dim, level);
+        let rule = build_smolyak(dim, level)?;
         Ok(Self { rule, level })
     }
 
@@ -139,7 +155,28 @@ fn quantise(x: f64) -> i64 {
 ///   `Q_{q,d}` = Σ_{max(q-d+1,0) ≤ |l|-d ≤ q-d} (-1)^{q-|l|+d} C(d-1, q-|l|+d) · (`Q_{l_1}` ⊗ ... ⊗ `Q_{l_d}`)
 ///
 /// Using 0-indexed levels where `l_j` ≥ 0.
-fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
+/// Maximum sparse-grid dimension (keeps the binomial coefficients within usize).
+const MAX_DIM: usize = 32;
+/// Maximum sparse-grid level. The finest 1-D Clenshaw-Curtis rule has
+/// `2^level + 1` points and is built in O(n²), so this also bounds the
+/// precompute cost (`cc_order(12) = 4097`).
+const MAX_LEVEL: usize = 12;
+/// Maximum cumulative sparse-grid tensor-point count before construction is
+/// rejected. Bounds the total inner-loop work, the outer multi-index
+/// enumeration (each contributes >= 1), and the memory.
+const MAX_SPARSE_POINTS: usize = 1_000_000;
+
+fn build_smolyak(dim: usize, level: usize) -> Result<CubatureRule, QuadratureError> {
+    // Reject up front if the number of Smolyak multi-indices, C(level+dim, dim),
+    // already exceeds the budget: each contributes at least one point, so the
+    // grid is infeasible, and this avoids the expensive enumeration entirely.
+    // The level/dim caps keep this binomial within usize.
+    if binomial(level + dim, dim) > MAX_SPARSE_POINTS {
+        return Err(QuadratureError::InvalidInput(
+            "sparse grid point count exceeds the supported maximum (1e6)",
+        ));
+    }
+
     // Precompute CC rules for each level we'll need
     let max_level = level;
     let cc_rules: Vec<(Vec<f64>, Vec<f64>)> =
@@ -167,6 +204,11 @@ fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
     let sum_min = (q + 1).saturating_sub(dim);
     let sum_max = q;
 
+    // Running count of tensor-point insertions across all multi-indices. Bailing
+    // when it exceeds the budget bounds the inner-loop work, the number of
+    // multi-indices enumerated (each adds >= 1), and the memory in one check.
+    let mut cumulative = 0usize;
+
     for s in sum_min..=sum_max {
         let diff = q - s;
         // Binomial coefficients for sparse grid levels are small enough to fit in f64.
@@ -187,6 +229,17 @@ fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
             let orders: Vec<usize> = multi_idx.iter().map(|&l| cc_order(l)).collect();
             let total: usize = orders.iter().product();
 
+            // Reject infeasibly large grids (the curse of dimensionality) rather
+            // than hanging in the inner loop or exhausting memory. With level <=
+            // MAX_LEVEL each factor is small enough that `total` (and the running
+            // sum, which bails at <= 2*MAX_SPARSE_POINTS) cannot overflow.
+            cumulative += total;
+            if cumulative > MAX_SPARSE_POINTS {
+                return Err(QuadratureError::InvalidInput(
+                    "sparse grid point count exceeds the supported maximum (1e6)",
+                ));
+            }
+
             let mut indices = vec![0usize; dim];
             for _ in 0..total {
                 let mut w = coeff;
@@ -248,7 +301,7 @@ fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
         weights.push(w);
     }
 
-    CubatureRule::new(nodes_flat, weights, dim)
+    Ok(CubatureRule::new(nodes_flat, weights, dim))
 }
 
 /// Generate the next weak composition of `s` into `d` non-negative parts.
@@ -304,6 +357,22 @@ fn binomial(n: usize, k: usize) -> usize {
 mod tests {
     use super::*;
 
+    #[test]
+    fn extreme_level_and_dim_rejected() {
+        // Beyond the supported caps, construction must error rather than
+        // overflow `1 << level` or build a gigantic single 1-D rule.
+        assert!(SparseGrid::clenshaw_curtis(3, 64).is_err());
+        assert!(SparseGrid::clenshaw_curtis(100, 3).is_err());
+    }
+
+    #[test]
+    fn infeasible_point_count_rejected() {
+        // Within the level/dim caps, but the combined point count is far too
+        // large (curse of dimensionality). Must reject (quickly) rather than
+        // exhausting memory or hanging in the multi-index enumeration.
+        assert!(SparseGrid::clenshaw_curtis(32, 6).is_err());
+    }
+
     #[test]
     fn level_0_single_point() {
         let sg = SparseGrid::clenshaw_curtis(3, 0).unwrap();
diff --git a/src/cubature/tensor.rs b/src/cubature/tensor.rs
index 013a1ea..46c4783 100644
--- a/src/cubature/tensor.rs
+++ b/src/cubature/tensor.rs
@@ -51,9 +51,22 @@ impl TensorProductRule {
 
         let dim = rules_1d.len();
         let orders: Vec<usize> = rules_1d.iter().map(|r| r.order()).collect();
-        let total_points: usize = orders.iter().product();
-
-        let mut nodes_flat = Vec::with_capacity(total_points * dim);
+        // Checked product: a wraparound here would otherwise yield a wrong (often
+        // empty) rule returned as Ok — e.g. isotropic(order-2, 64) overflows to 0
+        // and would integrate everything to zero.
+        let total_points = orders
+            .iter()
+            .try_fold(1usize, |acc, &o| acc.checked_mul(o))
+            .ok_or(QuadratureError::InvalidInput(
+                "tensor-product point count overflows usize",
+            ))?;
+        let flat_len = total_points
+            .checked_mul(dim)
+            .ok_or(QuadratureError::InvalidInput(
+                "tensor-product node array length overflows usize",
+            ))?;
+
+        let mut nodes_flat = Vec::with_capacity(flat_len);
         let mut weights = Vec::with_capacity(total_points);
         let mut indices = vec![0usize; dim];
 
@@ -126,6 +139,15 @@ mod tests {
         assert_eq!(tp.dim(), 2);
     }
 
+    #[test]
+    fn point_count_overflow_rejected() {
+        // An order-2 rule in 64 dims is 2^64 points, which overflows usize.
+        // It must error rather than wrapping to 0 and returning an empty rule
+        // that integrates everything to zero.
+        let gl = GaussLegendre::new(2).unwrap();
+        assert!(TensorProductRule::isotropic(gl.rule(), 64).is_err());
+    }
+
     #[test]
     fn isotropic_point_count() {
         let gl = GaussLegendre::new(4).unwrap();
diff --git a/src/gauss_jacobi.rs b/src/gauss_jacobi.rs
index 0bce676..1512de0 100644
--- a/src/gauss_jacobi.rs
+++ b/src/gauss_jacobi.rs
@@ -53,7 +53,8 @@ impl GaussJacobi {
     ///
     /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1`, `beta <= -1`,
-    /// or either parameter is non-finite.
+    /// either parameter is non-finite, or the weight integral `μ₀` overflows
+    /// `f64` (very large `alpha`/`beta`).
     pub fn new(n: usize, alpha: f64, beta: f64) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -105,48 +106,56 @@ fn compute_jacobi(
 ) -> Result<(Vec<f64>, Vec<f64>), QuadratureError> {
     let ab = alpha + beta;
 
-    // Diagonal: α_k = (β²-α²) / ((2k+ab)(2k+ab+2))
-    // Special handling when 2k+ab is near zero (only k=0 and ab≈0)
+    // Diagonal: a_k = (β²-α²) / ((2k+ab)(2k+ab+2)).
+    // For k = 0 the identity β+α = ab makes this (β-α)·ab / (ab(ab+2)) =
+    // (β-α)/(ab+2): pole-free for valid α,β > -1 (ab+2 > 0) and free of the
+    // β²-α² cancellation. For k ≥ 1, 2k+ab ≥ 2+ab > 0 so the denominator is
+    // strictly positive; the factored numerator (β-α)(β+α) is cancellation-free.
     let diag: Vec<f64> = (0..n)
         .map(|k| {
-            let two_k_ab = 2.0 * k as f64 + ab;
-            let denom = two_k_ab * (two_k_ab + 2.0);
-            if denom.abs() < 1e-300 {
-                // For ab=0: (β²-α²) = (β-α)(β+α) = 0 when α=β
-                // In general, use L'Hôpital: (β-α)/(ab+2) when k=0, ab→0
-                if k == 0 {
-                    (beta - alpha) / (ab + 2.0)
-                } else {
-                    0.0
-                }
+            if k == 0 {
+                (beta - alpha) / (ab + 2.0)
             } else {
-                (beta * beta - alpha * alpha) / denom
+                let two_k_ab = 2.0 * k as f64 + ab;
+                let denom = two_k_ab * (two_k_ab + 2.0);
+                (beta - alpha) * (beta + alpha) / denom
             }
         })
         .collect();
 
-    // Off-diagonal squared: β_k for k = 1, ..., n-1
-    // Special case: when k+α+β=0 and 2k+α+β-1=0 (i.e., k=1, α+β=-1),
-    // both numerator and denominator vanish. The limit is 2(1+α)(1+β).
+    // Off-diagonal squared: b_k for k = 1, ..., n-1, with
+    //   b_k = 4k(k+α)(k+β)(k+ab) / ((2k+ab)²(2k+ab+1)(2k+ab-1)).
+    // For k = 1 the (2k+ab-1) = (1+ab) denominator factor cancels the
+    // (k+ab) = (1+ab) numerator factor, giving b_1 = 4(1+α)(1+β)/((2+ab)²(3+ab)),
+    // which is regular at the α+β = -1 resonance (denominator 1·2 there) with no
+    // cancellation. For k ≥ 2, 2k+ab-1 ≥ 3+ab > 1 > 0, so the general form has a
+    // strictly positive denominator.
     let off_diag_sq: Vec<f64> = (1..n)
         .map(|k| {
-            let k = k as f64;
-            let two_k_ab = 2.0 * k + ab;
-            let denom = two_k_ab * two_k_ab * (two_k_ab + 1.0) * (two_k_ab - 1.0);
-            if denom.abs() < 1e-300 {
-                // 0/0 case: k=1, α+β=-1. Limit = 2(1+α)(1+β)
-                2.0 * (1.0 + alpha) * (1.0 + beta)
+            let kf = k as f64;
+            let two_k_ab = 2.0 * kf + ab;
+            if k == 1 {
+                4.0 * (1.0 + alpha) * (1.0 + beta) / (two_k_ab * two_k_ab * (two_k_ab + 1.0))
             } else {
-                let numer = 4.0 * k * (k + alpha) * (k + beta) * (k + ab);
+                let denom = two_k_ab * two_k_ab * (two_k_ab + 1.0) * (two_k_ab - 1.0);
+                let numer = 4.0 * kf * (kf + alpha) * (kf + beta) * (kf + ab);
                 numer / denom
             }
         })
         .collect();
 
-    // μ₀ = 2^(ab+1) Γ(α+1)Γ(β+1) / Γ(ab+2)
-    let mu0 = ((ab + 1.0) * core::f64::consts::LN_2 + ln_gamma(alpha + 1.0) + ln_gamma(beta + 1.0)
-        - ln_gamma(ab + 2.0))
-    .exp();
+    // μ₀ = 2^(ab+1) Γ(α+1)Γ(β+1) / Γ(ab+2), formed in log-space.
+    let ln_mu0 =
+        (ab + 1.0) * core::f64::consts::LN_2 + ln_gamma(alpha + 1.0) + ln_gamma(beta + 1.0)
+            - ln_gamma(ab + 2.0);
+    // Reject overflow (and the inf+inf-inf = NaN case for enormous α, β) rather
+    // than returning all-infinite weights with no diagnostic.
+    if !ln_mu0.is_finite() || ln_mu0 > f64::MAX.ln() {
+        return Err(QuadratureError::InvalidInput(
+            "Jacobi weight integral mu0 overflows f64 (alpha/beta too large)",
+        ));
+    }
+    let mu0 = ln_mu0.exp();
 
     golub_welsch(&diag, &off_diag_sq, mu0)
 }
@@ -239,6 +248,47 @@ mod tests {
         log_val.exp()
     }
 
+    /// Near-degenerate parameters exercise the recurrence-coefficient corners.
+    /// Low moments of the Jacobi weight are independent closed-form oracles:
+    ///   m₁ = μ₀·a₀,  m₂ = μ₀·(b₁ + a₀²),  with
+    ///   a₀ = (β-α)/(α+β+2),  b₁ = 4(1+α)(1+β)/((α+β+2)²(α+β+3)).
+    /// The rule's Σwᵢxᵢ and Σwᵢxᵢ² come from the Golub-Welsch eigensolve, a
+    /// wholly separate computation. The former `(β²-α²)` diagonal and general
+    /// off-diagonal forms lost precision near α+β=0 and α+β=-1, perturbing the
+    /// nodes/weights away from these moments.
+    #[test]
+    fn near_degenerate_moments_accurate() {
+        for (alpha, beta) in [(0.5, -0.5 + 1e-10), (-0.3, -0.7 + 1e-7)] {
+            let gj = GaussJacobi::new(8, alpha, beta).unwrap();
+            let mu0 = jacobi_integral(alpha, beta);
+            let ab = alpha + beta;
+            let a0 = (beta - alpha) / (ab + 2.0);
+            let b1 = 4.0 * (1.0 + alpha) * (1.0 + beta) / ((ab + 2.0).powi(2) * (ab + 3.0));
+            let m1: f64 = gj
+                .nodes()
+                .iter()
+                .zip(gj.weights())
+                .map(|(&x, &w)| x * w)
+                .sum();
+            let m2: f64 = gj
+                .nodes()
+                .iter()
+                .zip(gj.weights())
+                .map(|(&x, &w)| x * x * w)
+                .sum();
+            assert!(
+                (m1 - mu0 * a0).abs() <= 1e-12 * mu0,
+                "m1={m1}, expected {} (a={alpha}, b={beta})",
+                mu0 * a0
+            );
+            assert!(
+                (m2 - mu0 * (b1 + a0 * a0)).abs() <= 1e-12 * mu0,
+                "m2={m2}, expected {} (a={alpha}, b={beta})",
+                mu0 * (b1 + a0 * a0)
+            );
+        }
+    }
+
     #[test]
     fn nodes_sorted_and_bounded() {
         let gj = GaussJacobi::new(20, 1.0, 0.5).unwrap();
diff --git a/src/gauss_kronrod.rs b/src/gauss_kronrod.rs
index 37b64da..e57df13 100644
--- a/src/gauss_kronrod.rs
+++ b/src/gauss_kronrod.rs
@@ -137,11 +137,11 @@ impl GaussKronrod {
         (detail.estimate, detail.error)
     }
 
-    /// Internal detailed integration returning full QUADPACK output.
+    /// Internal detailed Gauss-Kronrod evaluation of a single panel.
     ///
-    /// Returns estimate, error, integral of |f|, and integral of |f - mean|
-    /// over the interval. The adaptive integrator uses these for roundoff
-    /// detection.
+    /// Returns the Kronrod estimate and the QUADPACK error estimate, which
+    /// scales the Gauss-Kronrod difference against `resasc` (the integral of
+    /// `|f - mean|` over the panel).
     #[allow(clippy::too_many_lines)] // single cohesive QUADPACK algorithm, splitting would obscure the logic
     #[allow(clippy::float_cmp)] // exact comparison for degenerate zero-width interval
     pub(crate) fn integrate_detail<G>(&self, a: f64, b: f64, f: G) -> GKDetail
@@ -161,7 +161,14 @@ impl GaussKronrod {
 
         let mut kronrod_sum = 0.0_f64;
         let mut gauss_sum = 0.0_f64;
-        let mut abs_sum = 0.0_f64;
+
+        // Function values per abscissa, retained for the resasc second pass
+        // below: resasc needs the Kronrod mean, which is only known after the
+        // full sum. The Kronrod order is at most 51 (G25-K51), so xgk has at
+        // most 26 entries — a fixed stack buffer avoids a per-panel allocation.
+        debug_assert!(n <= 32);
+        let mut fv1 = [0.0_f64; 32];
+        let mut fv2 = [0.0_f64; 32];
 
         // QUADPACK convention: xgk is stored from outermost to center.
         // Gauss nodes are at xgk[1], xgk[3], xgk[5], ... (odd indices, 0-based).
@@ -177,8 +184,8 @@ impl GaussKronrod {
             if xk == 0.0 {
                 // Center node
                 let fx = f(midpoint);
+                fv1[i] = fx;
                 kronrod_sum += wk * fx;
-                abs_sum += wk * fx.abs();
 
                 // Center is a Gauss node for odd Gauss order (G7, G15, G25)
                 if self.pair.gauss_order() % 2 == 1 && gauss_idx < self.wg.len() {
@@ -188,10 +195,11 @@ impl GaussKronrod {
                 // Symmetric pair ±xk
                 let f_pos = f(midpoint + half_width * xk);
                 let f_neg = f(midpoint - half_width * xk);
+                fv1[i] = f_pos;
+                fv2[i] = f_neg;
                 let f_sum = f_pos + f_neg;
 
                 kronrod_sum += wk * f_sum;
-                abs_sum += wk * (f_pos.abs() + f_neg.abs());
 
                 // Gauss nodes are at odd indices (0-based): 1, 3, 5, ...
                 if i % 2 == 1 && gauss_idx < self.wg.len() {
@@ -203,17 +211,32 @@ impl GaussKronrod {
 
         let estimate = half_width * kronrod_sum;
         let gauss_result = half_width * gauss_sum;
-        abs_sum *= half_width.abs();
-
-        // QUADPACK error estimation heuristic.
-        // This is an intentional simplification of the full QUADPACK error
-        // formula (which also considers |f - I/(b-a)| for roundoff detection).
-        // The simplified version omits the roundoff term but retains the
-        // core (200·δ/S)^1.5 scaling that prevents error underestimation
-        // for smooth integrands. See Piessens et al. (1983), §2.2.
+
+        // resasc = ∫|f − f̄| over the panel, where the Kronrod mean f̄ = reskh
+        // is half the unscaled Kronrod sum (the reference weights sum to 2).
+        // 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.
+        let reskh = kronrod_sum * 0.5;
+        let mut resasc = 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();
+            } else {
+                resasc += wk * ((fv1[i] - reskh).abs() + (fv2[i] - reskh).abs());
+            }
+        }
+        resasc *= 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
+        // law, capped at resasc itself (dmin1(1.0, …)) so the estimate can never
+        // exceed the panel's total variation. A constant integrand gives
+        // resasc ≈ 0 and δ ≈ 0, so the error is left at the vanishing δ.
         let mut error = (estimate - gauss_result).abs();
-        if abs_sum > 0.0 && error > 0.0 {
-            error = abs_sum * (200.0 * error / abs_sum).min(1.0).powf(1.5);
+        if resasc != 0.0 && error > 0.0 {
+            error = resasc * (200.0 * error / resasc).powf(1.5).min(1.0);
         }
 
         GKDetail { estimate, error }
@@ -645,6 +668,45 @@ mod tests {
         assert!(err < 1e-10, "x^12 error estimate = {err}");
     }
 
+    /// The error estimate scales against resasc = ∫|f − mean|, so a large
+    /// constant offset (integrated exactly by both rules) must not change it.
+    /// With the previous resabs (∫|f|) scaling, the offset deflated the error
+    /// by ~√(∫|f|), so the adaptive integrator converged prematurely on
+    /// offset-heavy integrands.
+    #[test]
+    fn error_estimate_independent_of_dc_offset() {
+        let gk = GaussKronrod::new(GKPair::G7K15);
+        let varying = gk.integrate_detail(0.0, 1.0, |x| (20.0 * x).sin());
+        let offset = gk.integrate_detail(0.0, 1.0, |x| 1.0e5 + (20.0 * x).sin());
+        assert!(varying.error > 0.0, "expected a nonzero error estimate");
+        // Both error estimates reflect only the varying part; the 1e5 offset
+        // (deflation factor ~316 under the old scaling) must not shrink it.
+        assert!(
+            offset.error > 0.5 * varying.error && offset.error < 2.0 * varying.error,
+            "varying.error={}, offset.error={}",
+            varying.error,
+            offset.error
+        );
+    }
+
+    /// QUADPACK invariant (dqk15.f `dmin1(1.0, …)`): the panel error estimate
+    /// never exceeds resasc = ∫|f − mean| ≤ (b − a)·(max f − min f). On a panel
+    /// that G7-K15 resolves poorly, the (200·δ/resasc)^1.5 factor saturates its
+    /// cap; that cap must be 1·resasc, not 1000·resasc.
+    #[test]
+    fn error_estimate_bounded_by_variation_on_rough_panel() {
+        let gk = GaussKronrod::new(GKPair::G7K15);
+        // sin(50x) over [0,1] is badly under-resolved by 15 nodes, so the cap
+        // engages. |sin| ≤ 1 ⇒ resasc ≤ 1·2 = 2, so a correct estimate ≤ 2;
+        // with a 1000× cap it would be in the hundreds.
+        let d = gk.integrate_detail(0.0, 1.0, |x| (50.0 * x).sin());
+        assert!(
+            d.error <= 2.0,
+            "error estimate {} exceeds the panel's total variation",
+            d.error
+        );
+    }
+
     /// Degenerate interval gives zero.
     #[test]
     fn degenerate_interval() {
diff --git a/src/gauss_laguerre.rs b/src/gauss_laguerre.rs
index 409620e..9599231 100644
--- a/src/gauss_laguerre.rs
+++ b/src/gauss_laguerre.rs
@@ -45,8 +45,8 @@ impl GaussLaguerre {
     /// # Errors
     ///
     /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
-    /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1` or `alpha`
-    /// is non-finite.
+    /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1`, `alpha` is
+    /// non-finite, or `μ₀ = Γ(alpha+1)` overflows `f64` (`alpha` ≳ 170).
     pub fn new(n: usize, alpha: f64) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -89,7 +89,15 @@ fn compute_laguerre(n: usize, alpha: f64) -> Result<(Vec<f64>, Vec<f64>), Quadra
             k * (k + alpha)
         })
         .collect();
-    let mu0 = ln_gamma(alpha + 1.0).exp();
+    // μ₀ = Γ(α+1), formed in log-space. For α ≳ 170 this overflows f64; reject
+    // it rather than returning all-infinite weights with no diagnostic.
+    let ln_mu0 = ln_gamma(alpha + 1.0);
+    if !ln_mu0.is_finite() || ln_mu0 > f64::MAX.ln() {
+        return Err(QuadratureError::InvalidInput(
+            "Laguerre weight integral mu0 = Gamma(alpha+1) overflows f64 (alpha too large)",
+        ));
+    }
+    let mu0 = ln_mu0.exp();
 
     golub_welsch(&diag, &off_diag_sq, mu0)
 }
@@ -150,6 +158,15 @@ mod tests {
         assert!(GaussLaguerre::new(5, f64::NAN).is_err());
     }
 
+    #[test]
+    fn large_alpha_overflow_rejected() {
+        // μ₀ = Γ(alpha+1) overflows f64 for alpha ≳ 170; reject rather than
+        // returning all-infinite weights with no diagnostic.
+        assert!(GaussLaguerre::new(5, 200.0).is_err());
+        // A value comfortably under the overflow threshold still works.
+        assert!(GaussLaguerre::new(5, 150.0).is_ok());
+    }
+
     /// Polynomial exactness: integral of x^k e^(-x) over [0,inf) = k! = Gamma(k+1).
     #[test]
     fn polynomial_exactness() {
diff --git a/src/oscillatory.rs b/src/oscillatory.rs
index 73dc563..0f24156 100644
--- a/src/oscillatory.rs
+++ b/src/oscillatory.rs
@@ -97,7 +97,10 @@ impl OscillatoryIntegrator {
     /// # Errors
     ///
     /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
-    /// is non-finite. May also propagate errors from the adaptive fallback when
+    /// is non-finite. Returns [`QuadratureError::InvalidInput`] if `omega` is so
+    /// large that the Filon-Clenshaw-Curtis moment quadrature would need an
+    /// impractical number of nodes (`|omega·(b − a)/2|` above ~100,000). May
+    /// also propagate errors from the adaptive fallback when
     /// `|omega * (b - a) / 2|` is small.
     #[allow(clippy::many_single_char_names)] // a, b, f, n are conventional in quadrature
     #[allow(clippy::similar_names)] // sum_c_half / sum_s_half are intentionally parallel names
@@ -127,6 +130,14 @@ impl OscillatoryIntegrator {
         let mid = 0.5 * (a + b);
         let theta = self.omega * half;
 
+        // `b - a` can overflow to ±inf even when a and b are individually finite
+        // (an interval wider than f64::MAX); with omega == 0 that yields
+        // theta = 0 * inf = NaN, which would slip past both the small-theta
+        // fallback and the moment-node guard below into a silent NaN result.
+        if !theta.is_finite() {
+            return Err(QuadratureError::DegenerateInterval);
+        }
+
         // For small theta, fall back to adaptive integration
         if theta.abs() < 2.0 {
             return self.adaptive_fallback(a, b, &f);
@@ -135,6 +146,21 @@ impl OscillatoryIntegrator {
         // Filon-Clenshaw-Curtis method
         let n = self.order.max(4);
 
+        // The Chebyshev moments ∫T_j(x)·{cos,sin}(θx) dx are evaluated with a
+        // Gauss-Legendre rule whose node count must exceed the kernel bandwidth
+        // ~|θ| to resolve the oscillation. Bound that node count: an extreme ω
+        // would otherwise force a multi-gigabyte allocation, and the f64→usize
+        // cast of θ inside `modified_chebyshev_moments` could saturate and then
+        // overflow. The check is done in f64 before any cast, so an infinite or
+        // huge θ is rejected cleanly rather than wrapping to a tiny node count.
+        const MAX_MOMENT_NODES: usize = 100_000;
+        let m_required = n as f64 + theta.abs().ceil() + 32.0;
+        if m_required > MAX_MOMENT_NODES as f64 {
+            return Err(QuadratureError::InvalidInput(
+                "omega too large for Filon-Clenshaw-Curtis moments",
+            ));
+        }
+
         // Step 1: Evaluate f at Clenshaw-Curtis nodes on [-1, 1]
         let f_vals: Vec<f64> = (0..=n)
             .map(|k| {
@@ -336,7 +362,9 @@ fn modified_chebyshev_moments(theta: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
 /// # Errors
 ///
 /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
-/// is non-finite.
+/// is non-finite, or [`QuadratureError::InvalidInput`] if `omega` is too large
+/// for the Filon-Clenshaw-Curtis moment quadrature (see
+/// [`OscillatoryIntegrator::integrate`]).
 pub fn integrate_oscillatory_sin<G>(
     f: G,
     a: f64,
@@ -368,7 +396,9 @@ where
 /// # Errors
 ///
 /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
-/// is non-finite.
+/// is non-finite, or [`QuadratureError::InvalidInput`] if `omega` is too large
+/// for the Filon-Clenshaw-Curtis moment quadrature (see
+/// [`OscillatoryIntegrator::integrate`]).
 pub fn integrate_oscillatory_cos<G>(
     f: G,
     a: f64,
@@ -496,4 +526,17 @@ mod tests {
             result.value
         );
     }
+
+    #[test]
+    fn extreme_omega_errors_cleanly() {
+        // A frequency far beyond the Filon moment-node budget must return
+        // InvalidInput rather than attempting a multi-gigabyte allocation or
+        // overflowing the f64→usize node-count cast (which previously wrapped
+        // to a tiny node count and returned a silently-wrong result).
+        let r = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 1e9, 1e-8);
+        assert!(
+            matches!(r, Err(QuadratureError::InvalidInput(_))),
+            "expected InvalidInput, got {r:?}"
+        );
+    }
 }
diff --git a/src/tanh_sinh.rs b/src/tanh_sinh.rs
index 464fc63..dc392cc 100644
--- a/src/tanh_sinh.rs
+++ b/src/tanh_sinh.rs
@@ -107,6 +107,10 @@ impl TanhSinh {
 
         let mid = 0.5 * (a + b);
         let half = 0.5 * (b - a);
+        // Open-interval bounds for the singular-endpoint guard. Ordered so the
+        // membership test holds for reversed limits (a > b); the sign of the
+        // integral is carried by `half`, which is negative when a > b.
+        let (lo, hi) = (a.min(b), a.max(b));
         let pi_2 = core::f64::consts::FRAC_PI_2;
 
         let mut total_evals = 0usize;
@@ -179,7 +183,7 @@ impl TanhSinh {
                 // Positive t: x = mid + half * tanh(u)
                 let mut local_contrib = 0.0_f64;
                 let x_pos = mid + half * tanh_u;
-                if x_pos > a && x_pos < b {
+                if x_pos > lo && x_pos < hi {
                     let fval = f(x_pos);
                     if fval.is_finite() {
                         new_sum += weight * fval;
@@ -190,7 +194,7 @@ impl TanhSinh {
 
                 // Negative t: x = mid - half * tanh(u)
                 let x_neg = mid - half * tanh_u;
-                if x_neg > a && x_neg < b {
+                if x_neg > lo && x_neg < hi {
                     let fval = f(x_neg);
                     if fval.is_finite() {
                         new_sum += weight * fval;
@@ -201,7 +205,7 @@ impl TanhSinh {
 
                 // Check if actual contributions (including f(x)) are negligible
                 let threshold = f64::EPSILON * prev_estimate.abs().max(1e-300);
-                if local_contrib * half < threshold {
+                if local_contrib * half.abs() < threshold {
                     consecutive_tiny += 1;
                     if consecutive_tiny >= 3 {
                         break;
@@ -308,6 +312,30 @@ mod tests {
         assert!((result.value - 2.0).abs() < 1e-7, "value={}", result.value);
     }
 
+    #[test]
+    fn reversed_bounds_negate() {
+        // tanh-sinh must honour reversed limits like the rest of the library:
+        // ∫_b^a f = -∫_a^b f. Here the endpoint singularity at 0 becomes the
+        // upper bound. Before the fix, the membership test dropped every
+        // interior sample for a > b and returned a meaningless near-zero value.
+        let rev = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 1.0, 0.0, 1e-10).unwrap();
+        assert!(
+            rev.converged,
+            "did not converge: err={}",
+            rev.error_estimate
+        );
+        assert!((rev.value - (-2.0)).abs() < 1e-7, "value={}", rev.value);
+
+        // Equals the negation of the forward integral.
+        let fwd = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap();
+        assert!(
+            (rev.value + fwd.value).abs() < 1e-9,
+            "fwd={}, rev={}",
+            fwd.value,
+            rev.value
+        );
+    }
+
     #[test]
     fn log_singularity() {
         // ∫₀¹ ln(x) dx = -1
diff --git a/src/transforms.rs b/src/transforms.rs
index c3f5f9b..0de6ea9 100644
--- a/src/transforms.rs
+++ b/src/transforms.rs
@@ -69,7 +69,16 @@ where
         let one_minus_t = 1.0 - t;
         let x = a + t / one_minus_t;
         let jacobian = 1.0 / (one_minus_t * one_minus_t);
-        f(x) * jacobian
+        // Jacobian overflow just inside the endpoint can give 0*inf = NaN or
+        // finite*inf = inf; treat a non-finite contribution as a zero tail (the
+        // QUADPACK qagi guard) so one bad panel cannot poison the estimate. Note
+        // this also masks a genuine non-finite f at an interior mapped point.
+        let v = f(x) * jacobian;
+        if v.is_finite() {
+            v
+        } else {
+            0.0
+        }
     };
 
     AdaptiveIntegrator::default()
@@ -121,7 +130,16 @@ where
         let one_minus_t = 1.0 - t;
         let x = b - t / one_minus_t;
         let jacobian = 1.0 / (one_minus_t * one_minus_t);
-        f(x) * jacobian
+        // Jacobian overflow just inside the endpoint can give 0*inf = NaN or
+        // finite*inf = inf; treat a non-finite contribution as a zero tail (the
+        // QUADPACK qagi guard) so one bad panel cannot poison the estimate. Note
+        // this also masks a genuine non-finite f at an interior mapped point.
+        let v = f(x) * jacobian;
+        if v.is_finite() {
+            v
+        } else {
+            0.0
+        }
     };
 
     AdaptiveIntegrator::default()
@@ -163,7 +181,16 @@ where
         let one_minus_t2 = 1.0 - t2;
         let x = t / one_minus_t2;
         let jacobian = (1.0 + t2) / (one_minus_t2 * one_minus_t2);
-        f(x) * jacobian
+        // Jacobian overflow just inside the endpoint can give 0*inf = NaN or
+        // finite*inf = inf; treat a non-finite contribution as a zero tail (the
+        // QUADPACK qagi guard) so one bad panel cannot poison the estimate. Note
+        // this also masks a genuine non-finite f at an interior mapped point.
+        let v = f(x) * jacobian;
+        if v.is_finite() {
+            v
+        } else {
+            0.0
+        }
     };
 
     AdaptiveIntegrator::default()
@@ -179,6 +206,17 @@ mod tests {
     use super::*;
     use core::f64::consts::PI;
 
+    /// A non-finite integrand value at a reachable mapped node must not poison
+    /// the whole estimate. f is +inf beyond |x| = 50, which the outer nodes of
+    /// the very first panel already reach (|x| up to ~113); the finiteness guard
+    /// masks that contribution to 0 so the result stays finite instead of NaN.
+    #[test]
+    fn non_finite_integrand_contribution_does_not_poison() {
+        let f = |x: f64| if x.abs() > 50.0 { f64::INFINITY } else { 1.0 };
+        let r = integrate_infinite(f, 1e-8).unwrap();
+        assert!(r.value.is_finite(), "value={}", r.value);
+    }
+
     /// exp(-x) over [0, inf) = 1.
     #[test]
     fn exp_decay_upper() {