From 0ebb1eea13be934597d2c0d4245bf64bf1a3a1d0 Mon Sep 17 00:00:00 2001
From: Greg von Nessi <greg.vonnessi@entrolution.ai>
Date: Thu, 26 Feb 2026 18:30:33 +0000
Subject: [PATCH] Codebase cleanup: CI, docs, pedantic lints, dedup, tests

- CI: add feature matrix (default/all/no_std), remove orphaned codecov.yml
- Docs: populate CHANGELOG, add examples to 8 convenience fns + crate docs
- Lints: resolve all 613 pedantic clippy warnings (0 remaining)
- Dedup: extract impl_rule_accessors! macro (9 types), QmcSequence trait,
  shared CC weight computation
- Tests: 15 new tests for weighted domains, high-dim Sobol, reversed bounds
---
 .github/workflows/ci.yml    |  59 +----
 CHANGELOG.md                |  38 +++-
 CONTRIBUTING.md             |   2 +-
 codecov.yml                 |  18 --
 src/adaptive.rs             |  90 +++++++-
 src/cauchy_pv.rs            |  27 +++
 src/clenshaw_curtis.rs      |  35 +--
 src/cubature/adaptive.rs    |  38 +++-
 src/cubature/halton.rs      |  20 +-
 src/cubature/mod.rs         |  17 +-
 src/cubature/monte_carlo.rs | 229 ++++++++++----------
 src/cubature/sobol.rs       |  95 +++++++-
 src/cubature/sparse_grid.rs |  57 ++---
 src/cubature/tensor.rs      |  11 +
 src/gauss_chebyshev.rs      |  62 ++----
 src/gauss_hermite.rs        |  33 +--
 src/gauss_jacobi.rs         |  60 +++---
 src/gauss_kronrod.rs        | 418 ++++++++++++++++++------------------
 src/gauss_laguerre.rs       |  36 +---
 src/gauss_legendre.rs       | 257 +++++++++++-----------
 src/gauss_lobatto.rs        |  33 +--
 src/gauss_radau.rs          |  35 +--
 src/golub_welsch.rs         |   7 +-
 src/lib.rs                  |  98 +++++++++
 src/oscillatory.rs          |  58 ++++-
 src/rule.rs                 |  18 ++
 src/tanh_sinh.rs            |  25 +++
 src/weighted.rs             | 127 ++++++++++-
 28 files changed, 1196 insertions(+), 807 deletions(-)
 delete mode 100644 codecov.yml

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d105245..bc7c8c8 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -14,6 +14,9 @@ jobs:
   test:
     name: Test
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        features: ["", "--all-features", "--no-default-features"]
     steps:
       - uses: actions/checkout@v6
 
@@ -24,10 +27,10 @@ jobs:
         uses: Swatinem/rust-cache@v2
 
       - name: Build
-        run: cargo build
+        run: cargo build ${{ matrix.features }}
 
       - name: Run tests
-        run: cargo test
+        run: cargo test ${{ matrix.features }}
 
   lint:
     name: Lint
@@ -47,7 +50,7 @@ jobs:
         run: cargo fmt --all -- --check
 
       - name: Clippy
-        run: cargo clippy -- -D warnings
+        run: cargo clippy --all-features -- -D warnings
 
   docs:
     name: Documentation
@@ -83,53 +86,3 @@ jobs:
 
       - name: Test (default features)
         run: cargo test
-
-  coverage:
-    name: Coverage
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v6
-
-      - name: Install Rust
-        uses: dtolnay/rust-toolchain@stable
-
-      - name: Cache cargo
-        uses: Swatinem/rust-cache@v2
-
-      - name: Install cargo-tarpaulin
-        run: cargo install cargo-tarpaulin
-
-      - name: Generate coverage
-        run: cargo tarpaulin --out xml --timeout 300
-        continue-on-error: true
-
-      - name: Upload coverage
-        uses: codecov/codecov-action@v5
-        with:
-          files: cobertura.xml
-          token: ${{ secrets.CODECOV_TOKEN }}
-          fail_ci_if_error: false
-
-  security:
-    name: Security Audit
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v6
-
-      - name: Install Rust
-        uses: dtolnay/rust-toolchain@stable
-
-      - name: Cache cargo
-        uses: Swatinem/rust-cache@v2
-
-      - name: Install cargo-audit
-        run: cargo install cargo-audit --locked
-
-      - name: Audit dependencies
-        run: cargo audit
-
-      - name: Install cargo-deny
-        run: cargo install cargo-deny --locked
-
-      - name: Check dependency policies
-        run: cargo deny check
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1c15487..84d81bb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,40 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
+## [0.1.0] - Unreleased
 
