[boost] accurate sum accumulator (kahan)

Subject: [boost] accurate sum accumulator (kahan)
From: Gaetano Mendola (mendola_at_[hidden])
Date: 2010-07-25 11:23:54

Hi all,
I had to implement in terms of boost accumulators the Kahan
summation algorithm. I will paste at the end of this message
the accumulator, the feature and the extractor source code.

This accumulator reduces the numerical error obtained with
standard sequential sum, as example the sum of 1e6f elements
equal to 1e-6f should give 1.0f, this is the results obtained
with sum and sum_kahan:

tag::sum:       1.009038925
tag::sum_kahan: 1.000000000

I'm not sure if that sum_kahan is worth to be inserted in the
accumulators provided with boost, I'm not an expert of
boost:accumulators (I just followed the guide to make this new
accumulator) so not sure if you prefer to have a fast_sum
and an accurate_sum instead of sum and sum_kahan.

Regards
Gaetano Mendola

// File kahan_sum.hpp
#ifndef KAHAN_SUM
#define KAHAN_SUM

#include <boost/accumulators/framework/accumulator_base.hpp>
#include <boost/accumulators/framework/parameters/sample.hpp>
#include <boost/numeric/conversion/cast.hpp>

namespace boost {
namespace accumulators {
namespace impl {

template<typename Sample>
struct sum_kahan_accumulator : accumulator_base
{
     typedef Sample result_type;

     template<typename Args>
     sum_kahan_accumulator(Args const & args)
       : sum(args[sample | Sample()]),
         compensation(boost::numeric_cast<Sample>(0.0))
     {
     }

     template<typename Args>
     void operator ()(Args const & args)
     {
       const Sample myTemp1 = args[sample] - this->compensation;
       const Sample myTemp2 = this->sum + myTemp1;
       this->compensation = (myTemp2 - this->sum) - myTemp1;
       this->sum = myTemp2;
     }

     result_type result(dont_care) const
     {
       return this->sum;
     }
private:
     Sample sum;
     Sample compensation;
};

}}}

namespace boost { namespace accumulators { namespace tag {

struct sum_kahan : depends_on< >
{
     typedef impl::sum_kahan_accumulator< mpl::_1 > impl;
};

}}}

namespace boost {
namespace accumulators {
namespace extract {

extractor< tag::sum_kahan > const sum_kahan = {};

}
using extract::sum_kahan;
}}

#endif