HiveBrain v1.2.0
Get Started
← Back to all entries
patterncppMinor

Kahan summation

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
codereviewcppstackoverflowc++floating-point
summationkahanstackoverflow

Problem

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
#include
#include
#include
#include

template
typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) {
typedef typename std::iterator_traits::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 doub

Solution

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 mytype addend{1e-15};
    const counttype count{1e7};
    const mytype multiplier{1e19};

    std::vector d;

    d.push_back(initial);

    std::fill_n(std::back_inserter(d), count, addend);

    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 
typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) {
    typedef typename std::iterator_traits::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 destructio

Code Snippets

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;
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;
}
typedef float counttype;
typedef Goofy mytype;


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

    std::vector<mytype> d;

    d.push_back(initial);

    std::fill_n(std::back_inserter(d), count, addend);

    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__);
}
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;
}
static_assert(Default_is_zero<real>(), "real value must initialize to 0");

Context

StackExchange Code Review Q#56532, answer score: 6

Revisions (0)

No revisions yet.