Floating-point error

Deep dive

How it works

In: f64 (or f32) values to add, multiply-accumulate or compare · Out: a sum or dot product with a bounded error, or a robust comparison · State: running sum plus a correction term (compensated sums), none for pairwise · Cost: n adds for pairwise; about 4 flops per element for Kahan

Error model. fl(a ∘ b) = (a ∘ b)(1 + δ) with |δ| ≤ u for ∘ ∈ {+, −, ×, ÷}, as long as no overflow/underflow. Errors of a computation are bounded by backward error analysis: the result is the exact result for slightly perturbed inputs.

Cancellation. If a and b agree in their first k bits, a − b is exact (Sterbenz lemma when b/2 ≤ a ≤ 2b), but any error already in a and b is now relatively 2ᵏ times larger. Examples and fixes:

• 1 − cos x for small x → 2·sin²(x/2).
• √(x + 1) − √x → 1/(√(x + 1) + √x).
• Quadratic roots: compute q = −(b + sign(b)·√(b² − 4ac))/2, then x₁ = q/a, x₂ = c/q.
• ln(1 + x), eˣ − 1 → ln_1p, exp_m1.

Compensated summation. Kahan keeps a running correction c of the low-order bits lost when adding y to a large sum. Neumaier's variant also handles the case where the new term is larger than the sum. Error-free transformations generalize this: two_sum gives the exact rounding error of an addition, two_prod that of a product (via FMA), and chaining them gives doubled precision (double-double, "Sum2"/"Dot2" of Ogita-Rump-Oishi).

Pairwise summation. Recursively sum halves; depth log₂n means error O(u·log n) at the same cost as the naive loop. With a base case of 8-128 elements summed naively, speed equals the naive loop.

Code

/// Error-free sum: a + b = s + e exactly (Knuth's TwoSum, 6 flops, no branch).
#[inline] pub fn two_sum(a: f64, b: f64) -> (f64, f64) {
    let s = a + b;
    let bb = s - a;
    (s, (a - (s - bb)) + (b - bb))
}

/// Error-free product via fused multiply-add: a * b = p + e exactly (barring underflow).
#[inline] pub fn two_prod(a: f64, b: f64) -> (f64, f64) {
    let p = a * b;
    (p, a.mul_add(b, -p))
}

pub fn naive_sum(x: &[f64]) -> f64 { x.iter().sum() }

pub fn kahan_sum(x: &[f64]) -> f64 {
    let (mut s, mut c) = (0.0, 0.0);
    for &v in x {
        let y = v - c;
        let t = s + y;
        c = (t - s) - y; // what was lost in t = s + y (negated)
        s = t;
    }
    s
}

/// Neumaier (Kahan-Babuska): also correct when |v| > |s|.
pub fn neumaier_sum(x: &[f64]) -> f64 {
    let (mut s, mut c) = (0.0f64, 0.0f64);
    for &v in x {
        let t = s + v;
        c += if s.abs() >= v.abs() { (s - t) + v } else { (v - t) + s };
        s = t;
    }
    s + c
}

/// Pairwise summation with a naive base case (vectorizes).
pub fn pairwise_sum(x: &[f64]) -> f64 {
    if x.len() <= 64 { return x.iter().sum(); }
    let mid = x.len() / 2;
    pairwise_sum(&x[..mid]) + pairwise_sum(&x[mid..])
}

/// Compensated dot product (Ogita-Rump-Oishi Dot2): as accurate as computing in twice the precision.
pub fn dot2(x: &[f64], y: &[f64]) -> f64 {
    let (mut s, mut c) = (0.0, 0.0);
    for (&a, &b) in x.iter().zip(y) {
        let (p, ep) = two_prod(a, b);
        let (t, es) = two_sum(s, p);
        s = t;
        c += ep + es;
    }
    s + c
}

/// Map f64 to a monotone i64 so that adjacent floats differ by 1 (+0 and -0 both map to 0).
fn ordered(x: f64) -> i64 {
    let b = x.to_bits() as i64;
    if b < 0 { i64::MIN - b } else { b }
}

/// Number of representable f64 values between a and b (0 if equal). None for NaN.
pub fn ulp_distance(a: f64, b: f64) -> Option<u64> {
    if a.is_nan() || b.is_nan() { return None; }
    Some((ordered(a) as i128 - ordered(b) as i128).unsigned_abs() as u64)
}

/// Mixed tolerance: |a - b| <= atol + rtol * max(|a|, |b|).
pub fn approx_eq(a: f64, b: f64, rtol: f64, atol: f64) -> bool {
    a == b || (a - b).abs() <= atol + rtol * a.abs().max(b.abs())
}

/// Stable quadratic roots of a x^2 + b x + c (real roots only).
pub fn quadratic_roots(a: f64, b: f64, c: f64) -> Option<(f64, f64)> {
    let d = b.mul_add(b, -4.0 * a * c);
    if d < 0.0 { return None; }
    let q = -0.5 * (b + b.signum() * d.sqrt());
    Some((q / a, c / q))
}

Complexity and constants

