Kahan summation

Inspired by another question, I decided to implement Kahan summation in C++ (though that question implemented a different summation algorithm).

Since I was writing C++, I decided to make the code generic. Although it's a little difficult to imagine anybody bothering to use Kahan summation on single-precision floating point, I suppose it's possible--and while doubles are probably the most common type, using it on various container types is probably more common.

Anyway, I've included a quick test that attempts to show how much difference an accurate summation can make. It starts with a relatively large number (1e4), then adds a much smaller number (1e-15) to it many (1e7) times. Then it subtracts the initial starting value from that result, and multiplies what's left by 1e19.

With naive summation, the difference in magnitude prevents any of the additions from changing the result (at all). So, when we subtract the initial value, we get 0. Multiplying by 1e19 leaves that as 0. Along with the Kahan summation, I've provided a reference: instead of adding the small number to the large one many times, it multiplies the smaller number by the count, and adds that to the larger number.

At least in my testing, the version using Kahan summation matches the reference to twenty digits of precision, while the version using naive summation doesn't produce even a single digit correctly.

#include <iostream>
#include <iterator>
#include <iomanip>
#include <vector>
#include <numeric>

template <class InIt>
typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
typedef typename std::iterator_traits<InIt>::value_type real;
real sum = real();
real running_error = real();

for ( ; begin != end; ++begin) {
real difference = *begin - running_error;
real temp = sum + difference;
running_error = (temp - sum) - difference;
sum = temp;
}
return sum;
}

int main(){
const double initial = 1e3;
const double addend = 1e-15;
const double count = 1e7;
const double multiplier = 1e19;

std::vector<double> d;

d.push_back(initial);

double result = (std::accumulate(d.begin(), d.end(), 0.0) - initial) * multiplier;

double result2 = (accumulate(d.begin(), d.end()) - initial) * multiplier;

double reference = ((initial + count * addend) - initial) * multiplier;

std::cout << "   simple: " << std::setprecision(20) << result << "\n";
std::cout << "    Kahan: " << std::setprecision(20) << result2 << "\n";
std::cout << "Reference: " << std::setprecision(20) << reference << "\n";
}


While all the code is (of course) open to critique, I'm obviously much more interested in comments on the implementation of the summation algorithm than the accompanying test code.

• Your implementation looks sound, but there is no guarantee it is portable. If you are only targeting one compiler, there may be command-line switches available to ensure the optimiser doesn't make over-zealous assumptions about algebraic equivalence. There's a good article about it on MSDN. AFAIK, the only truly portable way would be to put the function in a separate .cpp file, and compile it separately, without optimisations (assuming you can set up a suitable build system). Tricky with templates. – linguamachina Jul 9 '14 at 9:27
• Consider returning a final value of running_error along with sum. In a context of numerical analysis the result without an error margin is more or less useless. – vnp Jul 9 '14 at 17:44
• Why are you assuming the type of the sum is the same as that of the element? I mean, that could be a default but I would make that a template argument. – einpoklum Jan 2 '17 at 20:11
• Should you be using a volatile somewhere? I am guessing that some optimizing compiler might optimize running_error and make it 0. (Disclaimer: this just a guess and not backed by any experiments/literature review from my side) – lakshayg Jun 13 '17 at 13:54

I've taken the liberty of freely using C++11 in my answer because about half of the improvements I propose require it.

An instrumentation class

I decided to do a little instrumentation to see how this template would do with an artificial type so I created my own Goofy math type.

class Goofy
{
private:
double num;
public:
Goofy(double n = 0) : num(n) { ++constructions; }
Goofy(const Goofy &g2) : num(g2.num) { ++copyconstructions; }
Goofy(const Goofy &&g2) : num(g2.num) { ++moves; }
~Goofy() { ++destructions; }
Goofy &operator=(const Goofy &g2) { num = g2.num; return *this; }
Goofy &operator-=(const Goofy &g2) { num -= g2.num; return *this; }
Goofy &operator+=(const Goofy &g2) { num += g2.num; return *this; }
// none of the code below is needed by the new version of the function
Goofy &operator*=(const Goofy &g2) { num *= g2.num; return *this; }
friend std::ostream &operator<<(std::ostream &out, const Goofy &g2) {
return out << g2.num;
}
static void report(int line) {
std::cout << "At line " << line
<< "\nconstructions = " << Goofy::constructions
<< "\n       copies = " << Goofy::copyconstructions
<< "\n        moves = " << Goofy::moves
<< "\n destructions = " << Goofy::destructions
<< "\n     existing = " << Goofy::constructions +
Goofy::copyconstructions + Goofy::moves -
Goofy::destructions
<< '\n';
}
static long constructions;
static long copyconstructions;
static long moves;
static long destructions;
};

