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 class KahanSum;
  16 namespace sc::op
  17 {
  18 // Checkout available optimization options.
  19 // Note that it turned out to be problematic to support CPU-specific code
  20 // that's not guaranteed to be available on that specific platform (see
  21 // git commit 2d36e7f5186ba5215f2b228b98c24520bd4f2882). SSE2 is guaranteed on
  22 // x86_64 and it is our baseline requirement for x86 on Windows, so SSE2 use is
  23 // hardcoded on those platforms.
  24 // Whenever we raise baseline to e.g. AVX, this may get
  25 // replaced with AVX code (get it from mentioned git commit).
  26 // Do it similarly with other platforms.
  27 #if defined(X86_64) || (defined(INTEL) && defined(_WIN32))
  28 #define SC_USE_SSE2 1
  29 KahanSum executeSSE2(size_t& i, size_t nSize, const double* pCurrent);
  30 #else
  31 #define SC_USE_SSE2 0
  32 #endif
  33 }
  34
  35 /**
  36   * This class provides LO with Kahan summation algorithm
  37   * About this algorithm: https://en.wikipedia.org/wiki/Kahan_summation_algorithm
  38   * For general purpose software we assume first order error is enough.
  39   *
  40   * Additionally queue and remember the last recent non-zero value and add it
  41   * similar to approxAdd() when obtaining the final result to further eliminate
  42   * accuracy errors. (e.g. for the dreaded 0.1 + 0.2 - 0.3 != 0.0)
  43   */
  44
  45 class KahanSum
  46 {
  47 public:
  48     constexpr KahanSum() = default;
  49
  50     constexpr KahanSum(double x_0)
  51         : m_fSum(x_0)
  52     {
  53     }
  54
  55     constexpr KahanSum(double x_0, double err_0)
  56         : m_fSum(x_0)
  57         , m_fError(err_0)
  58     {
  59     }
  60
  61     constexpr KahanSum(const KahanSum& fSum) = default;
  62
  63 public:
  64     /**
  65       * Performs one step of the Neumaier sum of doubles.
  66       * Overwrites the summand and error.
  67       * T could be double or long double.
  68       */
  69     template <typename T> static inline void sumNeumaierNormal(T& sum, T& err, const double& value)
  70     {
  71         T t = sum + value;
  72         if (std::abs(sum) >= std::abs(value))
  73             err += (sum - t) + value;
  74         else
  75             err += (value - t) + sum;
  76         sum = t;
  77     }
  78
  79     /**
  80       * Adds a value to the sum using Kahan summation.
  81       * @param x_i
  82       */
  83     void add(double x_i)
  84     {
  85         if (x_i == 0.0)
  86             return;
  87
  88         if (!m_fMem)
  89         {
  90             m_fMem = x_i;
  91             return;
  92         }
  93
  94         sumNeumaierNormal(m_fSum, m_fError, m_fMem);
  95         m_fMem = x_i;
  96     }
  97
  98     /**
  99       * Adds a value to the sum using Kahan summation.
 100       * @param fSum
 101       */
 102     inline void add(const KahanSum& fSum)
 103     {
 104 #if SC_USE_SSE2
 105         add(fSum.m_fSum + fSum.m_fError);
 106         add(fSum.m_fMem);
 107 #else
 108         // Without SSE2 the sum+compensation value fails badly. Continue
 109         // keeping the old though slightly off (see tdf#156985) explicit
 110         // addition of the compensation value.
 111         double sum = fSum.m_fSum;
 112         double err = fSum.m_fError;
 113         if (fSum.m_fMem != 0.0)
 114             sumNeumaierNormal(sum, err, fSum.m_fMem);
 115         add(sum);
 116         add(err);
 117 #endif
 118     }
 119
 120     /**
 121       * Substracts a value to the sum using Kahan summation.
 122       * @param fSum
 123       */
 124     inline void subtract(const KahanSum& fSum) { add(-fSum); }
 125
 126 public:
 127     constexpr KahanSum operator-() const
 128     {
 129         KahanSum fKahanSum;
 130         fKahanSum.m_fSum = -m_fSum;
 131         fKahanSum.m_fError = -m_fError;
 132         fKahanSum.m_fMem = -m_fMem;
 133         return fKahanSum;
 134     }
 135
 136     constexpr KahanSum& operator=(double fSum)
 137     {
 138         m_fSum = fSum;
 139         m_fError = 0;
 140         m_fMem = 0;
 141         return *this;
 142     }
 143
 144     constexpr KahanSum& operator=(const KahanSum& fSum) = default;
 145
 146     inline void operator+=(const KahanSum& fSum) { add(fSum); }
 147
 148     inline void operator+=(double fSum) { add(fSum); }
 149
 150     inline void operator-=(const KahanSum& fSum) { subtract(fSum); }
 151
 152     inline void operator-=(double fSum) { add(-fSum); }
 153
 154     inline KahanSum operator+(double fSum) const
 155     {
 156         KahanSum fNSum(*this);
 157         fNSum.add(fSum);
 158         return fNSum;
 159     }
 160
 161     inline KahanSum operator+(const KahanSum& fSum) const
 162     {
 163         KahanSum fNSum(*this);
 164         fNSum += fSum;
 165         return fNSum;
 166     }
 167
 168     inline KahanSum operator-(double fSum) const
 169     {
 170         KahanSum fNSum(*this);
 171         fNSum.add(-fSum);
 172         return fNSum;
 173     }
 174
 175     inline KahanSum operator-(const KahanSum& fSum) const
 176     {
 177         KahanSum fNSum(*this);
 178         fNSum -= fSum;
 179         return fNSum;
 180     }
 181
 182     /**
 183       * In some parts of the code of interpr_.cxx this may be used for
 184       * product instead of sum. This operator shall be used for that task.
 185       */
 186     constexpr void operator*=(double fTimes)
 187     {
 188         if (m_fMem)
 189         {
 190             m_fSum = get() * fTimes;
 191             m_fMem = 0.0;
 192         }
 193         else
 194         {
 195             m_fSum = (m_fSum + m_fError) * fTimes;
 196         }
 197         m_fError = 0.0;
 198     }
 199
 200     constexpr void operator/=(double fDivides)
 201     {
 202         if (m_fMem)
 203         {
 204             m_fSum = get() / fDivides;
 205             m_fMem = 0.0;
 206         }
 207         else
 208         {
 209             m_fSum = (m_fSum + m_fError) / fDivides;
 210         }
 211         m_fError = 0.0;
 212     }
 213
 214     inline KahanSum operator*(const KahanSum& fTimes) const { return get() * fTimes.get(); }
 215
 216     inline KahanSum operator*(double fTimes) const { return get() * fTimes; }
 217
 218     inline KahanSum operator/(const KahanSum& fDivides) const { return get() / fDivides.get(); }
 219
 220     inline KahanSum operator/(double fDivides) const { return get() / fDivides; }
 221
 222     inline bool operator<(const KahanSum& fSum) const { return get() < fSum.get(); }
 223
 224     inline bool operator<(double fSum) const { return get() < fSum; }
 225
 226     inline bool operator>(const KahanSum& fSum) const { return get() > fSum.get(); }
 227
 228     inline bool operator>(double fSum) const { return get() > fSum; }
 229
 230     inline bool operator<=(const KahanSum& fSum) const { return get() <= fSum.get(); }
 231
 232     inline bool operator<=(double fSum) const { return get() <= fSum; }
 233
 234     inline bool operator>=(const KahanSum& fSum) const { return get() >= fSum.get(); }
 235
 236     inline bool operator>=(double fSum) const { return get() >= fSum; }
 237
 238     inline bool operator==(const KahanSum& fSum) const { return get() == fSum.get(); }
 239
 240     inline bool operator!=(const KahanSum& fSum) const { return get() != fSum.get(); }
 241
 242 public:
 243     /**
 244       * Returns the final sum.
 245       * @return final sum
 246       */
 247     double get() const
 248     {
 249         const double fTotal = m_fSum + m_fError;
 250         if (!m_fMem)
 251             return fTotal;
 252
 253         // Check the same condition as rtl::math::approxAdd() and if true
 254         // return 0.0, if false use another Kahan summation adding m_fMem.
 255         if (((m_fMem < 0.0 && fTotal > 0.0) || (fTotal < 0.0 && m_fMem > 0.0))
 256             && rtl::math::approxEqual(m_fMem, -fTotal))
 257         {
 258             /* TODO: should we reset all values to zero here for further
 259              * summation, or is it better to keep them as they are? */
 260             return 0.0;
 261         }
 262
 263         // The actual argument passed to add() here does not matter as long as
 264         // it is not 0, m_fMem is not 0 and will be added anyway, see add().
 265         const_cast<KahanSum*>(this)->add(m_fMem);
 266         const_cast<KahanSum*>(this)->m_fMem = 0.0;
 267         return m_fSum + m_fError;
 268     }
 269
 270 private:
 271     double m_fSum = 0;
 272     double m_fError = 0;
 273     double m_fMem = 0;
 274 };
 275
 276 /* vim:set shiftwidth=4 softtabstop=4 expandtab cinoptions=b1,g0,N-s cinkeys+=0=break: */