Compensated Summation and How Summation Order Bounds Error

Summing an array looks like the one operation that cannot go wrong. It has no divisions, no subtractions of similar quantities, no transcendental functions. Yet the error of a naive sum grows linearly with the number of terms in the worst case, and for a million-element array that growth is large enough to matter. The three standard responses — reorder the additions, compensate for the discarded bits, or use an exact accumulator — sit at different points on the same trade-off between accuracy and the ability to use the machine’s full width.

The error of a left-to-right sum

Write the sum as s₁ = x₁, sₖ = sₖ₋₁ + xₖ, and let ŝₖ be the computed value. Each addition introduces a relative error of at most u, so the computed result satisfies

$$ |\hat{s}n - s_n| \le \gamma{n-1} \sum_{i=1}^{n} |x_i|, \qquad \gamma_k = \frac{k u}{1 - k u} $$

For small k, γ_k ≈ ku, so the bound is essentially (n−1)u Σ|xᵢ|. The bound is in terms of the sum of absolute values, not the sum itself. That distinction is the whole content of the result: if the terms all have the same sign, Σ|xᵢ| = |Σxᵢ| and the bound is a relative error of (n−1)u. If the terms cancel, the absolute error stays the same while the result shrinks, and the relative error grows without limit.

The growth with n is the part that reordering can fix. For binary64, u = 2^-53, so a million-term sum has a worst-case bound of about 1.1e-10 relative to Σ|xᵢ|, which would leave roughly ten of the sixteen digits correct on the unluckiest rounding path. The bound is worst case and is not usually attained; rounding errors of the same sign all accumulating is the unlucky path, not the typical one. But it is attainable, and a bound that can be reached is what a correctness argument needs.

Reordering: the logarithm in pairwise reduction

A left-to-right sum builds a chain of n − 1 additions, and error from the first addition is carried through every later one, so it is amplified n times. A balanced binary tree of additions has depth ⌈log₂ n⌉, and every partial sum is involved in at most that many roundings. The bound becomes γ_{⌈log₂ n⌉} Σ|xᵢ|, which for a million terms is about 20u instead of 10⁶u. The same depth argument governs how the error behaves once a compiler reshapes a reduction into lanes or a runtime splits it across threads, as described in Reassociation, Excess Precision, and the Limits of Bitwise Reproducibility.

The cost is memory access pattern rather than arithmetic. A recursive halving touches the array in a cache-unfriendly order at the upper levels; the practical compromise is blocked pairwise summation, which sums fixed-size blocks sequentially and then combines the block results as a tree. With a block size chosen to fit in cache, the depth is log₂(n/block) plus the intra-block chain, so the bound keeps the logarithmic character while the inner loop stays vectorisable.

The same idea appears in parallel reductions without being called pairwise summation. Each thread or lane accumulates a partial sum, and the partials are combined in a tree. The depth of the tree is the number of roundings a given input passes through, so a parallel reduction is usually more accurate than a serial one over the same data, not less.

Reordering is not free of a different cost: it changes the result. Two programs that sum the same array in different orders produce different bits, and if the order depends on thread scheduling, the difference is nondeterministic. That is why a library may document that its reduction result is not bit-reproducible across thread counts.

Kahan’s compensation

The alternative to reordering is to keep the bits that each addition discards. Kahan’s algorithm maintains a running correction c alongside the sum, and the update is arranged so that the low-order part of each addition is captured rather than lost.

/* Kahan compensated summation. The loop-carried dependency is the point:
   every iteration must see the previous correction. */
double kahan_sum(const double *a, size_t n)
{
    double s = 0.0, c = 0.0;
    for (size_t i = 0; i < n; i++) {
        double y = a[i] - c;      /* apply the correction from the last step */
        double t = s + y;         /* rounded sum */
        c = (t - s) - y;          /* exact residual, negated */
        s = t;
    }
    return s;
}

The line c = (t - s) - y is the mechanism. The quantity t − s is the amount actually added to s, computed exactly, because t and s are close and the subtraction is exact by Sterbenz. Subtracting y from that gives the part of y that was rounded away. The next iteration subtracts it from the incoming term, so the discarded amount is reinjected rather than dropped.

