sum.hpp 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178
  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> // std::abs
  11. #include <type_traits> // std::is_floating_point, std::common_type
  12. namespace boost {
  13. namespace histogram {
  14. namespace accumulators {
  15. /**
  16. Uses Neumaier algorithm to compute accurate sums of floats.
  17. The algorithm is an improved Kahan algorithm
  18. (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). The algorithm uses memory for
  19. two numbers and is three to five times slower compared to using a single number to
  20. accumulate a sum, but the relative error of the sum is at the level of the machine
  21. precision, independent of the number of samples.
  22. A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51.
  23. */
  24. template <class ValueType>
  25. class sum {
  26. static_assert(std::is_floating_point<ValueType>::value,
  27. "ValueType must be a floating point type");
  28. public:
  29. using value_type = ValueType;
  30. using const_reference = const value_type&;
  31. sum() = default;
  32. /// Initialize sum to value and allow implicit conversion
  33. sum(const_reference value) noexcept : sum(value, 0) {}
  34. /// Allow implicit conversion from sum<T>
  35. template <class T>
  36. sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {}
  37. /// Initialize sum explicitly with large and small parts
  38. sum(const_reference large, const_reference small) noexcept
  39. : large_(large), small_(small) {}
  40. /// Increment sum by one
  41. sum& operator++() noexcept { return operator+=(1); }
  42. /// Increment sum by value
  43. sum& operator+=(const_reference value) noexcept {
  44. // prevent compiler optimization from destroying the algorithm
  45. // when -ffast-math is enabled
  46. volatile value_type l;
  47. value_type s;
  48. if (std::abs(large_) >= std::abs(value)) {
  49. l = large_;
  50. s = value;
  51. } else {
  52. l = value;
  53. s = large_;
  54. }
  55. large_ += value;
  56. l = l - large_;
  57. l = l + s;
  58. small_ += l;
  59. return *this;
  60. }
  61. /// Add another sum
  62. sum& operator+=(const sum& s) noexcept {
  63. operator+=(s.large_);
  64. small_ += s.small_;
  65. return *this;
  66. }
  67. /// Scale by value
  68. sum& operator*=(const_reference value) noexcept {
  69. large_ *= value;
  70. small_ *= value;
  71. return *this;
  72. }
  73. bool operator==(const sum& rhs) const noexcept {
  74. return large_ + small_ == rhs.large_ + rhs.small_;
  75. }
  76. bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); }
  77. /// Return value of the sum.
  78. value_type value() const noexcept { return large_ + small_; }
  79. /// Return large part of the sum.
  80. const_reference large() const noexcept { return large_; }
  81. /// Return small part of the sum.
  82. const_reference small() const noexcept { return small_; }
  83. // lossy conversion to value type must be explicit
  84. explicit operator value_type() const noexcept { return value(); }
  85. template <class Archive>
  86. void serialize(Archive& ar, unsigned /* version */) {
  87. ar& make_nvp("large", large_);
  88. ar& make_nvp("small", small_);
  89. }
  90. // begin: extra operators to make sum behave like a regular number
  91. sum& operator*=(const sum& rhs) noexcept {
  92. const auto scale = static_cast<value_type>(rhs);
  93. large_ *= scale;
  94. small_ *= scale;
  95. return *this;
  96. }
  97. sum operator*(const sum& rhs) const noexcept {
  98. sum s = *this;
  99. s *= rhs;
  100. return s;
  101. }
  102. sum& operator/=(const sum& rhs) noexcept {
  103. const auto scale = 1.0 / static_cast<value_type>(rhs);
  104. large_ *= scale;
  105. small_ *= scale;
  106. return *this;
  107. }
  108. sum operator/(const sum& rhs) const noexcept {
  109. sum s = *this;
  110. s /= rhs;
  111. return s;
  112. }
  113. bool operator<(const sum& rhs) const noexcept {
  114. return operator value_type() < rhs.operator value_type();
  115. }
  116. bool operator>(const sum& rhs) const noexcept {
  117. return operator value_type() > rhs.operator value_type();
  118. }
  119. bool operator<=(const sum& rhs) const noexcept {
  120. return operator value_type() <= rhs.operator value_type();
  121. }
  122. bool operator>=(const sum& rhs) const noexcept {
  123. return operator value_type() >= rhs.operator value_type();
  124. }
  125. // end: extra operators
  126. private:
  127. value_type large_{};
  128. value_type small_{};
  129. };
  130. } // namespace accumulators
  131. } // namespace histogram
  132. } // namespace boost
  133. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  134. namespace std {
  135. template <class T, class U>
  136. struct common_type<boost::histogram::accumulators::sum<T>,
  137. boost::histogram::accumulators::sum<U>> {
  138. using type = boost::histogram::accumulators::sum<common_type_t<T, U>>;
  139. };
  140. } // namespace std
  141. #endif
  142. #endif