From 543782d78b6fb5ec57eaac3eeb0cd62cffe4d2c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= Date: Sun, 25 Dec 2016 22:50:14 +0100 Subject: [PATCH] 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(). --- src/Magnum/Math/Algorithms/CMakeLists.txt | 1 + src/Magnum/Math/Algorithms/KahanSum.h | 85 +++++++++ .../Math/Algorithms/Test/CMakeLists.txt | 1 + .../Math/Algorithms/Test/KahanSumTest.cpp | 166 ++++++++++++++++++ 4 files changed, 253 insertions(+) create mode 100644 src/Magnum/Math/Algorithms/KahanSum.h create mode 100644 src/Magnum/Math/Algorithms/Test/KahanSumTest.cpp diff --git a/src/Magnum/Math/Algorithms/CMakeLists.txt b/src/Magnum/Math/Algorithms/CMakeLists.txt index cfddc216e..af9b46620 100644 --- 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) diff --git a/src/Magnum/Math/Algorithms/KahanSum.h b/src/Magnum/Math/Algorithms/KahanSum.h new file mode 100644 index 000000000..d601d7825 --- /dev/null +++ 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š + + 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 + +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 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 pixels; +Float sum = 0.0f, c = 0.0f; +for(UnsignedByte pixel: pixels) { + Float value = Math::normalize(pixel); + sum = Math::Algorithms::kahanSum(&value, &value + 1, sum, &c); +} +@endcode +*/ +template())>::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 diff --git a/src/Magnum/Math/Algorithms/Test/CMakeLists.txt b/src/Magnum/Math/Algorithms/Test/CMakeLists.txt index 926d21613..02b948ee1 100644 --- 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) diff --git a/src/Magnum/Math/Algorithms/Test/KahanSumTest.cpp b/src/Magnum/Math/Algorithms/Test/KahanSumTest.cpp new file mode 100644 index 000000000..d5ea05ab9 --- /dev/null +++ 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š + + 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 +#include + +#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 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 + #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 begin{1.0f}; + Iterator 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 begin{1}; + Iterator end{1, 100000}; + + CORRADE_COMPARE(kahanSum(begin, end), 100000); + CORRADE_COMPARE(std::accumulate(begin, end, 0.0f), 100000); +} + +void KahanSumTest::iterative() { + Iterator begin{1.0f}; + Iterator 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 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 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)