Stable Summation Algorithms

Numerical Methods

When summing floating-point numbers, naive summation can suffer from catastrophic cancellation. This page demonstrates three summation algorithms: naive, Kahan, and Neumaier, which provide increasing levels of numerical stability.

The Problem

Consider summing the vector . Naive summation can lose precision due to floating-point arithmetic limitations. The exact sum is , but naive summation may give .

Kahan Summation

Kahan's algorithm maintains a running compensation term to track lost low-order bits:

Neumaier Summation

Neumaier's algorithm is an improvement that handles cases where the sum itself is large:

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