unlimited_storage.hpp 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924
  1. // Copyright 2015-2017 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_UNLIMTED_STORAGE_HPP
  7. #define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
  8. #include <algorithm>
  9. #include <boost/assert.hpp>
  10. #include <boost/config/workaround.hpp>
  11. #include <boost/cstdint.hpp>
  12. #include <boost/histogram/detail/meta.hpp>
  13. #include <boost/histogram/fwd.hpp>
  14. #include <boost/iterator/iterator_adaptor.hpp>
  15. #include <boost/mp11/algorithm.hpp>
  16. #include <boost/mp11/list.hpp>
  17. #include <cmath>
  18. #include <functional>
  19. #include <limits>
  20. #include <memory>
  21. #include <type_traits>
  22. namespace boost {
  23. namespace histogram {
  24. namespace detail {
  25. // version of std::equal_to<> which handles comparison of signed and unsigned
  26. struct equal {
  27. template <class T, class U>
  28. bool operator()(const T& t, const U& u) const noexcept {
  29. return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
  30. }
  31. template <class T, class U>
  32. bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
  33. return t == u;
  34. }
  35. template <class T, class U>
  36. bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
  37. return u >= 0 && t == make_unsigned(u);
  38. }
  39. template <class T, class U>
  40. bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
  41. return t >= 0 && make_unsigned(t) == u;
  42. }
  43. template <class T, class U>
  44. bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
  45. return t == u;
  46. }
  47. };
  48. // version of std::less<> which handles comparison of signed and unsigned
  49. struct less {
  50. template <class T, class U>
  51. bool operator()(const T& t, const U& u) const noexcept {
  52. return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
  53. }
  54. template <class T, class U>
  55. bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
  56. return t < u;
  57. }
  58. template <class T, class U>
  59. bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
  60. return u >= 0 && t < make_unsigned(u);
  61. }
  62. template <class T, class U>
  63. bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
  64. return t < 0 || make_unsigned(t) < u;
  65. }
  66. template <class T, class U>
  67. bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
  68. return t < u;
  69. }
  70. };
  71. // version of std::greater<> which handles comparison of signed and unsigned
  72. struct greater {
  73. template <class T, class U>
  74. bool operator()(const T& t, const U& u) const noexcept {
  75. return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
  76. }
  77. template <class T, class U>
  78. bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
  79. return t > u;
  80. }
  81. template <class T, class U>
  82. bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
  83. return u < 0 || t > make_unsigned(u);
  84. }
  85. template <class T, class U>
  86. bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
  87. return t >= 0 && make_unsigned(t) > u;
  88. }
  89. template <class T, class U>
  90. bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
  91. return t > u;
  92. }
  93. };
  94. template <class Allocator>
  95. struct mp_int;
  96. template <class T>
  97. struct is_unsigned_integral : mp11::mp_and<std::is_integral<T>, std::is_unsigned<T>> {};
  98. template <class T>
  99. bool safe_increment(T& t) {
  100. if (t < std::numeric_limits<T>::max()) {
  101. ++t;
  102. return true;
  103. }
  104. return false;
  105. }
  106. template <class T, class U>
  107. bool safe_radd(T& t, const U& u) {
  108. static_assert(is_unsigned_integral<T>::value, "T must be unsigned integral type");
  109. static_assert(is_unsigned_integral<U>::value, "T must be unsigned integral type");
  110. if (static_cast<T>(std::numeric_limits<T>::max() - t) >= u) {
  111. t += static_cast<T>(u); // static_cast to suppress conversion warning
  112. return true;
  113. }
  114. return false;
  115. }
  116. // use boost.multiprecision.cpp_int in your own code, it is much more sophisticated
  117. // than this implementation; we use it here to reduce coupling between boost libs
  118. template <class Allocator>
  119. struct mp_int {
  120. explicit mp_int(Allocator a = {}) : data(1, 0, std::move(a)) {}
  121. explicit mp_int(uint64_t v, Allocator a = {}) : data(1, v, std::move(a)) {}
  122. mp_int(const mp_int&) = default;
  123. mp_int& operator=(const mp_int&) = default;
  124. mp_int(mp_int&&) = default;
  125. mp_int& operator=(mp_int&&) = default;
  126. mp_int& operator=(uint64_t o) {
  127. data = decltype(data)(1, o);
  128. return *this;
  129. }
  130. mp_int& operator++() {
  131. BOOST_ASSERT(data.size() > 0u);
  132. std::size_t i = 0;
  133. while (!safe_increment(data[i])) {
  134. data[i] = 0;
  135. ++i;
  136. if (i == data.size()) {
  137. data.push_back(1);
  138. break;
  139. }
  140. }
  141. return *this;
  142. }
  143. mp_int& operator+=(const mp_int& o) {
  144. if (this == &o) {
  145. auto tmp{o};
  146. return operator+=(tmp);
  147. }
  148. bool carry = false;
  149. std::size_t i = 0;
  150. for (uint64_t oi : o.data) {
  151. auto& di = maybe_extend(i);
  152. if (carry) {
  153. if (safe_increment(oi))
  154. carry = false;
  155. else {
  156. ++i;
  157. continue;
  158. }
  159. }
  160. if (!safe_radd(di, oi)) {
  161. add_remainder(di, oi);
  162. carry = true;
  163. }
  164. ++i;
  165. }
  166. while (carry) {
  167. auto& di = maybe_extend(i);
  168. if (safe_increment(di)) break;
  169. di = 0;
  170. ++i;
  171. }
  172. return *this;
  173. }
  174. mp_int& operator+=(uint64_t o) {
  175. BOOST_ASSERT(data.size() > 0u);
  176. if (safe_radd(data[0], o)) return *this;
  177. add_remainder(data[0], o);
  178. // carry the one, data may grow several times
  179. std::size_t i = 1;
  180. while (true) {
  181. auto& di = maybe_extend(i);
  182. if (safe_increment(di)) break;
  183. di = 0;
  184. ++i;
  185. }
  186. return *this;
  187. }
  188. operator double() const noexcept {
  189. BOOST_ASSERT(data.size() > 0u);
  190. double result = static_cast<double>(data[0]);
  191. std::size_t i = 0;
  192. while (++i < data.size())
  193. result += static_cast<double>(data[i]) * std::pow(2.0, i * 64);
  194. return result;
  195. }
  196. // total ordering for mp_int, mp_int
  197. bool operator<(const mp_int& o) const noexcept {
  198. BOOST_ASSERT(data.size() > 0u);
  199. BOOST_ASSERT(o.data.size() > 0u);
  200. // no leading zeros allowed
  201. BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
  202. BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
  203. if (data.size() < o.data.size()) return true;
  204. if (data.size() > o.data.size()) return false;
  205. auto s = data.size();
  206. while (s > 0u) {
  207. --s;
  208. if (data[s] < o.data[s]) return true;
  209. if (data[s] > o.data[s]) return false;
  210. }
  211. return false; // args are equal
  212. }
  213. bool operator==(const mp_int& o) const noexcept {
  214. BOOST_ASSERT(data.size() > 0u);
  215. BOOST_ASSERT(o.data.size() > 0u);
  216. // no leading zeros allowed
  217. BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
  218. BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
  219. if (data.size() != o.data.size()) return false;
  220. return std::equal(data.begin(), data.end(), o.data.begin());
  221. }
  222. // copied from boost/operators.hpp
  223. friend bool operator>(const mp_int& x, const mp_int& y) { return y < x; }
  224. friend bool operator<=(const mp_int& x, const mp_int& y) { return !(y < x); }
  225. friend bool operator>=(const mp_int& x, const mp_int& y) { return !(x < y); }
  226. friend bool operator!=(const mp_int& x, const mp_int& y) { return !(x == y); }
  227. // total ordering for mp_int, uint64; partial ordering for mp_int, double
  228. template <class U>
  229. bool operator<(const U& o) const noexcept {
  230. BOOST_ASSERT(data.size() > 0u);
  231. return static_if<is_unsigned_integral<U>>(
  232. [this](uint64_t o) { return data.size() == 1 && data[0] < o; },
  233. [this](double o) { return operator double() < o; }, o);
  234. }
  235. template <class U>
  236. bool operator>(const U& o) const noexcept {
  237. BOOST_ASSERT(data.size() > 0u);
  238. BOOST_ASSERT(data.back() > 0u); // no leading zeros allowed
  239. return static_if<is_unsigned_integral<U>>(
  240. [this](uint64_t o) { return data.size() > 1 || data[0] > o; },
  241. [this](double o) { return operator double() > o; }, o);
  242. }
  243. template <class U>
  244. bool operator==(const U& o) const noexcept {
  245. BOOST_ASSERT(data.size() > 0u);
  246. return static_if<is_unsigned_integral<U>>(
  247. [this](uint64_t o) { return data.size() == 1 && data[0] == o; },
  248. [this](double o) { return operator double() == o; }, o);
  249. }
  250. // adapted copy from boost/operators.hpp
  251. template <class U>
  252. friend bool operator<=(const mp_int& x, const U& y) {
  253. if (is_unsigned_integral<U>::value) return !(x > y);
  254. return (x < y) || (x == y);
  255. }
  256. template <class U>
  257. friend bool operator>=(const mp_int& x, const U& y) {
  258. if (is_unsigned_integral<U>::value) return !(x < y);
  259. return (x > y) || (x == y);
  260. }
  261. template <class U>
  262. friend bool operator>(const U& x, const mp_int& y) {
  263. if (is_unsigned_integral<U>::value) return y < x;
  264. return y < x;
  265. }
  266. template <class U>
  267. friend bool operator<(const U& x, const mp_int& y) {
  268. if (is_unsigned_integral<U>::value) return y > x;
  269. return y > x;
  270. }
  271. template <class U>
  272. friend bool operator<=(const U& x, const mp_int& y) {
  273. if (is_unsigned_integral<U>::value) return !(y < x);
  274. return (y > x) || (y == x);
  275. }
  276. template <class U>
  277. friend bool operator>=(const U& x, const mp_int& y) {
  278. if (is_unsigned_integral<U>::value) return !(y > x);
  279. return (y < x) || (y == x);
  280. }
  281. template <class U>
  282. friend bool operator==(const U& y, const mp_int& x) {
  283. return x == y;
  284. }
  285. template <class U>
  286. friend bool operator!=(const U& y, const mp_int& x) {
  287. return !(x == y);
  288. }
  289. template <class U>
  290. friend bool operator!=(const mp_int& y, const U& x) {
  291. return !(y == x);
  292. }
  293. uint64_t& maybe_extend(std::size_t i) {
  294. while (i >= data.size()) data.push_back(0);
  295. return data[i];
  296. }
  297. static void add_remainder(uint64_t& d, const uint64_t o) noexcept {
  298. BOOST_ASSERT(d > 0u);
  299. // in decimal system it would look like this:
  300. // 8 + 8 = 6 = 8 - (9 - 8) - 1
  301. // 9 + 1 = 0 = 9 - (9 - 1) - 1
  302. auto tmp = std::numeric_limits<uint64_t>::max();
  303. tmp -= o;
  304. --d -= tmp;
  305. }
  306. std::vector<uint64_t, Allocator> data;
  307. };
  308. template <class Allocator>
  309. auto create_buffer(Allocator& a, std::size_t n) {
  310. using AT = std::allocator_traits<Allocator>;
  311. auto ptr = AT::allocate(a, n); // may throw
  312. static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
  313. "ptr must be trivially copyable");
  314. auto it = ptr;
  315. const auto end = ptr + n;
  316. try {
  317. // this loop may throw
  318. while (it != end) AT::construct(a, it++, typename AT::value_type{});
  319. } catch (...) {
  320. // release resources that were already acquired before rethrowing
  321. while (it != ptr) AT::destroy(a, --it);
  322. AT::deallocate(a, ptr, n);
  323. throw;
  324. }
  325. return ptr;
  326. }
  327. template <class Allocator, class Iterator>
  328. auto create_buffer(Allocator& a, std::size_t n, Iterator iter) {
  329. BOOST_ASSERT(n > 0u);
  330. using AT = std::allocator_traits<Allocator>;
  331. auto ptr = AT::allocate(a, n); // may throw
  332. static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
  333. "ptr must be trivially copyable");
  334. auto it = ptr;
  335. const auto end = ptr + n;
  336. try {
  337. // this loop may throw
  338. while (it != end) AT::construct(a, it++, *iter++);
  339. } catch (...) {
  340. // release resources that were already acquired before rethrowing
  341. while (it != ptr) AT::destroy(a, --it);
  342. AT::deallocate(a, ptr, n);
  343. throw;
  344. }
  345. return ptr;
  346. }
  347. template <class Allocator>
  348. void destroy_buffer(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
  349. std::size_t n) {
  350. BOOST_ASSERT(p);
  351. BOOST_ASSERT(n > 0u);
  352. using AT = std::allocator_traits<Allocator>;
  353. auto it = p + n;
  354. while (it != p) AT::destroy(a, --it);
  355. AT::deallocate(a, p, n);
  356. }
  357. } // namespace detail
  358. /**
  359. Memory-efficient storage for integral counters which cannot overflow.
  360. This storage provides a no-overflow-guarantee if it is filled with integral weights
  361. only. This storage implementation keeps a contiguous array of elemental counters, one
  362. for each cell. If an operation is requested, which would overflow a counter, the whole
  363. array is replaced with another of a wider integral type, then the operation is executed.
  364. The storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
  365. integral type, similar to those in
  366. [Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
  367. A scaling operation or adding a floating point number turns the elements into doubles,
  368. which voids the no-overflow-guarantee.
  369. */
  370. template <class Allocator>
  371. class unlimited_storage {
  372. static_assert(
  373. std::is_same<typename std::allocator_traits<Allocator>::pointer,
  374. typename std::allocator_traits<Allocator>::value_type*>::value,
  375. "unlimited_storage requires allocator with trivial pointer type");
  376. public:
  377. using allocator_type = Allocator;
  378. using value_type = double;
  379. using mp_int = detail::mp_int<
  380. typename std::allocator_traits<allocator_type>::template rebind_alloc<uint64_t>>;
  381. private:
  382. using types = mp11::mp_list<uint8_t, uint16_t, uint32_t, uint64_t, mp_int, double>;
  383. template <class T>
  384. static constexpr char type_index() noexcept {
  385. return static_cast<char>(mp11::mp_find<types, T>::value);
  386. }
  387. struct buffer_type {
  388. allocator_type alloc;
  389. std::size_t size = 0;
  390. char type = 0;
  391. void* ptr = nullptr;
  392. template <class F, class... Ts>
  393. decltype(auto) apply(F&& f, Ts&&... ts) const {
  394. // this is intentionally not a switch, the if-chain is faster in benchmarks
  395. if (type == type_index<uint8_t>())
  396. return f(static_cast<uint8_t*>(ptr), std::forward<Ts>(ts)...);
  397. if (type == type_index<uint16_t>())
  398. return f(static_cast<uint16_t*>(ptr), std::forward<Ts>(ts)...);
  399. if (type == type_index<uint32_t>())
  400. return f(static_cast<uint32_t*>(ptr), std::forward<Ts>(ts)...);
  401. if (type == type_index<uint64_t>())
  402. return f(static_cast<uint64_t*>(ptr), std::forward<Ts>(ts)...);
  403. if (type == type_index<mp_int>())
  404. return f(static_cast<mp_int*>(ptr), std::forward<Ts>(ts)...);
  405. return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
  406. }
  407. buffer_type(allocator_type a = {}) : alloc(std::move(a)) {}
  408. buffer_type(buffer_type&& o) noexcept
  409. : alloc(std::move(o.alloc)), size(o.size), type(o.type), ptr(o.ptr) {
  410. o.size = 0;
  411. o.type = 0;
  412. o.ptr = nullptr;
  413. }
  414. buffer_type& operator=(buffer_type&& o) noexcept {
  415. if (this != &o) {
  416. using std::swap;
  417. swap(alloc, o.alloc);
  418. swap(size, o.size);
  419. swap(type, o.type);
  420. swap(ptr, o.ptr);
  421. }
  422. return *this;
  423. }
  424. buffer_type(const buffer_type& o) : alloc(o.alloc) {
  425. o.apply([this, &o](auto* otp) {
  426. using T = detail::remove_cvref_t<decltype(*otp)>;
  427. this->template make<T>(o.size, otp);
  428. });
  429. }
  430. buffer_type& operator=(const buffer_type& o) {
  431. *this = buffer_type(o);
  432. return *this;
  433. }
  434. ~buffer_type() noexcept { destroy(); }
  435. void destroy() noexcept {
  436. BOOST_ASSERT((ptr == nullptr) == (size == 0));
  437. if (ptr == nullptr) return;
  438. apply([this](auto* tp) {
  439. using T = detail::remove_cvref_t<decltype(*tp)>;
  440. using alloc_type =
  441. typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
  442. alloc_type a(alloc); // rebind allocator
  443. detail::destroy_buffer(a, tp, size);
  444. });
  445. size = 0;
  446. type = 0;
  447. ptr = nullptr;
  448. }
  449. #if defined(_MSC_VER)
  450. #pragma warning(push)
  451. #pragma warning(disable : 4244) // possible loss of data
  452. #endif
  453. template <class T>
  454. void make(std::size_t n) {
  455. // note: order of commands is to not leave buffer in invalid state upon throw
  456. destroy();
  457. if (n > 0) {
  458. // rebind allocator
  459. using alloc_type =
  460. typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
  461. alloc_type a(alloc);
  462. ptr = detail::create_buffer(a, n); // may throw
  463. }
  464. size = n;
  465. type = type_index<T>();
  466. }
  467. template <class T, class U>
  468. void make(std::size_t n, U iter) {
  469. // note: iter may be current ptr, so create new buffer before deleting old buffer
  470. T* new_ptr = nullptr;
  471. const auto new_type = type_index<T>();
  472. if (n > 0) {
  473. // rebind allocator
  474. using alloc_type =
  475. typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
  476. alloc_type a(alloc);
  477. new_ptr = detail::create_buffer(a, n, iter); // may throw
  478. }
  479. destroy();
  480. size = n;
  481. type = new_type;
  482. ptr = new_ptr;
  483. }
  484. #if defined(_MSC_VER)
  485. #pragma warning(pop)
  486. #endif
  487. };
  488. template <class Buffer>
  489. class reference_t {
  490. public:
  491. reference_t(Buffer* b, std::size_t i) : buffer_(b), idx_(i) {}
  492. reference_t(const reference_t&) = default;
  493. reference_t& operator=(const reference_t&) = delete; // references do not rebind
  494. reference_t& operator=(reference_t&&) = delete; // references do not rebind
  495. // minimal operators for partial ordering
  496. bool operator<(reference_t rhs) const { return op<detail::less>(rhs); }
  497. bool operator>(reference_t rhs) const { return op<detail::greater>(rhs); }
  498. bool operator==(reference_t rhs) const { return op<detail::equal>(rhs); }
  499. // adapted copy from boost/operators.hpp for partial ordering
  500. friend bool operator<=(reference_t x, reference_t y) { return !(y < x); }
  501. friend bool operator>=(reference_t x, reference_t y) { return !(y > x); }
  502. friend bool operator!=(reference_t y, reference_t x) { return !(x == y); }
  503. template <class U>
  504. bool operator<(const U& rhs) const {
  505. return op<detail::less>(rhs);
  506. }
  507. template <class U>
  508. bool operator>(const U& rhs) const {
  509. return op<detail::greater>(rhs);
  510. }
  511. template <class U>
  512. bool operator==(const U& rhs) const {
  513. return op<detail::equal>(rhs);
  514. }
  515. // adapted copy from boost/operators.hpp
  516. template <class U>
  517. friend bool operator<=(reference_t x, const U& y) {
  518. if (detail::is_unsigned_integral<U>::value) return !(x > y);
  519. return (x < y) || (x == y);
  520. }
  521. template <class U>
  522. friend bool operator>=(reference_t x, const U& y) {
  523. if (detail::is_unsigned_integral<U>::value) return !(x < y);
  524. return (x > y) || (x == y);
  525. }
  526. template <class U>
  527. friend bool operator>(const U& x, reference_t y) {
  528. return y < x;
  529. }
  530. template <class U>
  531. friend bool operator<(const U& x, reference_t y) {
  532. return y > x;
  533. }
  534. template <class U>
  535. friend bool operator<=(const U& x, reference_t y) {
  536. if (detail::is_unsigned_integral<U>::value) return !(y < x);
  537. return (y > x) || (y == x);
  538. }
  539. template <class U>
  540. friend bool operator>=(const U& x, reference_t y) {
  541. if (detail::is_unsigned_integral<U>::value) return !(y > x);
  542. return (y < x) || (y == x);
  543. }
  544. template <class U>
  545. friend bool operator==(const U& y, reference_t x) {
  546. return x == y;
  547. }
  548. template <class U>
  549. friend bool operator!=(const U& y, reference_t x) {
  550. return !(x == y);
  551. }
  552. template <class U>
  553. friend bool operator!=(reference_t y, const U& x) {
  554. return !(y == x);
  555. }
  556. operator double() const {
  557. return buffer_->apply(
  558. [this](const auto* tp) { return static_cast<double>(tp[idx_]); });
  559. }
  560. protected:
  561. template <class Binary, class U>
  562. bool op(const reference_t<U>& rhs) const {
  563. const auto i = idx_;
  564. const auto j = rhs.idx_;
  565. return buffer_->apply([i, j, &rhs](const auto* ptr) {
  566. const auto& pi = ptr[i];
  567. return rhs.buffer_->apply([&pi, j](const auto* q) { return Binary()(pi, q[j]); });
  568. });
  569. }
  570. template <class Binary, class U>
  571. bool op(const U& rhs) const {
  572. const auto i = idx_;
  573. return buffer_->apply([i, &rhs](const auto* tp) { return Binary()(tp[i], rhs); });
  574. }
  575. template <class U>
  576. friend class reference_t;
  577. Buffer* buffer_;
  578. std::size_t idx_;
  579. };
  580. public:
  581. using const_reference = reference_t<const buffer_type>;
  582. class reference : public reference_t<buffer_type> {
  583. using base_type = reference_t<buffer_type>;
  584. public:
  585. using base_type::base_type;
  586. reference operator=(reference t) {
  587. t.buffer_->apply([this, &t](const auto* otp) { *this = otp[t.idx_]; });
  588. return *this;
  589. }
  590. reference operator=(const_reference t) {
  591. t.buffer_->apply([this, &t](const auto* otp) { *this = otp[t.idx_]; });
  592. return *this;
  593. }
  594. template <class U>
  595. reference operator=(const U& t) {
  596. base_type::buffer_->apply([this, &t](auto* tp) {
  597. tp[this->idx_] = 0;
  598. adder()(tp, *(this->buffer_), this->idx_, t);
  599. });
  600. return *this;
  601. }
  602. template <class U>
  603. reference operator+=(const U& t) {
  604. base_type::buffer_->apply(adder(), *base_type::buffer_, base_type::idx_, t);
  605. return *this;
  606. }
  607. template <class U>
  608. reference operator*=(const U& t) {
  609. base_type::buffer_->apply(multiplier(), *base_type::buffer_, base_type::idx_, t);
  610. return *this;
  611. }
  612. template <class U>
  613. reference operator-=(const U& t) {
  614. return operator+=(-t);
  615. }
  616. template <class U>
  617. reference operator/=(const U& t) {
  618. return operator*=(1.0 / static_cast<double>(t));
  619. }
  620. reference operator++() {
  621. base_type::buffer_->apply(incrementor(), *base_type::buffer_, base_type::idx_);
  622. return *this;
  623. }
  624. // minimal operators for partial ordering
  625. bool operator<(reference rhs) const { return base_type::operator<(rhs); }
  626. bool operator>(reference rhs) const { return base_type::operator>(rhs); }
  627. bool operator==(reference rhs) const { return base_type::operator==(rhs); }
  628. // adapted copy from boost/operators.hpp for partial ordering
  629. friend bool operator<=(reference x, reference y) { return !(y < x); }
  630. friend bool operator>=(reference x, reference y) { return !(y > x); }
  631. friend bool operator!=(reference y, reference x) { return !(x == y); }
  632. };
  633. private:
  634. template <class Value, class Reference, class Buffer>
  635. class iterator_t
  636. : public boost::iterator_adaptor<iterator_t<Value, Reference, Buffer>, std::size_t,
  637. Value, std::random_access_iterator_tag, Reference,
  638. std::ptrdiff_t> {
  639. public:
  640. iterator_t() = default;
  641. template <class V, class R, class B>
  642. iterator_t(const iterator_t<V, R, B>& it)
  643. : iterator_t::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
  644. iterator_t(Buffer* b, std::size_t i) noexcept
  645. : iterator_t::iterator_adaptor_(i), buffer_(b) {}
  646. protected:
  647. template <class V, class R, class B>
  648. bool equal(const iterator_t<V, R, B>& rhs) const noexcept {
  649. return buffer_ == rhs.buffer_ && this->base() == rhs.base();
  650. }
  651. Reference dereference() const { return {buffer_, this->base()}; }
  652. friend class ::boost::iterator_core_access;
  653. template <class V, class R, class B>
  654. friend class iterator_t;
  655. private:
  656. Buffer* buffer_ = nullptr;
  657. };
  658. public:
  659. using const_iterator = iterator_t<const value_type, const_reference, const buffer_type>;
  660. using iterator = iterator_t<value_type, reference, buffer_type>;
  661. explicit unlimited_storage(allocator_type a = {}) : buffer(std::move(a)) {}
  662. unlimited_storage(const unlimited_storage&) = default;
  663. unlimited_storage& operator=(const unlimited_storage&) = default;
  664. unlimited_storage(unlimited_storage&&) = default;
  665. unlimited_storage& operator=(unlimited_storage&&) = default;
  666. template <class T>
  667. unlimited_storage(const storage_adaptor<T>& s) {
  668. using V = detail::remove_cvref_t<decltype(s[0])>;
  669. constexpr auto ti = type_index<V>();
  670. detail::static_if_c<(ti < mp11::mp_size<types>::value)>(
  671. [&](auto) { buffer.template make<V>(s.size(), s.begin()); },
  672. [&](auto) { buffer.template make<double>(s.size(), s.begin()); }, 0);
  673. }
  674. template <class Iterable, class = detail::requires_iterable<Iterable>>
  675. unlimited_storage& operator=(const Iterable& s) {
  676. *this = unlimited_storage(s);
  677. return *this;
  678. }
  679. allocator_type get_allocator() const { return buffer.alloc; }
  680. void reset(std::size_t s) { buffer.template make<uint8_t>(s); }
  681. std::size_t size() const noexcept { return buffer.size; }
  682. reference operator[](std::size_t i) noexcept { return {&buffer, i}; }
  683. const_reference operator[](std::size_t i) const noexcept { return {&buffer, i}; }
  684. bool operator==(const unlimited_storage& o) const noexcept {
  685. if (size() != o.size()) return false;
  686. return buffer.apply([&o](const auto* ptr) {
  687. return o.buffer.apply([ptr, &o](const auto* optr) {
  688. return std::equal(ptr, ptr + o.size(), optr, detail::equal{});
  689. });
  690. });
  691. }
  692. template <class T>
  693. bool operator==(const T& o) const {
  694. if (size() != o.size()) return false;
  695. return buffer.apply([&o](const auto* ptr) {
  696. return std::equal(ptr, ptr + o.size(), std::begin(o), detail::equal{});
  697. });
  698. }
  699. unlimited_storage& operator*=(const double x) {
  700. buffer.apply(multiplier(), buffer, x);
  701. return *this;
  702. }
  703. iterator begin() noexcept { return {&buffer, 0}; }
  704. iterator end() noexcept { return {&buffer, size()}; }
  705. const_iterator begin() const noexcept { return {&buffer, 0}; }
  706. const_iterator end() const noexcept { return {&buffer, size()}; }
  707. /// @private used by unit tests, not part of generic storage interface
  708. template <class T>
  709. unlimited_storage(std::size_t s, const T* p, allocator_type a = {})
  710. : buffer(std::move(a)) {
  711. buffer.template make<T>(s, p);
  712. }
  713. template <class Archive>
  714. void serialize(Archive&, unsigned);
  715. private:
  716. struct incrementor {
  717. template <class T, class Buffer>
  718. void operator()(T* tp, Buffer& b, std::size_t i) {
  719. if (!detail::safe_increment(tp[i])) {
  720. using U = mp11::mp_at_c<types, (type_index<T>() + 1)>;
  721. b.template make<U>(b.size, tp);
  722. ++static_cast<U*>(b.ptr)[i];
  723. }
  724. }
  725. template <class Buffer>
  726. void operator()(mp_int* tp, Buffer&, std::size_t i) {
  727. ++tp[i];
  728. }
  729. template <class Buffer>
  730. void operator()(double* tp, Buffer&, std::size_t i) {
  731. ++tp[i];
  732. }
  733. };
  734. struct adder {
  735. template <class Buffer, class U>
  736. void operator()(double* tp, Buffer&, std::size_t i, const U& x) {
  737. tp[i] += static_cast<double>(x);
  738. }
  739. template <class T, class Buffer, class U>
  740. void operator()(T* tp, Buffer& b, std::size_t i, const U& x) {
  741. U_is_integral(std::is_integral<U>{}, tp, b, i, x);
  742. }
  743. template <class T, class Buffer, class U>
  744. void U_is_integral(std::false_type, T* tp, Buffer& b, std::size_t i, const U& x) {
  745. b.template make<double>(b.size, tp);
  746. operator()(static_cast<double*>(b.ptr), b, i, x);
  747. }
  748. template <class T, class Buffer, class U>
  749. void U_is_integral(std::true_type, T* tp, Buffer& b, std::size_t i, const U& x) {
  750. U_is_unsigned_integral(std::is_unsigned<U>{}, tp, b, i, x);
  751. }
  752. template <class T, class Buffer, class U>
  753. void U_is_unsigned_integral(std::false_type, T* tp, Buffer& b, std::size_t i,
  754. const U& x) {
  755. if (x >= 0)
  756. U_is_unsigned_integral(std::true_type{}, tp, b, i, detail::make_unsigned(x));
  757. else
  758. U_is_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
  759. }
  760. template <class Buffer, class U>
  761. void U_is_unsigned_integral(std::true_type, mp_int* tp, Buffer&, std::size_t i,
  762. const U& x) {
  763. tp[i] += x;
  764. }
  765. template <class T, class Buffer, class U>
  766. void U_is_unsigned_integral(std::true_type, T* tp, Buffer& b, std::size_t i,
  767. const U& x) {
  768. if (detail::safe_radd(tp[i], x)) return;
  769. using V = mp11::mp_at_c<types, (type_index<T>() + 1)>;
  770. b.template make<V>(b.size, tp);
  771. U_is_unsigned_integral(std::true_type{}, static_cast<V*>(b.ptr), b, i, x);
  772. }
  773. };
  774. struct multiplier {
  775. template <class T, class Buffer>
  776. void operator()(T* tp, Buffer& b, const double x) {
  777. // potential lossy conversion that cannot be avoided
  778. b.template make<double>(b.size, tp);
  779. operator()(static_cast<double*>(b.ptr), b, x);
  780. }
  781. template <class Buffer>
  782. void operator()(double* tp, Buffer& b, const double x) {
  783. for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
  784. }
  785. template <class T, class Buffer, class U>
  786. void operator()(T* tp, Buffer& b, std::size_t i, const U& x) {
  787. b.template make<double>(b.size, tp);
  788. operator()(static_cast<double*>(b.ptr), b, i, x);
  789. }
  790. template <class Buffer, class U>
  791. void operator()(double* tp, Buffer&, std::size_t i, const U& x) {
  792. tp[i] *= static_cast<double>(x);
  793. }
  794. };
  795. buffer_type buffer;
  796. };
  797. } // namespace histogram
  798. } // namespace boost
  799. #endif