123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178 |
- // Copyright 2018 Hans Dembinski
- //
- // Distributed under the Boost Software License, version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
- #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
- #include <boost/core/nvp.hpp>
- #include <boost/histogram/fwd.hpp> // for sum<>
- #include <cmath> // std::abs
- #include <type_traits> // std::is_floating_point, std::common_type
- namespace boost {
- namespace histogram {
- namespace accumulators {
- /**
- Uses Neumaier algorithm to compute accurate sums of floats.
- The algorithm is an improved Kahan algorithm
- (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). The algorithm uses memory for
- two numbers and is three to five times slower compared to using a single number to
- accumulate a sum, but the relative error of the sum is at the level of the machine
- precision, independent of the number of samples.
- A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51.
- */
- template <class ValueType>
- class sum {
- static_assert(std::is_floating_point<ValueType>::value,
- "ValueType must be a floating point type");
- public:
- using value_type = ValueType;
- using const_reference = const value_type&;
- sum() = default;
- /// Initialize sum to value and allow implicit conversion
- sum(const_reference value) noexcept : sum(value, 0) {}
- /// Allow implicit conversion from sum<T>
- template <class T>
- sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {}
- /// Initialize sum explicitly with large and small parts
- sum(const_reference large, const_reference small) noexcept
- : large_(large), small_(small) {}
- /// Increment sum by one
- sum& operator++() noexcept { return operator+=(1); }
- /// Increment sum by value
- sum& operator+=(const_reference value) noexcept {
- // prevent compiler optimization from destroying the algorithm
- // when -ffast-math is enabled
- volatile value_type l;
- value_type s;
- if (std::abs(large_) >= std::abs(value)) {
- l = large_;
- s = value;
- } else {
- l = value;
- s = large_;
- }
- large_ += value;
- l = l - large_;
- l = l + s;
- small_ += l;
- return *this;
- }
- /// Add another sum
- sum& operator+=(const sum& s) noexcept {
- operator+=(s.large_);
- small_ += s.small_;
- return *this;
- }
- /// Scale by value
- sum& operator*=(const_reference value) noexcept {
- large_ *= value;
- small_ *= value;
- return *this;
- }
- bool operator==(const sum& rhs) const noexcept {
- return large_ + small_ == rhs.large_ + rhs.small_;
- }
- bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); }
- /// Return value of the sum.
- value_type value() const noexcept { return large_ + small_; }
- /// Return large part of the sum.
- const_reference large() const noexcept { return large_; }
- /// Return small part of the sum.
- const_reference small() const noexcept { return small_; }
- // lossy conversion to value type must be explicit
- explicit operator value_type() const noexcept { return value(); }
- template <class Archive>
- void serialize(Archive& ar, unsigned /* version */) {
- ar& make_nvp("large", large_);
- ar& make_nvp("small", small_);
- }
- // begin: extra operators to make sum behave like a regular number
- sum& operator*=(const sum& rhs) noexcept {
- const auto scale = static_cast<value_type>(rhs);
- large_ *= scale;
- small_ *= scale;
- return *this;
- }
- sum operator*(const sum& rhs) const noexcept {
- sum s = *this;
- s *= rhs;
- return s;
- }
- sum& operator/=(const sum& rhs) noexcept {
- const auto scale = 1.0 / static_cast<value_type>(rhs);
- large_ *= scale;
- small_ *= scale;
- return *this;
- }
- sum operator/(const sum& rhs) const noexcept {
- sum s = *this;
- s /= rhs;
- return s;
- }
- bool operator<(const sum& rhs) const noexcept {
- return operator value_type() < rhs.operator value_type();
- }
- bool operator>(const sum& rhs) const noexcept {
- return operator value_type() > rhs.operator value_type();
- }
- bool operator<=(const sum& rhs) const noexcept {
- return operator value_type() <= rhs.operator value_type();
- }
- bool operator>=(const sum& rhs) const noexcept {
- return operator value_type() >= rhs.operator value_type();
- }
- // end: extra operators
- private:
- value_type large_{};
- value_type small_{};
- };
- } // namespace accumulators
- } // namespace histogram
- } // namespace boost
- #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
- namespace std {
- template <class T, class U>
- struct common_type<boost::histogram::accumulators::sum<T>,
- boost::histogram::accumulators::sum<U>> {
- using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;
- };
- } // namespace std
- #endif
- #endif
|