Kahan summation algorithm

In numerical analysis, the Kahan summation algorithm, also known as compensated summation,^[1] significantly reduces the numerical error in the total obtained by adding a sequence of finite-precision floating-point numbers, compared to the naive approach. This is done by keeping a separate running compensation (a variable to accumulate small errors), in effect extending the precision of the sum by the precision of the compensation variable.

In particular, simply summing ${\displaystyle n}$ numbers in sequence has a worst-case error that grows proportional to ${\displaystyle n}$, and a root mean square error that grows as ${\displaystyle {\sqrt {n}}}$ for random inputs (the roundoff errors form a random walk).^[2] With compensated summation, using a compensation variable with sufficiently high precision the worst-case error bound is effectively independent of ${\displaystyle n}$, so a large number of values can be summed with an error that only depends on the floating-point precision of the result.^[2]

The algorithm is attributed to William Kahan;^[3] Ivo Babuška seems to have come up with a similar algorithm independently (hence Kahan–Babuška summation).^[4] Similar, earlier techniques are, for example, Bresenham's line algorithm, keeping track of the accumulated error in integer operations (although first documented around the same time^[5]) and the delta-sigma modulation.^[6]

In pseudocode, the algorithm will be:

function KahanSum(input)
    // Prepare the accumulator.
    var sum = 0.0
    // A running compensation for lost low-order bits.
    var c = 0.0
    // The array input has elements indexed input[1] to input[input.length].
    for i = 1 to input.length do
        // c is zero the first time around.
        var y = input[i] - c
        // Alas, sum is big, y small, so low-order digits of y are lost.         
        var t = sum + y
        // (t - sum) cancels the high-order part of y;
        // subtracting y recovers negative (low part of y)
        c = (t - sum) - y
        // Algebraically, c should always be zero. Beware
        // overly-aggressive optimizing compilers!
        sum = t
    // Next time around, the lost low part will be added to y in a fresh attempt.
    next i

    return sum

This algorithm can also be rewritten to use the Fast2Sum algorithm:^[7]

function KahanSum2(input)
    // Prepare the accumulator.
    var sum = 0.0
    // A running compensation for lost low-order bits.
    var c = 0.0
    // The array input has elements indexed 
    for i = 1 to input.length do
        // c is zero the first time around.
        var y = input[i] + c
        // sum + c is an approximation to the exact sum.
        (sum,c) = Fast2Sum(sum,y)
    // Next time around, the lost low part will be added to y in a fresh attempt.
    next i
    return sum

The algorithm does not mandate any specific choice of radix, only for the arithmetic to "normalize floating-point sums before rounding or truncating".^[3] Computers typically use binary arithmetic, but to make the example easier to read, it will be given in decimal. Suppose we are using six-digit decimal floating-point arithmetic, sum has attained the value 10000.0, and the next two values of input[i] are 3.14159 and 2.71828. The exact result is 10005.85987, which rounds to 10005.9. With a plain summation, each incoming value would be aligned with sum, and many low-order digits would be lost (by truncation or rounding). The first result, after rounding, would be 10003.1. The second result would be 10005.81828 before rounding and 10005.8 after rounding. This is not correct.

However, with compensated summation, we get the correctly rounded result of 10005.9.

Assume that c has the initial value zero. Trailing zeros shown where they are significant for the six-digit floating-point number.

  y = 3.14159 - 0.00000             y = input[i] - c
  t = 10000.0 + 3.14159             t = sum + y
    = 10003.14159                   Normalization done, next round off to six digits.
    = 10003.1                       Few digits from input[i] met those of sum. Many digits have been lost!
  c = (10003.1 - 10000.0) - 3.14159 c = (t - sum) - y  (Note: Parenthesis must be evaluated first!)
    = 3.10000 - 3.14159             The assimilated part of y minus the original full y.
    = -0.0415900                    Because c is close to zero, normalization retains many digits after the floating point. 
sum = 10003.1                       sum = t

The sum is so large that only the high-order digits of the input numbers are being accumulated. But on the next step, c, an approximation of the running error, counteracts the problem.

  y = 2.71828 - (-0.0415900)        Most digits meet, since c is of a size similar to y.
    = 2.75987                       The shortfall (low-order digits lost) of previous iteration successfully reinstated.
  t = 10003.1 + 2.75987             But still only few meet the digits of sum.
    = 10005.85987                   Normalization done, next round to six digits.
    = 10005.9                       Again, many digits have been lost, but c helped nudge the round-off.
  c = (10005.9 - 10003.1) - 2.75987 Estimate the accumulated error, based on the adjusted y.
    = 2.80000 - 2.75987             As expected, the low-order parts can be retained in c with no or minor round-off effects.
    = 0.0401300                     In this iteration, t was a bit too high, the excess will be subtracted off in next iteration.
sum = 10005.9                       Exact result is 10005.85987, sum is correct, rounded to 6 digits.

The algorithm performs summation with two accumulators: sum holds the sum, and c accumulates the parts not assimilated into sum, to nudge the low-order part of sum the next time around. Thus the summation proceeds with "guard digits" in c, which is better than not having any, but is not as good as performing the calculations with double the precision of the input. However, simply increasing the precision of the calculations is not practical in general; if input is already in double precision, few systems supply quadruple precision, and if they did, input could then be in quadruple precision.