Math: added Algorithms::kahanSum().

The test is running pretty slow on debug builds, though (sorry). I see no way how to make it faster while still exposing the issues with std::accumulate().
9 years ago · 543782d78b
4 changed files with 253 additions and 0 deletions
--- a/src/Magnum/Math/Algorithms/CMakeLists.txt
+++ b/src/Magnum/Math/Algorithms/CMakeLists.txt
@ -26,6 +26,7 @@
 set(MagnumMathAlgorithms_HEADERS
    GaussJordan.h
    GramSchmidt.h
+    KahanSum.h
    Qr.h
    Svd.h)

--- a/src/Magnum/Math/Algorithms/KahanSum.h
+++ b/src/Magnum/Math/Algorithms/KahanSum.h
@ -0,0 +1,85 @@
+#ifndef Magnum_Math_Algorithms_KahanSum_h
+#define Magnum_Math_Algorithms_KahanSum_h
+/*
+    This file is part of Magnum.
+
+    Copyright © 2010, 2011, 2012, 2013, 2014, 2015, 2016
+              Vladimír Vondruš <mosra@centrum.cz>
+
+    Permission is hereby granted, free of charge, to any person obtaining a
+    copy of this software and associated documentation files (the "Software"),
+    to deal in the Software without restriction, including without limitation
+    the rights to use, copy, modify, merge, publish, distribute, sublicense,
+    and/or sell copies of the Software, and to permit persons to whom the
+    Software is furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included
+    in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+    THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+    DEALINGS IN THE SOFTWARE.
+*/
+
+/** @file
+ * @brief Function @ref Magnum::Math::Algorithms::kahanSum()
+ */
+
+#include <utility>
+
+namespace Magnum { namespace Math { namespace Algorithms {
+
+/**
+@brief Kahan summation algorithm
+@param[in] begin            Range begin
+@param[in] end              Range end
+@param[in] sum              Initial value for the sum
+@param[in,out] compensation Floating-point roundoff error compensation value.
+    If non-`nullptr`, value behind the pointer is used as initial compensation
+    value and the resulting value is
+
+Calculates a sum of a large range of floating-point numbers with roundoff error
+compensation. Compared to for example `std::accumulate()` the algorithm
+significantly reduces numerical error in the total. See Wikipedia for an
+in-depth explanation: https://en.wikipedia.org/wiki/Kahan_summation_algorithm
+
+Example with summation of a hundred million ones:
+@code
+std::vector<Float> data(100000000, 1.0f);
+Float a = std::accumulate(data.begin(), data.end());            // 1.667e7f
+Float b = Math::Algorithms::kahanSum(data.begin(), data.end()); // 1.000e8f
+@endcode
+
+If required, it is also possible to use this algorithm on non-contiguous ranges
+or single values (for example when calculating sum of pixel values in an image
+with some row padding or when the inputs are generated / converted from other
+values):
+@code
+Containers::ArrayView<UnsignedByte> pixels;
+Float sum = 0.0f, c = 0.0f;
+for(UnsignedByte pixel: pixels) {
+    Float value = Math::normalize<Float>(pixel);
+    sum = Math::Algorithms::kahanSum(&value, &value + 1, sum, &c);
+}
+@endcode
+*/
+template<class Iterator, class T = typename std::decay<decltype(*std::declval<Iterator>())>::type> T kahanSum(Iterator begin, Iterator end, T sum = T(0), T* compensation = nullptr) {
+    T c = compensation ? *compensation : T(0);
+    for(Iterator it = begin; it != end; ++it) {
+        const T y = *it - c;
+        const T t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+    }
+
+    if(compensation) *compensation = c;
+    return sum;
+}
+
+}}}
+
+#endif
--- a/src/Magnum/Math/Algorithms/Test/CMakeLists.txt
+++ b/src/Magnum/Math/Algorithms/Test/CMakeLists.txt
@ -25,5 +25,6 @@

 corrade_add_test(MathAlgorithmsGaussJordanTest GaussJordanTest.cpp LIBRARIES MagnumMathTestLib)
 corrade_add_test(MathAlgorithmsGramSchmidtTest GramSchmidtTest.cpp LIBRARIES MagnumMathTestLib)
+corrade_add_test(MathAlgorithmsKahanSumTest KahanSumTest.cpp LIBRARIES MagnumMathTestLib)
 corrade_add_test(MathAlgorithmsQrTest QrTest.cpp LIBRARIES MagnumMathTestLib)
 corrade_add_test(MathAlgorithmsSvdTest SvdTest.cpp LIBRARIES MagnumMathTestLib)
