]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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<> | |
f67539c2 TL |
12 | #include <cmath> // std::abs |
13 | #include <type_traits> // std::is_floating_point, std::common_type | |
92f5a8d4 TL |
14 | |
15 | namespace boost { | |
16 | namespace histogram { | |
17 | namespace accumulators { | |
18 | ||
19 | /** | |
f67539c2 | 20 | Uses Neumaier algorithm to compute accurate sums of floats. |
92f5a8d4 | 21 | |
f67539c2 TL |
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. | |
92f5a8d4 | 27 | |
f67539c2 | 28 | A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51. |
92f5a8d4 | 29 | */ |
f67539c2 | 30 | template <class ValueType> |
92f5a8d4 | 31 | class sum { |
f67539c2 TL |
32 | static_assert(std::is_floating_point<ValueType>::value, |
33 | "ValueType must be a floating point type"); | |
34 | ||
92f5a8d4 | 35 | public: |
f67539c2 TL |
36 | using value_type = ValueType; |
37 | using const_reference = const value_type&; | |
38 | ||
92f5a8d4 TL |
39 | sum() = default; |
40 | ||
f67539c2 TL |
41 | /// Initialize sum to value and allow implicit conversion |
42 | sum(const_reference value) noexcept : sum(value, 0) {} | |
92f5a8d4 | 43 | |
f67539c2 TL |
44 | /// Allow implicit conversion from sum<T> |
45 | template <class T> | |
46 | sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {} | |
47 | ||
48 | /// Initialize sum explicitly with large and small parts | |
49 | sum(const_reference large, const_reference small) noexcept | |
50 | : large_(large), small_(small) {} | |
92f5a8d4 TL |
51 | |
52 | /// Increment sum by one | |
f67539c2 | 53 | sum& operator++() noexcept { return operator+=(1); } |
92f5a8d4 TL |
54 | |
55 | /// Increment sum by value | |
f67539c2 TL |
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 | |
76 | sum& operator+=(const sum& s) noexcept { | |
77 | operator+=(s.large_); | |
78 | small_ += s.small_; | |
92f5a8d4 TL |
79 | return *this; |
80 | } | |
81 | ||
82 | /// Scale by value | |
f67539c2 | 83 | sum& operator*=(const_reference value) noexcept { |
92f5a8d4 TL |
84 | large_ *= value; |
85 | small_ *= value; | |
86 | return *this; | |
87 | } | |
88 | ||
f67539c2 TL |
89 | bool operator==(const sum& rhs) const noexcept { |
90 | return large_ + small_ == rhs.large_ + rhs.small_; | |
92f5a8d4 TL |
91 | } |
92 | ||
f67539c2 TL |
93 | bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); } |
94 | ||
95 | /// Return value of the sum. | |
96 | value_type value() const noexcept { return large_ + small_; } | |
92f5a8d4 TL |
97 | |
98 | /// Return large part of the sum. | |
f67539c2 | 99 | const_reference large() const noexcept { return large_; } |
92f5a8d4 TL |
100 | |
101 | /// Return small part of the sum. | |
f67539c2 | 102 | const_reference small() const noexcept { return small_; } |
92f5a8d4 | 103 | |
f67539c2 TL |
104 | // lossy conversion to value type must be explicit |
105 | explicit operator value_type() const noexcept { return value(); } | |
92f5a8d4 TL |
106 | |
107 | template <class Archive> | |
108 | void serialize(Archive& ar, unsigned /* version */) { | |
109 | ar& make_nvp("large", large_); | |
110 | ar& make_nvp("small", small_); | |
111 | } | |
112 | ||
f67539c2 TL |
113 | // begin: extra operators to make sum behave like a regular number |
114 | ||
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 | ||
122 | sum operator*(const sum& rhs) const noexcept { | |
123 | sum s = *this; | |
124 | s *= rhs; | |
125 | return s; | |
126 | } | |
127 | ||
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 | ||
135 | sum operator/(const sum& rhs) const noexcept { | |
136 | sum s = *this; | |
137 | s /= rhs; | |
138 | return s; | |
139 | } | |
140 | ||
141 | bool operator<(const sum& rhs) const noexcept { | |
142 | return operator value_type() < rhs.operator value_type(); | |
143 | } | |
144 | ||
145 | bool operator>(const sum& rhs) const noexcept { | |
146 | return operator value_type() > rhs.operator value_type(); | |
147 | } | |
148 | ||
149 | bool operator<=(const sum& rhs) const noexcept { | |
150 | return operator value_type() <= rhs.operator value_type(); | |
151 | } | |
152 | ||
153 | bool operator>=(const sum& rhs) const noexcept { | |
154 | return operator value_type() >= rhs.operator value_type(); | |
155 | } | |
156 | ||
157 | // end: extra operators | |
158 | ||
92f5a8d4 | 159 | private: |
f67539c2 TL |
160 | value_type large_{}; |
161 | value_type small_{}; | |
92f5a8d4 TL |
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 |