histogram/accumulators/sum.hpp

0001 // Copyright 2018 Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
0008 #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
0009
0010 #include <boost/core/nvp.hpp>
0011 #include <boost/histogram/fwd.hpp> // for sum<>
0012 #include <cmath>                   // std::abs
0013 #include <type_traits>             // std::is_floating_point, std::common_type
0014
0015 namespace boost {
0016 namespace histogram {
0017 namespace accumulators {
0018
0019 /**
0020   Uses Neumaier algorithm to compute accurate sums of floats.
0021
0022   The algorithm is an improved Kahan algorithm
0023   (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). The algorithm uses memory for
0024   two numbers and is three to five times slower compared to using a single number to
0025   accumulate a sum, but the relative error of the sum is at the level of the machine
0026   precision, independent of the number of samples.
0027
0028   A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51.
0029 */
0030 template <class ValueType>
0031 class sum {
0032   static_assert(std::is_floating_point<ValueType>::value,
0033                 "ValueType must be a floating point type");
0034
0035 public:
0036   using value_type = ValueType;
0037   using const_reference = const value_type&;
0038
0039   sum() = default;
0040
0041   /// Initialize sum to value and allow implicit conversion
0042   sum(const_reference value) noexcept : sum(value, 0) {}
0043
0044   /// Allow implicit conversion from sum<T>
0045   template <class T>
0046   sum(const sum<T>& s) noexcept : sum(s.large_part(), s.small_part()) {}
0047
0048   /// Initialize sum explicitly with large and small parts
0049   sum(const_reference large_part, const_reference small_part) noexcept
0050       : large_(large_part), small_(small_part) {}
0051
0052   /// Increment sum by one
0053   sum& operator++() noexcept { return operator+=(1); }
0054
0055   /// Increment sum by value
0056   sum& operator+=(const_reference value) noexcept {
0057     // prevent compiler optimization from destroying the algorithm
0058     // when -ffast-math is enabled
0059     volatile value_type l;
0060     value_type s;
0061     if (std::abs(large_) >= std::abs(value)) {
0062       l = large_;
0063       s = value;
0064     } else {
0065       l = value;
0066       s = large_;
0067     }
0068     large_ += value;
0069     l = l - large_;
0070     l = l + s;
0071     small_ += l;
0072     return *this;
0073   }
0074
0075   /// Add another sum
0076   sum& operator+=(const sum& s) noexcept {
0077     operator+=(s.large_);
0078     small_ += s.small_;
0079     return *this;
0080   }
0081
0082   /// Scale by value
0083   sum& operator*=(const_reference value) noexcept {
0084     large_ *= value;
0085     small_ *= value;
0086     return *this;
0087   }
0088
0089   bool operator==(const sum& rhs) const noexcept {
0090     return large_ + small_ == rhs.large_ + rhs.small_;
0091   }
0092
0093   bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); }
0094
0095   /// Return value of the sum.
0096   value_type value() const noexcept { return large_ + small_; }
0097
0098   /// Return large part of the sum.
0099   const_reference large_part() const noexcept { return large_; }
0100
0101   /// Return small part of the sum.
0102   const_reference small_part() const noexcept { return small_; }
0103   // note: windows.h illegially uses `#define small char`, cannot use method "small"
0104
0105   // lossy conversion to value type must be explicit
0106   explicit operator value_type() const noexcept { return value(); }
0107
0108   template <class Archive>
0109   void serialize(Archive& ar, unsigned /* version */) {
0110     ar& make_nvp("large", large_);
0111     ar& make_nvp("small", small_);
0112   }
0113
0114   // begin: extra operators to make sum behave like a regular number
0115
0116   sum& operator*=(const sum& rhs) noexcept {
0117     const auto scale = static_cast<value_type>(rhs);
0118     large_ *= scale;
0119     small_ *= scale;
0120     return *this;
0121   }
0122
0123   sum operator*(const sum& rhs) const noexcept {
0124     sum s = *this;
0125     s *= rhs;
0126     return s;
0127   }
0128
0129   sum& operator/=(const sum& rhs) noexcept {
0130     const auto scale = 1.0 / static_cast<value_type>(rhs);
0131     large_ *= scale;
0132     small_ *= scale;
0133     return *this;
0134   }
0135
0136   sum operator/(const sum& rhs) const noexcept {
0137     sum s = *this;
0138     s /= rhs;
0139     return s;
0140   }
0141
0142   bool operator<(const sum& rhs) const noexcept {
0143     return operator value_type() < rhs.operator value_type();
0144   }
0145
0146   bool operator>(const sum& rhs) const noexcept {
0147     return operator value_type() > rhs.operator value_type();
0148   }
0149
0150   bool operator<=(const sum& rhs) const noexcept {
0151     return operator value_type() <= rhs.operator value_type();
0152   }
0153
0154   bool operator>=(const sum& rhs) const noexcept {
0155     return operator value_type() >= rhs.operator value_type();
0156   }
0157
0158   // end: extra operators
0159
0160 private:
0161   value_type large_{};
0162   value_type small_{};
0163 };
0164
0165 } // namespace accumulators
0166 } // namespace histogram
0167 } // namespace boost
0168
0169 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
0170 namespace std {
0171 template <class T, class U>
0172 struct common_type<boost::histogram::accumulators::sum<T>,
0173                    boost::histogram::accumulators::sum<U>> {
0174   using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;
0175 };
0176 } // namespace std
0177 #endif
0178
0179 #endif