Kahan summation

A million purchases at 149.90 tenge is exactly 149,900,000 tenge. Add them up with an ordinary loop over float numbers and you get 149,900,000.0028. Almost nothing, but some extra tiyn have appeared from somewhere, and with every million there are more of them. In float32, which graphics cards and many models compute with, the same sum comes up more than a hundred and twenty thousand tenge short. Below, the sum grows right in your browser in three ways, and the chart shows how far each has drifted from the truth.

The blue line, float32, goes up fastest of all: it has only about seven significant digits. The orange one, float64, grows more slowly, but it grows too. The green one is the same float64, but with Kahan’s correction: its error stays put however many terms there are. Switch to the harmonic series: in float32 it stops growing at the 2,097,152nd term and freezes forever at 15.403683, although in mathematics it grows to infinity. Why 0.1 + 0.2 is not 0.3 and how floating-point numbers are built is told in the course chapters on variables and on bits; here is what happens when there are millions of such additions.

The idea

Each addition rounds the result to the nearest number that fits the format and loses the low bits. The loss is at most half the step between neighbouring numbers near the sum, but the step grows with the sum, and the losses pile up. Kahan proposed keeping a second variable beside the sum, what was lost in the previous addition, and adding it to the next term. What was lost does not vanish but comes back into the sum at the next step.

Here is the whole algorithm. The variable c is the correction, zero at first.

y = x - c          # the term together with what was lost last time
t = s + y          # the new sum: the low bits of y are lost here
c = (t - s) - y    # (t - s) is how much was really added; minus y is how much was lost
s = t

The third line looks like zero: (s + y − s) − y. In exact arithmetic it is, but in machine arithmetic t − s is what was really added to the sum, and its difference from y is exactly the lost bits, and this difference itself is computed without error. William Kahan published this technique in 1965. For a plain loop the error of a sum of n terms is proportional to n in the worst case; with Kahan’s correction it is bounded by two units in the last place times the sum of the absolute values of the terms, and hardly depends on n.

Where Kahan is powerless

The correction is one number, and it is good when small numbers are added to a large sum many times. But if a small number is squeezed between huge ones of different sign, 10¹⁰⁰, then 1, then −10¹⁰⁰, the correction itself is lost against the huge sum, and the result is zero instead of one. In 1974 Arnold Neumaier proposed a variant that compares the absolute values of the sum and the term and keeps the low bits of whichever is smaller; it copes with this example, and it is exactly what works inside the built-in sum() in Python since version 3.12.

There are two more ways. Pairwise summation adds the numbers as a tree: half with half, quarter with quarter. Each number goes through not n additions but log₂ n, and the error grows as a logarithm; this is how numpy.sum adds. Exact summation does not round at all until the end: math.fsum keeps the sum as several float numbers that together represent it without loss (Jonathan Shewchuk’s algorithm), and in Go the same can be done with a long big.Float.

And money is not counted with fractional numbers at all. 149.90 tenge is 14,990 tiyn, an integer, and a sum of integers is exact as long as it fits in 64 bits, that is, up to 92 quadrillion tenge. Databases keep money in the DECIMAL type, and Python has the decimal module.

Code

Five methods on two experiments, and the harmonic series in float32:

Kahan’s and Neumaier’s corrections give the exact answer in these experiments. Pairwise summation of 149.90 was off by one unit in the last place, much better than the plain loop, but not exact. The exact sum is the one printed in the decimal line, with one reservation: in the first experiment the term is not one tenth but the nearest binary number to it, slightly larger, and the exact sum of a million such numbers rounds to 100,000 only because the difference is less than half a step.

What it costs

Kahan and Neumaier are O(n) and O(1) memory, but four arithmetic operations per term instead of one. In practice the time difference is often smaller: a long loop is limited by reading memory, not by additions. Pairwise is the same O(n) additions plus recursion, which implementations unroll into blocks. Exact is O(n) on average, but with long numbers in Go noticeably slower.

When to take what. The sum of a hundred numbers: a plain loop. Millions of terms in a simulation, accumulating time or coordinates step by step: Kahan’s or Neumaier’s correction. A total that must match to the last bit: the exact sum. Money: integers in the smallest units.

Problems

The problem is to compute both the sum an ordinary loop gives and the exact one. Among the tests are the daily rainfall in Almaty over 85 years and numbers that cancel each other.