1 // Copyright 2018 Hans Dembinski 2 // 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) 6 7 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP 8 #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP 9 10 #include <boost/core/nvp.hpp> 11 #include <boost/histogram/fwd.hpp> // for sum<> 12 #include <cmath> // std::abs 13 #include <type_traits> // std::is_floating_point, std::common_type 14 15 namespace boost { 16 namespace histogram { 17 namespace accumulators { 18 19 /** 20 Uses Neumaier algorithm to compute accurate sums of floats. 21 22 The algorithm is an improved Kahan algorithm 23 (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). The algorithm uses memory for 24 two numbers and is three to five times slower compared to using a single number to 25 accumulate a sum, but the relative error of the sum is at the level of the machine 26 precision, independent of the number of samples. 27 28 A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51. 29 */ 30 template <class ValueType> 31 class sum { 32 static_assert(std::is_floating_point<ValueType>::value, 33 "ValueType must be a floating point type"); 34 35 public: 36 using value_type = ValueType; 37 using const_reference = const value_type&; 38 39 sum() = default; 40 41 /// Initialize sum to value and allow implicit conversion sum(const_reference value)42 sum(const_reference value) noexcept : sum(value, 0) {} 43 44 /// Allow implicit conversion from sum<T> 45 template <class T> sum(const sum<T> & s)46 sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {} 47 48 /// Initialize sum explicitly with large and small parts sum(const_reference large,const_reference small)49 sum(const_reference large, const_reference small) noexcept 50 : large_(large), small_(small) {} 51 52 /// Increment sum by one operator ++()53 sum& operator++() noexcept { return operator+=(1); } 54 55 /// Increment sum by value operator +=(const_reference value)56 sum& operator+=(const_reference value) noexcept { 57 // prevent compiler optimization from destroying the algorithm 58 // when -ffast-math is enabled 59 volatile value_type l; 60 value_type s; 61 if (std::abs(large_) >= std::abs(value)) { 62 l = large_; 63 s = value; 64 } else { 65 l = value; 66 s = large_; 67 } 68 large_ += value; 69 l -= large_; 70 l += s; 71 small_ += l; 72 return *this; 73 } 74 75 /// Add another sum operator +=(const sum & s)76 sum& operator+=(const sum& s) noexcept { 77 operator+=(s.large_); 78 small_ += s.small_; 79 return *this; 80 } 81 82 /// Scale by value operator *=(const_reference value)83 sum& operator*=(const_reference value) noexcept { 84 large_ *= value; 85 small_ *= value; 86 return *this; 87 } 88 operator ==(const sum & rhs) const89 bool operator==(const sum& rhs) const noexcept { 90 return large_ + small_ == rhs.large_ + rhs.small_; 91 } 92 operator !=(const sum & rhs) const93 bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); } 94 95 /// Return value of the sum. value() const96 value_type value() const noexcept { return large_ + small_; } 97 98 /// Return large part of the sum. large() const99 const_reference large() const noexcept { return large_; } 100 101 /// Return small part of the sum. small() const102 const_reference small() const noexcept { return small_; } 103 104 // lossy conversion to value type must be explicit operator value_type() const105 explicit operator value_type() const noexcept { return value(); } 106 107 template <class Archive> serialize(Archive & ar,unsigned)108 void serialize(Archive& ar, unsigned /* version */) { 109 ar& make_nvp("large", large_); 110 ar& make_nvp("small", small_); 111 } 112 113 // begin: extra operators to make sum behave like a regular number 114 operator *=(const sum & rhs)115 sum& operator*=(const sum& rhs) noexcept { 116 const auto scale = static_cast<value_type>(rhs); 117 large_ *= scale; 118 small_ *= scale; 119 return *this; 120 } 121 operator *(const sum & rhs) const122 sum operator*(const sum& rhs) const noexcept { 123 sum s = *this; 124 s *= rhs; 125 return s; 126 } 127 operator /=(const sum & rhs)128 sum& operator/=(const sum& rhs) noexcept { 129 const auto scale = 1.0 / static_cast<value_type>(rhs); 130 large_ *= scale; 131 small_ *= scale; 132 return *this; 133 } 134 operator /(const sum & rhs) const135 sum operator/(const sum& rhs) const noexcept { 136 sum s = *this; 137 s /= rhs; 138 return s; 139 } 140 operator <(const sum & rhs) const141 bool operator<(const sum& rhs) const noexcept { 142 return operator value_type() < rhs.operator value_type(); 143 } 144 operator >(const sum & rhs) const145 bool operator>(const sum& rhs) const noexcept { 146 return operator value_type() > rhs.operator value_type(); 147 } 148 operator <=(const sum & rhs) const149 bool operator<=(const sum& rhs) const noexcept { 150 return operator value_type() <= rhs.operator value_type(); 151 } 152 operator >=(const sum & rhs) const153 bool operator>=(const sum& rhs) const noexcept { 154 return operator value_type() >= rhs.operator value_type(); 155 } 156 157 // end: extra operators 158 159 private: 160 value_type large_{}; 161 value_type small_{}; 162 }; 163 164 } // namespace accumulators 165 } // namespace histogram 166 } // namespace boost 167 168 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED 169 namespace std { 170 template <class T, class U> 171 struct common_type<boost::histogram::accumulators::sum<T>, 172 boost::histogram::accumulators::sum<U>> { 173 using type = boost::histogram::accumulators::sum<common_type_t<T, U>>; 174 }; 175 } // namespace std 176 #endif 177 178 #endif 179