“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

Stable Summation Algorithms

Numerical Methods

When you add numbers in floating point, they don't always add up to what you expect. A double carries about 16 significant digits, so when a small number is added to a much larger sum, its low digits — sometimes all of its digits — are rounded away. Add a million small numbers to a big one naively and the losses compound. This page builds up three summers: the naive loop, Kahan's fix, and Neumaier's fix of Kahan's blind spot.

The Problem

The brutal demonstration: sum . Step one computes , and since 1 is far below the last representable digit of , the answer is exactly — the 1 is absorbed without a trace. Step two subtracts and leaves 0. The exact answer is 1; every digit of it was lost. (Run on this machine: naive gives 0.0 — and so does Kahan. Only Neumaier returns 1.0. Keep reading for why.)

Kahan Summation

Kahan's idea: after every addition, measure the rounding error you just made and feed it back into the next addition. The compensation variable c holds exactly the piece that was rounded away — computable in floating point itself, because (t − sum) − y recovers what t = sum + y lost. On the everyday case (a big sum absorbing many small addends) this works beautifully: summing 1 plus a million copies of , the naive loop returns exactly 1.0 — every addend absorbed — while Kahan returns the correct 1.0000000001.

Neumaier Summation

Kahan's blind spot: it assumes the loss happens on the addend's side. When an incoming value is larger than the running sum — as when arrives after the 1 — the digits destroyed belong to the sum, and Kahan's correction measures the wrong quantity (that is why it also returns 0 on the demonstration above). Neumaier's fix is one comparison: check which of the two operands is smaller in magnitude, and compensate for that one:

The final result is .

C++ Implementation

The following code implements all three methods:

#include <vector>
#include <cmath>
#include <iostream>

enum class SumType {
    Naive,
    Kahan,
    Neumaier
};

class Summation {
public:
    explicit Summation(SumType t)
        : type(t), sum(0.0), c(0.0) {}

    void add(double x) {
        switch (type) {

        case SumType::Naive:
            sum += x;
            break;

        case SumType::Kahan:
        {
            double y = x - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
        break;

        case SumType::Neumaier:
        {
            double t = sum + x;
            if (std::fabs(sum) >= std::fabs(x))
                c += (sum - t) + x;   // sum is bigger
            else
                c += (x - t) + sum;   // x is bigger
            sum = t;
        }
        break;
        }
    }

    double value() const {
        return sum + c;   // for naive: c = 0
    }

    // convenient overload to sum a vector
    double operator()(const std::vector<double>& v) {
        for (double x : v) add(x);
        return value();
    }

    void reset() {
        sum = 0.0;
        c = 0.0;
    }

private:
    SumType type;
    double sum;   // main sum
    double c;     // correction
};

// Example usage
int main() {
    std::vector<double> v = {1e100, 1.0, -1e100};

    Summation naive(SumType::Naive);
    Summation kahan(SumType::Kahan);
    Summation neumaier(SumType::Neumaier);

    std::cout << "Naive:    " << naive(v) << std::endl;
    std::cout << "Kahan:    " << kahan(v) << std::endl;
    std::cout << "Neumaier: " << neumaier(v) << std::endl;
    
    return 0;
}

When to Use