File indexing completed on 2025-12-16 09:53:11
0001
0002
0003
0004
0005
0006
0007
0008 #ifndef BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
0009 #define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
0010
0011 #include <algorithm>
0012 #include <boost/core/alloc_construct.hpp>
0013 #include <boost/core/exchange.hpp>
0014 #include <boost/core/nvp.hpp>
0015 #include <boost/histogram/detail/array_wrapper.hpp>
0016 #include <boost/histogram/detail/iterator_adaptor.hpp>
0017 #include <boost/histogram/detail/large_int.hpp>
0018 #include <boost/histogram/detail/operators.hpp>
0019 #include <boost/histogram/detail/safe_comparison.hpp>
0020 #include <boost/histogram/fwd.hpp>
0021 #include <boost/mp11/algorithm.hpp>
0022 #include <boost/mp11/list.hpp>
0023 #include <boost/mp11/utility.hpp>
0024 #include <cassert>
0025 #include <cmath>
0026 #include <cstdint>
0027 #include <functional>
0028 #include <iterator>
0029 #include <memory>
0030 #include <type_traits>
0031
0032 namespace boost {
0033 namespace histogram {
0034 namespace detail {
0035
0036 template <class T>
0037 struct is_large_int : std::false_type {};
0038
0039 template <class A>
0040 struct is_large_int<large_int<A>> : std::true_type {};
0041
0042 template <class T, class ReturnType>
0043 using if_arithmetic_or_large_int =
0044 std::enable_if_t<(std::is_arithmetic<T>::value || is_large_int<T>::value),
0045 ReturnType>;
0046
0047 template <class L, class T>
0048 using next_type = mp11::mp_at_c<L, (mp11::mp_find<L, T>::value + 1)>;
0049
0050 template <class Allocator>
0051 class construct_guard {
0052 public:
0053 using pointer = typename std::allocator_traits<Allocator>::pointer;
0054
0055 construct_guard(Allocator& a, pointer p, std::size_t n) noexcept
0056 : a_(a), p_(p), n_(n) {}
0057
0058 ~construct_guard() {
0059 if (p_) { a_.deallocate(p_, n_); }
0060 }
0061
0062 void release() { p_ = pointer(); }
0063
0064 construct_guard(const construct_guard&) = delete;
0065 construct_guard& operator=(const construct_guard&) = delete;
0066
0067 private:
0068 Allocator& a_;
0069 pointer p_;
0070 std::size_t n_;
0071 };
0072
0073 template <class Allocator>
0074 void* buffer_create(Allocator& a, std::size_t n) {
0075 auto ptr = a.allocate(n);
0076 static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
0077 "ptr must be trivially copyable");
0078 construct_guard<Allocator> guard(a, ptr, n);
0079 boost::alloc_construct_n(a, ptr, n);
0080 guard.release();
0081 return static_cast<void*>(ptr);
0082 }
0083
0084 template <class Allocator, class Iterator>
0085 auto buffer_create(Allocator& a, std::size_t n, Iterator iter) {
0086 assert(n > 0u);
0087 auto ptr = a.allocate(n);
0088 static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
0089 "ptr must be trivially copyable");
0090 construct_guard<Allocator> guard(a, ptr, n);
0091 using T = typename std::allocator_traits<Allocator>::value_type;
0092 struct casting_iterator {
0093 void operator++() noexcept { ++iter_; }
0094 T operator*() noexcept {
0095 return static_cast<T>(*iter_);
0096 }
0097 Iterator iter_;
0098 };
0099 boost::alloc_construct_n(a, ptr, n, casting_iterator{iter});
0100 guard.release();
0101 return ptr;
0102 }
0103
0104 template <class Allocator>
0105 void buffer_destroy(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
0106 std::size_t n) {
0107 assert(p);
0108 assert(n > 0u);
0109 boost::alloc_destroy_n(a, p, n);
0110 a.deallocate(p, n);
0111 }
0112
0113 }
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129 template <class Allocator>
0130 class unlimited_storage {
0131 static_assert(
0132 std::is_same<typename std::allocator_traits<Allocator>::pointer,
0133 typename std::allocator_traits<Allocator>::value_type*>::value,
0134 "unlimited_storage requires allocator with trivial pointer type");
0135 using U8 = std::uint8_t;
0136 using U16 = std::uint16_t;
0137 using U32 = std::uint32_t;
0138 using U64 = std::uint64_t;
0139
0140 public:
0141 static constexpr bool has_threading_support = false;
0142
0143 using allocator_type = Allocator;
0144 using value_type = double;
0145 using large_int = detail::large_int<
0146 typename std::allocator_traits<allocator_type>::template rebind_alloc<U64>>;
0147
0148 struct buffer_type {
0149
0150 using types = mp11::mp_list<U8, U16, U32, U64, large_int, double>;
0151
0152 template <class T>
0153 static constexpr unsigned type_index() noexcept {
0154 return static_cast<unsigned>(mp11::mp_find<types, T>::value);
0155 }
0156
0157 template <class F, class... Ts>
0158 decltype(auto) visit(F&& f, Ts&&... ts) const {
0159
0160 if (type == type_index<U8>())
0161 return f(static_cast<U8*>(ptr), std::forward<Ts>(ts)...);
0162 if (type == type_index<U16>())
0163 return f(static_cast<U16*>(ptr), std::forward<Ts>(ts)...);
0164 if (type == type_index<U32>())
0165 return f(static_cast<U32*>(ptr), std::forward<Ts>(ts)...);
0166 if (type == type_index<U64>())
0167 return f(static_cast<U64*>(ptr), std::forward<Ts>(ts)...);
0168 if (type == type_index<large_int>())
0169 return f(static_cast<large_int*>(ptr), std::forward<Ts>(ts)...);
0170 return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
0171 }
0172
0173 buffer_type(const allocator_type& a = {}) : alloc(a) {}
0174
0175 buffer_type(buffer_type&& o) noexcept
0176 : alloc(std::move(o.alloc))
0177 , size(boost::exchange(o.size, 0))
0178 , type(boost::exchange(o.type, 0))
0179 , ptr(boost::exchange(o.ptr, nullptr)) {}
0180
0181 buffer_type& operator=(buffer_type&& o) noexcept {
0182 using std::swap;
0183 swap(alloc, o.alloc);
0184 swap(size, o.size);
0185 swap(type, o.type);
0186 swap(ptr, o.ptr);
0187 return *this;
0188 }
0189
0190 buffer_type(const buffer_type& x) : alloc(x.alloc) {
0191 x.visit([this, n = x.size](const auto* xp) {
0192 using T = std::decay_t<decltype(*xp)>;
0193 this->template make<T>(n, xp);
0194 });
0195 }
0196
0197 buffer_type& operator=(const buffer_type& o) {
0198 *this = buffer_type(o);
0199 return *this;
0200 }
0201
0202 ~buffer_type() noexcept { destroy(); }
0203
0204 void destroy() noexcept {
0205 assert((ptr == nullptr) == (size == 0));
0206 if (ptr == nullptr) return;
0207 visit([this](auto* p) {
0208 using T = std::decay_t<decltype(*p)>;
0209 using alloc_type =
0210 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0211 alloc_type a(alloc);
0212 detail::buffer_destroy(a, p, this->size);
0213 });
0214 size = 0;
0215 type = 0;
0216 ptr = nullptr;
0217 }
0218
0219 template <class T>
0220 void make(std::size_t n) {
0221
0222 destroy();
0223 if (n > 0) {
0224
0225 using alloc_type =
0226 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0227 alloc_type a(alloc);
0228 ptr = detail::buffer_create(a, n);
0229 }
0230 size = n;
0231 type = type_index<T>();
0232 }
0233
0234 template <class T, class U>
0235 void make(std::size_t n, U iter) {
0236
0237 void* new_ptr = nullptr;
0238 const auto new_type = type_index<T>();
0239 if (n > 0) {
0240
0241 using alloc_type =
0242 typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
0243 alloc_type a(alloc);
0244 new_ptr = detail::buffer_create(a, n, iter);
0245 }
0246 destroy();
0247 size = n;
0248 type = new_type;
0249 ptr = new_ptr;
0250 }
0251
0252 allocator_type alloc;
0253 std::size_t size = 0;
0254 unsigned type = 0;
0255 mutable void* ptr = nullptr;
0256 };
0257
0258 class reference;
0259
0260
0261 class const_reference
0262 : detail::partially_ordered<const_reference, const_reference, void> {
0263 public:
0264 const_reference(buffer_type& b, std::size_t i) noexcept : bref_(b), idx_(i) {
0265 assert(idx_ < bref_.size);
0266 }
0267
0268 const_reference(const const_reference&) noexcept = default;
0269
0270
0271 const_reference& operator=(const const_reference&) = delete;
0272 const_reference& operator=(const_reference&&) = delete;
0273
0274 operator double() const noexcept {
0275 return bref_.visit(
0276 [this](const auto* p) { return static_cast<double>(p[this->idx_]); });
0277 }
0278
0279 bool operator<(const const_reference& o) const noexcept {
0280 return apply_binary<detail::safe_less>(o);
0281 }
0282
0283 bool operator==(const const_reference& o) const noexcept {
0284 return apply_binary<detail::safe_equal>(o);
0285 }
0286
0287 template <class U>
0288 detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
0289 return apply_binary<detail::safe_less>(o);
0290 }
0291
0292 template <class U>
0293 detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
0294 return apply_binary<detail::safe_greater>(o);
0295 }
0296
0297 template <class U>
0298 detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
0299 return apply_binary<detail::safe_equal>(o);
0300 }
0301
0302 private:
0303 template <class Binary>
0304 bool apply_binary(const const_reference& x) const noexcept {
0305 return x.bref_.visit([this, ix = x.idx_](const auto* xp) {
0306 return this->apply_binary<Binary>(xp[ix]);
0307 });
0308 }
0309
0310 template <class Binary, class U>
0311 bool apply_binary(const U& x) const noexcept {
0312 return bref_.visit([i = idx_, &x](const auto* p) { return Binary()(p[i], x); });
0313 }
0314
0315 protected:
0316 buffer_type& bref_;
0317 std::size_t idx_;
0318 friend class reference;
0319 };
0320
0321
0322 class reference : public const_reference,
0323 public detail::partially_ordered<reference, reference, void> {
0324 public:
0325 reference(buffer_type& b, std::size_t i) noexcept : const_reference(b, i) {}
0326
0327
0328 reference(const reference& x) noexcept = default;
0329
0330
0331 reference& operator=(const reference& x) {
0332 return operator=(static_cast<const_reference>(x));
0333 }
0334
0335
0336 reference& operator=(const const_reference& x) {
0337
0338 x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator=(xp[ix]); });
0339 return *this;
0340 }
0341
0342 template <class U>
0343 detail::if_arithmetic_or_large_int<U, reference&> operator=(const U& x) {
0344 this->bref_.visit([this, &x](auto* p) {
0345
0346
0347 adder()((p[this->idx_] = 0, p), this->bref_, this->idx_, x);
0348 });
0349 return *this;
0350 }
0351
0352 bool operator<(const reference& o) const noexcept {
0353 return const_reference::operator<(o);
0354 }
0355
0356 bool operator==(const reference& o) const noexcept {
0357 return const_reference::operator==(o);
0358 }
0359
0360 template <class U>
0361 detail::if_arithmetic_or_large_int<U, bool> operator<(const U& o) const noexcept {
0362 return const_reference::operator<(o);
0363 }
0364
0365 template <class U>
0366 detail::if_arithmetic_or_large_int<U, bool> operator>(const U& o) const noexcept {
0367 return const_reference::operator>(o);
0368 }
0369
0370 template <class U>
0371 detail::if_arithmetic_or_large_int<U, bool> operator==(const U& o) const noexcept {
0372 return const_reference::operator==(o);
0373 }
0374
0375 reference& operator+=(const const_reference& x) {
0376 x.bref_.visit([this, ix = x.idx_](const auto* xp) { this->operator+=(xp[ix]); });
0377 return *this;
0378 }
0379
0380 template <class U>
0381 detail::if_arithmetic_or_large_int<U, reference&> operator+=(const U& x) {
0382 this->bref_.visit(adder(), this->bref_, this->idx_, x);
0383 return *this;
0384 }
0385
0386 reference& operator-=(const double x) { return operator+=(-x); }
0387
0388 reference& operator*=(const double x) {
0389 this->bref_.visit(multiplier(), this->bref_, this->idx_, x);
0390 return *this;
0391 }
0392
0393 reference& operator/=(const double x) { return operator*=(1.0 / x); }
0394
0395 reference& operator++() {
0396 this->bref_.visit(incrementor(), this->bref_, this->idx_);
0397 return *this;
0398 }
0399 };
0400
0401 private:
0402 template <class Value, class Reference>
0403 class iterator_impl : public detail::iterator_adaptor<iterator_impl<Value, Reference>,
0404 std::size_t, Reference, Value> {
0405 public:
0406 iterator_impl() = default;
0407 template <class V, class R>
0408 iterator_impl(const iterator_impl<V, R>& it)
0409 : iterator_impl::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
0410 iterator_impl(buffer_type* b, std::size_t i) noexcept
0411 : iterator_impl::iterator_adaptor_(i), buffer_(b) {}
0412
0413 Reference operator*() const noexcept { return {*buffer_, this->base()}; }
0414
0415 template <class V, class R>
0416 friend class iterator_impl;
0417
0418 private:
0419 mutable buffer_type* buffer_ = nullptr;
0420 };
0421
0422 public:
0423 using const_iterator = iterator_impl<const value_type, const_reference>;
0424 using iterator = iterator_impl<value_type, reference>;
0425
0426 explicit unlimited_storage(const allocator_type& a = {}) : buffer_(a) {}
0427 unlimited_storage(const unlimited_storage&) = default;
0428 unlimited_storage& operator=(const unlimited_storage&) = default;
0429 unlimited_storage(unlimited_storage&&) = default;
0430 unlimited_storage& operator=(unlimited_storage&&) = default;
0431
0432
0433
0434
0435
0436 template <class Iterable, class = detail::requires_iterable<Iterable>>
0437 explicit unlimited_storage(const Iterable& s) {
0438 using std::begin;
0439 using std::end;
0440 auto s_begin = begin(s);
0441 auto s_end = end(s);
0442 using V = typename std::iterator_traits<decltype(begin(s))>::value_type;
0443
0444 auto ti = buffer_type::template type_index<V>();
0445 auto nt = mp11::mp_size<typename buffer_type::types>::value;
0446 const std::size_t size = static_cast<std::size_t>(std::distance(s_begin, s_end));
0447 if (ti < nt)
0448 buffer_.template make<V>(size, s_begin);
0449 else
0450 buffer_.template make<double>(size, s_begin);
0451 }
0452
0453 template <class Iterable, class = detail::requires_iterable<Iterable>>
0454 unlimited_storage& operator=(const Iterable& s) {
0455 *this = unlimited_storage(s);
0456 return *this;
0457 }
0458
0459 allocator_type get_allocator() const { return buffer_.alloc; }
0460
0461 void reset(std::size_t n) { buffer_.template make<U8>(n); }
0462
0463 std::size_t size() const noexcept { return buffer_.size; }
0464
0465 reference operator[](std::size_t i) noexcept { return {buffer_, i}; }
0466 const_reference operator[](std::size_t i) const noexcept { return {buffer_, i}; }
0467
0468 bool operator==(const unlimited_storage& x) const noexcept {
0469 if (size() != x.size()) return false;
0470 return buffer_.visit([&x](const auto* p) {
0471 return x.buffer_.visit([p, n = x.size()](const auto* xp) {
0472 return std::equal(p, p + n, xp, detail::safe_equal{});
0473 });
0474 });
0475 }
0476
0477 template <class Iterable>
0478 bool operator==(const Iterable& iterable) const {
0479 if (size() != iterable.size()) return false;
0480 return buffer_.visit([&iterable](const auto* p) {
0481 return std::equal(p, p + iterable.size(), std::begin(iterable),
0482 detail::safe_equal{});
0483 });
0484 }
0485
0486 unlimited_storage& operator*=(const double x) {
0487 buffer_.visit(multiplier(), buffer_, x);
0488 return *this;
0489 }
0490
0491 iterator begin() noexcept { return {&buffer_, 0}; }
0492 iterator end() noexcept { return {&buffer_, size()}; }
0493 const_iterator begin() const noexcept { return {&buffer_, 0}; }
0494 const_iterator end() const noexcept { return {&buffer_, size()}; }
0495
0496
0497 template <class T>
0498 unlimited_storage(std::size_t s, const T* p, const allocator_type& a = {})
0499 : buffer_(std::move(a)) {
0500 buffer_.template make<T>(s, p);
0501 }
0502
0503 template <class Archive>
0504 void serialize(Archive& ar, unsigned ) {
0505 if (Archive::is_loading::value) {
0506 buffer_type tmp(buffer_.alloc);
0507 std::size_t size;
0508 ar& make_nvp("type", tmp.type);
0509 ar& make_nvp("size", size);
0510 tmp.visit([this, size](auto* tp) {
0511 assert(tp == nullptr);
0512 using T = std::decay_t<decltype(*tp)>;
0513 buffer_.template make<T>(size);
0514 });
0515 } else {
0516 ar& make_nvp("type", buffer_.type);
0517 ar& make_nvp("size", buffer_.size);
0518 }
0519 buffer_.visit([this, &ar](auto* tp) {
0520 auto w = detail::make_array_wrapper(tp, this->buffer_.size);
0521 ar& make_nvp("buffer", w);
0522 });
0523 }
0524
0525 private:
0526 struct incrementor {
0527 template <class T>
0528 void operator()(T* tp, buffer_type& b, std::size_t i) {
0529 assert(tp && i < b.size);
0530 if (!detail::safe_increment(tp[i])) {
0531 using U = detail::next_type<typename buffer_type::types, T>;
0532 b.template make<U>(b.size, tp);
0533 ++static_cast<U*>(b.ptr)[i];
0534 }
0535 }
0536
0537 void operator()(large_int* tp, buffer_type&, std::size_t i) { ++tp[i]; }
0538
0539 void operator()(double* tp, buffer_type&, std::size_t i) { ++tp[i]; }
0540 };
0541
0542 struct adder {
0543 template <class U>
0544 void operator()(double* tp, buffer_type&, std::size_t i, const U& x) {
0545 tp[i] += static_cast<double>(x);
0546 }
0547
0548 void operator()(large_int* tp, buffer_type&, std::size_t i, const large_int& x) {
0549 tp[i] += x;
0550 }
0551
0552 template <class T, class U>
0553 void operator()(T* tp, buffer_type& b, std::size_t i, const U& x) {
0554 is_x_integral(std::is_integral<U>{}, tp, b, i, x);
0555 }
0556
0557 template <class T, class U>
0558 void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
0559 const U& x) {
0560
0561 const auto v = static_cast<double>(x);
0562 b.template make<double>(b.size, tp);
0563 operator()(static_cast<double*>(b.ptr), b, i, v);
0564 }
0565
0566 template <class T>
0567 void is_x_integral(std::false_type, T* tp, buffer_type& b, std::size_t i,
0568 const large_int& x) {
0569
0570 const auto v = static_cast<large_int>(x);
0571 b.template make<large_int>(b.size, tp);
0572 operator()(static_cast<large_int*>(b.ptr), b, i, v);
0573 }
0574
0575 template <class T, class U>
0576 void is_x_integral(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
0577 is_x_unsigned(std::is_unsigned<U>{}, tp, b, i, x);
0578 }
0579
0580 template <class T, class U>
0581 void is_x_unsigned(std::false_type, T* tp, buffer_type& b, std::size_t i,
0582 const U& x) {
0583 if (x >= 0)
0584 is_x_unsigned(std::true_type{}, tp, b, i, detail::make_unsigned(x));
0585 else
0586 is_x_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
0587 }
0588
0589 template <class T, class U>
0590 void is_x_unsigned(std::true_type, T* tp, buffer_type& b, std::size_t i, const U& x) {
0591 if (detail::safe_radd(tp[i], x)) return;
0592
0593 const auto y = x;
0594 using TN = detail::next_type<typename buffer_type::types, T>;
0595 b.template make<TN>(b.size, tp);
0596 is_x_unsigned(std::true_type{}, static_cast<TN*>(b.ptr), b, i, y);
0597 }
0598
0599 template <class U>
0600 void is_x_unsigned(std::true_type, large_int* tp, buffer_type&, std::size_t i,
0601 const U& x) {
0602 tp[i] += x;
0603 }
0604 };
0605
0606 struct multiplier {
0607 template <class T>
0608 void operator()(T* tp, buffer_type& b, const double x) {
0609
0610 b.template make<double>(b.size, tp);
0611 operator()(static_cast<double*>(b.ptr), b, x);
0612 }
0613
0614 void operator()(double* tp, buffer_type& b, const double x) {
0615 for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
0616 }
0617
0618 template <class T>
0619 void operator()(T* tp, buffer_type& b, std::size_t i, const double x) {
0620 b.template make<double>(b.size, tp);
0621 operator()(static_cast<double*>(b.ptr), b, i, x);
0622 }
0623
0624 void operator()(double* tp, buffer_type&, std::size_t i, const double x) {
0625 tp[i] *= static_cast<double>(x);
0626 }
0627 };
0628
0629 mutable buffer_type buffer_;
0630 friend struct unsafe_access;
0631 };
0632
0633 }
0634 }
0635
0636 #endif