viva_math/precision
High-precision numerical primitives.
IEEE-754 doubles lose precision when summing values that differ widely in magnitude. This module provides compensated and exact summation algorithms plus Pébay’s online higher-order moment accumulators.
Algorithms shipped
| Function | Worst-case error | Cost | When to use |
|---|---|---|---|
neumaier_sum | O(ε) | ~3x | Default — fast + safe |
kahan_sum | O(ε) | ~4x | Reference / didactic |
pairwise_sum | O(log n · ε) | ~1x | Reductions on long lists |
fsum | round-once exact | ~10x | Critical accuracy |
neumaier_sum is what CPython 3.12 sum() uses by default and matches
the Gonum and Julia recommendations. Catastrophic example:
import viva_math/precision
precision.neumaier_sum([1.0, 1.0e100, 1.0, -1.0e100])
// -> 2.0 (naive sum would give 0.0)
References
- Kahan (1965) “Pracniques: further remarks on reducing truncation errors”
- Neumaier (1974) “Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen”
- Shewchuk (1997) “Adaptive Precision Floating-Point Arithmetic”
- Pébay (2008) “Formulas for robust, one-pass parallel computation of covariances and arbitrary-order statistical moments” (Sandia)
- Pébay, Terriberry et al. (2016) “Numerically Stable, Scalable Formulas for Parallel and Online Computation of Higher-Order Multivariate Central Moments with Arbitrary Weights”
- CPython issue #100425 (Python 3.12 switched to Neumaier)
Types
Online accumulator for moments up to fourth order.
Updates Pébay-style central moments without storing the data. Numerically stable through cancellation tricks; computes variance, skewness, and kurtosis from streaming input.
pub type Moments {
Moments(
count: Int,
mean: Float,
m2: Float,
m3: Float,
m4: Float,
)
}Constructors
-
Moments(count: Int, mean: Float, m2: Float, m3: Float, m4: Float)
Values
pub fn cmp_abs_desc(a: Float, b: Float) -> order.OrderStable cmp by absolute value (descending). Useful for sorting before summation in worst-case scenarios.
pub fn fsum(xs: List(Float)) -> FloatShewchuk’s fsum: maintains a list of non-overlapping partial sums.
Round-to-nearest exact result up to a final rounding. Matches Python’s
math.fsum. The most accurate option available, at ~10x the cost of
naive summation.
pub fn kahan_sum(xs: List(Float)) -> FloatClassical Kahan compensated sum.
Slightly less accurate than neumaier_sum (fails on the pathological
[1, 1e100, 1, -1e100] example) but easier to reason about. Useful as a
reference implementation.
pub fn moments_combine(a: Moments, b: Moments) -> MomentsCombine two Pébay accumulators computed in parallel (Chan formula).
Useful for distributed / parallel reductions: split a stream across workers, then merge.
pub fn moments_excess_kurtosis(m: Moments) -> Result(Float, Nil)Excess kurtosis γ₂ = n · M₄ / M₂² - 3.
pub fn moments_from_list(xs: List(Float)) -> MomentsBuild accumulator from a list.
pub fn moments_sample_variance(m: Moments) -> Result(Float, Nil)Sample variance s² = M₂ / (n-1).
pub fn moments_skewness(m: Moments) -> Result(Float, Nil)Population skewness γ₁ = (M₃/n) / (M₂/n)(3/2) = √n · M₃ / M₂(3/2).
pub fn moments_update(m: Moments, x: Float) -> MomentsUpdate with a single sample (Pébay 2008 recurrence).
Let n be the new count, δ = x - mean, δ_n = δ/n, δ_n² and term1 = δ·δ_n·(n-1). Then: M₂ ← M₂ + term1 M₃ ← M₃ + term1·δ_n·(n-2) - 3·δ_n·M₂ M₄ ← M₄ + term1·δ_n²·(n²-3n+3) + 6·δ_n²·M₂ - 4·δ_n·M₃
Updates the mean last to preserve the previous-iteration moments above.
pub fn moments_variance(m: Moments) -> Result(Float, Nil)Population variance σ² = M₂ / n.
pub fn neumaier_mean(xs: List(Float)) -> Result(Float, Nil)Mean using Neumaier-compensated sum. Errors on empty input.
pub fn neumaier_sum(xs: List(Float)) -> FloatNeumaier (Kahan-Babuška improved) compensated sum.
Maintains a separate compensation accumulator and swaps the role of the running sum / next term depending on magnitude. Recovers about 16 extra bits of precision over naive summation at ~3x the cost.
Example
precision.neumaier_sum([1.0, 1.0e100, 1.0, -1.0e100])
// -> 2.0 (naive sum returns 0.0)
pub fn pairwise_sum(xs: List(Float)) -> FloatPairwise sum: recursively splits and combines.
Error grows as O(log n · ε) instead of O(n · ε) for naive sum, at the same asymptotic cost. NumPy’s default. Good for very long lists where the constant factor of Neumaier matters.