-[Unreleased]: https://github.com/Entrolution/bilby/compare/v0.1.0...HEAD
+Initial release of bilby.
+
+### Added
+
+- **Gauss-Legendre** quadrature with Bogaert (2014) asymptotic expansion for large n and Newton iteration for small n
+- **Gauss-Kronrod** embedded pairs (G7K15, G10K21, G15K31, G20K41, G25K51, G30K61) for error estimation
+- **Classical families** via Golub-Welsch eigenvalue solver:
+  - Gauss-Jacobi (arbitrary alpha, beta)
+  - Gauss-Hermite (physicists' convention)
+  - Gauss-Laguerre (generalised, arbitrary alpha)
+  - Gauss-Chebyshev Types I and II (closed-form)
+  - Gauss-Lobatto (both endpoints)
+  - Gauss-Radau (left or right endpoint)
+  - Clenshaw-Curtis (nested, Chebyshev extrema)
+- **Adaptive integration** -- QUADPACK-style global subdivision with configurable GK pair, break points, absolute/relative tolerance
+- **Domain transforms** -- `integrate_infinite`, `integrate_semi_infinite_lower`, `integrate_semi_infinite_upper` for unbounded domains
+- **Tanh-sinh** (double-exponential) quadrature for endpoint singularities
+- **Cauchy principal value** integration via the subtraction technique
+- **Oscillatory integration** via Filon-Clenshaw-Curtis with adaptive fallback for small omega
+- **Weighted integration** -- unified API for product integration with Jacobi, Laguerre, Hermite, Chebyshev, and log weight functions
+- **Multi-dimensional cubature**:
+  - Tensor product rules (isotropic and anisotropic)
+  - Smolyak sparse grids with Clenshaw-Curtis basis
+  - Genz-Malik adaptive cubature (degree-7/degree-5 embedded rule)
+  - Monte Carlo (plain pseudo-random with Welford variance estimation)
+  - Quasi-Monte Carlo via Sobol and Halton low-discrepancy sequences
+- **Generic over `F: Float`** via num-traits for core rule types
+- **`no_std` compatible** -- works with `alloc` only (disable default `std` feature)
+- **`parallel` feature** -- rayon-based parallel variants:
+  - `GaussLegendre::new_par` for parallel node generation
+  - `integrate_composite_par` for parallel composite rules
+  - `integrate_par` / `integrate_box_par` for cubature rules
+  - `MonteCarloIntegrator::integrate_par` for parallel MC/QMC
+- **Precomputed rule cache** (`cache` module, requires `std`) with `GL5`, `GL10`, `GL15`, `GL20`, `GL50`, `GL100` lazy singletons
+- **Benchmarks** via criterion for node generation, 1D integration, and cubature
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 8dac9a3..bb65511 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -56,7 +56,7 @@ Before submitting a PR, ensure:
 cargo fmt
 
 # Run clippy
-cargo clippy -- -D warnings
+cargo clippy --all-features -- -D warnings
 
 # Check documentation
 cargo doc --no-deps
diff --git a/codecov.yml b/codecov.yml
deleted file mode 100644
index 79baf69..0000000
--- a/codecov.yml
+++ /dev/null
@@ -1,18 +0,0 @@
-coverage:
-  precision: 2
-  round: down
-  range: "70...100"
-
-status:
-  project:
-    default:
-      target: 85%
-      threshold: 2%
-  patch:
-    default:
-      target: 80%
-
-comment:
-  layout: "reach,diff,flags,files"
-  behavior: default
-  require_changes: true
diff --git a/src/adaptive.rs b/src/adaptive.rs
index 124258d..f0f7515 100644
--- a/src/adaptive.rs
+++ b/src/adaptive.rs
@@ -97,24 +97,28 @@ impl Default for AdaptiveIntegrator {
 
 impl AdaptiveIntegrator {
     /// Set the Gauss-Kronrod pair used for panel evaluation.
+    #[must_use]
     pub fn with_pair(mut self, pair: GKPair) -> Self {
         self.pair = pair;
         self
     }
 
     /// Set the absolute error tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set the relative error tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Set the maximum number of function evaluations.
+    #[must_use]
     pub fn with_max_evals(mut self, n: usize) -> Self {
         self.max_evals = n;
         self
@@ -125,6 +129,11 @@ impl AdaptiveIntegrator {
     /// Returns a [`QuadratureResult`] with the best estimate and error bound.
     /// If the tolerance could not be achieved within the evaluation budget,
     /// `converged` will be `false` but the best estimate is still returned.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if both tolerances are non-positive
+    /// or if `max_evals` is less than one Gauss-Kronrod panel evaluation.
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -142,6 +151,11 @@ impl AdaptiveIntegrator {
     /// Break points split the interval at known discontinuities or singularities.
     /// The adaptive algorithm then refines each sub-interval independently in a
     /// global priority queue.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if any break point is outside
+    /// `(a, b)`, is NaN, or if both tolerances are non-positive.
     pub fn integrate_with_breaks<G>(
         &self,
         a: f64,
@@ -189,6 +203,9 @@ impl AdaptiveIntegrator {
 
         // Seed the priority queue with initial intervals
         for &(ia, ib) in intervals {
+            // Exact comparison is intentional: degenerate zero-width
+            // intervals are an edge case, not a floating-point tolerance issue.
+            #[allow(clippy::float_cmp)]
             if ia == ib {
                 continue;
             }
@@ -210,10 +227,7 @@ impl AdaptiveIntegrator {
         while !self.tolerance_met(total_estimate, total_error)
             && num_evals + 2 * evals_per_call <= self.max_evals
         {
-            let worst = match heap.pop() {
-                Some(s) => s,
-                None => break,
-            };
+            let Some(worst) = heap.pop() else { break };
 
             // Bisect the worst interval
             let mid = 0.5 * (worst.a + worst.b);
@@ -340,6 +354,16 @@ where
 /// Break points indicate locations of discontinuities or singularities.
 /// The interval is split at these points before adaptive refinement begins.
 ///
+/// # Example
+///
+/// ```
+/// use bilby::adaptive_integrate_with_breaks;
+///
+/// // Integrate |x| over [-1, 1] with a break at the cusp
+/// let result = adaptive_integrate_with_breaks(|x: f64| x.abs(), -1.0, 1.0, &[0.0], 1e-12).unwrap();
+/// assert!((result.value - 1.0).abs() < 1e-12);
+/// ```
+///
 /// # Errors
 ///
 /// Returns [`QuadratureError::InvalidInput`] if any break point is outside `(a, b)` or NaN.
@@ -507,4 +531,62 @@ mod tests {
         assert!(r.num_evals > 0);
         assert!(r.num_evals <= 200, "num_evals={}", r.num_evals);
     }
+
+    /// Reversed bounds for exp(x): ∫₁⁰ exp(x) dx = -(e - 1).
+    #[test]
+    fn reversed_bounds_exp() {
+        let forward = adaptive_integrate(f64::exp, 0.0, 1.0, 1e-12).unwrap();
+        let reverse = adaptive_integrate(f64::exp, 1.0, 0.0, 1e-12).unwrap();
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
+
+    /// Reversed bounds for a polynomial: ∫₁⁰ x² dx = -1/3.
+    #[test]
+    fn reversed_bounds_polynomial() {
+        let r = adaptive_integrate(|x| x * x, 1.0, 0.0, 1e-12).unwrap();
+        assert!(r.converged, "err={}", r.error_estimate);
+        assert!(
+            (r.value + 1.0 / 3.0).abs() < 1e-12,
+            "value={}, expected -1/3",
+            r.value
+        );
+    }
+
+    /// Reversed bounds with break points: ∫₁⁻¹ |x| dx = -1.
+    #[test]
+    fn reversed_bounds_with_breaks() {
+        let forward =
+            adaptive_integrate_with_breaks(|x: f64| x.abs(), -1.0, 1.0, &[0.0], 1e-12).unwrap();
+        let reverse =
+            adaptive_integrate_with_breaks(|x: f64| x.abs(), 1.0, -1.0, &[0.0], 1e-12).unwrap();
+        assert!(forward.converged, "forward err={}", forward.error_estimate);
+        assert!(reverse.converged, "reverse err={}", reverse.error_estimate);
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
+
+    /// Reversed bounds with break points for step function.
+    #[test]
+    fn reversed_bounds_step_with_breaks() {
+        let step = |x: f64| if x < 2.0 { 1.0 } else { 3.0 };
+        let forward = adaptive_integrate_with_breaks(step, 0.0, 4.0, &[2.0], 1e-12).unwrap();
+        let reverse = adaptive_integrate_with_breaks(step, 4.0, 0.0, &[2.0], 1e-12).unwrap();
+        assert!(forward.converged, "forward err={}", forward.error_estimate);
+        assert!(reverse.converged, "reverse err={}", reverse.error_estimate);
+        assert!(
+            (forward.value + reverse.value).abs() < 1e-12,
+            "forward={}, reverse={}",
+            forward.value,
+            reverse.value
+        );
+    }
 }
diff --git a/src/cauchy_pv.rs b/src/cauchy_pv.rs
index d6a6e80..4c3f712 100644
--- a/src/cauchy_pv.rs
+++ b/src/cauchy_pv.rs
@@ -48,18 +48,21 @@ impl Default for CauchyPV {
 
 impl CauchyPV {
     /// Set absolute tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set relative tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Set maximum number of function evaluations.
+    #[must_use]
     pub fn with_max_evals(mut self, n: usize) -> Self {
         self.max_evals = n;
         self
@@ -68,6 +71,13 @@ impl CauchyPV {
     /// Compute PV ∫ₐᵇ f(x)/(x-c) dx.
     ///
     /// The singularity c must be strictly inside (a, b).
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN.
+    /// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside
+    /// `(a, b)` or if `f(c)` is not finite.
+    #[allow(clippy::many_single_char_names)] // a, b, c, f, g are conventional in quadrature
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -128,6 +138,23 @@ impl CauchyPV {
 }
 
 /// Convenience: Cauchy principal value integration with default settings.
+///
+/// # Example
+///
+/// ```
+/// use bilby::pv_integrate;
+///
+/// // PV integral of x^2/(x - 0.3) over [0, 1]
+/// let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln();
+/// let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap();
+/// assert!((result.value - exact).abs() < 1e-7);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::DegenerateInterval`] if any bound or `c` is NaN.
+/// Returns [`QuadratureError::InvalidInput`] if `c` is not strictly inside
+/// `(a, b)` or if `f(c)` is not finite.
 pub fn pv_integrate<G>(
     f: G,
     a: f64,
diff --git a/src/clenshaw_curtis.rs b/src/clenshaw_curtis.rs
index f6a5971..e8a2053 100644
--- a/src/clenshaw_curtis.rs
+++ b/src/clenshaw_curtis.rs
@@ -38,6 +38,10 @@ impl ClenshawCurtis {
     /// Create a new n-point Clenshaw-Curtis rule.
     ///
     /// Requires `n >= 1`. For `n == 1`, returns the midpoint rule.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -48,37 +52,16 @@ impl ClenshawCurtis {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(ClenshawCurtis, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 /// Compute n-point Clenshaw-Curtis nodes and weights.
 ///
 /// For n == 1: midpoint rule (x=0, w=2).
-/// For n >= 2: nodes at x_k = cos(k π / (n-1)), weights via explicit formula.
-fn compute_clenshaw_curtis(n: usize) -> (Vec<f64>, Vec<f64>) {
+/// For n >= 2: nodes at `x_k` = cos(k π / (n-1)), weights via explicit formula.
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
+pub(crate) fn compute_clenshaw_curtis(n: usize) -> (Vec<f64>, Vec<f64>) {
     use core::f64::consts::PI;
 
     if n == 1 {
diff --git a/src/cubature/adaptive.rs b/src/cubature/adaptive.rs
index ddd4549..32c9809 100644
--- a/src/cubature/adaptive.rs
+++ b/src/cubature/adaptive.rs
@@ -53,24 +53,33 @@ impl Default for AdaptiveCubature {
 
 impl AdaptiveCubature {
     /// Set absolute tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set relative tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Set maximum number of function evaluations.
+    #[must_use]
     pub fn with_max_evals(mut self, n: usize) -> Self {
         self.max_evals = n;
         self
     }
 
     /// Adaptively integrate `f` over the hyperrectangle \[lower, upper\].
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
+    /// different lengths or are empty. Returns [`QuadratureError::DegenerateInterval`]
+    /// if any bound is NaN.
     pub fn integrate<G>(
         &self,
         lower: &[f64],
@@ -123,10 +132,7 @@ impl AdaptiveCubature {
         while global_error > self.abs_tol.max(self.rel_tol * global_estimate.abs())
             && total_evals < self.max_evals
         {
-            let worst = match heap.pop() {
-                Some(r) => r,
-                None => break,
-            };
+            let Some(worst) = heap.pop() else { break };
 
             global_estimate -= worst.estimate;
             global_error -= worst.error;
@@ -176,6 +182,22 @@ impl AdaptiveCubature {
 }
 
 /// Convenience: adaptive cubature with default settings.
+///
+/// # Example
+///
+/// ```
+/// use bilby::cubature::adaptive_cubature;
+///
+/// // Integral of x*y over [0,1]^2 = 1/4
+/// let result = adaptive_cubature(|x| x[0] * x[1], &[0.0, 0.0], &[1.0, 1.0], 1e-10).unwrap();
+/// assert!((result.value - 0.25).abs() < 1e-10);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
+/// different lengths or are empty. Returns [`QuadratureError::DegenerateInterval`]
+/// if any bound is NaN.
 pub fn adaptive_cubature<G>(
     f: G,
     lower: &[f64],
@@ -254,6 +276,8 @@ where
     let lambda4 = (9.0_f64 / 10.0).sqrt();
     let lambda5 = (9.0_f64 / 19.0).sqrt();
 
+    // Dimension is bounded by practical limits (d <= ~7), so this cast is safe.
+    #[allow(clippy::cast_precision_loss)]
     let d_f = d as f64;
 
     // Weights for degree-7 rule (scaled for [-1,1]^d)
@@ -261,7 +285,8 @@ where
     let w2 = 980.0 / 6561.0;
     let w3 = (1820.0 - 400.0 * d_f) / 19683.0;
     let w4 = 200.0 / 19683.0;
-    let w5 = 6859.0 / 19683.0 / (1 << d) as f64;
+    // d is bounded by practical limits (d <= ~7), so (1 << d) fits in i32.
+    let w5 = 6859.0 / 19683.0 / f64::from(1i32 << d);
 
     // Weights for degree-5 rule
     let wp1 = (729.0 - 950.0 * d_f + 50.0 * d_f * d_f) / 729.0;
@@ -355,8 +380,7 @@ where
         .iter()
         .enumerate()
         .max_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(Ordering::Equal))
-        .map(|(i, _)| i)
-        .unwrap_or(0);
+        .map_or(0, |(i, _)| i);
 
     GMDetail {
         estimate: est7,
diff --git a/src/cubature/halton.rs b/src/cubature/halton.rs
index 38b20d4..f88a2aa 100644
--- a/src/cubature/halton.rs
+++ b/src/cubature/halton.rs
@@ -43,6 +43,10 @@ impl HaltonSequence {
     /// Create a new Halton sequence generator for `dim` dimensions.
     ///
     /// Supports up to 100 dimensions.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero or exceeds 100.
     pub fn new(dim: usize) -> Result<Self, QuadratureError> {
         if dim == 0 {
             return Err(QuadratureError::InvalidInput("dimension must be >= 1"));
@@ -60,6 +64,10 @@ impl HaltonSequence {
     }
 
     /// Generate the next point in \[0, 1)^d.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `point.len()` is less than the sequence dimension.
     pub fn next_point(&mut self, point: &mut [f64]) {
         assert!(point.len() >= self.dim);
         self.index += 1;
@@ -69,11 +77,13 @@ impl HaltonSequence {
     }
 
     /// Current index (number of points generated so far).
+    #[must_use]
     pub fn index(&self) -> u64 {
         self.index
     }
 
     /// Spatial dimension.
+    #[must_use]
     pub fn dim(&self) -> usize {
         self.dim
     }
@@ -81,12 +91,16 @@ impl HaltonSequence {
 
 /// Radical-inverse function: reverse the digits of `n` in the given `base`.
 fn radical_inverse(mut n: u64, base: u32) -> f64 {
-    let base_f = base as f64;
+    let base_f = f64::from(base);
+    let base_u64 = u64::from(base);
     let mut result = 0.0;
     let mut factor = 1.0 / base_f;
     while n > 0 {
-        result += (n % base as u64) as f64 * factor;
-        n /= base as u64;
+        // n % base_u64 is always < base (a small u32), so fits exactly in f64.
+        #[allow(clippy::cast_precision_loss)]
+        let digit = (n % base_u64) as f64;
+        result += digit * factor;
+        n /= base_u64;
         factor /= base_f;
     }
     result
diff --git a/src/cubature/mod.rs b/src/cubature/mod.rs
index 2d583c2..d08b7f6 100644
--- a/src/cubature/mod.rs
+++ b/src/cubature/mod.rs
@@ -52,7 +52,7 @@ use alloc::{vec, vec::Vec};
 pub struct CubatureRule {
     /// Nodes stored flat: node i at indices [i*dim .. (i+1)*dim].
     nodes: Vec<f64>,
-    /// Weights (length = num_points).
+    /// Weights (length = `num_points`).
     weights: Vec<f64>,
     /// Spatial dimension.
     dim: usize,
@@ -62,6 +62,7 @@ impl CubatureRule {
     /// Create a new cubature rule from flat node data and weights.
     ///
     /// `nodes_flat` must have length `weights.len() * dim`.
+    #[must_use]
     pub fn new(nodes_flat: Vec<f64>, weights: Vec<f64>, dim: usize) -> Self {
         debug_assert_eq!(nodes_flat.len(), weights.len() * dim);
         Self {
@@ -73,24 +74,28 @@ impl CubatureRule {
 
     /// Number of cubature points.
     #[inline]
+    #[must_use]
     pub fn num_points(&self) -> usize {
         self.weights.len()
     }
 
     /// Spatial dimension.
     #[inline]
+    #[must_use]
     pub fn dim(&self) -> usize {
         self.dim
     }
 
     /// Access the i-th node as a slice of length `dim`.
     #[inline]
+    #[must_use]
     pub fn node(&self, i: usize) -> &[f64] {
         &self.nodes[i * self.dim..(i + 1) * self.dim]
     }
 
     /// The weights.
     #[inline]
+    #[must_use]
     pub fn weights(&self) -> &[f64] {
         &self.weights
     }
@@ -111,6 +116,11 @@ impl CubatureRule {
     /// Integrate `f` over the hyperrectangle \[lower, upper\].
     ///
     /// Applies an affine transform from the reference domain \[-1, 1\]^d.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `lower.len()` or `upper.len()` does not equal the rule's
+    /// spatial dimension.
     pub fn integrate_box<G>(&self, lower: &[f64], upper: &[f64], f: G) -> f64
     where
         G: Fn(&[f64]) -> f64,
@@ -155,6 +165,11 @@ impl CubatureRule {
     /// Parallel integration over the hyperrectangle \[lower, upper\].
     ///
     /// Identical to [`integrate_box`](Self::integrate_box) but evaluates points in parallel.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `lower.len()` or `upper.len()` does not equal the rule's
+    /// spatial dimension.
     #[cfg(feature = "parallel")]
     pub fn integrate_box_par<G>(&self, lower: &[f64], upper: &[f64], f: G) -> f64
     where
diff --git a/src/cubature/monte_carlo.rs b/src/cubature/monte_carlo.rs
index adca729..09ef58e 100644
--- a/src/cubature/monte_carlo.rs
+++ b/src/cubature/monte_carlo.rs
@@ -13,6 +13,34 @@ use alloc::{vec, vec::Vec};
 #[cfg(not(feature = "std"))]
 use num_traits::Float as _;
 
+/// Trait abstracting a low-discrepancy sequence for QMC integration.
+///
+/// Implemented by [`SobolSequence`] and [`HaltonSequence`].
+trait QmcSequence {
+    fn new(dim: usize) -> Result<Self, QuadratureError>
+    where
+        Self: Sized;
+    fn next_point(&mut self, point: &mut [f64]);
+}
+
+impl QmcSequence for SobolSequence {
+    fn new(dim: usize) -> Result<Self, QuadratureError> {
+        SobolSequence::new(dim)
+    }
+    fn next_point(&mut self, point: &mut [f64]) {
+        self.next_point(point);
+    }
+}
+
+impl QmcSequence for HaltonSequence {
+    fn new(dim: usize) -> Result<Self, QuadratureError> {
+        HaltonSequence::new(dim)
+    }
+    fn next_point(&mut self, point: &mut [f64]) {
+        self.next_point(point);
+    }
+}
+
 /// Monte Carlo integration method.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum MCMethod {
@@ -54,24 +82,34 @@ impl Default for MonteCarloIntegrator {
 
 impl MonteCarloIntegrator {
     /// Set the MC method.
+    #[must_use]
     pub fn with_method(mut self, method: MCMethod) -> Self {
         self.method = method;
         self
     }
 
     /// Set the number of samples.
+    #[must_use]
     pub fn with_samples(mut self, n: usize) -> Self {
         self.num_samples = n;
         self
     }
 
     /// Set the random seed (only used for `MCMethod::Plain`).
+    #[must_use]
     pub fn with_seed(mut self, seed: u64) -> Self {
         self.seed = seed;
         self
     }
 
     /// Integrate `f` over the hyperrectangle \[lower, upper\].
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
+    /// different lengths, are empty, or if the number of samples is zero.
+    /// 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>(
         &self,
         lower: &[f64],
@@ -99,11 +137,19 @@ impl MonteCarloIntegrator {
 
         match self.method {
             MCMethod::Plain => self.integrate_plain(d, lower, &widths, volume, n, &f),
-            MCMethod::Sobol => self.integrate_qmc_sobol(d, lower, &widths, volume, n, &f),
-            MCMethod::Halton => self.integrate_qmc_halton(d, lower, &widths, volume, n, &f),
+            MCMethod::Sobol => {
+                self.integrate_qmc::<SobolSequence, _>(d, lower, &widths, volume, n, &f)
+            }
+            MCMethod::Halton => {
+                self.integrate_qmc::<HaltonSequence, _>(d, lower, &widths, volume, n, &f)
+            }
         }
     }
 
+    // Returns Result for API consistency with the other integrate_* methods,
+    // which may fail (e.g., Sobol dimension limit).
+    #[allow(clippy::unnecessary_wraps)]
+    #[allow(clippy::cast_precision_loss)] // usize-to-f64 casts are for sample counts, always small enough
     fn integrate_plain<G>(
         &self,
         d: usize,
@@ -145,7 +191,10 @@ impl MonteCarloIntegrator {
         })
     }
 
-    fn integrate_qmc_sobol<G>(
+    #[allow(clippy::unused_self)] // &self for API consistency with integrate_plain
+    #[allow(clippy::many_single_char_names)] // d, n, f, x, u are conventional in MC integration
+    #[allow(clippy::cast_precision_loss)] // usize-to-f64 casts are for sample counts
+    fn integrate_qmc<S, G>(
         &self,
         d: usize,
         lower: &[f64],
@@ -155,15 +204,16 @@ impl MonteCarloIntegrator {
         f: &G,
     ) -> Result<QuadratureResult<f64>, QuadratureError>
     where
+        S: QmcSequence,
         G: Fn(&[f64]) -> f64,
     {
-        let mut sob = SobolSequence::new(d)?;
+        let mut seq = S::new(d)?;
         let mut x = vec![0.0; d];
         let mut u = vec![0.0; d];
         let mut sum = 0.0;
 
         for _ in 0..n {
-            sob.next_point(&mut u);
+            seq.next_point(&mut u);
             for j in 0..d {
                 x[j] = lower[j] + widths[j] * u[j];
             }
@@ -174,63 +224,11 @@ impl MonteCarloIntegrator {
 
         // Heuristic error estimate: compare N/2 and N estimates
         let half_sum = {
-            let mut s2 = SobolSequence::new(d)?;
-            let half = n / 2;
-            let mut sm = 0.0;
-            for _ in 0..half {
-                s2.next_point(&mut u);
-                for j in 0..d {
-                    x[j] = lower[j] + widths[j] * u[j];
-                }
-                sm += f(&x);
-            }
-            volume * sm / half as f64
-        };
-
-        let error = (estimate - half_sum).abs();
-
-        Ok(QuadratureResult {
-            value: estimate,
-            error_estimate: error,
-            num_evals: n + n / 2,
-            converged: true,
-        })
-    }
-
-    fn integrate_qmc_halton<G>(
-        &self,
-        d: usize,
-        lower: &[f64],
-        widths: &[f64],
-        volume: f64,
-        n: usize,
-        f: &G,
-    ) -> Result<QuadratureResult<f64>, QuadratureError>
-    where
-        G: Fn(&[f64]) -> f64,
-    {
-        let mut hal = HaltonSequence::new(d)?;
-        let mut x = vec![0.0; d];
-        let mut u = vec![0.0; d];
-        let mut sum = 0.0;
-
-        for _ in 0..n {
-            hal.next_point(&mut u);
-            for j in 0..d {
-                x[j] = lower[j] + widths[j] * u[j];
-            }
-            sum += f(&x);
-        }
-
-        let estimate = volume * sum / n as f64;
-
-        // Heuristic error: compare N/2 and N
-        let half_sum = {
-            let mut h2 = HaltonSequence::new(d)?;
+            let mut seq2 = S::new(d)?;
             let half = n / 2;
             let mut sm = 0.0;
             for _ in 0..half {
-                h2.next_point(&mut u);
+                seq2.next_point(&mut u);
                 for j in 0..d {
                     x[j] = lower[j] + widths[j] * u[j];
                 }
@@ -256,6 +254,13 @@ impl MonteCarloIntegrator {
     /// evaluations are parallelised. For Plain MC, the sample space is split
     /// into independent chunks with separate PRNGs, giving statistically
     /// equivalent but numerically different results compared to sequential.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
+    /// different lengths, are empty, or if the number of samples is zero.
+    /// 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")]
     pub fn integrate_par<G>(
         &self,
@@ -284,12 +289,19 @@ impl MonteCarloIntegrator {
 
         match self.method {
             MCMethod::Plain => self.integrate_plain_par(d, lower, &widths, volume, n, &f),
-            MCMethod::Sobol => self.integrate_qmc_par_sobol(d, lower, &widths, volume, n, &f),
-            MCMethod::Halton => self.integrate_qmc_par_halton(d, lower, &widths, volume, n, &f),
+            MCMethod::Sobol => {
+                self.integrate_qmc_par::<SobolSequence, _>(d, lower, &widths, volume, n, &f)
+            }
+            MCMethod::Halton => {
+                self.integrate_qmc_par::<HaltonSequence, _>(d, lower, &widths, volume, n, &f)
+            }
         }
     }
 
+    // Returns Result for API consistency with the other integrate_*_par methods.
     #[cfg(feature = "parallel")]
+    #[allow(clippy::unnecessary_wraps)]
+    #[allow(clippy::cast_precision_loss)] // usize-to-f64 casts are for sample counts
     fn integrate_plain_par<G>(
         &self,
         d: usize,
@@ -312,7 +324,7 @@ impl MonteCarloIntegrator {
         let chunk_results: Vec<(f64, f64, usize)> = (0..num_chunks)
             .into_par_iter()
             .map(|chunk_id| {
-                let my_n = chunk_size + if chunk_id < remainder { 1 } else { 0 };
+                let my_n = chunk_size + usize::from(chunk_id < remainder);
                 if my_n == 0 {
                     return (0.0, 0.0, 0);
                 }
@@ -322,6 +334,7 @@ impl MonteCarloIntegrator {
                 let mut mean = 0.0;
                 let mut m2 = 0.0;
 
+                #[allow(clippy::cast_precision_loss)]
                 for i in 1..=my_n {
                     for j in 0..d {
                         x[j] = lower[j] + widths[j] * rng.next_f64();
@@ -333,7 +346,9 @@ impl MonteCarloIntegrator {
                     m2 += delta * delta2;
                 }
 
-                (mean * my_n as f64, m2, my_n)
+                #[allow(clippy::cast_precision_loss)]
+                let sum = mean * my_n as f64;
+                (sum, m2, my_n)
             })
             .collect();
 
@@ -347,12 +362,15 @@ impl MonteCarloIntegrator {
             total_n += count;
         }
 
+        #[allow(clippy::cast_precision_loss)]
         let mean = total_sum / total_n as f64;
+        #[allow(clippy::cast_precision_loss)]
         let variance = if total_n > 1 {
             total_m2 / (total_n - 1) as f64
         } else {
             0.0
         };
+        #[allow(clippy::cast_precision_loss)]
         let std_error = (variance / total_n as f64).sqrt() * volume;
 
         Ok(QuadratureResult {
@@ -364,7 +382,9 @@ impl MonteCarloIntegrator {
     }
 
     #[cfg(feature = "parallel")]
-    fn integrate_qmc_par_sobol<G>(
+    #[allow(clippy::unused_self)] // &self for API consistency with integrate_plain_par
+    #[allow(clippy::cast_precision_loss)] // usize-to-f64 casts are for sample counts
+    fn integrate_qmc_par<S, G>(
         &self,
         d: usize,
         lower: &[f64],
@@ -374,16 +394,17 @@ impl MonteCarloIntegrator {
         f: &G,
     ) -> Result<QuadratureResult<f64>, QuadratureError>
     where
+        S: QmcSequence,
         G: Fn(&[f64]) -> f64 + Sync,
     {
         use rayon::prelude::*;
 
-        // Pre-generate all Sobol points (sequential — gray-code is inherently serial)
-        let mut sob = SobolSequence::new(d)?;
+        // Pre-generate all points (sequential — sequences are inherently serial)
+        let mut seq = S::new(d)?;
         let mut points = vec![0.0; n * d];
         let mut u = vec![0.0; d];
         for i in 0..n {
-            sob.next_point(&mut u);
+            seq.next_point(&mut u);
             for j in 0..d {
                 points[i * d + j] = lower[j] + widths[j] * u[j];
             }
@@ -413,60 +434,25 @@ impl MonteCarloIntegrator {
             converged: true,
         })
     }
-
-    #[cfg(feature = "parallel")]
-    fn integrate_qmc_par_halton<G>(
-        &self,
-        d: usize,
-        lower: &[f64],
-        widths: &[f64],
-        volume: f64,
-        n: usize,
-        f: &G,
-    ) -> Result<QuadratureResult<f64>, QuadratureError>
-    where
-        G: Fn(&[f64]) -> f64 + Sync,
-    {
-        use rayon::prelude::*;
-
-        // Pre-generate all Halton points
-        let mut hal = HaltonSequence::new(d)?;
-        let mut points = vec![0.0; n * d];
-        let mut u = vec![0.0; d];
-        for i in 0..n {
-            hal.next_point(&mut u);
-            for j in 0..d {
-                points[i * d + j] = lower[j] + widths[j] * u[j];
-            }
-        }
-
-        // Parallel evaluation
-        let sum: f64 = (0..n)
-            .into_par_iter()
-            .map(|i| f(&points[i * d..(i + 1) * d]))
-            .sum();
-
-        let estimate = volume * sum / n as f64;
-
-        // Heuristic error
-        let half = n / 2;
-        let half_sum: f64 = (0..half)
-            .into_par_iter()
-            .map(|i| f(&points[i * d..(i + 1) * d]))
-            .sum();
-        let half_estimate = volume * half_sum / half as f64;
-        let error = (estimate - half_estimate).abs();
-
-        Ok(QuadratureResult {
-            value: estimate,
-            error_estimate: error,
-            num_evals: n + half,
-            converged: true,
-        })
-    }
 }
 
 /// Convenience: quasi-Monte Carlo integration using Sobol sequences.
+///
+/// # Example
+///
+/// ```
+/// use bilby::cubature::monte_carlo_integrate;
+///
+/// // Integral of 1 over [0,1]^3 = 1
+/// let result = monte_carlo_integrate(|_| 1.0, &[0.0; 3], &[1.0; 3], 10000).unwrap();
+/// assert!((result.value - 1.0).abs() < 0.01);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::InvalidInput`] if `lower` and `upper` have
+/// different lengths, are empty, if `num_samples` is zero, or if the
+/// dimension exceeds the Sobol sequence limit of 40.
 pub fn monte_carlo_integrate<G>(
     f: G,
     lower: &[f64],
@@ -494,9 +480,9 @@ impl Xoshiro256pp {
         let mut s = [0u64; 4];
         let mut z = seed;
         for slot in &mut s {
-            z = z.wrapping_add(0x9e3779b97f4a7c15);
-            z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
-            z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
+            z = z.wrapping_add(0x9e37_79b9_7f4a_7c15);
+            z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
+            z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
             *slot = z ^ (z >> 31);
         }
         Self { s }
@@ -516,6 +502,7 @@ impl Xoshiro256pp {
         result
     }
 
+    #[allow(clippy::cast_precision_loss)] // Shifting right by 11 keeps 53 bits, which fit exactly in f64
     fn next_f64(&mut self) -> f64 {
         // Generate a double in [0, 1) using the upper 53 bits
         (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
diff --git a/src/cubature/sobol.rs b/src/cubature/sobol.rs
index 3a2c551..a2b2d7e 100644
--- a/src/cubature/sobol.rs
+++ b/src/cubature/sobol.rs
@@ -44,6 +44,10 @@ impl SobolSequence {
     /// Create a new Sobol sequence generator for `dim` dimensions.
     ///
     /// Supports up to 40 dimensions (expandable with more direction numbers).
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero or exceeds 40.
     pub fn new(dim: usize) -> Result<Self, QuadratureError> {
         if dim == 0 {
             return Err(QuadratureError::InvalidInput("dimension must be >= 1"));
@@ -59,7 +63,10 @@ impl SobolSequence {
 
         // Dimension 0: Van der Corput (base 2)
         for (i, d) in direction.iter_mut().enumerate().take(b) {
-            *d = 1u32 << (BITS - 1 - i as u32);
+            // i < BITS (32), so i as u32 cannot truncate.
+            #[allow(clippy::cast_possible_truncation)]
+            let i_u32 = i as u32;
+            *d = 1u32 << (BITS - 1 - i_u32);
         }
 
         // Higher dimensions: from direction number table
@@ -70,7 +77,10 @@ impl SobolSequence {
 
             // Initial direction numbers from the table
             for i in 0..s as usize {
-                direction[j * b + i] = entry.m[i] << (BITS - 1 - i as u32);
+                // i < s <= 8, so i as u32 cannot truncate.
+                #[allow(clippy::cast_possible_truncation)]
+                let i_u32 = i as u32;
+                direction[j * b + i] = entry.m[i] << (BITS - 1 - i_u32);
             }
 
             // Generate remaining direction numbers via recurrence
@@ -95,6 +105,10 @@ impl SobolSequence {
     }
 
     /// Generate the next point in \[0, 1)^d.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `point.len()` is less than the sequence dimension.
     pub fn next_point(&mut self, point: &mut [f64]) {
         assert!(point.len() >= self.dim);
         self.index += 1;
@@ -102,11 +116,13 @@ impl SobolSequence {
         // Find the rightmost zero bit of (index - 1) (gray-code position)
         let c = (self.index - 1).trailing_ones() as usize;
         let b = BITS as usize;
+        // 1u64 << 32 = 4294967296, fits exactly in f64.
+        #[allow(clippy::cast_precision_loss)]
         let norm = 1.0 / (1u64 << BITS) as f64;
 
         for (j, p) in point.iter_mut().enumerate().take(self.dim) {
             self.state[j] ^= self.direction[j * b + c];
-            *p = self.state[j] as f64 * norm;
+            *p = f64::from(self.state[j]) * norm;
         }
     }
 
@@ -133,11 +149,13 @@ impl SobolSequence {
     }
 
     /// Current index.
+    #[must_use]
     pub fn index(&self) -> u64 {
         self.index
     }
 
     /// Spatial dimension.
+    #[must_use]
     pub fn dim(&self) -> usize {
         self.dim
     }
@@ -152,14 +170,14 @@ struct SobolEntry {
     degree: u32,
     /// Coefficients of the primitive polynomial (excluding leading and trailing 1).
     coeffs: u32,
-    /// Initial direction numbers m_1, ..., m_s.
+    /// Initial direction numbers `m_1`, ..., `m_s`.
     m: [u32; 8],
 }
 
 /// Direction numbers from Joe & Kuo (2010) for dimensions 2..40.
 /// Each entry: (degree s, polynomial coefficients a, initial m-values).
 ///
-/// Source: https://web.maths.unsw.edu.au/~fkuo/sobol/joe-kuo-old.1111
+/// Source: <https://web.maths.unsw.edu.au/~fkuo/sobol/joe-kuo-old.1111>
 static SOBOL_TABLE: [SobolEntry; 39] = [
     SobolEntry {
         degree: 1,
@@ -412,4 +430,71 @@ mod tests {
         assert!(SobolSequence::new(0).is_err());
         assert!(SobolSequence::new(41).is_err());
     }
+
+    #[test]
+    fn high_dim_points_in_unit_cube() {
+        let dim = 20;
+        let mut sob = SobolSequence::new(dim).unwrap();
+        let mut pt = vec![0.0; dim];
+        for i in 0..200 {
+            sob.next_point(&mut pt);
+            for (j, &x) in pt.iter().enumerate() {
+                assert!(
+                    x >= 0.0 && x < 1.0,
+                    "point {i}, dim {j}: x={x} out of [0,1)"
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn high_dim_30_points_in_unit_cube() {
+        let dim = 30;
+        let mut sob = SobolSequence::new(dim).unwrap();
+        let mut pt = vec![0.0; dim];
+        for i in 0..100 {
+            sob.next_point(&mut pt);
+            for (j, &x) in pt.iter().enumerate() {
+                assert!(
+                    x >= 0.0 && x < 1.0,
+                    "point {i}, dim {j}: x={x} out of [0,1)"
+                );
+            }
+        }
+    }
+
+    #[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;
+        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
+        }
+        let estimate = sum / n as f64;
+        assert!(
+            (estimate - 1.0).abs() < 1e-14,
+            "estimate={estimate}, expected 1.0"
+        );
+    }
+
+    #[test]
+    fn max_dim_construction() {
+        // Verify we can construct the maximum supported dimension (40).
+        let sob = SobolSequence::new(40);
+        assert!(sob.is_ok());
+        let mut sob = sob.unwrap();
+        assert_eq!(sob.dim(), 40);
+
+        let mut pt = vec![0.0; 40];
+        sob.next_point(&mut pt);
+        for (j, &x) in pt.iter().enumerate() {
+            assert!(x >= 0.0 && x < 1.0, "dim {j}: x={x} out of [0,1)");
+        }
+    }
 }
diff --git a/src/cubature/sparse_grid.rs b/src/cubature/sparse_grid.rs
index a168df2..801cc50 100644
--- a/src/cubature/sparse_grid.rs
+++ b/src/cubature/sparse_grid.rs
@@ -52,6 +52,10 @@ impl SparseGrid {
     /// Level 0 uses a single center point. Higher levels add more points
     /// for greater accuracy. The rule is exact for polynomials of increasing
     /// total degree as the level increases.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero.
     pub fn new(dim: usize, level: usize, _basis: SparseGridBasis) -> Result<Self, QuadratureError> {
         if dim == 0 {
             return Err(QuadratureError::InvalidInput("dimension must be >= 1"));
@@ -62,30 +66,38 @@ impl SparseGrid {
     }
 
     /// Construct using Clenshaw-Curtis basis (convenience).
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero.
     pub fn clenshaw_curtis(dim: usize, level: usize) -> Result<Self, QuadratureError> {
         Self::new(dim, level, SparseGridBasis::ClenshawCurtis)
     }
 
     /// Returns a reference to the underlying cubature rule.
     #[inline]
+    #[must_use]
     pub fn rule(&self) -> &CubatureRule {
         &self.rule
     }
 
     /// Number of cubature points.
     #[inline]
+    #[must_use]
     pub fn num_points(&self) -> usize {
         self.rule.num_points()
     }
 
     /// Spatial dimension.
     #[inline]
+    #[must_use]
     pub fn dim(&self) -> usize {
         self.rule.dim()
     }
 
     /// Smolyak level.
     #[inline]
+    #[must_use]
     pub fn level(&self) -> usize {
         self.level
     }
@@ -105,42 +117,9 @@ fn cc_order(level: usize) -> usize {
 
 /// Compute Clenshaw-Curtis nodes and weights for a given number of points.
 ///
-/// Returns nodes on [-1, 1] sorted ascending.
+/// Delegates to the shared implementation in [`crate::clenshaw_curtis`].
 fn cc_rule(n: usize) -> (Vec<f64>, Vec<f64>) {
-    if n == 1 {
-        return (vec![0.0], vec![2.0]);
-    }
-
-    let nm1 = n - 1;
-    let nm1_f = nm1 as f64;
-    let pi = core::f64::consts::PI;
-
-    // Nodes: cos(k*pi/(n-1)) for k=0..n-1, reversed to ascending order
-    let nodes: Vec<f64> = (0..n)
-        .map(|k| -(((nm1 - k) as f64) * pi / nm1_f).cos())
-        .collect();
-
-    // Weights via explicit formula
-    let m = nm1 / 2;
-    let mut weights = vec![0.0; n];
-    for (i, w) in weights.iter_mut().enumerate() {
-        let k = nm1 - i;
-        let mut sum = 0.0;
-        for j in 0..=m {
-            let j_f = j as f64;
-            let b_j = if j == 0 || (nm1.is_multiple_of(2) && j == m) {
-                1.0
-            } else {
-                2.0
-            };
-            let denom = 1.0 - 4.0 * j_f * j_f;
-            sum += b_j / denom * (2.0 * j_f * k as f64 * pi / nm1_f).cos();
-        }
-        let c_k = if k == 0 || k == nm1 { 1.0 } else { 2.0 };
-        *w = c_k * sum / nm1_f;
-    }
-
-    (nodes, weights)
+    crate::clenshaw_curtis::compute_clenshaw_curtis(n)
 }
 
 /// Quantise a float to an integer key for exact point merging.
@@ -148,6 +127,8 @@ fn cc_rule(n: usize) -> (Vec<f64>, Vec<f64>) {
 /// CC nodes are cos(k*pi/n) which are algebraic numbers. Using 48-bit
 /// quantisation avoids floating-point comparison issues when merging
 /// duplicate points from different tensor products.
+// x is in [-1, 1], so x * 2^48 fits in i64. The casts are intentional.
+#[allow(clippy::cast_precision_loss, clippy::cast_possible_truncation)]
 fn quantise(x: f64) -> i64 {
     (x * (1i64 << 48) as f64).round() as i64
 }
@@ -155,9 +136,9 @@ fn quantise(x: f64) -> i64 {
 /// Build a Smolyak sparse grid using the combination technique.
 ///
 /// The Smolyak formula:
-///   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})
+///   `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.
+/// Using 0-indexed levels where `l_j` ≥ 0.
 fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
     // Precompute CC rules for each level we'll need
     let max_level = level;
@@ -188,6 +169,8 @@ fn build_smolyak(dim: usize, level: usize) -> CubatureRule {
 
     for s in sum_min..=sum_max {
         let diff = q - s;
+        // Binomial coefficients for sparse grid levels are small enough to fit in f64.
+        #[allow(clippy::cast_precision_loss)]
         let coeff =
             if diff.is_multiple_of(2) { 1.0 } else { -1.0 } * binomial(dim - 1, diff) as f64;
 
diff --git a/src/cubature/tensor.rs b/src/cubature/tensor.rs
index ceceb39..013a1ea 100644
--- a/src/cubature/tensor.rs
+++ b/src/cubature/tensor.rs
@@ -38,6 +38,10 @@ impl TensorProductRule {
     /// Construct from a slice of 1D rules (one per dimension).
     ///
     /// Total points = product of all 1D orders.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `rules_1d` is empty.
     pub fn new(rules_1d: &[&QuadratureRule<f64>]) -> Result<Self, QuadratureError> {
         if rules_1d.is_empty() {
             return Err(QuadratureError::InvalidInput(
@@ -77,6 +81,10 @@ impl TensorProductRule {
     }
 
     /// Construct from a single 1D rule applied to all dimensions.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `dim` is zero.
     pub fn isotropic(rule_1d: &QuadratureRule<f64>, dim: usize) -> Result<Self, QuadratureError> {
         let refs: Vec<&QuadratureRule<f64>> = (0..dim).map(|_| rule_1d).collect();
         Self::new(&refs)
@@ -84,18 +92,21 @@ impl TensorProductRule {
 
     /// Returns a reference to the underlying cubature rule.
     #[inline]
+    #[must_use]
     pub fn rule(&self) -> &CubatureRule {
         &self.rule
     }
 
     /// Number of cubature points.
     #[inline]
+    #[must_use]
     pub fn num_points(&self) -> usize {
         self.rule.num_points()
     }
 
     /// Spatial dimension.
     #[inline]
+    #[must_use]
     pub fn dim(&self) -> usize {
         self.rule.dim()
     }
diff --git a/src/gauss_chebyshev.rs b/src/gauss_chebyshev.rs
index 3f5b2ad..e5bef3b 100644
--- a/src/gauss_chebyshev.rs
+++ b/src/gauss_chebyshev.rs
@@ -45,6 +45,11 @@ pub struct GaussChebyshevFirstKind {
 
 impl GaussChebyshevFirstKind {
     /// Create a new n-point Gauss-Chebyshev Type I rule.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -68,32 +73,10 @@ impl GaussChebyshevFirstKind {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussChebyshevFirstKind, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 /// Gauss-Chebyshev Type II quadrature rule.
 ///
 /// Weight function `w(x) = sqrt(1 - x²)` on \[-1, 1\].
@@ -113,6 +96,11 @@ pub struct GaussChebyshevSecondKind {
 
 impl GaussChebyshevSecondKind {
     /// Create a new n-point Gauss-Chebyshev Type II rule.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    #[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -137,32 +125,10 @@ impl GaussChebyshevSecondKind {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussChebyshevSecondKind, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 #[cfg(test)]
 mod tests {
     use super::*;
diff --git a/src/gauss_hermite.rs b/src/gauss_hermite.rs
index 62eceff..8525c7b 100644
--- a/src/gauss_hermite.rs
+++ b/src/gauss_hermite.rs
@@ -35,6 +35,10 @@ pub struct GaussHermite {
 
 impl GaussHermite {
     /// Create a new n-point Gauss-Hermite rule.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -45,39 +49,18 @@ impl GaussHermite {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on (-∞, ∞).
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussHermite, nodes_doc: "Returns the nodes on (-∞, ∞).");
+
 /// Compute n Gauss-Hermite nodes and weights via Golub-Welsch.
 ///
 /// Monic physicists' Hermite recurrence:
-///   x h̃_k = h̃_{k+1} + 0·h̃_k + (k/2)·h̃_{k-1}
+///   x `h̃_k` = `h̃_{k+1}` + 0·`h̃_k` + (k/2)·`h̃_{k-1}`
 ///
 /// Jacobi matrix: diagonal = 0, off-diagonal² = (k+1)/2 for k=0..n-2.
 /// μ₀ = √π.
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_hermite(n: usize) -> (Vec<f64>, Vec<f64>) {
     let diag = vec![0.0; n];
     let off_diag_sq: Vec<f64> = (1..n).map(|k| k as f64 / 2.0).collect();
diff --git a/src/gauss_jacobi.rs b/src/gauss_jacobi.rs
index aa3d33c..1e53132 100644
--- a/src/gauss_jacobi.rs
+++ b/src/gauss_jacobi.rs
@@ -48,6 +48,12 @@ impl GaussJacobi {
     /// Create a new n-point Gauss-Jacobi rule with parameters α and β.
     ///
     /// Requires `n >= 1`, `α > -1`, `β > -1`.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1`, `beta <= -1`,
+    /// or either parameter is NaN.
     pub fn new(n: usize, alpha: f64, beta: f64) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -67,50 +73,31 @@ impl GaussJacobi {
     }
 
     /// Returns the α parameter.
+    #[must_use]
     pub fn alpha(&self) -> f64 {
         self.alpha
     }
 
     /// Returns the β parameter.
+    #[must_use]
     pub fn beta(&self) -> f64 {
         self.beta
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussJacobi, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 /// Compute n Gauss-Jacobi nodes and weights via the Golub-Welsch algorithm.
 ///
 /// The monic Jacobi polynomial recurrence:
-///   x p_k = p_{k+1} + α_k p_k + β_k p_{k-1}
+///   x `p_k` = `p_{k+1}` + `α_k` `p_k` + `β_k` `p_{k-1}`
 ///
 /// with:
-///   α_k = (β²-α²) / ((2k+α+β)(2k+α+β+2))
-///   β_k = 4k(k+α)(k+β)(k+α+β) / ((2k+α+β)²(2k+α+β+1)(2k+α+β-1))   for k ≥ 1
+///   `α_k` = (β²-α²) / ((2k+α+β)(2k+α+β+2))
+///   `β_k` = 4k(k+α)(k+β)(k+α+β) / ((2k+α+β)²(2k+α+β+1)(2k+α+β-1))   for k ≥ 1
 ///
 /// μ₀ = 2^(α+β+1) Γ(α+1)Γ(β+1) / Γ(α+β+2)
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> (Vec<f64>, Vec<f64>) {
     let ab = alpha + beta;
 
@@ -161,6 +148,7 @@ fn compute_jacobi(n: usize, alpha: f64, beta: f64) -> (Vec<f64>, Vec<f64>) {
 }
 
 #[allow(clippy::excessive_precision)]
+#[allow(clippy::cast_precision_loss)] // Loop index i is at most 8, fits exactly in f64
 /// Log-gamma function using the Lanczos approximation.
 ///
 /// Accurate to ~15 digits for positive arguments.
@@ -168,15 +156,15 @@ pub(crate) fn ln_gamma(x: f64) -> f64 {
     // Lanczos approximation with g=7, n=9
     // Coefficients from Numerical Recipes
     const COEFF: [f64; 9] = [
-        0.99999999999980993,
-        676.5203681218851,
-        -1259.1392167224028,
-        771.32342877765313,
-        -176.61502916214059,
-        12.507343278686905,
-        -0.13857109526572012,
-        9.9843695780195716e-6,
-        1.5056327351493116e-7,
+        0.999_999_999_999_809_93,
+        676.520_368_121_885_1,
+        -1_259.139_216_722_402_8,
+        771.323_428_777_653_13,
+        -176.615_029_162_140_59,
+        12.507_343_278_686_905,
+        -0.138_571_095_265_720_12,
+        9.984_369_578_019_571_6e-6,
+        1.505_632_735_149_311_6e-7,
     ];
 
     if x < 0.5 {
diff --git a/src/gauss_kronrod.rs b/src/gauss_kronrod.rs
index 181a0f6..8ab8018 100644
--- a/src/gauss_kronrod.rs
+++ b/src/gauss_kronrod.rs
@@ -46,6 +46,7 @@ pub enum GKPair {
 impl GKPair {
     /// Number of Kronrod nodes.
     #[inline]
+    #[must_use]
     pub fn kronrod_order(self) -> usize {
         match self {
             Self::G7K15 => 15,
@@ -58,6 +59,7 @@ impl GKPair {
 
     /// Number of Gauss nodes embedded in the Kronrod rule.
     #[inline]
+    #[must_use]
     pub fn gauss_order(self) -> usize {
         match self {
             Self::G7K15 => 7,
@@ -103,6 +105,7 @@ pub struct GaussKronrod {
 
 impl GaussKronrod {
     /// Create a new Gauss-Kronrod rule for the given pair.
+    #[must_use]
     pub fn new(pair: GKPair) -> Self {
         let (xgk, wgk, wg) = match pair {
             GKPair::G7K15 => (&GK15_XGK[..], &GK15_WGK[..], &GK15_WG[..]),
@@ -116,11 +119,12 @@ impl GaussKronrod {
 
     /// Returns the pair type.
     #[inline]
+    #[must_use]
     pub fn pair(&self) -> GKPair {
         self.pair
     }
 
-    /// Integrate `f` over \[a, b\], returning (estimate, error_estimate).
+    /// Integrate `f` over \[a, b\], returning (estimate, `error_estimate`).
     ///
     /// The estimate uses the full Kronrod rule. The error estimate is based
     /// on the difference between the Kronrod and Gauss results, using the
@@ -138,6 +142,8 @@ impl GaussKronrod {
     /// Returns estimate, error, integral of |f|, and integral of |f - mean|
     /// over the interval. The adaptive integrator uses these for roundoff
     /// detection.
+    #[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
     where
         G: Fn(f64) -> f64,
@@ -209,11 +215,13 @@ impl GaussKronrod {
     }
 
     /// Returns the Kronrod order.
+    #[must_use]
     pub fn kronrod_order(&self) -> usize {
         self.pair.kronrod_order()
     }
 
     /// Returns the embedded Gauss order.
+    #[must_use]
     pub fn gauss_order(&self) -> usize {
         self.pair.gauss_order()
     }
@@ -238,275 +246,275 @@ impl GaussKronrod {
 // G7-K15: 8 abscissae, 8 Kronrod weights, 4 Gauss weights
 #[rustfmt::skip]
 static GK15_XGK: [f64; 8] = [
-    0.9914553711208126392068546975263e0,
-    0.9491079123427585245261896840479e0,
-    0.8648644233597690727897127886093e0,
-    0.7415311855993944398638647732808e0,
-    0.5860872354676911302941448382587e0,
-    0.4058451513773971669066064120770e0,
-    0.2077849550078984676006894037733e0,
-    0.0000000000000000000000000000000e0,
+    0.991_455_371_120_812_639_206_854_697_526_3e0,
+    0.949_107_912_342_758_524_526_189_684_047_9e0,
+    0.864_864_423_359_769_072_789_712_788_609_3e0,
+    0.741_531_185_599_394_439_863_864_773_280_8e0,
+    0.586_087_235_467_691_130_294_144_838_258_7e0,
+    0.405_845_151_377_397_166_906_606_412_077_0e0,
+    0.207_784_955_007_898_467_600_689_403_773_3e0,
+    0.000_000_000_000_000_000_000_000_000_000_0e0,
 ];
 
 #[rustfmt::skip]
 static GK15_WGK: [f64; 8] = [
-    0.0229353220105292249637320080590e0,
-    0.0630920926299785532907006631892e0,
-    0.1047900103222501838398763225415e0,
-    0.1406532597155259187451895951024e0,
-    0.1690047266392679028265834265986e0,
-    0.1903505780647854099132564024211e0,
-    0.2044329400752988924141619992347e0,
-    0.2094821410847278280129991748917e0,
+    0.022_935_322_010_529_224_963_732_008_059_0e0,
+    0.063_092_092_629_978_553_290_700_663_189_2e0,
+    0.104_790_010_322_250_183_839_876_322_541_5e0,
+    0.140_653_259_715_525_918_745_189_595_102_4e0,
+    0.169_004_726_639_267_902_826_583_426_598_6e0,
+    0.190_350_578_064_785_409_913_256_402_421_1e0,
+    0.204_432_940_075_298_892_414_161_999_234_7e0,
+    0.209_482_141_084_727_828_012_999_174_891_7e0,
 ];
 
 // G7 weights: correspond to XGK indices 1, 3, 5, 7 (center)
 #[rustfmt::skip]
 static GK15_WG: [f64; 4] = [
-    0.1294849661688696932706114326791e0,
-    0.2797053914892766679014677714238e0,
-    0.3818300505051189449503697754890e0,
-    0.4179591836734693877551020408163e0,
+    0.129_484_966_168_869_693_270_611_432_679_1e0,
+    0.279_705_391_489_276_667_901_467_771_423_8e0,
+    0.381_830_050_505_118_944_950_369_775_489_0e0,
+    0.417_959_183_673_469_387_755_102_040_816_3e0,
 ];
 
 // G10-K21: 11 abscissae, 11 Kronrod weights, 5 Gauss weights
 #[rustfmt::skip]
 static GK21_XGK: [f64; 11] = [
-    0.9956571630258080807355272806890e0,
-    0.9739065285171717200779640120845e0,
-    0.9301574913557082260012071800595e0,
-    0.8650633666889845107320966884235e0,
-    0.7808177265864168970637175783450e0,
-    0.6794095682990244062343271146987e0,
-    0.5627571346686046833390009927269e0,
-    0.4333953941292471907992659431657e0,
-    0.2943928627014601981311260310387e0,
-    0.1488743389816312108848260011297e0,
-    0.0000000000000000000000000000000e0,
+    0.995_657_163_025_808_080_735_527_280_689_0e0,
+    0.973_906_528_517_171_720_077_964_012_084_5e0,
+    0.930_157_491_355_708_226_001_207_180_059_5e0,
+    0.865_063_366_688_984_510_732_096_688_423_5e0,
+    0.780_817_726_586_416_897_063_717_578_345_0e0,
+    0.679_409_568_299_024_406_234_327_114_698_7e0,
+    0.562_757_134_668_604_683_339_000_992_726_9e0,
+    0.433_395_394_129_247_190_799_265_943_165_7e0,
+    0.294_392_862_701_460_198_131_126_031_038_7e0,
+    0.148_874_338_981_631_210_884_826_001_129_7e0,
+    0.000_000_000_000_000_000_000_000_000_000_0e0,
 ];
 
 #[rustfmt::skip]
 static GK21_WGK: [f64; 11] = [
-    0.0116946388673718742780640396062e0,
-    0.0325581623079647274788189724594e0,
-    0.0547558965743519960313813002446e0,
-    0.0750396748109199527670431409619e0,
-    0.0931254545836976055350654650834e0,
-    0.1093871588022976418992105903258e0,
-    0.1234919762620658510779581098311e0,
-    0.1347092173114733259280540017717e0,
-    0.1427759385770600807970942731387e0,
-    0.1477391049013384913748415159721e0,
-    0.1494455540029169056649366438982e0,
+    0.011_694_638_867_371_874_278_064_039_606_2e0,
+    0.032_558_162_307_964_727_478_818_972_459_4e0,
+    0.054_755_896_574_351_996_031_381_300_244_6e0,
+    0.075_039_674_810_919_952_767_043_140_961_9e0,
+    0.093_125_454_583_697_605_535_065_465_083_4e0,
+    0.109_387_158_802_297_641_899_210_590_325_8e0,
+    0.123_491_976_262_065_851_077_958_109_831_1e0,
+    0.134_709_217_311_473_325_928_054_001_771_7e0,
+    0.142_775_938_577_060_080_797_094_273_138_7e0,
+    0.147_739_104_901_338_491_374_841_515_972_1e0,
+    0.149_445_554_002_916_905_664_936_643_898_2e0,
 ];
 
 // G10 weights: correspond to XGK indices 1, 3, 5, 7, 9
 #[rustfmt::skip]
 static GK21_WG: [f64; 5] = [
-    0.0666713443086881375935688098933e0,
-    0.1494513491505805931457763396577e0,
-    0.2190863625159820439955349342816e0,
-    0.2692667193099963550912269216947e0,
-    0.2955242247147528701738929465134e0,
+    0.066_671_344_308_688_137_593_568_809_893_3e0,
+    0.149_451_349_150_580_593_145_776_339_657_7e0,
+    0.219_086_362_515_982_043_995_534_934_281_6e0,
+    0.269_266_719_309_996_355_091_226_921_694_7e0,
+    0.295_524_224_714_752_870_173_892_946_513_4e0,
 ];
 
 // G15-K31: 16 abscissae, 16 Kronrod weights, 8 Gauss weights
 #[rustfmt::skip]
 static GK31_XGK: [f64; 16] = [
-    0.9980022986933970602851728401523e0,
-    0.9879925180204854284956518651868e0,
-    0.9677390756791391342573479878434e0,
-    0.9372733924007059043077834927821e0,
-    0.8972645323440819008825096564545e0,
-    0.8482065834104272162006482077822e0,
-    0.7904185014424659329766494817995e0,
-    0.7244177313601700474161606535169e0,
-    0.6509967412974169705337358953133e0,
-    0.5709721726085388475372263895391e0,
-    0.4850818636402396806936554402325e0,
-    0.3941513470775633689720720398105e0,
-    0.2991800071531688121668028224627e0,
-    0.2011940939974345223006283034596e0,
-    0.1011420669187174990270741231439e0,
-    0.0000000000000000000000000000000e0,
+    0.998_002_298_693_397_060_285_172_840_152_3e0,
+    0.987_992_518_020_485_428_495_651_865_186_8e0,
+    0.967_739_075_679_139_134_257_347_987_843_4e0,
+    0.937_273_392_400_705_904_307_783_492_782_1e0,
+    0.897_264_532_344_081_900_882_509_656_454_5e0,
+    0.848_206_583_410_427_216_200_648_207_782_2e0,
+    0.790_418_501_442_465_932_976_649_481_799_5e0,
+    0.724_417_731_360_170_047_416_160_653_516_9e0,
+    0.650_996_741_297_416_970_533_735_895_313_3e0,
+    0.570_972_172_608_538_847_537_226_389_539_1e0,
+    0.485_081_863_640_239_680_693_655_440_232_5e0,
+    0.394_151_347_077_563_368_972_072_039_810_5e0,
+    0.299_180_007_153_168_812_166_802_822_462_7e0,
+    0.201_194_093_997_434_522_300_628_303_459_6e0,
+    0.101_142_066_918_717_499_027_074_123_143_9e0,
+    0.000_000_000_000_000_000_000_000_000_000_0e0,
 ];
 
 #[rustfmt::skip]
 static GK31_WGK: [f64; 16] = [
-    0.0053774798729233489879207051413e0,
-    0.0150079473293161225383747630758e0,
-    0.0254608473267153201868740010197e0,
-    0.0353463607913758462220397848436e0,
-    0.0445897513247648766082272993728e0,
-    0.0534815246909280872653431472394e0,
-    0.0620095678006706402851392306080e0,
-    0.0698541213187282587095007799915e0,
-    0.0768496807577203788943277248266e0,
-    0.0830805028231330210382892472861e0,
-    0.0885644430562117706472754436938e0,
-    0.0931265981708253212252486872735e0,
-    0.0966427269836236785051790627259e0,
-    0.0991735987217919593323931739846e0,
-    0.1007698455238755950449266617570e0,
-    0.1013300070147915490173747926749e0,
+    0.005_377_479_872_923_348_987_920_705_141_3e0,
+    0.015_007_947_329_316_122_538_374_763_075_8e0,
+    0.025_460_847_326_715_320_186_874_001_019_7e0,
+    0.035_346_360_791_375_846_222_039_784_843_6e0,
+    0.044_589_751_324_764_876_608_227_299_372_8e0,
+    0.053_481_524_690_928_087_265_343_147_239_4e0,
+    0.062_009_567_800_670_640_285_139_230_608_0e0,
+    0.069_854_121_318_728_258_709_500_779_991_5e0,
+    0.076_849_680_757_720_378_894_327_724_826_6e0,
+    0.083_080_502_823_133_021_038_289_247_286_1e0,
+    0.088_564_443_056_211_770_647_275_443_693_8e0,
+    0.093_126_598_170_825_321_225_248_687_273_5e0,
+    0.096_642_726_983_623_678_505_179_062_725_9e0,
+    0.099_173_598_721_791_959_332_393_173_984_6e0,
+    0.100_769_845_523_875_595_044_926_661_757_0e0,
+    0.101_330_007_014_791_549_017_374_792_674_9e0,
 ];
 
 // G15 weights: correspond to XGK indices 1, 3, 5, 7, 9, 11, 13, 15 (center)
 #[rustfmt::skip]
 static GK31_WG: [f64; 8] = [
-    0.0307532419961172683546283935772e0,
-    0.0703660474881081247092674164507e0,
-    0.1071592204671719350118695468569e0,
-    0.1395706779261543144480479451103e0,
-    0.1662692058169939355320086048121e0,
-    0.1861610000155622110268005618642e0,
-    0.1984314853271115764561183264384e0,
-    0.2025782419255612728806019996752e0,
+    0.030_753_241_996_117_268_354_628_393_577_2e0,
+    0.070_366_047_488_108_124_709_267_416_450_7e0,
+    0.107_159_220_467_171_935_011_869_546_856_9e0,
+    0.139_570_677_926_154_314_448_047_945_110_3e0,
+    0.166_269_205_816_993_935_532_008_604_812_1e0,
+    0.186_161_000_015_562_211_026_800_561_864_2e0,
+    0.198_431_485_327_111_576_456_118_326_438_4e0,
+    0.202_578_241_925_561_272_880_601_999_675_2e0,
 ];
 
 // G20-K41: 21 abscissae, 21 Kronrod weights, 10 Gauss weights
 #[rustfmt::skip]
 static GK41_XGK: [f64; 21] = [
-    0.9988590315882776638383156530069e0,
-    0.9931285991850949247861223884713e0,
-    0.9815078774502502591933430720022e0,
-    0.9639719272779137912676661319728e0,
-    0.9408226338317547535199442212444e0,
-    0.9122344282513259058677524412033e0,
-    0.8782768112522819760774911307811e0,
-    0.8391169718222188823394529061702e0,
-    0.7950414288375511983506388272789e0,
-    0.7463319064601507926143050703556e0,
-    0.6932376563347513848054590571459e0,
-    0.6360536807265150254528369622629e0,
-    0.5751404468197103153429586642543e0,
-    0.5108670019508270980043640955253e0,
-    0.4435931752387251031999922134926e0,
-    0.3737060887154195606725481770493e0,
-    0.3016278681149130043205555356886e0,
-    0.2277858511416450780804961968576e0,
-    0.1526054652409226755052201216068e0,
-    0.0765265211334973337546409839884e0,
-    0.0000000000000000000000000000000e0,
+    0.998_859_031_588_277_663_838_315_653_006_9e0,
+    0.993_128_599_185_094_924_786_122_388_471_3e0,
+    0.981_507_877_450_250_259_193_343_072_002_2e0,
+    0.963_971_927_277_913_791_267_666_131_972_8e0,
+    0.940_822_633_831_754_753_519_944_221_244_4e0,
+    0.912_234_428_251_325_905_867_752_441_203_3e0,
+    0.878_276_811_252_281_976_077_491_130_781_1e0,
+    0.839_116_971_822_218_882_339_452_906_170_2e0,
+    0.795_041_428_837_551_198_350_638_827_278_9e0,
+    0.746_331_906_460_150_792_614_305_070_355_6e0,
+    0.693_237_656_334_751_384_805_459_057_145_9e0,
+    0.636_053_680_726_515_025_452_836_962_262_9e0,
+    0.575_140_446_819_710_315_342_958_664_254_3e0,
+    0.510_867_001_950_827_098_004_364_095_525_3e0,
+    0.443_593_175_238_725_103_199_992_213_492_6e0,
+    0.373_706_088_715_419_560_672_548_177_049_3e0,
+    0.301_627_868_114_913_004_320_555_535_688_6e0,
+    0.227_785_851_141_645_078_080_496_196_857_6e0,
+    0.152_605_465_240_922_675_505_220_121_606_8e0,
+    0.076_526_521_133_497_333_754_640_983_988_4e0,
+    0.000_000_000_000_000_000_000_000_000_000_0e0,
 ];
 
 #[rustfmt::skip]
 static GK41_WGK: [f64; 21] = [
-    0.0030735837185205315012182932460e0,
-    0.0086002698556429419866178795010e0,
-    0.0146261692569712529883787960309e0,
-    0.0203883734612665235980102314328e0,
-    0.0258821336049511583450506709615e0,
-    0.0312873067703279895895432233801e0,
-    0.0366001697582007980305572407072e0,
-    0.0416688733279736862637880593689e0,
-    0.0464348218674976747202318092611e0,
-    0.0509445739237286919327070705035e0,
-    0.0551951053482859947448329724198e0,
-    0.0591114008806395723749672206459e0,
-    0.0626532375547811680258701221255e0,
-    0.0658345971336184211156356956940e0,
-    0.0686486729285216193456234118537e0,
-    0.0710544235534440683057903617232e0,
-    0.0730306903327866674951894176891e0,
-    0.0745828754004991889865814183625e0,
-    0.0757044976845566746595427638462e0,
-    0.0763778676720807367050283538061e0,
-    0.0766007119179996564450499015301e0,
+    0.003_073_583_718_520_531_501_218_293_246_0e0,
+    0.008_600_269_855_642_941_986_617_879_501_0e0,
+    0.014_626_169_256_971_252_988_378_796_030_9e0,
+    0.020_388_373_461_266_523_598_010_231_432_8e0,
+    0.025_882_133_604_951_158_345_050_670_961_5e0,
+    0.031_287_306_770_327_989_589_543_223_380_1e0,
+    0.036_600_169_758_200_798_030_557_240_707_2e0,
+    0.041_668_873_327_973_686_263_788_059_368_9e0,
+    0.046_434_821_867_497_674_720_231_809_261_1e0,
+    0.050_944_573_923_728_691_932_707_070_503_5e0,
+    0.055_195_105_348_285_994_744_832_972_419_8e0,
+    0.059_111_400_880_639_572_374_967_220_645_9e0,
+    0.062_653_237_554_781_168_025_870_122_125_5e0,
+    0.065_834_597_133_618_421_115_635_695_694_0e0,
+    0.068_648_672_928_521_619_345_623_411_853_7e0,
+    0.071_054_423_553_444_068_305_790_361_723_2e0,
+    0.073_030_690_332_786_667_495_189_417_689_1e0,
+    0.074_582_875_400_499_188_986_581_418_362_5e0,
+    0.075_704_497_684_556_674_659_542_763_846_2e0,
+    0.076_377_867_672_080_736_705_028_353_806_1e0,
+    0.076_600_711_917_999_656_445_049_901_530_1e0,
 ];
 
 // G20 weights: correspond to XGK indices 1, 3, 5, 7, 9, 11, 13, 15, 17, 19
 #[rustfmt::skip]
 static GK41_WG: [f64; 10] = [
-    0.0176140071391521183118619623519e0,
-    0.0406014298003869413310399524974e0,
-    0.0626720483341090635695065351870e0,
-    0.0832767415767047487247581432205e0,
-    0.1019301198172404350367501354804e0,
-    0.1181945319615184173123773771138e0,
-    0.1316886384491766268989849749816e0,
-    0.1420961093183820513293006702817e0,
-    0.1491729864726037467878287370097e0,
-    0.1527533871307258506980433419510e0,
+    0.017_614_007_139_152_118_311_861_962_351_9e0,
+    0.040_601_429_800_386_941_331_039_952_497_4e0,
+    0.062_672_048_334_109_063_569_506_535_187_0e0,
+    0.083_276_741_576_704_748_724_758_143_220_5e0,
+    0.101_930_119_817_240_435_036_750_135_480_4e0,
+    0.118_194_531_961_518_417_312_377_377_113_8e0,
+    0.131_688_638_449_176_626_898_984_974_981_6e0,
+    0.142_096_109_318_382_051_329_300_670_281_7e0,
+    0.149_172_986_472_603_746_787_828_737_009_7e0,
+    0.152_753_387_130_725_850_698_043_341_951_0e0,
 ];
 
 // G25-K51: 26 abscissae, 26 Kronrod weights, 13 Gauss weights
 #[rustfmt::skip]
 static GK51_XGK: [f64; 26] = [
-    0.9992621049926098341934575486034e0,
-    0.9955569697904980990878487946894e0,
-    0.9880357945340772476385127157740e0,
-    0.9766639214595175114983153864598e0,
-    0.9616149864258425124181303366017e0,
-    0.9429745712289743394140116958647e0,
-    0.9207471152817015617463446613163e0,
-    0.8949919978782753688510420067828e0,
-    0.8658470652932755954499695958834e0,
-    0.8334426287608340014210069367057e0,
-    0.7978737979985000594104499430431e0,
-    0.7592592630373576305772826520436e0,
-    0.7177664068130843881866549077330e0,
-    0.6735663684734683644851203463762e0,
-    0.6268100990103174127881266182452e0,
-    0.5776629302412229677236898461265e0,
-    0.5263252843347191825996367815801e0,
-    0.4730027314457149605221821150919e0,
-    0.4178853821930377488518139459457e0,
-    0.3611723058093878377358182127264e0,
-    0.3030895389311078301674789098034e0,
-    0.2438668837209884320452190362797e0,
-    0.1837189394210488920159698889753e0,
-    0.1228646926107103963873598881804e0,
-    0.0615444830056850788865469263892e0,
-    0.0000000000000000000000000000000e0,
+    0.999_262_104_992_609_834_193_457_548_603_4e0,
+    0.995_556_969_790_498_099_087_848_794_689_4e0,
+    0.988_035_794_534_077_247_638_512_715_774_0e0,
+    0.976_663_921_459_517_511_498_315_386_459_8e0,
+    0.961_614_986_425_842_512_418_130_336_601_7e0,
+    0.942_974_571_228_974_339_414_011_695_864_7e0,
+    0.920_747_115_281_701_561_746_344_661_316_3e0,
+    0.894_991_997_878_275_368_851_042_006_782_8e0,
+    0.865_847_065_293_275_595_449_969_595_883_4e0,
+    0.833_442_628_760_834_001_421_006_936_705_7e0,
+    0.797_873_797_998_500_059_410_449_943_043_1e0,
+    0.759_259_263_037_357_630_577_282_652_043_6e0,
+    0.717_766_406_813_084_388_186_654_907_733_0e0,
+    0.673_566_368_473_468_364_485_120_346_376_2e0,
+    0.626_810_099_010_317_412_788_126_618_245_2e0,
+    0.577_662_930_241_222_967_723_689_846_126_5e0,
+    0.526_325_284_334_719_182_599_636_781_580_1e0,
+    0.473_002_731_445_714_960_522_182_115_091_9e0,
+    0.417_885_382_193_037_748_851_813_945_945_7e0,
+    0.361_172_305_809_387_837_735_818_212_726_4e0,
+    0.303_089_538_931_107_830_167_478_909_803_4e0,
+    0.243_866_883_720_988_432_045_219_036_279_7e0,
+    0.183_718_939_421_048_892_015_969_888_975_3e0,
+    0.122_864_692_610_710_396_387_359_888_180_4e0,
+    0.061_544_483_005_685_078_886_546_926_389_2e0,
+    0.000_000_000_000_000_000_000_000_000_000_0e0,
 ];
 
 #[rustfmt::skip]
 static GK51_WGK: [f64; 26] = [
-    0.0019873838923303159265078518288e0,
-    0.0055619321353567137580285023668e0,
-    0.0094739733861741516072077130124e0,
-    0.0132362291955716748137365405847e0,
-    0.0168478177091282982316675653634e0,
-    0.0204353711458828354568292235594e0,
-    0.0240099456069532162209289164881e0,
-    0.0274753175878517378094945551781e0,
-    0.0307923001673874888911090201523e0,
-    0.0340021302743293378367487529552e0,
-    0.0371162714834155435606030625368e0,
-    0.0400838255040323820748394846708e0,
-    0.0428728450201700494768957924395e0,
-    0.0455029130499217889098705847527e0,
-    0.0479825371388367139063925674920e0,
-    0.0502776790807156719633252593344e0,
-    0.0523628858064074758643661371279e0,
-    0.0542511298885454901454370598788e0,
-    0.0559508112204123173082406863275e0,
-    0.0574371163615678328535826939394e0,
-    0.0586896800223942079619745617888e0,
-    0.0597203403241740599790929193256e0,
-    0.0605394553760458629453602675157e0,
-    0.0611285097170530483058603013429e0,
-    0.0614711898714253166154431596526e0,
-    0.0615808180678329350759824240067e0,
+    0.001_987_383_892_330_315_926_507_851_828_8e0,
+    0.005_561_932_135_356_713_758_028_502_366_8e0,
+    0.009_473_973_386_174_151_607_207_713_012_4e0,
+    0.013_236_229_195_571_674_813_736_540_584_7e0,
+    0.016_847_817_709_128_298_231_667_565_363_4e0,
+    0.020_435_371_145_882_835_456_829_223_559_4e0,
+    0.024_009_945_606_953_216_220_928_916_488_1e0,
+    0.027_475_317_587_851_737_809_494_555_178_1e0,
+    0.030_792_300_167_387_488_891_109_020_152_3e0,
+    0.034_002_130_274_329_337_836_748_752_955_2e0,
+    0.037_116_271_483_415_543_560_603_062_536_8e0,
+    0.040_083_825_504_032_382_074_839_484_670_8e0,
+    0.042_872_845_020_170_049_476_895_792_439_5e0,
+    0.045_502_913_049_921_788_909_870_584_752_7e0,
+    0.047_982_537_138_836_713_906_392_567_492_0e0,
+    0.050_277_679_080_715_671_963_325_259_334_4e0,
+    0.052_362_885_806_407_475_864_366_137_127_9e0,
+    0.054_251_129_888_545_490_145_437_059_878_8e0,
+    0.055_950_811_220_412_317_308_240_686_327_5e0,
+    0.057_437_116_361_567_832_853_582_693_939_4e0,
+    0.058_689_680_022_394_207_961_974_561_788_8e0,
+    0.059_720_340_324_174_059_979_092_919_325_6e0,
+    0.060_539_455_376_045_862_945_360_267_515_7e0,
+    0.061_128_509_717_053_048_305_860_301_342_9e0,
+    0.061_471_189_871_425_316_615_443_159_652_6e0,
+    0.061_580_818_067_832_935_075_982_424_006_7e0,
 ];
 
 // G25 weights: correspond to XGK indices 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25 (center)
 #[rustfmt::skip]
 static GK51_WG: [f64; 13] = [
-    0.0113937985010262879479029641132e0,
-    0.0263549866150321372619018152953e0,
-    0.0409391567013063126556234877116e0,
-    0.0549046959758351919254159136154e0,
-    0.0680383338123569172071871856566e0,
-    0.0801407003350010180132349596696e0,
-    0.0910282619829636498114972207028e0,
-    0.1005359490670506442022068613869e0,
-    0.1085196244742636531160929150051e0,
-    0.1148582591457116483393255458696e0,
-    0.1194557635357847722281781265290e0,
-    0.1222424429903100416889595984585e0,
-    0.1231760537267154512039028730905e0,
+    0.011_393_798_501_026_287_947_902_964_113_2e0,
+    0.026_354_986_615_032_137_261_901_815_295_3e0,
+    0.040_939_156_701_306_312_655_623_487_711_6e0,
+    0.054_904_695_975_835_191_925_415_913_615_4e0,
+    0.068_038_333_812_356_917_207_187_185_656_6e0,
+    0.080_140_700_335_001_018_013_234_959_669_6e0,
+    0.091_028_261_982_963_649_811_497_220_702_8e0,
+    0.100_535_949_067_050_644_202_206_861_386_9e0,
+    0.108_519_624_474_263_653_116_092_915_005_1e0,
+    0.114_858_259_145_711_648_339_325_545_869_6e0,
+    0.119_455_763_535_784_772_228_178_126_529_0e0,
+    0.122_242_442_990_310_041_688_959_598_458_5e0,
+    0.123_176_053_726_715_451_203_902_873_090_5e0,
 ];
 
 #[cfg(test)]
diff --git a/src/gauss_laguerre.rs b/src/gauss_laguerre.rs
index df81eb0..58913e1 100644
--- a/src/gauss_laguerre.rs
+++ b/src/gauss_laguerre.rs
@@ -41,6 +41,12 @@ impl GaussLaguerre {
     /// Create a new n-point Gauss-Laguerre rule with parameter α.
     ///
     /// Requires `n >= 1`, `α > -1`.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
+    /// Returns [`QuadratureError::InvalidInput`] if `alpha <= -1` or `alpha`
+    /// is NaN.
     pub fn new(n: usize, alpha: f64) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -59,42 +65,22 @@ impl GaussLaguerre {
     }
 
     /// Returns the α parameter.
+    #[must_use]
     pub fn alpha(&self) -> f64 {
         self.alpha
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[0, ∞).
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussLaguerre, nodes_doc: "Returns the nodes on \\[0, ∞).");
+
 /// Compute n Gauss-Laguerre nodes and weights via Golub-Welsch.
 ///
 /// Monic generalised Laguerre recurrence:
-///   x L̃_k = L̃_{k+1} + (2k+1+α) L̃_k + k(k+α) L̃_{k-1}
+///   x `L̃_k` = `L̃_{k+1}` + (2k+1+α) `L̃_k` + k(k+α) `L̃_{k-1}`
 ///
 /// Jacobi matrix: diagonal = 2k+1+α, off-diagonal² = (k+1)(k+1+α).
 /// μ₀ = Γ(α+1).
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_laguerre(n: usize, alpha: f64) -> (Vec<f64>, Vec<f64>) {
     let diag: Vec<f64> = (0..n).map(|k| 2.0 * k as f64 + 1.0 + alpha).collect();
     let off_diag_sq: Vec<f64> = (1..n)
diff --git a/src/gauss_legendre.rs b/src/gauss_legendre.rs
index 15799da..6a5fcdf 100644
--- a/src/gauss_legendre.rs
+++ b/src/gauss_legendre.rs
@@ -41,6 +41,10 @@ impl GaussLegendre {
     /// Create a new n-point Gauss-Legendre rule.
     ///
     /// Returns an error if `n == 0`.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -52,35 +56,15 @@ impl GaussLegendre {
         })
     }
 
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
-
     /// Create a new n-point Gauss-Legendre rule with parallel node generation.
     ///
     /// Only beneficial for large n (> ~1000) where the Bogaert asymptotic
     /// path is used. For small n, the overhead of thread scheduling exceeds
     /// the computation time.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     #[cfg(feature = "parallel")]
     pub fn new_par(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
@@ -94,6 +78,8 @@ impl GaussLegendre {
     }
 }
 
+impl_rule_accessors!(GaussLegendre, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 // ---------------------------------------------------------------------------
 // Threshold: for n <= this value, use Newton refinement on the three-term
 // recurrence. For n > this value, use Bogaert's asymptotic expansion.
@@ -126,7 +112,9 @@ fn compute_gl_pair(n: usize) -> (Vec<f64>, Vec<f64>) {
 /// Compute GL nodes/weights via Newton iteration on the Legendre recurrence.
 ///
 /// The initial guess for each root uses the Tricomi approximation:
-///   theta_k ≈ pi * (4k - 1) / (4n + 2)
+///   `theta_k` ≈ pi * (4k - 1) / (4n + 2)
+#[allow(clippy::many_single_char_names)] // n, m, k, x, w are conventional in quadrature
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_newton(n: usize, m: usize, nodes: &mut [f64], weights: &mut [f64]) {
     let nf = n as f64;
 
@@ -160,8 +148,9 @@ fn compute_newton(n: usize, m: usize, nodes: &mut [f64], weights: &mut [f64]) {
     }
 }
 
-/// Evaluate the Legendre polynomial P_n(x) and its derivative P_n'(x)
+/// Evaluate the Legendre polynomial `P_n(x)` and its derivative `P_n'(x)`
 /// using the three-term recurrence.
+#[allow(clippy::cast_precision_loss)] // n and k are quadrature orders, always small enough for exact f64
 pub(crate) fn legendre_eval(n: usize, x: f64) -> (f64, f64) {
     let mut p_prev = 1.0; // P_0(x)
     let mut p_curr = x; // P_1(x)
@@ -191,35 +180,36 @@ pub(crate) fn legendre_eval(n: usize, x: f64) -> (f64, f64) {
 // Tabulated zeros of J_0(x) for k = 1..20
 #[rustfmt::skip]
 const BESSEL_J0_ZEROS: [f64; 20] = [
-     2.40482555769577276862163187933,   5.52007811028631064959660411281,
-     8.65372791291101221695419871266,  11.7915344390142816137430449119,
-    14.9309177084877859477625939974,  18.0710639679109225431478829756,
-    21.2116366298792589590783933505,  24.3524715307493027370579447632,
-    27.4934791320402547958772882346,  30.6346064684319751175495789269,
-    33.7758202135735686842385463467,  36.9170983536640439797694930633,
-    40.0584257646282392947993073740,  43.1997917131767303575240727287,
-    46.3411883716618140186857888791,  49.4826098973978171736027615332,
-    52.6240518411149960292512853804,  55.7655107550199793116834927735,
-    58.9069839260809421328344066346,  62.0484691902271698828525002646,
+     2.404_825_557_695_772_768_621_631_879_33,   5.520_078_110_286_310_649_596_604_112_81,
+     8.653_727_912_911_012_216_954_198_712_66,  11.791_534_439_014_281_613_743_044_911_9,
+    14.930_917_708_487_785_947_762_593_997_4,  18.071_063_967_910_922_543_147_882_975_6,
+    21.211_636_629_879_258_959_078_393_350_5,  24.352_471_530_749_302_737_057_944_763_2,
+    27.493_479_132_040_254_795_877_288_234_6,  30.634_606_468_431_975_117_549_578_926_9,
+    33.775_820_213_573_568_684_238_546_346_7,  36.917_098_353_664_043_979_769_493_063_3,
+    40.058_425_764_628_239_294_799_307_374_0,  43.199_791_713_176_730_357_524_072_728_7,
+    46.341_188_371_661_814_018_685_788_879_1,  49.482_609_897_397_817_173_602_761_533_2,
+    52.624_051_841_114_996_029_251_285_380_4,  55.765_510_755_019_979_311_683_492_773_5,
+    58.906_983_926_080_942_132_834_406_634_6,  62.048_469_190_227_169_882_852_500_264_6,
 ];
 
 // Tabulated J_1(j_{0,k})^2 for k = 1..21
 #[rustfmt::skip]
 const BESSEL_J1_SQUARED: [f64; 21] = [
-    0.269514123941916926139021992911, 0.115780138582203695807812836182,
-    0.0736863511364082151406476811985, 0.0540375731981162820417749182758,
-    0.0426614290172430912655106063495, 0.0352421034909961013587473033648,
-    0.0300210701030546726750888157688, 0.0261473914953080885904584675399,
-    0.0231591218246913922652676382178, 0.0207838291222678576039808057297,
-    0.0188504506693176678161056800214, 0.0172461575696650082995240053542,
-    0.0158935181059235978027065594287, 0.0147376260964721895895742982592,
-    0.0137384651453871179182880484134, 0.0128661817376151328791406637228,
-    0.0120980515486267975471075438497, 0.0114164712244916085168627222986,
-    0.0108075927911802040115547286830, 0.0102603729262807628110423992790,
-    0.00976589713979105054059846736696,
+    0.269_514_123_941_916_926_139_021_992_911, 0.115_780_138_582_203_695_807_812_836_182,
+    0.073_686_351_136_408_215_140_647_681_198_5, 0.054_037_573_198_116_282_041_774_918_275_8,
+    0.042_661_429_017_243_091_265_510_606_349_5, 0.035_242_103_490_996_101_358_747_303_364_8,
+    0.030_021_070_103_054_672_675_088_815_768_8, 0.026_147_391_495_308_088_590_458_467_539_9,
+    0.023_159_121_824_691_392_265_267_638_217_8, 0.020_783_829_122_267_857_603_980_805_729_7,
+    0.018_850_450_669_317_667_816_105_680_021_4, 0.017_246_157_569_665_008_299_524_005_354_2,
+    0.015_893_518_105_923_597_802_706_559_428_7, 0.014_737_626_096_472_189_589_574_298_259_2,
+    0.013_738_465_145_387_117_918_288_048_413_4, 0.012_866_181_737_615_132_879_140_663_722_8,
+    0.012_098_051_548_626_797_547_107_543_849_7, 0.011_416_471_224_491_608_516_862_722_298_6,
+    0.010_807_592_791_180_204_011_554_728_683_0, 0.010_260_372_926_280_762_811_042_399_279_0,
+    0.009_765_897_139_791_050_540_598_467_366_96,
 ];
 
-/// Compute the k-th zero of J_0(x), k >= 1.
+/// Compute the k-th zero of `J_0(x)`, k >= 1.
+#[allow(clippy::cast_precision_loss)] // k is a node index, always small enough for exact f64
 fn bessel_j0_zero(k: usize) -> f64 {
     if k <= 20 {
         BESSEL_J0_ZEROS[k - 1]
@@ -230,18 +220,20 @@ fn bessel_j0_zero(k: usize) -> f64 {
         let r2 = r * r;
         z + r
             * (0.125
-                + r2 * (-0.807291666666666666666666666667e-1
-                    + r2 * (0.246028645833333333333333333333
-                        + r2 * (-1.82443876720610119047619047619
-                            + r2 * (25.3364147973439050099206349206
-                                + r2 * (-567.644412135183381139802038240
-                                    + r2 * (18690.4765282320653831636345064
-                                        + r2 * (-8.49353580299148769921876983660e5
-                                            + 5.09225462402226769498681286758e7 * r2))))))))
+                + r2 * (-0.807_291_666_666_666_666_666_666_666_667e-1
+                    + r2 * (0.246_028_645_833_333_333_333_333_333_333
+                        + r2 * (-1.824_438_767_206_101_190_476_190_476_19
+                            + r2 * (25.336_414_797_343_905_009_920_634_920_6
+                                + r2 * (-567.644_412_135_183_381_139_802_038_240
+                                    + r2 * (18_690.476_528_232_065_383_163_634_506_4
+                                        + r2 * (-8.493_535_802_991_487_699_218_769_836_60e5
+                                            + 5.092_254_624_022_267_694_986_812_867_58e7
+                                                * r2))))))))
     }
 }
 
-/// Compute J_1(j_{0,k})^2, k >= 1.
+/// Compute `J_1(j_{0,k})^2`, k >= 1.
+#[allow(clippy::cast_precision_loss)] // k is a node index, always small enough for exact f64
 fn bessel_j1_squared(k: usize) -> f64 {
     if k <= 21 {
         BESSEL_J1_SQUARED[k - 1]
@@ -249,21 +241,22 @@ fn bessel_j1_squared(k: usize) -> f64 {
         // Asymptotic expansion
         let x = 1.0 / (k as f64 - 0.25);
         let x2 = x * x;
-        x * (0.202642367284675542887758926420
+        x * (0.202_642_367_284_675_542_887_758_926_420
             + x2 * x2
-                * (-0.303380429711290253026202643516e-3
-                    + x2 * (0.198924364245969295201137972743e-3
-                        + x2 * (-0.228969902772111653038747229723e-3
-                            + x2 * (0.433710719130746277915572905025e-3
-                                + x2 * (-0.123632349727175414724737657367e-2
-                                    + x2 * (0.496101423268883102872271417616e-2
-                                        + x2 * (-0.266837393702323757700998557826e-1
-                                            + 0.185395398206345628711318848386 * x2))))))))
+                * (-0.303_380_429_711_290_253_026_202_643_516e-3
+                    + x2 * (0.198_924_364_245_969_295_201_137_972_743e-3
+                        + x2 * (-0.228_969_902_772_111_653_038_747_229_723e-3
+                            + x2 * (0.433_710_719_130_746_277_915_572_905_025e-3
+                                + x2 * (-0.123_632_349_727_175_414_724_737_657_367e-2
+                                    + x2 * (0.496_101_423_268_883_102_872_271_417_616e-2
+                                        + x2 * (-0.266_837_393_702_323_757_700_998_557_826e-1
+                                            + 0.185_395_398_206_345_628_711_318_848_386 * x2))))))))
     }
 }
 
 /// Compute a single (theta, weight) pair using Bogaert's asymptotic expansion.
 /// k is 1-indexed, referring to the k-th node from the left (k <= ceil(n/2)).
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn bogaert_pair(n: usize, k: usize) -> (f64, f64) {
     let w = 1.0 / (n as f64 + 0.5);
     let nu = bessel_j0_zero(k);
@@ -271,96 +264,98 @@ fn bogaert_pair(n: usize, k: usize) -> (f64, f64) {
     let x = theta * theta;
 
     // Chebyshev interpolants for node correction
-    let sf1t = (((((-1.29052996274280508473467968379e-12 * x
-        + 2.40724685864330121825976175184e-10)
-        * x
-        - 3.13148654635992041468855740012e-8)
-        * x
-        + 0.275573168962061235623801563453e-5)
-        * x
-        - 0.148809523713909147898955880165e-3)
-        * x
-        + 0.416666666665193394525296923981e-2)
-        * x
-        - 0.416666666666662959639712457549e-1;
-
-    let sf2t =
-        (((((2.20639421781871003734786884322e-9 * x - 7.53036771373769326811030753538e-8) * x
-            + 0.161969259453836261731700382098e-5)
-            * x
-            - 0.253300326008232025914059965302e-4)
-            * x
-            + 0.282116886057560434805998583817e-3)
-            * x
-            - 0.209022248387852902722635654229e-2)
-            * x
-            + 0.815972221772932265640401128517e-2;
-
-    let sf3t =
-        (((((-2.97058225375526229899781956673e-8 * x + 5.55845330223796209655886325712e-7) * x
-            - 0.567797841356833081642185432056e-5)
-            * x
-            + 0.418498100329504574443885193835e-4)
-            * x
-            - 0.251395293283965914823026348764e-3)
-            * x
-            + 0.128654198542845137196151147483e-2)
-            * x
-            - 0.416012165620204364833694266818e-2;
+    let sf1t = (((((-1.290_529_962_742_805_084_734_679_683_79e-12 * x
+        + 2.407_246_858_643_301_218_259_761_751_84e-10)
+        * x
+        - 3.131_486_546_359_920_414_688_557_400_12e-8)
+        * x
+        + 0.275_573_168_962_061_235_623_801_563_453e-5)
+        * x
+        - 0.148_809_523_713_909_147_898_955_880_165e-3)
+        * x
+        + 0.416_666_666_665_193_394_525_296_923_981e-2)
+        * x
+        - 0.416_666_666_666_662_959_639_712_457_549e-1;
+
+    let sf2t = (((((2.206_394_217_818_710_037_347_868_843_22e-9 * x
+        - 7.530_367_713_737_693_268_110_307_535_38e-8)
+        * x
+        + 0.161_969_259_453_836_261_731_700_382_098e-5)
+        * x
+        - 0.253_300_326_008_232_025_914_059_965_302e-4)
+        * x
+        + 0.282_116_886_057_560_434_805_998_583_817e-3)
+        * x
+        - 0.209_022_248_387_852_902_722_635_654_229e-2)
+        * x
+        + 0.815_972_221_772_932_265_640_401_128_517e-2;
+
+    let sf3t = (((((-2.970_582_253_755_262_298_997_819_566_73e-8 * x
+        + 5.558_453_302_237_962_096_558_863_257_12e-7)
+        * x
+        - 0.567_797_841_356_833_081_642_185_432_056e-5)
+        * x
+        + 0.418_498_100_329_504_574_443_885_193_835e-4)
+        * x
+        - 0.251_395_293_283_965_914_823_026_348_764e-3)
+        * x
+        + 0.128_654_198_542_845_137_196_151_147_483e-2)
+        * x
+        - 0.416_012_165_620_204_364_833_694_266_818e-2;
 
     // Chebyshev interpolants for weight correction
-    let wsf1t = ((((((((-2.20902861044616638398573427475e-14 * x
-        + 2.30365726860377376873232578871e-12)
+    let wsf1t = ((((((((-2.209_028_610_446_166_383_985_734_274_75e-14 * x
+        + 2.303_657_268_603_773_768_732_325_788_71e-12)
         * x
-        - 1.75257700735423807659851042318e-10)
+        - 1.752_577_007_354_238_076_598_510_423_18e-10)
         * x
-        + 1.03756066927916795821098009353e-8)
+        + 1.037_560_669_279_167_958_210_980_093_53e-8)
         * x
-        - 4.63968647553221331251529631098e-7)
+        - 4.639_686_475_532_213_312_515_296_310_98e-7)
         * x
-        + 0.149644593625028648361395938176e-4)
+        + 0.149_644_593_625_028_648_361_395_938_176e-4)
         * x
-        - 0.326278659594412170300449074873e-3)
+        - 0.326_278_659_594_412_170_300_449_074_873e-3)
         * x
-        + 0.436507936507598105249726413120e-2)
+        + 0.436_507_936_507_598_105_249_726_413_120e-2)
         * x
-        - 0.305555555555553028279487898503e-1)
+        - 0.305_555_555_555_553_028_279_487_898_503e-1)
         * x
-        + 0.833333333333333302184063103900e-1;
+        + 0.833_333_333_333_333_302_184_063_103_900e-1;
 
-    let wsf2t = (((((((3.63117412152654783455929483029e-12 * x
-        + 7.67643545069893130779501844323e-11)
+    let wsf2t = (((((((3.631_174_121_526_547_834_559_294_830_29e-12 * x
+        + 7.676_435_450_698_931_307_795_018_443_23e-11)
         * x
-        - 7.12912857233642220650643150625e-9)
+        - 7.129_128_572_336_422_206_506_431_506_25e-9)
         * x
-        + 2.11483880685947151466370130277e-7)
+        + 2.114_838_806_859_471_514_663_701_302_77e-7)
         * x
-        - 0.381817918680045468483009307090e-5)
+        - 0.381_817_918_680_045_468_483_009_307_090e-5)
         * x
-        + 0.465969530694968391417927388162e-4)
+        + 0.465_969_530_694_968_391_417_927_388_162e-4)
         * x
-        - 0.407297185611335764191683161117e-3)
+        - 0.407_297_185_611_335_764_191_683_161_117e-3)
         * x
-        + 0.268959435694729660779984493795e-2)
+        + 0.268_959_435_694_729_660_779_984_493_795e-2)
         * x
-        - 0.111111111111214923138249347172e-1;
+        - 0.111_111_111_111_214_923_138_249_347_172e-1;
 
-    let wsf3t = (((((((2.01826791256703301806643264922e-9 * x
-        - 4.38647122520206649251063212545e-8)
+    let wsf3t = (((((((2.018_267_912_567_033_018_066_432_649_22e-9 * x
+        - 4.386_471_225_202_066_492_510_632_125_45e-8)
         * x
-        + 5.08898347288671653137451093208e-7)
+        + 5.088_983_472_886_716_531_374_510_932_08e-7)
         * x
-        - 0.397933316519135275712977531366e-5)
+        - 0.397_933_316_519_135_275_712_977_531_366e-5)
         * x
-        + 0.200559326396458326778521795392e-4)
+        + 0.200_559_326_396_458_326_778_521_795_392e-4)
         * x
-        - 0.422888059282921161626339411388e-4)
+        - 0.422_888_059_282_921_161_626_339_411_388e-4)
         * x
-        - 0.105646050254076140548678457002e-3)
+        - 0.105_646_050_254_076_140_548_678_457_002e-3)
         * x
-        - 0.947969308958577323145923317955e-4)
+        - 0.947_969_308_958_577_323_145_923_317_955e-4)
         * x
-        + 0.656966489926484797412985260842e-2;
+        + 0.656_966_489_926_484_797_412_985_260_842e-2;
 
     // Apply the expansions
     let nu_over_sin = nu / theta.sin();
@@ -380,6 +375,8 @@ fn bogaert_pair(n: usize, k: usize) -> (f64, f64) {
 /// Uses the same Newton/Bogaert split as the sequential version, but
 /// parallelises the per-node computation using rayon.
 #[cfg(feature = "parallel")]
+#[allow(clippy::many_single_char_names)] // n, m, k, x, w are conventional in quadrature
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_gl_pair_par(n: usize) -> (Vec<f64>, Vec<f64>) {
     use rayon::prelude::*;
 
diff --git a/src/gauss_lobatto.rs b/src/gauss_lobatto.rs
index c11a929..4f15645 100644
--- a/src/gauss_lobatto.rs
+++ b/src/gauss_lobatto.rs
@@ -40,6 +40,10 @@ impl GaussLobatto {
     /// Create a new n-point Gauss-Lobatto rule.
     ///
     /// Requires `n >= 2` (must include both endpoints).
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::InvalidInput`] if `n` is less than 2.
     pub fn new(n: usize) -> Result<Self, QuadratureError> {
         if n < 2 {
             return Err(QuadratureError::InvalidInput(
@@ -52,36 +56,15 @@ impl GaussLobatto {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussLobatto, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 /// Compute n-point Gauss-Lobatto nodes and weights.
 ///
 /// Nodes include -1 and 1. Interior nodes are zeros of P'_{n-1}(x).
-/// Weights: w_k = 2 / (n(n-1) [P_{n-1}(x_k)]^2).
+/// Weights: `w_k` = 2 / (n(n-1) \[`P_{n-1}(x_k)`\]^2).
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_lobatto(n: usize) -> (Vec<f64>, Vec<f64>) {
     let n_f = n as f64;
     let nm1 = n - 1;
diff --git a/src/gauss_radau.rs b/src/gauss_radau.rs
index b22c2e1..d89db18 100644
--- a/src/gauss_radau.rs
+++ b/src/gauss_radau.rs
@@ -40,6 +40,10 @@ impl GaussRadau {
     /// Create an n-point Gauss-Radau rule including the left endpoint -1.
     ///
     /// Requires `n >= 1`.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     pub fn left(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -53,6 +57,10 @@ impl GaussRadau {
     /// Create an n-point Gauss-Radau rule including the right endpoint +1.
     ///
     /// Requires `n >= 1`.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `n` is zero.
     pub fn right(n: usize) -> Result<Self, QuadratureError> {
         if n == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -65,37 +73,16 @@ impl GaussRadau {
             rule: QuadratureRule { nodes, weights },
         })
     }
-
-    /// Returns a reference to the underlying quadrature rule.
-    #[inline]
-    pub fn rule(&self) -> &QuadratureRule<f64> {
-        &self.rule
-    }
-
-    /// Returns the number of quadrature points.
-    #[inline]
-    pub fn order(&self) -> usize {
-        self.rule.order()
-    }
-
-    /// Returns the nodes on \[-1, 1\].
-    #[inline]
-    pub fn nodes(&self) -> &[f64] {
-        &self.rule.nodes
-    }
-
-    /// Returns the weights.
-    #[inline]
-    pub fn weights(&self) -> &[f64] {
-        &self.rule.weights
-    }
 }
 
+impl_rule_accessors!(GaussRadau, nodes_doc: "Returns the nodes on \\[-1, 1\\].");
+
 /// Compute n-point left Gauss-Radau nodes and weights (includes -1).
 ///
 /// Uses the Golub-Welsch algorithm with a Radau modification of the
 /// Legendre Jacobi matrix. The last diagonal element is modified so
 /// that -1 is an eigenvalue of the tridiagonal matrix.
+#[allow(clippy::cast_precision_loss)] // n is a quadrature order, always small enough for exact f64
 fn compute_radau_left(n: usize) -> (Vec<f64>, Vec<f64>) {
     if n == 1 {
         return (vec![-1.0], vec![2.0]);
diff --git a/src/golub_welsch.rs b/src/golub_welsch.rs
index d11f0b9..1df1631 100644
--- a/src/golub_welsch.rs
+++ b/src/golub_welsch.rs
@@ -1,9 +1,9 @@
 //! Golub-Welsch algorithm for computing Gaussian quadrature rules
 //! from three-term recurrence coefficients.
 //!
-//! Given the monic recurrence: x p_k = p_{k+1} + α_k p_k + β_k p_{k-1},
-//! the Jacobi matrix has diagonal α_k and off-diagonal √β_{k+1}.
-//! Eigenvalues = quadrature nodes, weights = μ₀ · z_k² where z_k is the
+//! Given the monic recurrence: x `p_k` = `p_{k+1}` + `α_k` `p_k` + `β_k` `p_{k-1}`,
+//! the Jacobi matrix has diagonal `α_k` and off-diagonal √`β_{k+1}`.
+//! Eigenvalues = quadrature nodes, weights = μ₀ · `z_k²` where `z_k` is the
 //! first component of the k-th eigenvector and μ₀ is the zeroth moment.
 //!
 //! Reference: Golub & Welsch (1969), "Calculation of Gauss Quadrature Rules",
@@ -93,6 +93,7 @@ pub(crate) fn radau_modify(diag: &mut [f64], off_diag_sq: &[f64], x0: f64) {
 ///
 /// On entry: `d` = diagonal, `e` = off-diagonal (length n-1).
 /// On exit: `d` = eigenvalues (unsorted), `z` = first row of eigenvector matrix.
+#[allow(clippy::many_single_char_names)] // d, e, z, m, l, g, r, s, c, p, f, b are standard names in tridiagonal eigenvalue algorithms
 fn symmetric_tridiag_eig(d: &mut [f64], e: &mut [f64], z: &mut [f64]) {
     let n = d.len();
     if n <= 1 {
diff --git a/src/lib.rs b/src/lib.rs
index c2948d3..ca73adc 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -93,10 +93,108 @@
 //! assert!((result - 1.0 / 3.0).abs() < 1e-14);
 //! # }
 //! ```
+//!
+//! ## Tanh-Sinh (Double Exponential)
+//!
+//! ```
+//! use bilby::tanh_sinh_integrate;
+//!
+//! // Endpoint singularity: integral of 1/sqrt(x) over [0, 1] = 2
+//! let result = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap();
+//! assert!((result.value - 2.0).abs() < 1e-7);
+//! ```
+//!
+//! ## Cauchy Principal Value
+//!
+//! ```
+//! use bilby::pv_integrate;
+//!
+//! // PV integral of x^2/(x - 0.3) over [0, 1]
+//! let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln();
+//! let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap();
+//! assert!((result.value - exact).abs() < 1e-7);
+//! ```
+//!
+//! ## Oscillatory Integration
+//!
+//! ```
+//! use bilby::integrate_oscillatory_sin;
+//!
+//! // Integral of sin(100x) over [0, 1]
+//! let exact = (1.0 - 100.0_f64.cos()) / 100.0;
+//! let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 100.0, 1e-10).unwrap();
+//! assert!((result.value - exact).abs() < 1e-8);
+//! ```
+//!
+//! ## Weighted Integration
+//!
+//! ```
+//! use bilby::weighted::{weighted_integrate, WeightFunction};
+//!
+//! // Gauss-Hermite: integral of e^(-x^2) over (-inf, inf) = sqrt(pi)
+//! let result = weighted_integrate(|_| 1.0, WeightFunction::Hermite, 20).unwrap();
+//! assert!((result - core::f64::consts::PI.sqrt()).abs() < 1e-12);
+//! ```
+//!
+//! ## Sparse Grids
+//!
+//! ```
+//! use bilby::cubature::SparseGrid;
+//!
+//! let sg = SparseGrid::clenshaw_curtis(3, 3).unwrap();
+//! let result = sg.rule().integrate_box(
+//!     &[0.0, 0.0, 0.0], &[1.0, 1.0, 1.0],
+//!     |x| x[0] * x[1] * x[2],
+//! );
+//! assert!((result - 0.125).abs() < 1e-12);
+//! ```
 
 #[cfg(not(feature = "std"))]
 extern crate alloc;
 
+/// Implement the standard quadrature rule accessor methods.
+///
+/// Most 1-D rule types wrap a `QuadratureRule<f64>` in a field called `rule`
+/// and expose the same four accessors: `rule()`, `order()`, `nodes()`, and
+/// `weights()`. This macro generates all four with the correct attributes
+/// and documentation.
+///
+/// `$nodes_doc` customises the doc comment on `nodes()` to describe the
+/// domain (e.g. "\\[-1, 1\\]", "\\[0, ∞)", "(-∞, ∞)").
+macro_rules! impl_rule_accessors {
+    ($ty:ty, nodes_doc: $nodes_doc:expr) => {
+        impl $ty {
+            /// Returns a reference to the underlying quadrature rule.
+            #[inline]
+            #[must_use]
+            pub fn rule(&self) -> &$crate::rule::QuadratureRule<f64> {
+                &self.rule
+            }
+
+            /// Returns the number of quadrature points.
+            #[inline]
+            #[must_use]
+            pub fn order(&self) -> usize {
+                self.rule.order()
+            }
+
+            #[doc = $nodes_doc]
+            #[inline]
+            #[must_use]
+            pub fn nodes(&self) -> &[f64] {
+                &self.rule.nodes
+            }
+
+            /// Returns the weights.
+            #[inline]
+            #[must_use]
+            pub fn weights(&self) -> &[f64] {
+                &self.rule.weights
+            }
+        }
+    };
+}
+
 pub mod adaptive;
 #[cfg(feature = "std")]
 pub mod cache;
diff --git a/src/oscillatory.rs b/src/oscillatory.rs
index bcb8617..36ab26b 100644
--- a/src/oscillatory.rs
+++ b/src/oscillatory.rs
@@ -60,6 +60,7 @@ pub struct OscillatoryIntegrator {
 
 impl OscillatoryIntegrator {
     /// Create a new oscillatory integrator.
+    #[must_use]
     pub fn new(kernel: OscillatoryKernel, omega: f64) -> Self {
         Self {
             kernel,
@@ -71,24 +72,36 @@ impl OscillatoryIntegrator {
     }
 
     /// Set the Chebyshev expansion order.
+    #[must_use]
     pub fn with_order(mut self, n: usize) -> Self {
         self.order = n;
         self
     }
 
     /// Set absolute tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set relative tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Integrate f(x)·kernel(ω·x) over \[a, b\].
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
+    /// is NaN. 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
+    #[allow(clippy::cast_precision_loss)] // k, n are small quadrature indices, always exact in f64
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -220,9 +233,10 @@ impl OscillatoryIntegrator {
 
 /// Compute Chebyshev coefficients from function values at CC nodes.
 ///
-/// Given f_k = f(cos(kπ/n)) for k = 0, ..., n, computes the coefficients
-/// c_j such that f(t) ≈ Σ c_j T_j(t) (with the convention that c_0 and
-/// c_n are halved).
+/// Given `f_k` = f(cos(kπ/n)) for k = 0, ..., n, computes the coefficients
+/// `c_j` such that f(t) ≈ Σ `c_j` `T_j(t)` (with the convention that `c_0` and
+/// `c_n` are halved).
+#[allow(clippy::cast_precision_loss)] // j, k, n are small Chebyshev indices, always exact in f64
 fn chebyshev_coefficients(f_vals: &[f64], n: usize) -> Vec<f64> {
     let mut coeffs = vec![0.0; n + 1];
     let pi_n = core::f64::consts::PI / n as f64;
@@ -245,13 +259,15 @@ fn chebyshev_coefficients(f_vals: &[f64], n: usize) -> Vec<f64> {
 
 /// Compute modified Chebyshev moments via Gauss-Legendre quadrature.
 ///
-/// Returns (μ_j^c, μ_j^s) for j = 0, ..., n where:
-///   μ_j^c = ∫₋₁¹ T_j(x) cos(θx) dx
-///   μ_j^s = ∫₋₁¹ T_j(x) sin(θx) dx
+/// Returns (`μ_j^c`, `μ_j^s`) for j = 0, ..., n where:
+///   `μ_j^c` = ∫₋₁¹ `T_j(x)` cos(θx) dx
+///   `μ_j^s` = ∫₋₁¹ `T_j(x)` sin(θx) dx
 ///
 /// Uses GL quadrature with enough nodes to resolve both the Chebyshev
 /// polynomial (degree j ≤ n) and the oscillatory kernel (effective
 /// frequency ~θ). Converges exponentially for these smooth integrands.
+#[allow(clippy::cast_precision_loss)] // j, n are small Chebyshev indices, always exact in f64
+#[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] // theta.abs().ceil() is a small positive integer, safe to cast to usize
 fn modified_chebyshev_moments(theta: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
     let mut mc = vec![0.0; n + 1];
     let mut ms = vec![0.0; n + 1];
@@ -306,6 +322,21 @@ fn modified_chebyshev_moments(theta: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
 }
 
 /// Convenience: integrate f(x)·sin(ωx) over \[a, b\].
+///
+/// # Example
+///
+/// ```
+/// use bilby::integrate_oscillatory_sin;
+///
+/// let exact = (1.0 - 100.0_f64.cos()) / 100.0;
+/// let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 100.0, 1e-10).unwrap();
+/// assert!((result.value - exact).abs() < 1e-8);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
+/// is NaN.
 pub fn integrate_oscillatory_sin<G>(
     f: G,
     a: f64,
@@ -323,6 +354,21 @@ where
 }
 
 /// Convenience: integrate f(x)·cos(ωx) over \[a, b\].
+///
+/// # Example
+///
+/// ```
+/// use bilby::integrate_oscillatory_cos;
+///
+/// let exact = 50.0_f64.sin() / 50.0;
+/// let result = integrate_oscillatory_cos(|_| 1.0, 0.0, 1.0, 50.0, 1e-10).unwrap();
+/// assert!((result.value - exact).abs() < 1e-8);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::DegenerateInterval`] if `a`, `b`, or `omega`
+/// is NaN.
 pub fn integrate_oscillatory_cos<G>(
     f: G,
     a: f64,
diff --git a/src/rule.rs b/src/rule.rs
index fe873d3..b0aeaba 100644
--- a/src/rule.rs
+++ b/src/rule.rs
@@ -20,6 +20,7 @@ pub struct QuadratureRule<F> {
 impl<F: Float> QuadratureRule<F> {
     /// Returns the number of points in this rule.
     #[inline]
+    #[must_use]
     pub fn order(&self) -> usize {
         self.nodes.len()
     }
@@ -30,6 +31,11 @@ impl<F: Float> QuadratureRule<F> {
     ///   x = (b - a) / 2 * t + (a + b) / 2
     ///
     /// The result is scaled by (b - a) / 2.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the conversion of `2.0` to `F` via `F::from` fails, which
+    /// should not occur for standard floating-point types.
     #[inline]
     pub fn integrate<G>(&self, a: F, b: F, f: G) -> F
     where
@@ -51,6 +57,12 @@ impl<F: Float> QuadratureRule<F> {
     ///
     /// Each panel applies this quadrature rule independently, then sums the results.
     /// Total function evaluations = `self.order() * n_panels`.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `n_panels` or the panel indices cannot be converted to `F`
+    /// via `F::from`, which should not occur for standard floating-point types
+    /// with reasonable panel counts.
     pub fn integrate_composite<G>(&self, a: F, b: F, n_panels: usize, f: G) -> F
     where
         G: Fn(F) -> F,
@@ -71,6 +83,12 @@ impl<F: Float> QuadratureRule<F> {
     ///
     /// Identical to [`integrate_composite`](Self::integrate_composite) but evaluates
     /// panels in parallel using rayon. Requires `F: Send + Sync` and `f: Fn(F) -> F + Sync`.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `n_panels` or the panel indices cannot be converted to `F`
+    /// via `F::from`, which should not occur for standard floating-point types
+    /// with reasonable panel counts.
     #[cfg(feature = "parallel")]
     pub fn integrate_composite_par<G>(&self, a: F, b: F, n_panels: usize, f: G) -> F
     where
diff --git a/src/tanh_sinh.rs b/src/tanh_sinh.rs
index 66433b7..c844eb6 100644
--- a/src/tanh_sinh.rs
+++ b/src/tanh_sinh.rs
@@ -55,24 +55,35 @@ impl Default for TanhSinh {
 
 impl TanhSinh {
     /// Set the maximum number of refinement levels.
+    #[must_use]
     pub fn with_max_levels(mut self, levels: usize) -> Self {
         self.max_levels = levels;
         self
     }
 
     /// Set absolute tolerance.
+    #[must_use]
     pub fn with_abs_tol(mut self, tol: f64) -> Self {
         self.abs_tol = tol;
         self
     }
 
     /// Set relative tolerance.
+    #[must_use]
     pub fn with_rel_tol(mut self, tol: f64) -> Self {
         self.rel_tol = tol;
         self
     }
 
     /// Integrate `f` over [a, b] using the tanh-sinh transform.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is NaN.
+    #[allow(clippy::many_single_char_names)] // a, b, f, h, k, t, u are conventional in tanh-sinh quadrature
+    #[allow(clippy::cast_precision_loss)] // k is a small loop index, always exact in f64
+    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)] // (7.0/h).ceil() is a small positive value, safe to cast to usize
+    #[allow(clippy::too_many_lines)] // single cohesive level-doubling loop, splitting would obscure the algorithm
     pub fn integrate<G>(
         &self,
         a: f64,
@@ -245,6 +256,20 @@ impl TanhSinh {
 }
 
 /// Convenience: tanh-sinh integration with default settings.
+///
+/// # Example
+///
+/// ```
+/// use bilby::tanh_sinh_integrate;
+///
+/// // Integral of 1/sqrt(x) over [0, 1] = 2 (endpoint singularity)
+/// let result = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap();
+/// assert!((result.value - 2.0).abs() < 1e-7);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::DegenerateInterval`] if `a` or `b` is NaN.
 pub fn tanh_sinh_integrate<G>(
     f: G,
     a: f64,
diff --git a/src/weighted.rs b/src/weighted.rs
index 4cd1a9a..610007a 100644
--- a/src/weighted.rs
+++ b/src/weighted.rs
@@ -59,6 +59,10 @@ pub struct WeightedIntegrator {
 
 impl WeightedIntegrator {
     /// Create a new weighted integrator.
+    ///
+    /// # Errors
+    ///
+    /// Returns [`QuadratureError::ZeroOrder`] if `order` is zero.
     pub fn new(weight: WeightFunction, order: usize) -> Result<Self, QuadratureError> {
         if order == 0 {
             return Err(QuadratureError::ZeroOrder);
@@ -67,6 +71,7 @@ impl WeightedIntegrator {
     }
 
     /// Set the quadrature order.
+    #[must_use]
     pub fn with_order(mut self, n: usize) -> Self {
         self.order = n;
         self
@@ -78,7 +83,14 @@ impl WeightedIntegrator {
     /// - Laguerre: \[0, ∞)
     /// - Hermite: (-∞, ∞)
     /// - ChebyshevI/II: \[-1, 1\]
-    /// - LogWeight: (0, 1\]
+    /// - `LogWeight`: (0, 1\]
+    ///
+    /// # 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.
     pub fn integrate<G>(&self, f: G) -> f64
     where
         G: Fn(f64) -> f64,
@@ -147,7 +159,7 @@ impl WeightedIntegrator {
     /// Integrate f(x) · w(x) over \[a, b\] via affine transform.
     ///
     /// Only applicable for finite-domain weights (Jacobi, ChebyshevI/II).
-    /// For LogWeight, maps (0, 1\] to (a, b\] with the log singularity at a.
+    /// For `LogWeight`, maps (0, 1\] to (a, b\] with the log singularity at a.
     /// For Laguerre/Hermite, a and b are ignored (uses natural domain).
     pub fn integrate_over<G>(&self, a: f64, b: f64, f: G) -> f64
     where
@@ -182,6 +194,20 @@ impl WeightedIntegrator {
 }
 
 /// Convenience: integrate f(x) · w(x) over the natural domain.
+///
+/// # Example
+///
+/// ```
+/// use bilby::weighted::{weighted_integrate, WeightFunction};
+///
+/// // Integral of 1/sqrt(1-x^2) over [-1,1] = pi
+/// let result = weighted_integrate(|_| 1.0, WeightFunction::ChebyshevI, 20).unwrap();
+/// assert!((result - core::f64::consts::PI).abs() < 1e-12);
+/// ```
+///
+/// # Errors
+///
+/// Returns [`QuadratureError::ZeroOrder`] if `order` is zero.
 pub fn weighted_integrate<G>(
     f: G,
     weight: WeightFunction,
@@ -294,4 +320,101 @@ mod tests {
     fn zero_order() {
         assert!(weighted_integrate(|_| 1.0, WeightFunction::Hermite, 0).is_err());
     }
+
+    #[test]
+    fn integrate_over_jacobi_legendre_on_0_2() {
+        // Jacobi(0,0) is Legendre weight w(x)=1.
+        // integrate_over maps [-1,1] to [0,2] via affine transform with Jacobian (b-a)/2 = 1.
+        // ∫₀² x² dx = 8/3
+        let wi = WeightedIntegrator::new(
+            WeightFunction::Jacobi {
+                alpha: 0.0,
+                beta: 0.0,
+            },
+            20,
+        )
+        .unwrap();
+        let result = wi.integrate_over(0.0, 2.0, |x| x * x);
+        assert!(
+            (result - 8.0 / 3.0).abs() < 1e-10,
+            "result={result}, expected={}",
+            8.0 / 3.0
+        );
+    }
+
+    #[test]
+    fn integrate_over_jacobi_legendre_on_1_3() {
+        // ∫₁³ x dx = (9 - 1)/2 = 4
+        let wi = WeightedIntegrator::new(
+            WeightFunction::Jacobi {
+                alpha: 0.0,
+                beta: 0.0,
+            },
+            10,
+        )
+        .unwrap();
+        let result = wi.integrate_over(1.0, 3.0, |x| x);
+        assert!((result - 4.0).abs() < 1e-10, "result={result}");
+    }
+
+    #[test]
+    fn integrate_over_chebyshev_i_on_0_2() {
+        // ChebyshevI weight: w(t) = 1/√(1-t²) on [-1,1].
+        // integrate_over(0, 2, f) maps t -> x = t + 1, so x ∈ [0,2].
+        // Result = half * ∫₋₁¹ f(half*t + mid) / √(1-t²) dt
+        //   where half = 1, mid = 1.
+        // With f(x) = 1: result = 1 * ∫₋₁¹ 1/√(1-t²) dt = π
+        let wi = WeightedIntegrator::new(WeightFunction::ChebyshevI, 20).unwrap();
+        let result = wi.integrate_over(0.0, 2.0, |_| 1.0);
+        assert!(
+            (result - core::f64::consts::PI).abs() < 1e-10,
+            "result={result}"
+        );
+    }
+
+    #[test]
+    fn integrate_over_chebyshev_ii_on_0_2() {
+        // ChebyshevII weight: w(t) = √(1-t²) on [-1,1].
+        // integrate_over(0, 2, f=1) = half * ∫₋₁¹ √(1-t²) dt = 1 * π/2
+        let wi = WeightedIntegrator::new(WeightFunction::ChebyshevII, 20).unwrap();
+        let result = wi.integrate_over(0.0, 2.0, |_| 1.0);
+        assert!(
+            (result - core::f64::consts::FRAC_PI_2).abs() < 1e-10,
+            "result={result}"
+        );
+    }
+
+    #[test]
+    fn integrate_over_laguerre_ignores_bounds() {
+        // Laguerre uses its natural domain [0,∞); bounds are ignored.
+        let wi = WeightedIntegrator::new(WeightFunction::Laguerre { alpha: 0.0 }, 10).unwrap();
+        let natural = wi.integrate(|_| 1.0);
+        let remapped = wi.integrate_over(5.0, 10.0, |_| 1.0);
+        assert!(
+            (natural - remapped).abs() < 1e-14,
+            "natural={natural}, remapped={remapped}"
+        );
+    }
+
+    #[test]
+    fn integrate_over_hermite_ignores_bounds() {
+        // Hermite uses its natural domain (-∞,∞); bounds are ignored.
+        let wi = WeightedIntegrator::new(WeightFunction::Hermite, 10).unwrap();
+        let natural = wi.integrate(|_| 1.0);
+        let remapped = wi.integrate_over(5.0, 10.0, |_| 1.0);
+        assert!(
+            (natural - remapped).abs() < 1e-14,
+            "natural={natural}, remapped={remapped}"
+        );
+    }
+
+    #[test]
+    fn integrate_over_log_weight_on_0_2() {
+        // LogWeight: integrate_over(0, 2, f) = width * integrate(|u| f(a + width*u))
+        //   = 2 * ∫₀¹ -ln(u) · f(2u) du
+        // With f(x) = 1: result = 2 * ∫₀¹ -ln(u) du = 2 * 1 = 2
+        let wi = WeightedIntegrator::new(WeightFunction::LogWeight, 20).unwrap();
+        let result = wi.integrate_over(0.0, 2.0, |_| 1.0);
+        assert!((result - 2.0).abs() < 1e-10, "result={result}");
+    }
 }