• Bounds (Higham): naive recursive sum |error| ≤ (n − 1)·u·Σ|xᵢ| + O(u²); pairwise ≤ ⌈log₂n⌉·u·Σ|xᵢ|; Kahan ≤ (2u + O(n·u²))·Σ|xᵢ|.
• Measured with the code above, 10⁷ copies of 0.1 (exact binary sum known): naive error 1.6·10⁻⁴, pairwise 4.7·10⁻¹⁰, Kahan and Neumaier 0. On 10⁷ uniform [0, 1) values: naive relative error 6.5·10⁻¹⁵, pairwise 1.9·10⁻¹⁶, Kahan/Neumaier 0 (correctly rounded). On [1, 10¹⁰⁰, 1, −10¹⁰⁰]: naive and Kahan return 0, Neumaier returns the correct 2. Dot2 on (10¹⁶, 1, −10¹⁶)·(1, 1, 1) returns 1, naive 0.
• Kahan's loop has a dependency chain of 4 operations per element and cannot be reassociated (the compiler must not apply -ffast-math); measured on Apple M4 Pro: naive iter().sum() 0.55 ns per element (one dependent add chain, not vectorized since Rust keeps the order), pairwise with a 64-element base 0.13 ns (the base case vectorizes), Neumaier 0.72 ns, Kahan 2.2 ns.
• f32 accumulators: u ≈ 6·10⁻⁸; summing 2·10⁷ ones in f32 stops at 16,777,216 = 2²⁴ (measured), after which adding 1 changes nothing.

Efficiency

• Cost: measured on Apple M4 Pro with the code above: naive iter().sum() 0.55 ns per element, pairwise with a 64-element base 0.13 ns, Neumaier 0.72 ns, Kahan 2.2 ns. Accuracy on 10⁷ uniform [0, 1) values: naive relative error 6.5·10⁻¹⁵, pairwise 1.9·10⁻¹⁶, Kahan/Neumaier correctly rounded. Memory is one accumulator (two with compensation). Bound: the dependent add chain (latency), not bandwidth.
• Levers:
  • pairwise summation with a 64-element base case — error O(u·log n) at 0.13 ns per element, faster than naive (measured above) — any long reduction
  • Neumaier instead of Kahan — 3× faster and also correct when a term exceeds the sum (measured above) — compensated sums
  • f64 accumulator for f32 data — avoids stalling at 2²⁴ (measured above) — f32 and half-precision inputs
  • FMA for two_prod — 2 flops instead of a long split sequence (Ogita-Rump-Oishi 2005) — compensated dot products
  • fixed reduction tree or superaccumulator — bitwise reproducible results (Neal's xsum) — parallel and distributed sums
• Naive → standard → state of the art:
  • recursive loop in f32 · stalls at 16,777,216 when summing ones (measured above)
  • pairwise in f64 · 0.13 ns per element, 1.9·10⁻¹⁶ error
  • Neumaier or Sum2/Dot2, or an exact superaccumulator · 0.72 ns per element, correctly rounded

Implementation notes

• Never compile with fast-math (-ffast-math, -Ofast) around compensated code: it may simplify (t − s) − y to 0. Rust does not reassociate floats by default.
• Reproducibility: parallel reductions change the order and therefore the result; fix the reduction tree, or use compensated / exact (superaccumulator) summation for bitwise-reproducible results.
• Accumulate in higher precision: f32 data summed into an f64 accumulator is often the cheapest fix (mixed-precision training accumulates in f32 for f16/bf16 data).
• Comparisons: use a relative tolerance (rtol ≈ 10⁻⁹-10⁻¹² in f64) plus an absolute floor for values near zero; ULP distance for testing math libraries.
• Subnormals (below 2.2·10⁻³⁰⁸) take a slow microcode path on many x86 CPUs; flush-to-zero (FTZ/DAZ) trades accuracy for speed.
• Variance: use Welford's update (Monte Carlo), not E[x²] − E[x]², which cancels catastrophically when the mean is large.
• Integers in floats: f64 represents every integer up to 2⁵³ = 9,007,199,254,740,992 exactly.

Pitfalls

• Equality tests on computed floats (x == 0.3), or loop counters in floats (while t < 1.0 { t += 0.1 } runs 11 times, not 10).
• Subtracting nearly equal quantities computed separately (finite differences with tiny h, variance by the textbook formula, 1 − cos x).
• Summing small terms into a large running total (time accumulators in long simulations: t += dt drifts; compute t = n·dt instead).
• Assuming the same result across compilers/CPUs: x87 extended precision, FMA contraction and libm differences change last bits.
• NaN propagation: NaN != NaN, max(NaN, x) semantics vary; use total_cmp for sorting.

Variants and alternatives

In production

Python math.fsum (exact Shewchuk-style partial sums), NumPy sum (pairwise with a block size of 128), Julia sum (pairwise), Shewchuk's robust geometric predicates (adaptive exact arithmetic) used in Triangle and CGAL, the QD library, the CORE-MATH and CRlibm correctly rounded libms, ReproBLAS for reproducible reductions. Machine learning: mixed-precision training keeps f32 master weights and accumulates in f32 (Micikevicius et al.).

Implementations

Sources

• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, chapters 2 and 4. book
• David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys 23(1), 1991. paper
• William Kahan, Further Remarks on Reducing Truncation Errors, Communications of the ACM 8(1), 1965. paper
• Takeshi Ogita, Siegfried M. Rump, Shin'ichi Oishi, Accurate Sum and Dot Product, SIAM Journal on Scientific Computing 26(6), 2005. paper
• Jonathan R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete & Computational Geometry 18(3), 1997. paper
• IEEE Std 754-2019, IEEE Standard for Floating-Point Arithmetic.
• Paulius Micikevicius et al., Mixed Precision Training, ICLR 2018. paper