--- a/src/Magnum/Math/Algorithms/Test/KahanSumTest.cpp
+++ b/src/Magnum/Math/Algorithms/Test/KahanSumTest.cpp
@ -0,0 +1,166 @@
+/*
+    This file is part of Magnum.
+
+    Copyright © 2010, 2011, 2012, 2013, 2014, 2015, 2016
+              Vladimír Vondruš <mosra@centrum.cz>
+
+    Permission is hereby granted, free of charge, to any person obtaining a
+    copy of this software and associated documentation files (the "Software"),
+    to deal in the Software without restriction, including without limitation
+    the rights to use, copy, modify, merge, publish, distribute, sublicense,
+    and/or sell copies of the Software, and to permit persons to whom the
+    Software is furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included
+    in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+    THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+    DEALINGS IN THE SOFTWARE.
+*/
+
+#include <numeric>
+#include <Corrade/TestSuite/Tester.h>
+
+#include "Magnum/Magnum.h"
+#include "Magnum/Math/Algorithms/KahanSum.h"
+
+namespace Magnum { namespace Math { namespace Algorithms { namespace Test {
+
+struct KahanSumTest: TestSuite::Tester {
+    explicit KahanSumTest();
+
+    void floats();
+    void integers();
+    void iterative();
+
+    void accumulate100k();
+    void kahan100k();
+
+};
+
+KahanSumTest::KahanSumTest() {
+    addTests({&KahanSumTest::floats,
+              &KahanSumTest::integers,
+              &KahanSumTest::iterative});
+
+    addBenchmarks({&KahanSumTest::accumulate100k,
+                   &KahanSumTest::kahan100k}, 50);
+}
+
+namespace {
+
+/* Custom iterator class to avoid allocating half a gigabyte for hundred
+   million values */
+template<class T> struct Iterator
+    #ifdef CORRADE_MSVC2015_COMPATIBILITY
+    /* This is needed on MSVC in order to have std::accumulate() eat these
+       iterators properly. However, std::iterator is being deprecated in C++17
+       and the go-to way of typedef'ing all the things doesn't work (the
+       compiler complains about iterator_category not being defined in
+       std::iterator_traits) */
+    : std::iterator<std::forward_iterator_tag, T>
+    #endif
+{
+    explicit Iterator(T value, std::size_t i = 0): _value{value}, _i{i} {}
+
+    bool operator!=(const Iterator& other) const { return _i != other._i; }
+
+    Iterator& operator++() {
+        ++_i;
+        return *this;
+    }
+
+    T operator*() const { return _value; }
+
+    private:
+        T _value{};
+        std::size_t _i{};
+};
+
+}
+
+void KahanSumTest::floats() {
+    Iterator<Float> begin{1.0f};
+    Iterator<Float> end{1.0f, 100000000};
+
+    CORRADE_COMPARE(kahanSum(begin, end), 1.0e8f);
+    {
+        #ifdef CORRADE_TARGET_EMSCRIPTEN
+        /* Last tested: emscripten 1.37.1 */
+        CORRADE_EXPECT_FAIL("Gives the same result as kahanSum(), might be "
+            "because doubles are used internally.");
+        #endif
+        CORRADE_COMPARE(std::accumulate(begin, end, 0.0f), 1.6777216e7f);
+    }
+}
+
+void KahanSumTest::integers() {
+    Iterator<Int> begin{1};
+    Iterator<Int> end{1, 100000};
+
+    CORRADE_COMPARE(kahanSum(begin, end), 100000);
+    CORRADE_COMPARE(std::accumulate(begin, end, 0.0f), 100000);
+}
+
+void KahanSumTest::iterative() {
+    Iterator<Float> begin{1.0f};
+    Iterator<Float> end{1.0f, 100000000};
+
+    Float sum{}, sumKahan{}, c{};
+    for(auto it = begin; it != end; ++it) {
+        Float d = *it;
+        sum += d;
+        sumKahan = kahanSum(&d, &d + 1, sumKahan, &c);
+    }
+
+    /* This is very weird on Emscripten. */
+    {
+        #ifdef CORRADE_TARGET_EMSCRIPTEN
+        /* Last tested: emscripten 1.37.1 */
+        CORRADE_EXPECT_FAIL_IF(sumKahan == 1.6777216e7f,
+            "Probably because of some misoptimization, happens only on -O1 "
+            "and not on larger optimization levels.");
+        #endif
+        CORRADE_COMPARE(sumKahan, 1.0e8f);
+    } {
+        #ifdef CORRADE_TARGET_EMSCRIPTEN
+        /* Last tested: emscripten 1.37.1 */
+        CORRADE_EXPECT_FAIL_IF(sumKahan == 1.0e8f,
+            "Gives the same result as kahanSum(), might be because doubles are "
+            "used internally. Happens only on larger optimization levels, not "
+            "on -O1.");
+        #endif
+        CORRADE_COMPARE(sum, 1.6777216e7f);
+    }
+}
+
+void KahanSumTest::accumulate100k() {
+    std::vector<Float> data(100000, 1.0f);
+
+    volatile Float a; /* to avoid optimizing the loop out */
+    CORRADE_BENCHMARK(10) {
+        a = std::accumulate(data.begin(), data.end(), 0.0f);
+    }
+
+    CORRADE_COMPARE(Float(a), 100000.0f);
+}
+
+void KahanSumTest::kahan100k() {
+    std::vector<Float> data(100000, 1.0f);
+
+    volatile Float a; /* to avoid optimizing the loop out */
+    CORRADE_BENCHMARK(10) {
+        a = kahanSum(data.begin(), data.end());
+    }
+
+    CORRADE_COMPARE(Float(a), 100000.0f);
+}
+
+}}}}
+
+CORRADE_TEST_MAIN(Magnum::Math::Algorithms::Test::KahanSumTest)