UT_Accumulator.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * COMMENTS:
 *    Accumulators for scalar and vector types. These are used to collect the
 *    sum of incrementally accumulated quantities. Using the Neumaier's
 *    extension to Kahan's summation method, these should vastly limit the
 *    truncation error of floating-points.
 */

#ifndef __UT_Accumulator__
#define __UT_Accumulator__

#include "UT_FixedVector.h"
#include <SYS/SYS_Types.h>


namespace UT
{

template <typename T>
class Accumulator
{
public:
    explicit             Accumulator(const T &init_val = T(0.0)) :
                                mySum(init_val), myCompensation(T(0.0)) { }

    SYS_FORCE_INLINE
    T                   &operator+=(const T& rhs)
    {

        T  t = mySum + rhs;
        if (SYSabs(mySum) >= SYSabs(rhs))
            myCompensation += (mySum - t) + rhs;
        else
            myCompensation += (rhs - t) + mySum;
        mySum = t;
        return mySum;
    }

    T                    sum() const { return mySum + myCompensation; }

private:
    T                    mySum;
    T                    myCompensation;
};

template <typename T, int S>
class Accumulator<UT_FixedVector<T, S>>
{
public:
    using V = UT_FixedVector<T, S>;
    
    explicit             Accumulator(const V &init_val = V(T(0.0))) :
                                mySum(init_val), myCompensation(T(0.0)) { }

    SYS_FORCE_INLINE
    V                   &operator+=(const V& rhs)
    {
        for (int i = 0, ie = S; i < ie; i++)
        {
            T t = mySum(i) + rhs(i);
            if (SYSabs(mySum(i)) >= SYSabs(rhs(i)))
                myCompensation(i) += (mySum(i) - t) + rhs(i);
            else
                myCompensation(i) += (rhs(i) - t) + mySum(i);
            mySum(i) = t;
        }
        return mySum;
    }

    V                    sum() const { return mySum + myCompensation; }
        
private:
    V                    mySum;
    V                    myCompensation;
};


template <typename T>
class Accumulator<std::complex<T>>
{
public:
    using C = std::complex<T>;
    
    explicit             Accumulator(const C &init_val = C(0.0)) :
                                mySum(init_val), myCompensation(C(0.0)) { }

    SYS_FORCE_INLINE
    C                   &operator+=(const C& rhs)
    {
        T t = mySum.real() + rhs.real();
        if (SYSabs(mySum.real()) >= SYSabs(rhs.real()))
            myCompensation.real(myCompensation.real()
                                + (mySum.real() - t) + rhs.real());
        else
            myCompensation.real(myCompensation.real()
                                    + (rhs.real() - t) + mySum.real());
        mySum.real(t);
    
        t = mySum.imag() + rhs.imag();
        if (SYSabs(mySum.imag()) >= SYSabs(rhs.imag()))
            myCompensation.imag(myCompensation.imag()
                                    + (mySum.imag() - t) + rhs.imag());
        else
            myCompensation.imag(myCompensation.imag()
                                    + (rhs.imag() - t) + mySum.imag());
        mySum.imag(t);
        
        return mySum;
    }

    C                    sum() const { return mySum + myCompensation; }
        
private:
    C                    mySum;
    C                    myCompensation;
};


} // namespace UT

template <typename T>
using UT_Accumulator = UT::Accumulator<T>;


#endif