sum.hpp 2.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111
  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. #ifndef BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
  7. #define BOOST_HISTOGRAM_ACCUMULATORS_SUM_HPP
  8. #include <boost/core/nvp.hpp>
  9. #include <boost/histogram/fwd.hpp> // for sum<>
  10. #include <cmath>
  11. #include <type_traits>
  12. namespace boost {
  13. namespace histogram {
  14. namespace accumulators {
  15. /**
  16. Uses Neumaier algorithm to compute accurate sums.
  17. The algorithm uses memory for two floats and is three to
  18. five times slower compared to a simple floating point
  19. number used to accumulate a sum, but the relative error
  20. of the sum is at the level of the machine precision,
  21. independent of the number of samples.
  22. A. Neumaier, Zeitschrift fuer Angewandte Mathematik
  23. und Mechanik 54 (1974) 39-51.
  24. */
  25. template <typename RealType>
  26. class sum {
  27. public:
  28. sum() = default;
  29. /// Initialize sum to value
  30. explicit sum(const RealType& value) noexcept : large_(value) {}
  31. /// Set sum to value
  32. sum& operator=(const RealType& value) noexcept {
  33. large_ = value;
  34. small_ = 0;
  35. return *this;
  36. }
  37. /// Increment sum by one
  38. sum& operator++() { return operator+=(1); }
  39. /// Increment sum by value
  40. sum& operator+=(const RealType& value) {
  41. auto temp = large_ + value; // prevent optimization
  42. if (std::abs(large_) >= std::abs(value))
  43. small_ += (large_ - temp) + value;
  44. else
  45. small_ += (value - temp) + large_;
  46. large_ = temp;
  47. return *this;
  48. }
  49. /// Scale by value
  50. sum& operator*=(const RealType& value) {
  51. large_ *= value;
  52. small_ *= value;
  53. return *this;
  54. }
  55. template <class T>
  56. bool operator==(const sum<T>& rhs) const noexcept {
  57. return large_ == rhs.large_ && small_ == rhs.small_;
  58. }
  59. template <class T>
  60. bool operator!=(const T& rhs) const noexcept {
  61. return !operator==(rhs);
  62. }
  63. /// Return large part of the sum.
  64. const RealType& large() const { return large_; }
  65. /// Return small part of the sum.
  66. const RealType& small() const { return small_; }
  67. // allow implicit conversion to RealType
  68. operator RealType() const { return large_ + small_; }
  69. template <class Archive>
  70. void serialize(Archive& ar, unsigned /* version */) {
  71. ar& make_nvp("large", large_);
  72. ar& make_nvp("small", small_);
  73. }
  74. private:
  75. RealType large_ = RealType();
  76. RealType small_ = RealType();
  77. };
  78. } // namespace accumulators
  79. } // namespace histogram
  80. } // namespace boost
  81. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  82. namespace std {
  83. template <class T, class U>
  84. struct common_type<boost::histogram::accumulators::sum<T>,
  85. boost::histogram::accumulators::sum<U>> {
  86. using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;
  87. };
  88. } // namespace std
  89. #endif
  90. #endif