1 // Copyright 2015-2019 Hans Dembinski
2 // Copyright 2019 Glen Joseph Fernandes (glenjofe@gmail.com)
4 // Distributed under the Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
8 #ifndef BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
9 #define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
12 #include <boost/assert.hpp>
13 #include <boost/config.hpp>
14 #include <boost/core/alloc_construct.hpp>
15 #include <boost/core/exchange.hpp>
16 #include <boost/core/nvp.hpp>
17 #include <boost/histogram/detail/array_wrapper.hpp>
18 #include <boost/histogram/detail/iterator_adaptor.hpp>
19 #include <boost/histogram/detail/large_int.hpp>
20 #include <boost/histogram/detail/operators.hpp>
21 #include <boost/histogram/detail/safe_comparison.hpp>
22 #include <boost/histogram/fwd.hpp>
23 #include <boost/mp11/algorithm.hpp>
24 #include <boost/mp11/list.hpp>
25 #include <boost/mp11/utility.hpp>
31 #include <type_traits>
38 struct is_large_int : std::false_type {};
41 struct is_large_int<large_int<A>> : std::true_type {};
43 template <class T, class ReturnType>
44 using if_arithmetic_or_large_int =
45 std::enable_if_t<(std::is_arithmetic<T>::value || is_large_int<T>::value),
48 template <class L, class T>
49 using next_type = mp11::mp_at_c<L, (mp11::mp_find<L, T>::value + 1)>;
51 template <class Allocator>
52 class construct_guard {
54 using pointer = typename std::allocator_traits<Allocator>::pointer;
56 construct_guard(Allocator& a, pointer p, std::size_t n) noexcept
57 : a_(a), p_(p), n_(n) {}
60 if (p_) { a_.deallocate(p_, n_); }
63 void release() { p_ = pointer(); }
65 construct_guard(const construct_guard&) = delete;
66 construct_guard& operator=(const construct_guard&) = delete;
74 template <class Allocator>
75 void* buffer_create(Allocator& a, std::size_t n) {
76 auto ptr = a.allocate(n); // may throw
77 static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
78 "ptr must be trivially copyable");
79 construct_guard<Allocator> guard(a, ptr, n);
80 boost::alloc_construct_n(a, ptr, n);
82 return static_cast<void*>(ptr);
85 template <class Allocator, class Iterator>
86 auto buffer_create(Allocator& a, std::size_t n, Iterator iter) {
88 auto ptr = a.allocate(n); // may throw
89 static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
90 "ptr must be trivially copyable");
91 construct_guard<Allocator> guard(a, ptr, n);
92 using T = typename std::allocator_traits<Allocator>::value_type;
93 struct casting_iterator {
94 void operator++() noexcept { ++iter_; }
95 T operator*() noexcept {
96 return static_cast<T>(*iter_);
97 } // silence conversion warnings
100 boost::alloc_construct_n(a, ptr, n, casting_iterator{iter});
105 template <class Allocator>
106 void buffer_destroy(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
109 BOOST_ASSERT(n > 0u);
110 boost::alloc_destroy_n(a, p, n);
114 } // namespace detail
117 Memory-efficient storage for integral counters which cannot overflow.
119 This storage provides a no-overflow-guarantee if the counters are incremented with
120 integer weights. It maintains a contiguous array of elemental counters, one for each
121 cell. If an operation is requested which would overflow a counter, the array is
122 replaced with another of a wider integral type, then the operation is executed. The
123 storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
124 integral type, similar to those in
125 [Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
127 A scaling operation or adding a floating point number triggers a conversion of the
128 elemental counters into doubles, which voids the no-overflow-guarantee.
130 template <class Allocator>
131 class unlimited_storage {
133 std::is_same<typename std::allocator_traits<Allocator>::pointer,
134 typename std::allocator_traits<Allocator>::value_type*>::value,
135 "unlimited_storage requires allocator with trivial pointer type");
136 using U8 = std::uint8_t;
137 using U16 = std::uint16_t;
138 using U32 = std::uint32_t;
139 using U64 = std::uint64_t;
142 static constexpr bool has_threading_support = false;
144 using allocator_type = Allocator;
145 using value_type = double;
146 using large_int = detail::large_int<
147 typename std::allocator_traits<allocator_type>::template rebind_alloc<U64>>;
150 // cannot be moved outside of scope of unlimited_storage, large_int is dependent type
151 using types = mp11::mp_list<U8, U16, U32, U64, large_int, double>;
154 static constexpr unsigned type_index() noexcept {
155 return static_cast<unsigned>(mp11::mp_find<types, T>::value);
158 template <class F, class... Ts>
159 decltype(auto) visit(F&& f, Ts&&... ts) const {
160 // this is intentionally not a switch, the if-chain is faster in benchmarks
161 if (type == type_index<U8>())
162 return f(static_cast<U8*>(ptr), std::forward<Ts>(ts)...);
163 if (type == type_index<U16>())
164 return f(static_cast<U16*>(ptr), std::forward<Ts>(ts)...);
165 if (type == type_index<U32>())
166 return f(static_cast<U32*>(ptr), std::forward<Ts>(ts)...);
167 if (type == type_index<U64>())
168 return f(static_cast<U64*>(ptr), std::forward<Ts>(ts)...);
169 if (type == type_index<large_int>())
170 return f(static_cast<large_int*>(ptr), std::forward<Ts>(ts)...);
171 return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
174 buffer_type(const allocator_type& a = {}) : alloc(a) {}
176 buffer_type(buffer_type&& o) noexcept
177 : alloc(std::move(o.alloc))
178 , size(boost::exchange(o.size, 0))
179 , type(boost::exchange(o.type, 0))
180 , ptr(boost::exchange(o.ptr, nullptr)) {}
182 buffer_type& operator=(buffer_type&& o) noexcept {
184 swap(alloc, o.alloc);
191 buffer_type(const buffer_type& x) : alloc(x.alloc) {
192 x.visit([this, n = x.size](const auto* xp) {
193 using T = std::decay_t<decltype(*xp)>;
194 this->template make<T>(n, xp);
198 buffer_type& operator=(const buffer_type& o) {
199 *this = buffer_type(o);
203 ~buffer_type() noexcept { destroy(); }
205 void destroy() noexcept {
206 BOOST_ASSERT((ptr == nullptr) == (size == 0));
207 if (ptr == nullptr) return;
208 visit([this](auto* p) {
209 using T = std::decay_t<decltype(*p)>;
211 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
212 alloc_type a(alloc); // rebind allocator
213 detail::buffer_destroy(a, p, this->size);
221 void make(std::size_t n) {
222 // note: order of commands is to not leave buffer in invalid state upon throw
227 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
229 ptr = detail::buffer_create(a, n); // may throw
232 type = type_index<T>();
235 template <class T, class U>
236 void make(std::size_t n, U iter) {
237 // note: iter may be current ptr, so create new buffer before deleting old buffer
238 void* new_ptr = nullptr;
239 const auto new_type = type_index<T>();
243 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
245 new_ptr = detail::buffer_create(a, n, iter); // may throw
253 allocator_type alloc;
254 std::size_t size = 0;
256 mutable void* ptr = nullptr;
259 class reference; // forward declare to make friend of const_reference
261 /// implementation detail
262 class const_reference
263 : detail::partially_ordered<const_reference, const_reference, void> {
265 const_reference(buffer_type& b, std::size_t i) noexcept : bref_(b), idx_(i) {
266 BOOST_ASSERT(idx_ < bref_.size);
269 const_reference(const const_reference&) noexcept = default;
271 // no assignment for const_references
272 const_reference& operator=(const const_reference&) = delete;
273 const_reference& operator=(const_reference&&) = delete;
275 operator double() const noexcept {
277 [this](const auto* p) { return static_cast<double>(p[this->idx_]); });
280 bool operator<(const const_reference& o) const noexcept {
281 return apply_binary<detail::safe_less>(o);
284 bool operator==(const const_reference& o) const noexcept {
285 return apply_binary<detail::safe_equal>(o);
289 detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
290 return apply_binary<detail::safe_less>(o);
294 detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
295 return apply_binary<detail::safe_greater>(o);
299 detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
300 return apply_binary<detail::safe_equal>(o);
304 template <class Binary>
305 bool apply_binary(const const_reference& x) const noexcept {
306 return x.bref_.visit([this, ix = x.idx_](const auto* xp) {
307 return this->apply_binary<Binary>(xp[ix]);
311 template <class Binary, class U>
312 bool apply_binary(const U& x) const noexcept {
313 return bref_.visit([i = idx_, &x](const auto* p) { return Binary()(p[i], x); });
319 friend class reference;
322 /// implementation detail
323 class reference : public const_reference,
324 public detail::partially_ordered<reference, reference, void> {
326 reference(buffer_type& b, std::size_t i) noexcept : const_reference(b, i) {}
328 // references do copy-construct
329 reference(const reference& x) noexcept = default;
331 // references do not rebind, assign through
332 reference& operator=(const reference& x) {
333 return operator=(static_cast<const_reference>(x));
336 // references do not rebind, assign through
337 reference& operator=(const const_reference& x) {
338 // safe for self-assignment, assigning matching type doesn't invalide buffer
339 x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator=(xp[ix]); });
344 detail::if_arithmetic_or_large_int<U, reference&> operator=(const U& x) {
345 this->bref_.visit([this, &x](auto* p) {
346 // gcc-8 optimizes the expression `p[this->idx_] = 0` away even at -O0,
347 // so we merge it into the next line which is properly counted
348 adder()((p[this->idx_] = 0, p), this->bref_, this->idx_, x);
353 bool operator<(const reference& o) const noexcept {
354 return const_reference::operator<(o);
357 bool operator==(const reference& o) const noexcept {
358 return const_reference::operator==(o);
362 detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
363 return const_reference::operator<(o);
367 detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
368 return const_reference::operator>(o);
372 detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
373 return const_reference::operator==(o);
376 reference& operator+=(const const_reference& x) {
377 x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator+=(xp[ix]); });
382 detail::if_arithmetic_or_large_int<U, reference&> operator+=(const U& x) {
383 this->bref_.visit(adder(), this->bref_, this->idx_, x);
387 reference& operator-=(const double x) { return operator+=(-x); }
389 reference& operator*=(const double x) {
390 this->bref_.visit(multiplier(), this->bref_, this->idx_, x);
394 reference& operator/=(const double x) { return operator*=(1.0 / x); }
396 reference& operator++() {
397 this->bref_.visit(incrementor(), this->bref_, this->idx_);
403 template <class Value, class Reference>
404 class iterator_impl : public detail::iterator_adaptor<iterator_impl<Value, Reference>,
405 std::size_t, Reference, Value> {
407 iterator_impl() = default;
408 template <class V, class R>
409 iterator_impl(const iterator_impl<V, R>& it)
410 : iterator_impl::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
411 iterator_impl(buffer_type* b, std::size_t i) noexcept
412 : iterator_impl::iterator_adaptor_(i), buffer_(b) {}
414 Reference operator*() const noexcept { return {*buffer_, this->base()}; }
416 template <class V, class R>
417 friend class iterator_impl;
420 mutable buffer_type* buffer_ = nullptr;
424 using const_iterator = iterator_impl<const value_type, const_reference>;
425 using iterator = iterator_impl<value_type, reference>;
427 explicit unlimited_storage(const allocator_type& a = {}) : buffer_(a) {}
428 unlimited_storage(const unlimited_storage&) = default;
429 unlimited_storage& operator=(const unlimited_storage&) = default;
430 unlimited_storage(unlimited_storage&&) = default;
431 unlimited_storage& operator=(unlimited_storage&&) = default;
434 // template <class Allocator>
435 // unlimited_storage(const unlimited_storage<Allocator>& s)
437 template <class Iterable, class = detail::requires_iterable<Iterable>>
438 explicit unlimited_storage(const Iterable& s) {
441 auto s_begin = begin(s);
443 using V = typename std::iterator_traits<decltype(begin(s))>::value_type;
444 constexpr auto ti = buffer_type::template type_index<V>();
445 constexpr auto nt = mp11::mp_size<typename buffer_type::types>::value;
446 const std::size_t size = static_cast<std::size_t>(std::distance(s_begin, s_end));
447 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
453 buffer_.template make<V>(size, s_begin);
455 buffer_.template make<double>(size, s_begin);
458 template <class Iterable, class = detail::requires_iterable<Iterable>>
459 unlimited_storage& operator=(const Iterable& s) {
460 *this = unlimited_storage(s);
464 allocator_type get_allocator() const { return buffer_.alloc; }
466 void reset(std::size_t n) { buffer_.template make<U8>(n); }
468 std::size_t size() const noexcept { return buffer_.size; }
470 reference operator[](std::size_t i) noexcept { return {buffer_, i}; }
471 const_reference operator[](std::size_t i) const noexcept { return {buffer_, i}; }
473 bool operator==(const unlimited_storage& x) const noexcept {
474 if (size() != x.size()) return false;
475 return buffer_.visit([&x](const auto* p) {
476 return x.buffer_.visit([p, n = x.size()](const auto* xp) {
477 return std::equal(p, p + n, xp, detail::safe_equal{});
482 template <class Iterable>
483 bool operator==(const Iterable& iterable) const {
484 if (size() != iterable.size()) return false;
485 return buffer_.visit([&iterable](const auto* p) {
486 return std::equal(p, p + iterable.size(), std::begin(iterable),
487 detail::safe_equal{});
491 unlimited_storage& operator*=(const double x) {
492 buffer_.visit(multiplier(), buffer_, x);
496 iterator begin() noexcept { return {&buffer_, 0}; }
497 iterator end() noexcept { return {&buffer_, size()}; }
498 const_iterator begin() const noexcept { return {&buffer_, 0}; }
499 const_iterator end() const noexcept { return {&buffer_, size()}; }
501 /// implementation detail; used by unit tests, not part of generic storage interface
503 unlimited_storage(std::size_t s, const T* p, const allocator_type& a = {})
504 : buffer_(std::move(a)) {
505 buffer_.template make<T>(s, p);
508 template <class Archive>
509 void serialize(Archive& ar, unsigned /* version */) {
510 if (Archive::is_loading::value) {
511 buffer_type tmp(buffer_.alloc);
513 ar& make_nvp("type", tmp.type);
514 ar& make_nvp("size", size);
515 tmp.visit([this, size](auto* tp) {
516 BOOST_ASSERT(tp == nullptr);
517 using T = std::decay_t<decltype(*tp)>;
518 buffer_.template make<T>(size);
521 ar& make_nvp("type", buffer_.type);
522 ar& make_nvp("size", buffer_.size);
524 buffer_.visit([this, &ar](auto* tp) {
525 auto w = detail::make_array_wrapper(tp, this->buffer_.size);
526 ar& make_nvp("buffer", w);
533 void operator()(T* tp, buffer_type& b, std::size_t i) {
534 BOOST_ASSERT(tp && i < b.size);
535 if (!detail::safe_increment(tp[i])) {
536 using U = detail::next_type<typename buffer_type::types, T>;
537 b.template make<U>(b.size, tp);
538 ++static_cast<U*>(b.ptr)[i];
542 void operator()(large_int* tp, buffer_type&, std::size_t i) { ++tp[i]; }
544 void operator()(double* tp, buffer_type&, std::size_t i) { ++tp[i]; }
549 void operator()(double* tp, buffer_type&, std::size_t i, const U& x) {
550 tp[i] += static_cast<double>(x);
553 void operator()(large_int* tp, buffer_type&, std::size_t i, const large_int& x) {
554 tp[i] += x; // potentially adding large_int to itself is safe
557 template <class T, class U>
558 void operator()(T* tp, buffer_type& b, std::size_t i, const U& x) {
559 is_x_integral(std::is_integral<U>{}, tp, b, i, x);
562 template <class T, class U>
563 void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
565 // x could be reference to buffer we manipulate, make copy before changing buffer
566 const auto v = static_cast<double>(x);
567 b.template make<double>(b.size, tp);
568 operator()(static_cast<double*>(b.ptr), b, i, v);
572 void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
573 const large_int& x) {
574 // x could be reference to buffer we manipulate, make copy before changing buffer
575 const auto v = static_cast<large_int>(x);
576 b.template make<large_int>(b.size, tp);
577 operator()(static_cast<large_int*>(b.ptr), b, i, v);
580 template <class T, class U>
581 void is_x_integral(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
582 is_x_unsigned(std::is_unsigned<U>{}, tp, b, i, x);
585 template <class T, class U>
586 void is_x_unsigned(std::false_type, T* tp, buffer_type& b, std::size_t i,
589 is_x_unsigned(std::true_type{}, tp, b, i, detail::make_unsigned(x));
591 is_x_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
594 template <class T, class U>
595 void is_x_unsigned(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
596 if (detail::safe_radd(tp[i], x)) return;
597 // x could be reference to buffer we manipulate, need to convert to value
599 using TN = detail::next_type<typename buffer_type::types, T>;
600 b.template make<TN>(b.size, tp);
601 is_x_unsigned(std::true_type{}, static_cast<TN*>(b.ptr), b, i, y);
605 void is_x_unsigned(std::true_type, large_int* tp, buffer_type&, std::size_t i,
613 void operator()(T* tp, buffer_type& b, const double x) {
614 // potential lossy conversion that cannot be avoided
615 b.template make<double>(b.size, tp);
616 operator()(static_cast<double*>(b.ptr), b, x);
619 void operator()(double* tp, buffer_type& b, const double x) {
620 for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
624 void operator()(T* tp, buffer_type& b, std::size_t i, const double x) {
625 b.template make<double>(b.size, tp);
626 operator()(static_cast<double*>(b.ptr), b, i, x);
629 void operator()(double* tp, buffer_type&, std::size_t i, const double x) {
630 tp[i] *= static_cast<double>(x);
634 mutable buffer_type buffer_;
635 friend struct unsafe_access;
638 } // namespace histogram