From 5cccb7277a66ef1beb65a5d83831875bb7bbb735 Mon Sep 17 00:00:00 2001 From: Robert Mongold Date: Sun, 12 Jul 2026 16:03:33 -0400 Subject: [PATCH] feature: Add opt-in Kahan compensated summation Note: Adds a public summation field on IterativeConfig (default Pairwise). Call sites that build the struct with a full literal need the new field; Default, from_parameters, and with_summation are unchanged in behavior for existing users. --- CHANGELOG.md | 2 + crates/multicalc/src/approximation/README.md | 1 + .../src/approximation/linear_approximation.rs | 24 ++- crates/multicalc/src/approximation/mod.rs | 12 +- .../approximation/quadratic_approximation.rs | 29 ++- .../src/numerical_integration/README.md | 1 + .../iterative_integration.rs | 176 +++++++++++++----- .../src/numerical_integration/mod.rs | 2 + crates/multicalc/src/utils/summation.rs | 109 ++++++++++- crates/multicalc/tests/approximation.rs | 26 +++ .../multicalc/tests/numerical_integration.rs | 24 +++ 11 files changed, 347 insertions(+), 59 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c668999..237e1d0 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/crates/multicalc/src/approximation/README.md b/crates/multicalc/src/approximation/README.md index 7bfa8b3..1bb658d 100644 --- a/crates/multicalc/src/approximation/README.md +++ b/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; diff --git a/crates/multicalc/src/approximation/linear_approximation.rs b/crates/multicalc/src/approximation/linear_approximation.rs index e7dc940..a428ced 100644 --- a/crates/multicalc/src/approximation/linear_approximation.rs +++ b/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 { 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 LinearApproximation { /// 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, const NUM_POINTS: usize>( @@ -70,6 +76,7 @@ impl LinearApproximation { points, &|x: &[T; NUM_VARS]| original_function.eval(x), NUM_VARS, // p = N linear coefficients + self.summation, ); LinearApproximationPredictionMetrics { @@ -86,12 +93,14 @@ impl LinearApproximation { /// ([`AutoDiffMulti`]); pass a finite-difference derivator explicitly to use that instead. pub struct LinearApproximator { derivator: D, + summation: SummationMethod, } impl Default for LinearApproximator { fn default() -> Self { LinearApproximator { derivator: D::default(), + summation: SummationMethod::Pairwise, } } } @@ -99,7 +108,19 @@ impl Default for LinearApproximator { impl LinearApproximator { /// 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 LinearApproximator { point: *point, value, gradient, + summation: self.summation, }) } } diff --git a/crates/multicalc/src/approximation/mod.rs b/crates/multicalc/src/approximation/mod.rs index 07bfdc2..f193608 100644 --- a/crates/multicalc/src/approximation/mod.rs +++ b/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, 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; diff --git a/crates/multicalc/src/approximation/quadratic_approximation.rs b/crates/multicalc/src/approximation/quadratic_approximation.rs index 267e729..97b81f8 100644 --- a/crates/multicalc/src/approximation/quadratic_approximation.rs +++ b/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 { 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,6 +59,10 @@ impl QuadraticApproximation { /// 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, const NUM_POINTS: usize>( @@ -64,14 +70,12 @@ impl QuadraticApproximation { points: &[[T; NUM_VARS]; NUM_POINTS], original_function: &O, ) -> QuadraticApproximationPredictionMetrics { - // 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,12 +92,14 @@ impl QuadraticApproximation { /// autodiff ([`AutoDiffMulti`]); pass a finite-difference derivator explicitly to use that instead. pub struct QuadraticApproximator { derivator: D, + summation: SummationMethod, } impl Default for QuadraticApproximator { fn default() -> Self { QuadraticApproximator { derivator: D::default(), + summation: SummationMethod::Pairwise, } } } @@ -101,7 +107,19 @@ impl Default for QuadraticApproximator { impl QuadraticApproximator { /// 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 QuadraticApproximator { value, gradient, hessian, + summation: self.summation, }) } } diff --git a/crates/multicalc/src/numerical_integration/README.md b/crates/multicalc/src/numerical_integration/README.md index eba8243..982d0f6 100644 --- a/crates/multicalc/src/numerical_integration/README.md +++ b/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` diff --git a/crates/multicalc/src/numerical_integration/iterative_integration.rs b/crates/multicalc/src/numerical_integration/iterative_integration.rs index 559dd41..84d7ecc 100644 --- a/crates/multicalc/src/numerical_integration/iterative_integration.rs +++ b/crates/multicalc/src/numerical_integration/iterative_integration.rs @@ -4,7 +4,7 @@ use crate::numerical_integration::integrator::*; use crate::numerical_integration::mode::IterativeMethod; use crate::scalar::Numeric; use crate::utils::error_codes::CalcError; -use crate::utils::summation::PairwiseSum; +use crate::utils::summation::{Acc, SummationMethod}; /// Default interval count. A multiple of 12 so Boole (needs a multiple of 4) and /// Simpson 3/8 (needs a multiple of 3) both align with the composite-rule weights. @@ -17,6 +17,8 @@ pub struct IterativeConfig { pub total_iterations: u64, /// The composite rule to use: Booles, Simpsons or Trapezoidal. pub integration_method: IterativeMethod, + /// Running sum algorithm for the composite-rule accumulators. Default: pairwise. + pub summation: SummationMethod, } impl Default for IterativeConfig { @@ -25,6 +27,7 @@ impl Default for IterativeConfig { IterativeConfig { total_iterations: DEFAULT_TOTAL_ITERATIONS, integration_method: IterativeMethod::Booles, + summation: SummationMethod::Pairwise, } } } @@ -35,6 +38,7 @@ impl IterativeConfig { IterativeConfig { total_iterations, integration_method, + summation: SummationMethod::Pairwise, } } @@ -66,20 +70,27 @@ fn integrate_rule T>( lo: T, hi: T, g: G, + summation: SummationMethod, ) -> T { match method { - IterativeMethod::Booles => booles(iterations, lo, hi, g), - IterativeMethod::Simpsons => simpsons(iterations, lo, hi, g), - IterativeMethod::Trapezoidal => trapezoidal(iterations, lo, hi, g), + IterativeMethod::Booles => booles(iterations, lo, hi, g, summation), + IterativeMethod::Simpsons => simpsons(iterations, lo, hi, g, summation), + IterativeMethod::Trapezoidal => trapezoidal(iterations, lo, hi, g, summation), } } /// Boole's composite rule over `[lo, hi]`. -fn booles T>(iterations: u64, lo: T, hi: T, mut g: G) -> T { +fn booles T>( + iterations: u64, + lo: T, + hi: T, + mut g: G, + summation: SummationMethod, +) -> T { let delta = (hi - lo) / T::from_u64(iterations); let mut point = lo; - let mut ans = PairwiseSum::new(); + let mut ans = Acc::new(summation); ans.add(T::from_f64(7.0) * g(point)); let mut multiplier = T::from_f64(32.0); @@ -102,11 +113,17 @@ fn booles T>(iterations: u64, lo: T, hi: T, mut g: G) } /// Simpson's 3/8 composite rule over `[lo, hi]`. -fn simpsons T>(iterations: u64, lo: T, hi: T, mut g: G) -> T { +fn simpsons T>( + iterations: u64, + lo: T, + hi: T, + mut g: G, + summation: SummationMethod, +) -> T { let delta = (hi - lo) / T::from_u64(iterations); let mut point = lo; - let mut ans = PairwiseSum::new(); + let mut ans = Acc::new(summation); ans.add(g(point)); let mut multiplier = T::from_f64(3.0); @@ -127,11 +144,17 @@ fn simpsons T>(iterations: u64, lo: T, hi: T, mut g: } /// Trapezoidal composite rule over `[lo, hi]`. -fn trapezoidal T>(iterations: u64, lo: T, hi: T, mut g: G) -> T { +fn trapezoidal T>( + iterations: u64, + lo: T, + hi: T, + mut g: G, + summation: SummationMethod, +) -> T { let delta = (hi - lo) / T::from_u64(iterations); let mut point = lo; - let mut ans = PairwiseSum::new(); + let mut ans = Acc::new(summation); ans.add(g(point)); for _ in 0..iterations - 1 { @@ -168,6 +191,15 @@ impl IterativeSingle { _marker: PhantomData, } } + + /// Selects the running-sum algorithm for the composite-rule accumulators. + /// + /// Default is [`SummationMethod::Pairwise`]. Pass [`SummationMethod::Kahan`] + /// for opt-in compensated summation. + pub fn with_kahan_summation(mut self) -> Self { + self.config.summation = SummationMethod::Kahan; + self + } } impl IterativeSingle { @@ -183,6 +215,7 @@ impl IterativeSingle { ) -> T { let method = self.config.integration_method; let iterations = self.config.total_iterations; + let summation = self.config.summation; let domain = match classify(&integration_limit[level - 1]) { Ok(d) => d, @@ -191,26 +224,40 @@ impl IterativeSingle { if level == 1 { return match domain { - Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, func), + Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, func, summation), _ => { let (lo, hi) = t_bounds(&domain); - integrate_rule(method, iterations, lo, hi, |t| { - let (x, jacobian) = map_sample(&domain, t); - func(x) * jacobian - }) + integrate_rule( + method, + iterations, + lo, + hi, + |t| { + let (x, jacobian) = map_sample(&domain, t); + func(x) * jacobian + }, + summation, + ) } }; } let inner = self.integrate(level - 1, func, integration_limit); match domain { - Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, |_| inner), + Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, |_| inner, summation), _ => { let (lo, hi) = t_bounds(&domain); - integrate_rule(method, iterations, lo, hi, |t| { - let (_, jacobian) = map_sample(&domain, t); - inner * jacobian - }) + integrate_rule( + method, + iterations, + lo, + hi, + |t| { + let (_, jacobian) = map_sample(&domain, t); + inner * jacobian + }, + summation, + ) } } } @@ -288,6 +335,12 @@ impl IterativeMulti { _marker: PhantomData, } } + + /// Opt in to Kahan compensated summation for the composite-rule accumulators. + pub fn with_kahan_summation(mut self) -> Self { + self.config.summation = SummationMethod::Kahan; + self + } } impl IterativeMulti { @@ -309,6 +362,7 @@ impl IterativeMulti { ) -> T { let method = self.config.integration_method; let iterations = self.config.total_iterations; + let summation = self.config.summation; let domain = match classify(&integration_limits[level - 1]) { Ok(d) => d, @@ -319,47 +373,75 @@ impl IterativeMulti { if level == 1 { let mut current = *point; return match domain { - Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, |x| { - current[var] = x; - func(¤t) - }), + Domain::Finite(a, b) => integrate_rule( + method, + iterations, + a, + b, + |x| { + current[var] = x; + func(¤t) + }, + summation, + ), _ => { let (lo, hi) = t_bounds(&domain); - integrate_rule(method, iterations, lo, hi, |t| { - let (x, jacobian) = map_sample(&domain, t); - current[var] = x; - func(¤t) * jacobian - }) + integrate_rule( + method, + iterations, + lo, + hi, + |t| { + let (x, jacobian) = map_sample(&domain, t); + current[var] = x; + func(¤t) * jacobian + }, + summation, + ) } }; } let mut current = *point; match domain { - Domain::Finite(a, b) => integrate_rule(method, iterations, a, b, |x| { - current[var] = x; - self.integrate( - level - 1, - idx_to_integrate, - func, - integration_limits, - ¤t, - ) - }), - _ => { - let (lo, hi) = t_bounds(&domain); - integrate_rule(method, iterations, lo, hi, |t| { - let (x, jacobian) = map_sample(&domain, t); + Domain::Finite(a, b) => integrate_rule( + method, + iterations, + a, + b, + |x| { current[var] = x; - let inner = self.integrate( + self.integrate( level - 1, idx_to_integrate, func, integration_limits, ¤t, - ); - inner * jacobian - }) + ) + }, + summation, + ), + _ => { + let (lo, hi) = t_bounds(&domain); + integrate_rule( + method, + iterations, + lo, + hi, + |t| { + let (x, jacobian) = map_sample(&domain, t); + current[var] = x; + let inner = self.integrate( + level - 1, + idx_to_integrate, + func, + integration_limits, + ¤t, + ); + inner * jacobian + }, + summation, + ) } } } diff --git a/crates/multicalc/src/numerical_integration/mod.rs b/crates/multicalc/src/numerical_integration/mod.rs index 750f05d..2fe91cd 100644 --- a/crates/multicalc/src/numerical_integration/mod.rs +++ b/crates/multicalc/src/numerical_integration/mod.rs @@ -2,3 +2,5 @@ pub mod gaussian_integration; pub mod integrator; pub mod iterative_integration; pub mod mode; + +pub use crate::utils::summation::SummationMethod; diff --git a/crates/multicalc/src/utils/summation.rs b/crates/multicalc/src/utils/summation.rs index ce29a1f..9d4acab 100644 --- a/crates/multicalc/src/utils/summation.rs +++ b/crates/multicalc/src/utils/summation.rs @@ -1,4 +1,9 @@ -//! Blocked pairwise (cascade) summation for the long running sums. +//! Running-sum accumulators for long numeric series. +//! +//! - [`PairwiseSum`] — blocked pairwise (cascade) summation; the default in iterative +//! integration and approximation metrics (O(log n · eps) error growth). +//! - [`KahanSum`] — classic compensated summation; opt-in via +//! `.with_kahan_summation()` on the iterative integrators / approximators. use crate::scalar::Numeric; @@ -81,6 +86,82 @@ impl PairwiseSum { } } +/// Which running-sum algorithm the iterative integrators / metrics use. +/// +/// [`SummationMethod::Pairwise`] is the default. Enable Kahan with +/// `.with_kahan_summation()` on the integrator or approximator. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum SummationMethod { + /// Blocked pairwise (cascade) summation - the historical default. + #[default] + Pairwise, + /// Classic Kahan compensated summation. + Kahan, +} + +/// Classic Kahan compensated summation. +/// +/// Tracks a running compensation term so each add recovers the low-order bits that +/// would otherwise be lost when adding a small value into a large partial sum. +pub(crate) struct KahanSum { + sum: T, + compensation: T, +} + +impl KahanSum { + #[inline] + pub(crate) fn new() -> Self { + Self { + sum: T::ZERO, + compensation: T::ZERO, + } + } + + #[inline] + pub(crate) fn add(&mut self, value: T) { + let y = value - self.compensation; + let t = self.sum + y; + self.compensation = (t - self.sum) - y; + self.sum = t; + } + + #[inline] + pub(crate) fn total(&self) -> T { + self.sum + } +} + +pub(crate) enum Acc { + Pairwise(PairwiseSum), + Kahan(KahanSum), +} + +impl Acc { + #[inline] + pub(crate) fn new(method: SummationMethod) -> Self { + match method { + SummationMethod::Pairwise => Acc::Pairwise(PairwiseSum::new()), + SummationMethod::Kahan => Acc::Kahan(KahanSum::new()), + } + } + + #[inline] + pub(crate) fn add(&mut self, v: T) { + match self { + Acc::Pairwise(a) => a.add(v), + Acc::Kahan(a) => a.add(v), + } + } + + #[inline] + pub(crate) fn total(self) -> T { + match self { + Acc::Pairwise(a) => a.total(), + Acc::Kahan(a) => a.total(), + } + } +} + #[cfg(test)] mod tests { use super::*; @@ -152,4 +233,30 @@ mod tests { "pairwise ({pairwise_err:e}) should be strictly closer than naive ({naive_err:e})" ); } + + #[test] + fn kahan_beats_naive_on_long_sum() { + // 1.0 then N tiny terms (half-ulp at 1.0). Naive stays at 1.0; Kahan recovers the mass. + let tiny = 2f64.powi(-53); + let n: u64 = 1 << 24; + let analytic = 1.0 + (n as f64) * tiny; //1.0 + 2^-29 + + let mut kahan = KahanSum::new(); + kahan.add(1.0); + let mut naive = 1.0f64; + for _ in 0..n { + kahan.add(tiny); + naive += tiny; + } + let got = kahan.total(); + + assert_eq!(naive, 1.0, "naive should lose every tiny term"); + let kahan_err = (got - analytic).abs(); + let naive_err = (naive - analytic).abs(); + assert!(kahan_err < 1e-12, "kahan error {kahan_err:e} too large"); + assert!( + kahan_err < naive_err, + "kahan ({kahan_err:e}) should be strictly closer than naive ({naive_err:e})" + ); + } } diff --git a/crates/multicalc/tests/approximation.rs b/crates/multicalc/tests/approximation.rs index 0f9cd4a..cd08989 100644 --- a/crates/multicalc/tests/approximation.rs +++ b/crates/multicalc/tests/approximation.rs @@ -104,6 +104,32 @@ fn test_linear_approximation_exact() { assert!(f64::abs(metrics.adjusted_r_squared - 1.0) < 1e-9); } +#[test] +fn kahan_metrics_exact_on_affine() { + // Same affine truth as the pairwise exactness test: Kahan path must also report + // a perfect fit so the opt-in wire-in is exercised without changing defaults. + let truth = scalar_fn!(|v: &[f64; 3]| c(5.0) + c(2.0) * v[0] + c(3.0) * v[1] - v[2]); + let point = [1.0, 2.0, 3.0]; + + let model = LinearApproximator::::default() + .with_kahan_summation() + .get(&truth, &point) + .unwrap(); + + let mut points = [[0.0; 3]; 10]; + for (i, p) in points.iter_mut().enumerate() { + let s = i as f64; + *p = [1.0 + s, 2.0 - s, 3.0 + 0.5 * s]; + } + + let metrics = model.get_prediction_metrics(&points, &truth); + + assert!(metrics.mean_absolute_error < 1e-9); + assert!(metrics.root_mean_squared_error < 1e-9); + assert!((metrics.r_squared - 1.0).abs() < 1e-9); + assert!((metrics.adjusted_r_squared - 1.0).abs() < 1e-9); +} + #[test] fn metrics_are_accurate_on_large_point_set() { //truth is x^2, so a linear approximation about `a` has the exact residual -(x - a)^2. diff --git a/crates/multicalc/tests/numerical_integration.rs b/crates/multicalc/tests/numerical_integration.rs index 178fad5..0454a42 100644 --- a/crates/multicalc/tests/numerical_integration.rs +++ b/crates/multicalc/tests/numerical_integration.rs @@ -592,3 +592,27 @@ fn test_gauss_legendre_integration_f32() { let val = integrator.get_single(&func, &[0.0, 2.0]).unwrap(); assert!(f32::abs(val - 8.0) < 1e-2, "got {val}"); } + +#[test] +fn kahan_integration_beats_naive_on_long_sum() { + let func = |x: f64| -> f64 { 1.0 / (1.0 + x * x) }; + let exact = core::f64::consts::PI / 4.0; + let iterations: u64 = 1 << 23; + + let integrator = iterative_integration::IterativeSingle::from_parameters( + iterations, + IterativeMethod::Trapezoidal, + ) + .with_kahan_summation(); + + let kahan = integrator.get_single(&func, &[0.0, 1.0]).unwrap(); + let naive = naive_trapezoidal(iterations, 0.0, 1.0, func); + + let kahan_err = f64::abs(kahan - exact); + let naive_err = f64::abs(naive - exact); + assert!(kahan_err < 1e-12, "kahan error {kahan_err:e} too large"); + assert!( + kahan_err < naive_err, + "kahan ({kahan_err:e}) should be closer than naive ({naive_err:e})" + ); +}