123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111 |
- // 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>
- #include <type_traits>
- namespace boost {
- namespace histogram {
- namespace accumulators {
- /**
- Uses Neumaier algorithm to compute accurate sums.
- The algorithm uses memory for two floats and is three to
- five times slower compared to a simple floating point
- number used 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 <typename RealType>
- class sum {
- public:
- sum() = default;
- /// Initialize sum to value
- explicit sum(const RealType& value) noexcept : large_(value) {}
- /// Set sum to value
- sum& operator=(const RealType& value) noexcept {
- large_ = value;
- small_ = 0;
- return *this;
- }
- /// Increment sum by one
- sum& operator++() { return operator+=(1); }
- /// Increment sum by value
- sum& operator+=(const RealType& value) {
- auto temp = large_ + value; // prevent optimization
- if (std::abs(large_) >= std::abs(value))
- small_ += (large_ - temp) + value;
- else
- small_ += (value - temp) + large_;
- large_ = temp;
- return *this;
- }
- /// Scale by value
- sum& operator*=(const RealType& value) {
- large_ *= value;
- small_ *= value;
- return *this;
- }
- template <class T>
- bool operator==(const sum<T>& rhs) const noexcept {
- return large_ == rhs.large_ && small_ == rhs.small_;
- }
- template <class T>
- bool operator!=(const T& rhs) const noexcept {
- return !operator==(rhs);
- }
- /// Return large part of the sum.
- const RealType& large() const { return large_; }
- /// Return small part of the sum.
- const RealType& small() const { return small_; }
- // allow implicit conversion to RealType
- operator RealType() const { return large_ + small_; }
- template <class Archive>
- void serialize(Archive& ar, unsigned /* version */) {
- ar& make_nvp("large", large_);
- ar& make_nvp("small", small_);
- }
- private:
- RealType large_ = RealType();
- RealType small_ = RealType();
- };
- } // 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
|