long Goofy::constructions = 0;
long Goofy::copyconstructions = 0;
long Goofy::moves = 0;
long Goofy::destructions = 0;


This class doesn't have quite everything necessary for the original code, however because it lacks three binary operators:

Goofy operator-(const Goofy &g1, const Goofy &g2)
{
Goofy result(g1);
result -= g2;
return result;
}

Goofy operator+(const Goofy &g1, const Goofy &g2)
{
Goofy result(g1);
result += g2;
return result;
}

Goofy operator*(const Goofy &g1, const Goofy &g2)
{
Goofy result(g1);
result *= g2;
return result;
}


Modified test harness

For convenience, I then modified your test code a bit:

typedef float counttype;
typedef Goofy mytype;

void test(){
const mytype initial{1e3};
const counttype count{1e7};
const mytype multiplier{1e19};

std::vector<mytype> d;

d.push_back(initial);

mytype result = (std::accumulate(d.begin(), d.end(), mytype()) - initial) * multiplier;

Goofy::report(__LINE__);
mytype result2 = (accumulate(d.cbegin(), d.cend()) - initial) * multiplier;
Goofy::report(__LINE__);
mytype result3 = (accumulateOriginal(d.cbegin(), d.cend()) - initial) * multiplier;
Goofy::report(__LINE__);

mytype reference = ((initial + count * addend) - initial) * multiplier;

std::cout << "   simple: " << std::setprecision(20) << result << "\n";
std::cout << "    Kahan: " << std::setprecision(20) << result2 << "\n";
std::cout << "origKahan: " << std::setprecision(20) << result3 << "\n";
std::cout << "Reference: " << std::setprecision(20) << reference << "\n";
}

int main(){
test();
Goofy::report(__LINE__);
}


As you can see, there are now two versions: accumulateOriginal is the code as posted, and accumulate is one I modified. Here's the modified version and below that is an explanation of what was done and why. I also chose to use const iterators just to verify that the original vector wasn't being modified. Unfortunately, I don't know of any standard way to indicate that in the templated function's code.

A modified template

template <class InIt>
typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
typedef typename std::iterator_traits<InIt>::value_type real;
real sum = real();
real running_error = real();
real temp;
real difference;

for (; begin != end; ++begin) {
difference = *begin;
difference -= running_error;
temp = sum;
temp += difference;
running_error = temp;
running_error -= sum;
running_error -= difference;
sum = std::move(temp);
}
return sum;
}


The first change was that real sum = real(); is exactly the same as real sum; so I chose the shorter form.(Edit: As @ruds points out in a comment, this isn't necessarily true for primitive types such as int or double.) I've also moved the declarations of temp and difference out of the loop. In this case, that saves 40000000 constructor calls and 40000000 destructor calls. In that same vein, I've used std::move to give the hint to the compiler that the value of temp doesn't need to be preserved. This can make a difference if the construction and/or destruction of these objects is expensive and hurts nothing if it's not, so I think it's an improvement. For example, on my machine, using std::complex<double> as the numeric type, it takes 479 millseconds for the modified version and 802 milliseconds for the original.

Additional changes are to use -= and += to avoid the requirement for a non-class member + and - binary operators.

It's worth noting that *= and the binary operator * are only required for the test code and not for the template itself. Since the binary operators are often defined by using the member operators, this may represent a slightly smaller requirement for the iterators used with this templated function. The other other part that's only needed by the test code is, of course, the operator<< code.

Finally, it's possible that there may be some use in providing a test to make sure that values initialize to zero. One way might be to use something like this just after the line defining real:

static_assert(Default_is_zero<real>(), "real value must initialize to 0");


The code for the assert is this:

template <typename T>
constexpr bool Default_is_zero() {
return T{} == T{0};
}


However, this adds two requirements not explictly needed otherwise, namely the ability to initialize a real type with an integer and the need for an operator==. Hence this works for std::complex<double> but fails with Goofy. Both are probably useful for a numeric type, but are beyond the bare minimum.

• real sum; is the same as real sum = real(); if real is Goofy, but not when real is double or float, which is the most common case. – ruds Jul 10 '14 at 8:55
• Are you sure you're not over-optimizing? I mean, do you really want to be planning for summation of classes with non-trivial constructors? I feel that's a bit overkill. – einpoklum Jan 2 '17 at 20:16