sc/inc/kahan.hxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; fill-column: 100 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  */
   9
  10 #pragma once
  11
  12 #include <rtl/math.hxx>
  13 #include <cmath>
  14
  15 /**
  16   * This class provides LO with Kahan summation algorithm
  17   * About this algorithm: https://en.wikipedia.org/wiki/Kahan_summation_algorithm
  18   * For general purpose software we assume first order error is enough.
  19   *
  20   * Additionally queue and remember the last recent non-zero value and add it
  21   * similar to approxAdd() when obtaining the final result to further eliminate
  22   * accuracy errors. (e.g. for the dreaded 0.1 + 0.2 - 0.3 != 0.0)
  23   */
  24
  25 class KahanSum
  26 {
  27 public:
  28     constexpr KahanSum() = default;
  29
  30     constexpr KahanSum(double x_0)
  31         : m_fSum(x_0)
  32     {
  33     }
  34
  35     constexpr KahanSum(double x_0, double err_0)
  36         : m_fSum(x_0)
  37         , m_fError(err_0)
  38     {
  39     }
  40
  41     constexpr KahanSum(const KahanSum& fSum) = default;
  42
  43 public:
  44     /**
  45       * Adds a value to the sum using Kahan summation.
  46       * @param x_i
  47       */
  48     void add(double x_i)
  49     {
  50         if (x_i == 0.0)
  51             return;
  52
  53         if (!m_fMem)
  54         {
  55             m_fMem = x_i;
  56             return;
  57         }
  58
  59         double t = m_fSum + m_fMem;
  60         if (std::abs(m_fSum) >= std::abs(m_fMem))
  61             m_fError += (m_fSum - t) + m_fMem;
  62         else
  63             m_fError += (m_fMem - t) + m_fSum;
  64         m_fSum = t;
  65         m_fMem = x_i;
  66     }
  67
  68     /**
  69       * Adds a value to the sum using Kahan summation.
  70       * @param fSum
  71       */
  72     inline void add(const KahanSum& fSum)
  73     {
  74         add(fSum.m_fSum);
  75         add(fSum.m_fError);
  76         add(fSum.m_fMem);
  77     }
  78
  79     /**
  80       * Substracts a value to the sum using Kahan summation.
  81       * @param fSum
  82       */
  83     inline void subtract(const KahanSum& fSum)
  84     {
  85         add(-fSum.m_fSum);
  86         add(-fSum.m_fError);
  87         add(-fSum.m_fMem);
  88     }
  89
  90 public:
  91     constexpr KahanSum operator-() const
  92     {
  93         KahanSum fKahanSum;
  94         fKahanSum.m_fSum = -m_fSum;
  95         fKahanSum.m_fError = -m_fError;
  96         fKahanSum.m_fMem = -m_fMem;
  97         return fKahanSum;
  98     }
  99
 100     constexpr KahanSum& operator=(double fSum)
 101     {
 102         m_fSum = fSum;
 103         m_fError = 0;
 104         m_fMem = 0;
 105         return *this;
 106     }
 107
 108     constexpr KahanSum& operator=(const KahanSum& fSum) = default;
 109
 110     inline void operator+=(const KahanSum& fSum) { add(fSum); }
 111
 112     inline void operator+=(double fSum) { add(fSum); }
 113
 114     inline void operator-=(const KahanSum& fSum) { subtract(fSum); }
 115
 116     inline void operator-=(double fSum) { add(-fSum); }
 117
 118     inline KahanSum operator+(double fSum) const
 119     {
 120         KahanSum fNSum(*this);
 121         fNSum.add(fSum);
 122         return fNSum;
 123     }
 124
 125     inline KahanSum operator+(const KahanSum& fSum) const
 126     {
 127         KahanSum fNSum(*this);
 128         fNSum += fSum;
 129         return fNSum;
 130     }
 131
 132     inline KahanSum operator-(double fSum) const
 133     {
 134         KahanSum fNSum(*this);
 135         fNSum.add(-fSum);
 136         return fNSum;
 137     }
 138
 139     inline KahanSum operator-(const KahanSum& fSum) const
 140     {
 141         KahanSum fNSum(*this);
 142         fNSum -= fSum;
 143         return fNSum;
 144     }
 145
 146     /**
 147       * In some parts of the code of interpr_.cxx this may be used for
 148       * product instead of sum. This operator shall be used for that task.
 149       */
 150     constexpr void operator*=(double fTimes)
 151     {
 152         if (m_fMem)
 153         {
 154             m_fSum = get() * fTimes;
 155             m_fMem = 0.0;
 156         }
 157         else
 158         {
 159             m_fSum = (m_fSum + m_fError) * fTimes;
 160         }
 161         m_fError = 0.0;
 162     }
 163
 164     constexpr void operator/=(double fDivides)
 165     {
 166         if (m_fMem)
 167         {
 168             m_fSum = get() / fDivides;
 169             m_fMem = 0.0;
 170         }
 171         else
 172         {
 173             m_fSum = (m_fSum + m_fError) / fDivides;
 174         }
 175         m_fError = 0.0;
 176     }
 177
 178     inline KahanSum operator*(const KahanSum& fTimes) const { return get() * fTimes.get(); }
 179
 180     inline KahanSum operator*(double fTimes) const { return get() * fTimes; }
 181
 182     inline KahanSum operator/(const KahanSum& fDivides) const { return get() / fDivides.get(); }
 183
 184     inline KahanSum operator/(double fDivides) const { return get() / fDivides; }
 185
 186     inline bool operator<(const KahanSum& fSum) const { return get() < fSum.get(); }
 187
 188     inline bool operator<(double fSum) const { return get() < fSum; }
 189
 190     inline bool operator>(const KahanSum& fSum) const { return get() > fSum.get(); }
 191
 192     inline bool operator>(double fSum) const { return get() > fSum; }
 193
 194     inline bool operator<=(const KahanSum& fSum) const { return get() <= fSum.get(); }
 195
 196     inline bool operator<=(double fSum) const { return get() <= fSum; }
 197
 198     inline bool operator>=(const KahanSum& fSum) const { return get() >= fSum.get(); }
 199
 200     inline bool operator>=(double fSum) const { return get() >= fSum; }
 201
 202     inline bool operator==(const KahanSum& fSum) const { return get() == fSum.get(); }
 203
 204     inline bool operator!=(const KahanSum& fSum) const { return get() != fSum.get(); }
 205
 206 public:
 207     /**
 208       * Returns the final sum.
 209       * @return final sum
 210       */
 211     double get() const
 212     {
 213         const double fTotal = m_fSum + m_fError;
 214         if (!m_fMem)
 215             return fTotal;
 216
 217         // Check the same condition as rtl::math::approxAdd() and if true
 218         // return 0.0, if false use another Kahan summation adding m_fMem.
 219         if (((m_fMem < 0.0 && fTotal > 0.0) || (fTotal < 0.0 && m_fMem > 0.0))
 220             && rtl::math::approxEqual(m_fMem, -fTotal))
 221         {
 222             /* TODO: should we reset all values to zero here for further
 223              * summation, or is it better to keep them as they are? */
 224             return 0.0;
 225         }
 226
 227         // The actual argument passed to add() here does not matter as long as
 228         // it is not 0, m_fMem is not 0 and will be added anyway, see add().
 229         const_cast<KahanSum*>(this)->add(m_fMem);
 230         const_cast<KahanSum*>(this)->m_fMem = 0.0;
 231         return m_fSum + m_fError;
 232     }
 233
 234 private:
 235     double m_fSum = 0;
 236     double m_fError = 0;
 237     double m_fMem = 0;
 238 };
 239
 240 /* vim:set shiftwidth=4 softtabstop=4 expandtab cinoptions=b1,g0,N-s cinkeys+=0=break: */