The error bound reflects the difference. Kahan’s sum satisfies |ŝ − s| ≤ (2u + O(nu²)) Σ|xᵢ|, with no first-order dependence on n. The linear growth is gone; what remains is a small constant times the absolute sum, plus a second-order term that only matters for very long sums.

The algorithm has a blind spot. The step c = (t − s) − y assumes that s is the larger operand, so that t − s is the small quantity that captures the rounding. When the incoming term is larger than the running sum — which happens at the start of a sum, or whenever a large value arrives late — the reasoning inverts and the correction is wrong. The failure is not a loss of the bound’s constant; it is a case where the compensation does not track the actual rounding at all.

Neumaier’s branch

Neumaier’s variant repairs the blind spot by choosing the correction according to which operand is larger.

double neumaier_sum(const double *a, size_t n)
{
    double s = 0.0, c = 0.0;
    for (size_t i = 0; i < n; i++) {
        double t = s + a[i];
        if (fabs(s) >= fabs(a[i]))
            c += (s - t) + a[i];   /* the term was the small one */
        else
            c += (a[i] - t) + s;   /* the sum was the small one */
        s = t;
    }
    return s + c;
}

Both branches compute the same thing in exact arithmetic — the residual of the addition — but they compute it with the operands ordered so that the subtraction is the exact one. The branch is predictable for monotone data and unpredictable for mixed magnitudes, and its presence blocks the vectorisation that a plain sum allows. A branch-free formulation exists using fmin/fmax to select, at the cost of more arithmetic per element.

The final s + c is a single rounding of the accumulated correction, which is where the last bit of accuracy is recovered. Omitting it, as some implementations do, leaves the correction unapplied and reduces the algorithm to a slightly better naive sum.

The cost, and when it is not worth paying

Compensation multiplies the arithmetic per element by roughly four and, more importantly, creates a loop-carried dependency: each iteration’s correction depends on the previous iteration’s sum. A plain sum can be split into independent lanes and vectorised; a compensated sum cannot, without restructuring it into per-lane compensated sums that are combined at the end. That restructuring is usually the right answer on a vector machine, and it is what a library implementing a compensated reduction will do internally.

The trade is accuracy against throughput, and the direction depends on the condition number. For data with a small condition number Σ|xᵢ|/|Σxᵢ|, the naive bound nu is already small and compensation buys nothing measurable; the extra work is pure cost. For data whose sum is a small residual of large cancelling terms, the naive result may have no correct digits at all and compensation is the difference between an answer and noise. Since the condition number can be computed alongside the sum at the cost of accumulating Σ|xᵢ|, a library can choose the algorithm adaptively.

There is a third option that removes the question. An exact accumulator represents the running sum as a short list of non-overlapping floating-point numbers, so no bits are discarded and the final result is the correctly rounded exact sum. The techniques come from the error-free transformations used in Error-Free Transformations and Double-Double Arithmetic, and they cost more per element than Kahan’s method while giving a guarantee that no rounding-based scheme can match. For sums where the exact result must be reproduced bit for bit, that guarantee is the requirement, and the cost is the price of determinism rather than of accuracy.

What the bounds do not cover

Every bound in this article assumes the model fl(x + y) = (x + y)(1 + δ) with |δ| ≤ u, which fails in three places. Results in the subnormal range have an absolute rather than relative error, so the Σ|xᵢ| bound has to be augmented with a term proportional to n times the smallest subnormal spacing. An overflow in the intermediate partial sum makes the bound vacuous even when the final result would be representable. And a compiler that reassociates the loop, or that contracts the additions into wider operations, has computed a different expression, so the bound applies to a program the compiler did not emit.

The last point is worth stating plainly: Kahan’s algorithm does not survive optimisation flags that permit reassociation. The expression (t - s) - y is an algebraic zero in exact arithmetic, and a compiler entitled to simplify it will delete the entire compensation.