/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; fill-column: 100 -*- */ /* * This file is part of the Collabora Office project. * * This Source Code Form is subject to the terms of the Mozilla Public * License, v. 2.0. If a copy of the MPL was not distributed with this * file, You can obtain one at http://mozilla.org/MPL/2.0/. */ #pragma once #include "arraysum.hxx" #include #include #include /** * This class provides LO with Kahan summation algorithm * About this algorithm: https://en.wikipedia.org/wiki/Kahan_summation_algorithm * For general purpose software we assume first order error is enough. * * Additionally queue and remember the last recent non-zero value and add it * similar to approxAdd() when obtaining the final result to further eliminate * accuracy errors. (e.g. for the dreaded 0.1 + 0.2 - 0.3 != 0.0) */ class KahanSum { public: constexpr KahanSum() = default; constexpr KahanSum(double x_0) : m_fSum(x_0) { } constexpr KahanSum(double x_0, double err_0) : m_fSum(x_0) , m_fError(err_0) { } constexpr KahanSum(const KahanSum& fSum) = default; public: /** * Performs one step of the Neumaier sum of doubles. * Overwrites the summand and error. * T could be double or long double. */ template static inline void sumNeumaierNormal(T& sum, T& err, const double& value) { T t = sum + value; if (std::abs(sum) >= std::abs(value)) err += (sum - t) + value; else err += (value - t) + sum; sum = t; } /** * Adds a value to the sum using Kahan summation. * @param x_i */ void add(double x_i) { if (x_i == 0.0) return; if (!m_fMem) { m_fMem = x_i; return; } sumNeumaierNormal(m_fSum, m_fError, m_fMem); m_fMem = x_i; } /** * Adds a value to the sum using Kahan summation. * @param fSum */ inline void add(const KahanSum& fSum) { #ifdef LO_X86_SIMD_AVAILABLE add(fSum.m_fSum + fSum.m_fError); add(fSum.m_fMem); #else // Without SSE2 the sum+compensation value fails badly. Continue // keeping the old though slightly off (see tdf#156985) explicit // addition of the compensation value. double sum = fSum.m_fSum; double err = fSum.m_fError; if (fSum.m_fMem != 0.0) sumNeumaierNormal(sum, err, fSum.m_fMem); add(sum); add(err); #endif } /** * Substracts a value to the sum using Kahan summation. * @param fSum */ inline void subtract(const KahanSum& fSum) { add(-fSum); } public: constexpr KahanSum operator-() const { KahanSum fKahanSum; fKahanSum.m_fSum = -m_fSum; fKahanSum.m_fError = -m_fError; fKahanSum.m_fMem = -m_fMem; return fKahanSum; } constexpr KahanSum& operator=(double fSum) { m_fSum = fSum; m_fError = 0; m_fMem = 0; return *this; } constexpr KahanSum& operator=(const KahanSum& fSum) = default; inline void operator+=(const KahanSum& fSum) { add(fSum); } inline void operator+=(double fSum) { add(fSum); } inline void operator-=(const KahanSum& fSum) { subtract(fSum); } inline void operator-=(double fSum) { add(-fSum); } inline KahanSum operator+(double fSum) const { KahanSum fNSum(*this); fNSum.add(fSum); return fNSum; } inline KahanSum operator+(const KahanSum& fSum) const { KahanSum fNSum(*this); fNSum += fSum; return fNSum; } inline KahanSum operator-(double fSum) const { KahanSum fNSum(*this); fNSum.add(-fSum); return fNSum; } inline KahanSum operator-(const KahanSum& fSum) const { KahanSum fNSum(*this); fNSum -= fSum; return fNSum; } /** * In some parts of the code of interpr_.cxx this may be used for * product instead of sum. This operator shall be used for that task. */ constexpr void operator*=(double fTimes) { if (m_fMem) { m_fSum = get() * fTimes; m_fMem = 0.0; } else { m_fSum = (m_fSum + m_fError) * fTimes; } m_fError = 0.0; } constexpr void operator/=(double fDivides) { if (m_fMem) { m_fSum = get() / fDivides; m_fMem = 0.0; } else { m_fSum = (m_fSum + m_fError) / fDivides; } m_fError = 0.0; } inline KahanSum operator*(const KahanSum& fTimes) const { return get() * fTimes.get(); } inline KahanSum operator*(double fTimes) const { return get() * fTimes; } inline KahanSum operator/(const KahanSum& fDivides) const { return o3tl::div_allow_zero(get(), fDivides.get()); } inline KahanSum operator/(double fDivides) const { return get() / fDivides; } inline bool operator<(const KahanSum& fSum) const { return get() < fSum.get(); } inline bool operator<(double fSum) const { return get() < fSum; } inline bool operator>(const KahanSum& fSum) const { return get() > fSum.get(); } inline bool operator>(double fSum) const { return get() > fSum; } inline bool operator<=(const KahanSum& fSum) const { return get() <= fSum.get(); } inline bool operator<=(double fSum) const { return get() <= fSum; } inline bool operator>=(const KahanSum& fSum) const { return get() >= fSum.get(); } inline bool operator>=(double fSum) const { return get() >= fSum; } inline bool operator==(const KahanSum& fSum) const { return get() == fSum.get(); } inline bool operator!=(const KahanSum& fSum) const { return get() != fSum.get(); } public: /** * Returns the final sum. * @return final sum */ double get() const { const double fTotal = m_fSum + m_fError; if (!m_fMem) return fTotal; // Check the same condition as rtl::math::approxAdd() and if true // return 0.0, if false use another Kahan summation adding m_fMem. if (((m_fMem < 0.0 && fTotal > 0.0) || (fTotal < 0.0 && m_fMem > 0.0)) && rtl::math::approxEqual(m_fMem, -fTotal)) { /* TODO: should we reset all values to zero here for further * summation, or is it better to keep them as they are? */ return 0.0; } // The actual argument passed to add() here does not matter as long as // it is not 0, m_fMem is not 0 and will be added anyway, see add(). const_cast(this)->add(m_fMem); const_cast(this)->m_fMem = 0.0; return m_fSum + m_fError; } private: double m_fSum = 0; double m_fError = 0; double m_fMem = 0; }; /* vim:set shiftwidth=4 softtabstop=4 expandtab cinoptions=b1,g0,N-s cinkeys+=0=break: */