1 // Copyright 2018 Hans Dembinski
3 // Distributed under the Boost Software License, version 1.0.
4 // (See accompanying file LICENSE_1_0.txt
5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
8 #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
10 #include <boost/core/nvp.hpp>
11 #include <boost/histogram/fwd.hpp> // for sum<>
13 #include <type_traits>
17 namespace accumulators {
20 Uses Neumaier algorithm to compute accurate sums.
22 The algorithm uses memory for two floats and is three to
23 five times slower compared to a simple floating point
24 number used to accumulate a sum, but the relative error
25 of the sum is at the level of the machine precision,
26 independent of the number of samples.
28 A. Neumaier, Zeitschrift fuer Angewandte Mathematik
29 und Mechanik 54 (1974) 39-51.
31 template <typename RealType>
36 /// Initialize sum to value
37 explicit sum(const RealType& value) noexcept : large_(value) {}
40 sum& operator=(const RealType& value) noexcept {
46 /// Increment sum by one
47 sum& operator++() { return operator+=(1); }
49 /// Increment sum by value
50 sum& operator+=(const RealType& value) {
51 auto temp = large_ + value; // prevent optimization
52 if (std::abs(large_) >= std::abs(value))
53 small_ += (large_ - temp) + value;
55 small_ += (value - temp) + large_;
61 sum& operator*=(const RealType& value) {
68 bool operator==(const sum<T>& rhs) const noexcept {
69 return large_ == rhs.large_ && small_ == rhs.small_;
73 bool operator!=(const T& rhs) const noexcept {
74 return !operator==(rhs);
77 /// Return large part of the sum.
78 const RealType& large() const { return large_; }
80 /// Return small part of the sum.
81 const RealType& small() const { return small_; }
83 // allow implicit conversion to RealType
84 operator RealType() const { return large_ + small_; }
86 template <class Archive>
87 void serialize(Archive& ar, unsigned /* version */) {
88 ar& make_nvp("large", large_);
89 ar& make_nvp("small", small_);
93 RealType large_ = RealType();
94 RealType small_ = RealType();
97 } // namespace accumulators
98 } // namespace histogram
101 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
103 template <class T, class U>
104 struct common_type<boost::histogram::accumulators::sum<T>,
105 boost::histogram::accumulators::sum<U>> {
106 using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;