feature: Add opt-in Kahan compensated summation

CHANGELOG.md

@@ -21,6 +21,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   matching-degree polynomials, metrics consistency). @rtmongold (#130)
 - Property tests that autodiff and central finite-difference derivatives agree on random
   polynomial compositions. @rtmongold (#129)
+- Opt-in Kahan compensated summation for iterative integration and approximation metrics
+  via `.with_kahan_summation()`; pairwise remains the default. @rtmongold (#134)
 
 ### Fixed
 

crates/multicalc/src/approximation/README.md

@@ -8,6 +8,7 @@ metrics.
   captures curvature.
 - `get` builds the model; `predict` evaluates it; `get_prediction_metrics` returns MAE, MSE, RMSE,
   R², and adjusted R² against sample points.
+- Metrics use pairwise summation by default; chain `.with_kahan_summation()` to opt into Kahan.
 
 ```rust
 use multicalc::approximation::linear_approximation::LinearApproximator;

crates/multicalc/src/approximation/linear_approximation.rs

@@ -3,6 +3,7 @@ use crate::numerical_derivative::autodiff::AutoDiffMulti;
 use crate::numerical_derivative::derivator::DerivatorMultiVariable;
 use crate::scalar::{Numeric, ScalarFnN};
 use crate::utils::error_codes::CalcError;
+use crate::utils::summation::SummationMethod;
 
 /// A first-order (linear) Taylor approximation of a function about a base point:
 /// `f(x) ≈ value + Σ gradient[i] * (x[i] - point[i])`.
@@ -11,6 +12,7 @@ pub struct LinearApproximation<const NUM_VARS: usize, T = f64> {
     point: [T; NUM_VARS],
     value: T,
     gradient: [T; NUM_VARS],
+    summation: SummationMethod,
 }
 
 /// Goodness-of-fit metrics for a [`LinearApproximation`] over a set of sample points.
@@ -58,6 +60,10 @@ impl<const NUM_VARS: usize, T: Numeric> LinearApproximation<NUM_VARS, T> {
 
     /// Computes goodness-of-fit metrics against `original_function` over `points`.
     ///
+    /// Uses the summation method chosen when the approximator was built
+    /// (pairwise by default; Kahan if [`LinearApproximator::with_kahan_summation`]
+    /// was used).
+    ///
     /// `r_squared` is `NaN` when the truth is constant over `points`;
     /// `adjusted_r_squared` is `NaN` when there are too few points.
     pub fn get_prediction_metrics<O: ScalarFnN<NUM_VARS>, const NUM_POINTS: usize>(
@@ -70,6 +76,7 @@ impl<const NUM_VARS: usize, T: Numeric> LinearApproximation<NUM_VARS, T> {
             points,
             &|x: &[T; NUM_VARS]| original_function.eval(x),
             NUM_VARS, // p = N linear coefficients
+            self.summation,
         );
 
         LinearApproximationPredictionMetrics {
@@ -86,20 +93,34 @@ impl<const NUM_VARS: usize, T: Numeric> LinearApproximation<NUM_VARS, T> {
 /// ([`AutoDiffMulti`]); pass a finite-difference derivator explicitly to use that instead.
 pub struct LinearApproximator<D: DerivatorMultiVariable = AutoDiffMulti> {
     derivator: D,
+    summation: SummationMethod,
 }
 
 impl<D: DerivatorMultiVariable + Default> Default for LinearApproximator<D> {
     fn default() -> Self {
         LinearApproximator {
             derivator: D::default(),
+            summation: SummationMethod::Pairwise,
         }
     }
 }
 
 impl<D: DerivatorMultiVariable> LinearApproximator<D> {
     /// Builds an approximator from an explicit derivator.
     pub fn from_derivator(derivator: D) -> Self {
-        LinearApproximator { derivator }
+        LinearApproximator {
+            derivator,
+            summation: SummationMethod::Pairwise,
+        }
+    }
+
+    /// Opt in to Kahan compensated summation for prediction metrics.
+    ///
+    /// Pairwise summation remains the default. Call this before [`Self::get`] so the
+    /// resulting [`LinearApproximation`] accumulates metrics with Kahan.
+    pub fn with_kahan_summation(mut self) -> Self {
+        self.summation = SummationMethod::Kahan;
+        self
     }
 
     /// Builds a linear (first-order Taylor) approximation of `function` about `point`.
@@ -139,6 +160,7 @@ impl<D: DerivatorMultiVariable> LinearApproximator<D> {
             point: *point,
             value,
             gradient,
+            summation: self.summation,
         })
     }
 }

crates/multicalc/src/approximation/mod.rs

@@ -1,5 +1,6 @@
 use crate::scalar::Numeric;
-use crate::utils::summation::PairwiseSum;
+use crate::utils::summation::Acc;
+pub use crate::utils::summation::SummationMethod;
 
 pub mod linear_approximation;
 pub mod quadratic_approximation;
@@ -15,6 +16,7 @@ pub(crate) fn compute_metrics<T, P, O, const NUM_VARS: usize, const NUM_POINTS:
     points: &[[T; NUM_VARS]; NUM_POINTS],
     original_function: &O,
     num_predictors: usize,
+    summation: SummationMethod,
 ) -> (T, T, T, T, T)
 where
     T: Numeric,
@@ -23,15 +25,15 @@ where
 {
     let n = T::from_usize(NUM_POINTS);
 
-    let mut mean_y = PairwiseSum::new();
+    let mut mean_y = Acc::new(summation);
     for point in points {
         mean_y.add(original_function(point));
     }
     let mean_y = mean_y.total() / n;
 
-    let mut sum_abs = PairwiseSum::new();
-    let mut ss_res = PairwiseSum::new();
-    let mut ss_tot = PairwiseSum::new();
+    let mut sum_abs = Acc::new(summation);
+    let mut ss_res = Acc::new(summation);
+    let mut ss_tot = Acc::new(summation);
     for point in points {
         let y = original_function(point);
         let residual = predict(point) - y;

crates/multicalc/src/approximation/quadratic_approximation.rs

@@ -3,6 +3,7 @@ use crate::numerical_derivative::autodiff::AutoDiffMulti;
 use crate::numerical_derivative::derivator::DerivatorMultiVariable;
 use crate::scalar::{Numeric, ScalarFnN};
 use crate::utils::error_codes::CalcError;
+use crate::utils::summation::SummationMethod;
 
 /// A second-order (quadratic) Taylor approximation of a function about a base point:
 /// `f(x) ≈ value + Σ gradient[i]·dx[i] + ½ Σ_i Σ_j hessian[i][j]·dx[i]·dx[j]`,
@@ -13,6 +14,7 @@ pub struct QuadraticApproximation<const NUM_VARS: usize, T = f64> {
     value: T,
     gradient: [T; NUM_VARS],
     hessian: [[T; NUM_VARS]; NUM_VARS],
+    summation: SummationMethod,
 }
 
 /// Goodness-of-fit metrics for a [`QuadraticApproximation`] over a set of sample points.
@@ -57,21 +59,23 @@ impl<const NUM_VARS: usize, T: Numeric> QuadraticApproximation<NUM_VARS, T> {
 
     /// Computes goodness-of-fit metrics against `original_function` over `points`.
     ///
+    /// Uses the summation method chosen when the approximator was built
+    /// (pairwise by default; Kahan if [`QuadraticApproximator::with_kahan_summation`]
+    /// was used).
+    ///
     /// `r_squared` is `NaN` when the truth is constant over `points`;
     /// `adjusted_r_squared` is `NaN` when there are too few points.
     pub fn get_prediction_metrics<O: ScalarFnN<NUM_VARS>, const NUM_POINTS: usize>(
         &self,
         points: &[[T; NUM_VARS]; NUM_POINTS],
         original_function: &O,
     ) -> QuadraticApproximationPredictionMetrics<T> {
-        // p = N gradient terms + N(N+1)/2 distinct (symmetric) Hessian terms
-        let num_predictors = NUM_VARS + NUM_VARS * (NUM_VARS + 1) / 2;
-
         let (mae, mse, rmse, r_squared, adjusted_r_squared) = crate::approximation::compute_metrics(
             |x| self.predict(x),
             points,
             &|x: &[T; NUM_VARS]| original_function.eval(x),
-            num_predictors,
+            NUM_VARS + NUM_VARS * (NUM_VARS + 1) / 2, // p = N gradient terms + N(N+1)/2 distinct (symmetric) Hessian terms
+            self.summation,
         );
 
         QuadraticApproximationPredictionMetrics {
@@ -88,20 +92,34 @@ impl<const NUM_VARS: usize, T: Numeric> QuadraticApproximation<NUM_VARS, T> {
 /// autodiff ([`AutoDiffMulti`]); pass a finite-difference derivator explicitly to use that instead.
 pub struct QuadraticApproximator<D: DerivatorMultiVariable = AutoDiffMulti> {
     derivator: D,
+    summation: SummationMethod,
 }
 
 impl<D: DerivatorMultiVariable + Default> Default for QuadraticApproximator<D> {
     fn default() -> Self {
         QuadraticApproximator {
             derivator: D::default(),
+            summation: SummationMethod::Pairwise,
         }
     }
 }
 
 impl<D: DerivatorMultiVariable> QuadraticApproximator<D> {
     /// Builds an approximator from an explicit derivator.
     pub fn from_derivator(derivator: D) -> Self {
-        QuadraticApproximator { derivator }
+        QuadraticApproximator {
+            derivator,
+            summation: SummationMethod::Pairwise,
+        }
+    }
+
+    /// Opt in to Kahan compensated summation for prediction metrics.
+    ///
+    /// Pairwise summation remains the default. Call this before [`Self::get`] so the
+    /// resulting [`QuadraticApproximation`] accumulates metrics with Kahan.
+    pub fn with_kahan_summation(mut self) -> Self {
+        self.summation = SummationMethod::Kahan;
+        self
     }
 
     /// Builds a quadratic (second-order Taylor) approximation of `function` about `point`.
@@ -158,6 +176,7 @@ impl<D: DerivatorMultiVariable> QuadraticApproximator<D> {
             value,
             gradient,
             hessian,
+            summation: self.summation,
         })
     }
 }

crates/multicalc/src/numerical_integration/README.md

@@ -5,6 +5,7 @@ semi-infinite, and infinite limits.
 
 - [`iterative_integration::IterativeSingle`](iterative_integration.rs) — Boole (default), Simpson, and
   Trapezoidal rules; pick the rule and interval count with `from_parameters`.
+- Pairwise summation is the default; chain `.with_kahan_summation()` to opt into Kahan.
 - [`gaussian_integration::GaussianSingle`](gaussian_integration.rs) — Gauss-Legendre, Gauss-Hermite,
   and Gauss-Laguerre. Pass the **bare** integrand; the weights already carry the weighting factor.
 - Both implement the [`integrator`](integrator.rs)`::IntegratorSingleVariable` / `MultiVariable`