FIX: Correct Gauss sum for even-order Gauss-Kronrod rules (#93)

The evaluation-reusing path (compute_gauss_kronrod_sums_stored,
introduced in faeaf2b) reconstructed the Gauss sum assuming an
odd-order Gauss rule: a center node at 0.0 plus pairs at Kronrod
offsets +/- 2j. Even-order rules (G10/G20/G30) have no node at 0.0;
their nodes sit at odd Kronrod offsets +/- (2j - 1). The old code
therefore added a spurious center term (a Kronrod-only node times an
ordinary pair weight), dropped one node pair, and read every pair
from Kronrod-only nodes, so the "Gauss" sum contained no Gauss node
at all. |G - K| then decayed only linearly with interval width, so
G10K21 / G20K41 / G30K61 (and R variants) never early-exited and
always subdivided to max_iter, even for a constant integrand.

Branch on the parity of g: odd keeps the previous mapping, even skips
the center term, sums all g/2 pairs at offsets +/- (2j - 1), and
starts pair weights at index g/2 - 1.

Tests:
* unit: reconstructed G matches direct Gauss-Legendre and |G - K| is
  round-off for a cubic, for all six rules (failed on G10K21 before)
* integration: a cubic finishes in exactly k evaluations for all
  twelve method variants, and the optimized path agrees with the
  reference gauss_kronrod_quadrature on a generic smooth integrand

Author: Axect

src/numerical/integral.rs

@@ -500,29 +500,46 @@ where
     let K = kronrod_sum * xh;
 
     // 6. Calculate Gauss sum (G) using stored evaluations and Gauss weights
-    // Assumes G(g) nodes are subset of K(k), g is odd, standard mapping: +/- xg_j maps to +/- xk_{2j}.
+    // The Gauss nodes are a subset of the Kronrod nodes, but their position
+    // relative to the Kronrod center node (0.0) depends on the parity of g:
+    // - odd g: 0.0 is a Gauss node; +/- xg_j maps to Kronrod indices center +/- 2j
+    // - even g: 0.0 is a Kronrod-only node; +/- xg_j maps to center +/- (2j - 1)
     let kronrod_center_idx = k / 2; // Index of 0.0 node in Kronrod/f_evals
-    let gauss_center_idx = g / 2; // Index of 0.0 weight in Gauss weights
-    let num_gauss_pairs = (g - 1) / 2;
 
-    // Contribution from center node (0.0)
-    let mut gauss_sum = f_evals[kronrod_center_idx].clone() * gauss_weights[gauss_center_idx];
+    let mut gauss_sum;
+    let num_gauss_pairs;
+    let node_offset_parity; // Kronrod offset of pair j is 2j - node_offset_parity
+    let gauss_weight_base; // Gauss weight index of pair j is gauss_weight_base + j
+    if g % 2 == 1 {
+        // Center node contribution with the true center Gauss weight
+        gauss_sum = f_evals[kronrod_center_idx].clone() * gauss_weights[g / 2];
+        num_gauss_pairs = (g - 1) / 2;
+        node_offset_parity = 0;
+        gauss_weight_base = g / 2;
+    } else {
+        // Even-order Gauss rules have no node at 0.0: every node is paired
+        gauss_sum = Y::ZERO;
+        num_gauss_pairs = g / 2;
+        node_offset_parity = 1;
+        gauss_weight_base = g / 2 - 1;
+    }
 
     // Contribution from paired Gauss nodes (+/- xg_j)
     for j in 1..=num_gauss_pairs {
-        if 2 * j > kronrod_center_idx {
+        let offset = 2 * j - node_offset_parity;
+        if offset > kronrod_center_idx {
             panic!(
                 "Kronrod node index underflow. k={}, center_idx={}, j={}",
                 k, kronrod_center_idx, j
             );
         }
 
-        // Indices in f_evals corresponding to Kronrod nodes +/- xk_{2j}
-        let kronrod_node_pos_idx = kronrod_center_idx + 2 * j;
-        let kronrod_node_neg_idx = kronrod_center_idx - 2 * j;
+        // Indices in f_evals corresponding to the Kronrod nodes hosting +/- xg_j
+        let kronrod_node_pos_idx = kronrod_center_idx + offset;
+        let kronrod_node_neg_idx = kronrod_center_idx - offset;
 
         // Index in gauss_weights for the symmetric weight of the pair +/- xg_j
-        let gauss_weight_pair_idx = gauss_center_idx + j;
+        let gauss_weight_pair_idx = gauss_weight_base + j;
 
         // Safety checks
         if kronrod_node_pos_idx >= k || kronrod_node_neg_idx >= k {
@@ -2192,3 +2209,56 @@ const KRONROD_WEIGHTS_61: [f64; 61] = [
     0.003890461127099884,
     0.001389013698677008,
 ];
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    const GK_ORDERS: [(usize, usize); 6] =
+        [(7, 15), (10, 21), (15, 31), (20, 41), (25, 51), (30, 61)];
+
+    #[test]
+    fn gauss_sum_reconstruction_is_exact_for_low_degree_polynomials() {
+        // Both G and K integrate a cubic exactly, so the reconstructed Gauss
+        // sum must match the Kronrod sum to round-off. Regression test for #93:
+        // the even-order reconstruction added a spurious center term and used
+        // wrong node offsets, so |G - K| stayed O(1) instead of vanishing.
+        let f = |x: f64| 1.0 + x + x * x + x * x * x;
+        let exact = 1.0 + 0.5 + 1.0 / 3.0 + 0.25;
+        for (g, k) in GK_ORDERS {
+            let (g_sum, k_sum, err): (f64, f64, f64) =
+                compute_gauss_kronrod_sums_stored(f, (0.0, 1.0), g, k);
+            assert!(
+                (g_sum - k_sum).abs() < 1e-13,
+                "G{}K{}: G = {:.17e}, K = {:.17e}",
+                g,
+                k,
+                g_sum,
+                k_sum
+            );
+            assert!((k_sum - exact).abs() < 1e-13);
+            assert!(err.abs() < 1e-13);
+        }
+    }
+
+    #[test]
+    fn gauss_sum_reconstruction_matches_direct_gauss_legendre() {
+        // For an arbitrary smooth integrand the Gauss sum reconstructed from
+        // the stored Kronrod evaluations must equal the directly computed
+        // Gauss-Legendre quadrature of the same order.
+        let f = |x: f64| (1.5 * x).sin() + (-x).exp();
+        for (g, k) in GK_ORDERS {
+            let (g_sum, _, _): (f64, f64, f64) =
+                compute_gauss_kronrod_sums_stored(f, (0.3, 2.1), g, k);
+            let direct: f64 = gauss_legendre_quadrature(f, g, (0.3, 2.1));
+            assert!(
+                (g_sum - direct).abs() < 1e-13,
+                "G{}K{}: reconstructed = {:.17e}, direct = {:.17e}",
+                g,
+                k,
+                g_sum,
+                direct
+            );
+        }
+    }
+}

tests/integral.rs

@@ -95,6 +95,65 @@ fn test_oscillatory() {
     }
 }
 
+#[test]
+fn test_early_exit_on_low_degree_polynomials() {
+    // A G_gK_k rule integrates polynomials of degree <= 2g - 1 exactly, so the
+    // very first |G - K| check must pass and the whole integration must finish
+    // with exactly k function evaluations (one interval, no subdivision).
+    // Regression test for #93: the even-order rules (G10K21, G20K41, G30K61
+    // and their R variants) never early-exited and subdivided to max_iter.
+    use std::cell::Cell;
+    let methods = [
+        (G7K15(1e-8, 20), 15),
+        (G10K21(1e-8, 20), 21),
+        (G15K31(1e-8, 20), 31),
+        (G20K41(1e-8, 20), 41),
+        (G25K51(1e-8, 20), 51),
+        (G30K61(1e-8, 20), 61),
+        (G7K15R(1e-8, 20), 15),
+        (G10K21R(1e-8, 20), 21),
+        (G15K31R(1e-8, 20), 31),
+        (G20K41R(1e-8, 20), 41),
+        (G25K51R(1e-8, 20), 51),
+        (G30K61R(1e-8, 20), 61),
+    ];
+    for (m, k) in methods {
+        let count = Cell::new(0u32);
+        let f = |x: f64| {
+            count.set(count.get() + 1);
+            x.powi(3)
+        };
+        let r: f64 = integrate(f, (0.0, 1.0), m);
+        assert_close(r, 0.25);
+        assert_eq!(
+            count.get(),
+            k,
+            "{:?} did not early-exit on a cubic: {} evaluations instead of {}",
+            m,
+            count.get(),
+            k
+        );
+    }
+}
+
+#[test]
+fn test_optimized_matches_reference_gauss_kronrod() {
+    // The evaluation-reusing path must agree with the straightforward
+    // independent-evaluation path on a generic smooth integrand.
+    let f = |x: f64| (2.0 * x).sin() * (-x).exp();
+    for m in gk_methods(1e-10) {
+        let opt: f64 = gauss_kronrod_quadrature_optimized(f, (0.0, 2.0), m);
+        let reference: f64 = gauss_kronrod_quadrature(f, (0.0, 2.0), m);
+        assert!(
+            (opt - reference).abs() < 1e-9,
+            "{:?}: optimized = {:.17e}, reference = {:.17e}",
+            m,
+            opt,
+            reference
+        );
+    }
+}
+
 #[test]
 fn test_endpoint_independence() {
     let m = || G10K21R(1e-10, 20);