924 lines
29 KiB
C++
924 lines
29 KiB
C++
// Copyright 2015-2017 Hans Dembinski
|
|
//
|
|
// Distributed under the Boost Software License, Version 1.0.
|
|
// (See accompanying file LICENSE_1_0.txt
|
|
// or copy at http://www.boost.org/LICENSE_1_0.txt)
|
|
|
|
#ifndef BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
|
|
#define BOOST_HISTOGRAM_UNLIMTED_STORAGE_HPP
|
|
|
|
#include <algorithm>
|
|
#include <boost/assert.hpp>
|
|
#include <boost/config/workaround.hpp>
|
|
#include <boost/cstdint.hpp>
|
|
#include <boost/histogram/detail/meta.hpp>
|
|
#include <boost/histogram/fwd.hpp>
|
|
#include <boost/iterator/iterator_adaptor.hpp>
|
|
#include <boost/mp11/algorithm.hpp>
|
|
#include <boost/mp11/list.hpp>
|
|
#include <cmath>
|
|
#include <functional>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <type_traits>
|
|
|
|
namespace boost {
|
|
namespace histogram {
|
|
|
|
namespace detail {
|
|
|
|
// version of std::equal_to<> which handles comparison of signed and unsigned
|
|
struct equal {
|
|
template <class T, class U>
|
|
bool operator()(const T& t, const U& u) const noexcept {
|
|
return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t == u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return u >= 0 && t == make_unsigned(u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t >= 0 && make_unsigned(t) == u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return t == u;
|
|
}
|
|
};
|
|
|
|
// version of std::less<> which handles comparison of signed and unsigned
|
|
struct less {
|
|
template <class T, class U>
|
|
bool operator()(const T& t, const U& u) const noexcept {
|
|
return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t < u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return u >= 0 && t < make_unsigned(u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t < 0 || make_unsigned(t) < u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return t < u;
|
|
}
|
|
};
|
|
|
|
// version of std::greater<> which handles comparison of signed and unsigned
|
|
struct greater {
|
|
template <class T, class U>
|
|
bool operator()(const T& t, const U& u) const noexcept {
|
|
return impl(std::is_signed<T>{}, std::is_signed<U>{}, t, u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t > u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::false_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return u < 0 || t > make_unsigned(u);
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::false_type, const T& t, const U& u) const noexcept {
|
|
return t >= 0 && make_unsigned(t) > u;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool impl(std::true_type, std::true_type, const T& t, const U& u) const noexcept {
|
|
return t > u;
|
|
}
|
|
};
|
|
|
|
template <class Allocator>
|
|
struct mp_int;
|
|
|
|
template <class T>
|
|
struct is_unsigned_integral : mp11::mp_and<std::is_integral<T>, std::is_unsigned<T>> {};
|
|
|
|
template <class T>
|
|
bool safe_increment(T& t) {
|
|
if (t < std::numeric_limits<T>::max()) {
|
|
++t;
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
template <class T, class U>
|
|
bool safe_radd(T& t, const U& u) {
|
|
static_assert(is_unsigned_integral<T>::value, "T must be unsigned integral type");
|
|
static_assert(is_unsigned_integral<U>::value, "T must be unsigned integral type");
|
|
if (static_cast<T>(std::numeric_limits<T>::max() - t) >= u) {
|
|
t += static_cast<T>(u); // static_cast to suppress conversion warning
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// use boost.multiprecision.cpp_int in your own code, it is much more sophisticated
|
|
// than this implementation; we use it here to reduce coupling between boost libs
|
|
template <class Allocator>
|
|
struct mp_int {
|
|
explicit mp_int(Allocator a = {}) : data(1, 0, std::move(a)) {}
|
|
explicit mp_int(uint64_t v, Allocator a = {}) : data(1, v, std::move(a)) {}
|
|
mp_int(const mp_int&) = default;
|
|
mp_int& operator=(const mp_int&) = default;
|
|
mp_int(mp_int&&) = default;
|
|
mp_int& operator=(mp_int&&) = default;
|
|
|
|
mp_int& operator=(uint64_t o) {
|
|
data = decltype(data)(1, o);
|
|
return *this;
|
|
}
|
|
|
|
mp_int& operator++() {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
std::size_t i = 0;
|
|
while (!safe_increment(data[i])) {
|
|
data[i] = 0;
|
|
++i;
|
|
if (i == data.size()) {
|
|
data.push_back(1);
|
|
break;
|
|
}
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
mp_int& operator+=(const mp_int& o) {
|
|
if (this == &o) {
|
|
auto tmp{o};
|
|
return operator+=(tmp);
|
|
}
|
|
bool carry = false;
|
|
std::size_t i = 0;
|
|
for (uint64_t oi : o.data) {
|
|
auto& di = maybe_extend(i);
|
|
if (carry) {
|
|
if (safe_increment(oi))
|
|
carry = false;
|
|
else {
|
|
++i;
|
|
continue;
|
|
}
|
|
}
|
|
if (!safe_radd(di, oi)) {
|
|
add_remainder(di, oi);
|
|
carry = true;
|
|
}
|
|
++i;
|
|
}
|
|
while (carry) {
|
|
auto& di = maybe_extend(i);
|
|
if (safe_increment(di)) break;
|
|
di = 0;
|
|
++i;
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
mp_int& operator+=(uint64_t o) {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
if (safe_radd(data[0], o)) return *this;
|
|
add_remainder(data[0], o);
|
|
// carry the one, data may grow several times
|
|
std::size_t i = 1;
|
|
while (true) {
|
|
auto& di = maybe_extend(i);
|
|
if (safe_increment(di)) break;
|
|
di = 0;
|
|
++i;
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
operator double() const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
double result = static_cast<double>(data[0]);
|
|
std::size_t i = 0;
|
|
while (++i < data.size())
|
|
result += static_cast<double>(data[i]) * std::pow(2.0, i * 64);
|
|
return result;
|
|
}
|
|
|
|
// total ordering for mp_int, mp_int
|
|
bool operator<(const mp_int& o) const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
BOOST_ASSERT(o.data.size() > 0u);
|
|
// no leading zeros allowed
|
|
BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
|
|
BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
|
|
if (data.size() < o.data.size()) return true;
|
|
if (data.size() > o.data.size()) return false;
|
|
auto s = data.size();
|
|
while (s > 0u) {
|
|
--s;
|
|
if (data[s] < o.data[s]) return true;
|
|
if (data[s] > o.data[s]) return false;
|
|
}
|
|
return false; // args are equal
|
|
}
|
|
|
|
bool operator==(const mp_int& o) const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
BOOST_ASSERT(o.data.size() > 0u);
|
|
// no leading zeros allowed
|
|
BOOST_ASSERT(data.size() == 1 || data.back() > 0u);
|
|
BOOST_ASSERT(o.data.size() == 1 || o.data.back() > 0u);
|
|
if (data.size() != o.data.size()) return false;
|
|
return std::equal(data.begin(), data.end(), o.data.begin());
|
|
}
|
|
|
|
// copied from boost/operators.hpp
|
|
friend bool operator>(const mp_int& x, const mp_int& y) { return y < x; }
|
|
friend bool operator<=(const mp_int& x, const mp_int& y) { return !(y < x); }
|
|
friend bool operator>=(const mp_int& x, const mp_int& y) { return !(x < y); }
|
|
friend bool operator!=(const mp_int& x, const mp_int& y) { return !(x == y); }
|
|
|
|
// total ordering for mp_int, uint64; partial ordering for mp_int, double
|
|
template <class U>
|
|
bool operator<(const U& o) const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
return static_if<is_unsigned_integral<U>>(
|
|
[this](uint64_t o) { return data.size() == 1 && data[0] < o; },
|
|
[this](double o) { return operator double() < o; }, o);
|
|
}
|
|
|
|
template <class U>
|
|
bool operator>(const U& o) const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
BOOST_ASSERT(data.back() > 0u); // no leading zeros allowed
|
|
return static_if<is_unsigned_integral<U>>(
|
|
[this](uint64_t o) { return data.size() > 1 || data[0] > o; },
|
|
[this](double o) { return operator double() > o; }, o);
|
|
}
|
|
|
|
template <class U>
|
|
bool operator==(const U& o) const noexcept {
|
|
BOOST_ASSERT(data.size() > 0u);
|
|
return static_if<is_unsigned_integral<U>>(
|
|
[this](uint64_t o) { return data.size() == 1 && data[0] == o; },
|
|
[this](double o) { return operator double() == o; }, o);
|
|
}
|
|
|
|
// adapted copy from boost/operators.hpp
|
|
template <class U>
|
|
friend bool operator<=(const mp_int& x, const U& y) {
|
|
if (is_unsigned_integral<U>::value) return !(x > y);
|
|
return (x < y) || (x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator>=(const mp_int& x, const U& y) {
|
|
if (is_unsigned_integral<U>::value) return !(x < y);
|
|
return (x > y) || (x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator>(const U& x, const mp_int& y) {
|
|
if (is_unsigned_integral<U>::value) return y < x;
|
|
return y < x;
|
|
}
|
|
template <class U>
|
|
friend bool operator<(const U& x, const mp_int& y) {
|
|
if (is_unsigned_integral<U>::value) return y > x;
|
|
return y > x;
|
|
}
|
|
template <class U>
|
|
friend bool operator<=(const U& x, const mp_int& y) {
|
|
if (is_unsigned_integral<U>::value) return !(y < x);
|
|
return (y > x) || (y == x);
|
|
}
|
|
template <class U>
|
|
friend bool operator>=(const U& x, const mp_int& y) {
|
|
if (is_unsigned_integral<U>::value) return !(y > x);
|
|
return (y < x) || (y == x);
|
|
}
|
|
template <class U>
|
|
friend bool operator==(const U& y, const mp_int& x) {
|
|
return x == y;
|
|
}
|
|
template <class U>
|
|
friend bool operator!=(const U& y, const mp_int& x) {
|
|
return !(x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator!=(const mp_int& y, const U& x) {
|
|
return !(y == x);
|
|
}
|
|
|
|
uint64_t& maybe_extend(std::size_t i) {
|
|
while (i >= data.size()) data.push_back(0);
|
|
return data[i];
|
|
}
|
|
|
|
static void add_remainder(uint64_t& d, const uint64_t o) noexcept {
|
|
BOOST_ASSERT(d > 0u);
|
|
// in decimal system it would look like this:
|
|
// 8 + 8 = 6 = 8 - (9 - 8) - 1
|
|
// 9 + 1 = 0 = 9 - (9 - 1) - 1
|
|
auto tmp = std::numeric_limits<uint64_t>::max();
|
|
tmp -= o;
|
|
--d -= tmp;
|
|
}
|
|
|
|
std::vector<uint64_t, Allocator> data;
|
|
};
|
|
|
|
template <class Allocator>
|
|
auto create_buffer(Allocator& a, std::size_t n) {
|
|
using AT = std::allocator_traits<Allocator>;
|
|
auto ptr = AT::allocate(a, n); // may throw
|
|
static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
|
|
"ptr must be trivially copyable");
|
|
auto it = ptr;
|
|
const auto end = ptr + n;
|
|
try {
|
|
// this loop may throw
|
|
while (it != end) AT::construct(a, it++, typename AT::value_type{});
|
|
} catch (...) {
|
|
// release resources that were already acquired before rethrowing
|
|
while (it != ptr) AT::destroy(a, --it);
|
|
AT::deallocate(a, ptr, n);
|
|
throw;
|
|
}
|
|
return ptr;
|
|
}
|
|
|
|
template <class Allocator, class Iterator>
|
|
auto create_buffer(Allocator& a, std::size_t n, Iterator iter) {
|
|
BOOST_ASSERT(n > 0u);
|
|
using AT = std::allocator_traits<Allocator>;
|
|
auto ptr = AT::allocate(a, n); // may throw
|
|
static_assert(std::is_trivially_copyable<decltype(ptr)>::value,
|
|
"ptr must be trivially copyable");
|
|
auto it = ptr;
|
|
const auto end = ptr + n;
|
|
try {
|
|
// this loop may throw
|
|
while (it != end) AT::construct(a, it++, *iter++);
|
|
} catch (...) {
|
|
// release resources that were already acquired before rethrowing
|
|
while (it != ptr) AT::destroy(a, --it);
|
|
AT::deallocate(a, ptr, n);
|
|
throw;
|
|
}
|
|
return ptr;
|
|
}
|
|
|
|
template <class Allocator>
|
|
void destroy_buffer(Allocator& a, typename std::allocator_traits<Allocator>::pointer p,
|
|
std::size_t n) {
|
|
BOOST_ASSERT(p);
|
|
BOOST_ASSERT(n > 0u);
|
|
using AT = std::allocator_traits<Allocator>;
|
|
auto it = p + n;
|
|
while (it != p) AT::destroy(a, --it);
|
|
AT::deallocate(a, p, n);
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
/**
|
|
Memory-efficient storage for integral counters which cannot overflow.
|
|
|
|
This storage provides a no-overflow-guarantee if it is filled with integral weights
|
|
only. This storage implementation keeps a contiguous array of elemental counters, one
|
|
for each cell. If an operation is requested, which would overflow a counter, the whole
|
|
array is replaced with another of a wider integral type, then the operation is executed.
|
|
The storage uses integers of 8, 16, 32, 64 bits, and then switches to a multiprecision
|
|
integral type, similar to those in
|
|
[Boost.Multiprecision](https://www.boost.org/doc/libs/develop/libs/multiprecision/doc/html/index.html).
|
|
|
|
A scaling operation or adding a floating point number turns the elements into doubles,
|
|
which voids the no-overflow-guarantee.
|
|
*/
|
|
template <class Allocator>
|
|
class unlimited_storage {
|
|
static_assert(
|
|
std::is_same<typename std::allocator_traits<Allocator>::pointer,
|
|
typename std::allocator_traits<Allocator>::value_type*>::value,
|
|
"unlimited_storage requires allocator with trivial pointer type");
|
|
|
|
public:
|
|
using allocator_type = Allocator;
|
|
using value_type = double;
|
|
using mp_int = detail::mp_int<
|
|
typename std::allocator_traits<allocator_type>::template rebind_alloc<uint64_t>>;
|
|
|
|
private:
|
|
using types = mp11::mp_list<uint8_t, uint16_t, uint32_t, uint64_t, mp_int, double>;
|
|
|
|
template <class T>
|
|
static constexpr char type_index() noexcept {
|
|
return static_cast<char>(mp11::mp_find<types, T>::value);
|
|
}
|
|
|
|
struct buffer_type {
|
|
allocator_type alloc;
|
|
std::size_t size = 0;
|
|
char type = 0;
|
|
void* ptr = nullptr;
|
|
|
|
template <class F, class... Ts>
|
|
decltype(auto) apply(F&& f, Ts&&... ts) const {
|
|
// this is intentionally not a switch, the if-chain is faster in benchmarks
|
|
if (type == type_index<uint8_t>())
|
|
return f(static_cast<uint8_t*>(ptr), std::forward<Ts>(ts)...);
|
|
if (type == type_index<uint16_t>())
|
|
return f(static_cast<uint16_t*>(ptr), std::forward<Ts>(ts)...);
|
|
if (type == type_index<uint32_t>())
|
|
return f(static_cast<uint32_t*>(ptr), std::forward<Ts>(ts)...);
|
|
if (type == type_index<uint64_t>())
|
|
return f(static_cast<uint64_t*>(ptr), std::forward<Ts>(ts)...);
|
|
if (type == type_index<mp_int>())
|
|
return f(static_cast<mp_int*>(ptr), std::forward<Ts>(ts)...);
|
|
return f(static_cast<double*>(ptr), std::forward<Ts>(ts)...);
|
|
}
|
|
|
|
buffer_type(allocator_type a = {}) : alloc(std::move(a)) {}
|
|
|
|
buffer_type(buffer_type&& o) noexcept
|
|
: alloc(std::move(o.alloc)), size(o.size), type(o.type), ptr(o.ptr) {
|
|
o.size = 0;
|
|
o.type = 0;
|
|
o.ptr = nullptr;
|
|
}
|
|
|
|
buffer_type& operator=(buffer_type&& o) noexcept {
|
|
if (this != &o) {
|
|
using std::swap;
|
|
swap(alloc, o.alloc);
|
|
swap(size, o.size);
|
|
swap(type, o.type);
|
|
swap(ptr, o.ptr);
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
buffer_type(const buffer_type& o) : alloc(o.alloc) {
|
|
o.apply([this, &o](auto* otp) {
|
|
using T = detail::remove_cvref_t<decltype(*otp)>;
|
|
this->template make<T>(o.size, otp);
|
|
});
|
|
}
|
|
|
|
buffer_type& operator=(const buffer_type& o) {
|
|
*this = buffer_type(o);
|
|
return *this;
|
|
}
|
|
|
|
~buffer_type() noexcept { destroy(); }
|
|
|
|
void destroy() noexcept {
|
|
BOOST_ASSERT((ptr == nullptr) == (size == 0));
|
|
if (ptr == nullptr) return;
|
|
apply([this](auto* tp) {
|
|
using T = detail::remove_cvref_t<decltype(*tp)>;
|
|
using alloc_type =
|
|
typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
|
|
alloc_type a(alloc); // rebind allocator
|
|
detail::destroy_buffer(a, tp, size);
|
|
});
|
|
size = 0;
|
|
type = 0;
|
|
ptr = nullptr;
|
|
}
|
|
|
|
#if defined(_MSC_VER)
|
|
#pragma warning(push)
|
|
#pragma warning(disable : 4244) // possible loss of data
|
|
#endif
|
|
|
|
template <class T>
|
|
void make(std::size_t n) {
|
|
// note: order of commands is to not leave buffer in invalid state upon throw
|
|
destroy();
|
|
if (n > 0) {
|
|
// rebind allocator
|
|
using alloc_type =
|
|
typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
|
|
alloc_type a(alloc);
|
|
ptr = detail::create_buffer(a, n); // may throw
|
|
}
|
|
size = n;
|
|
type = type_index<T>();
|
|
}
|
|
|
|
template <class T, class U>
|
|
void make(std::size_t n, U iter) {
|
|
// note: iter may be current ptr, so create new buffer before deleting old buffer
|
|
T* new_ptr = nullptr;
|
|
const auto new_type = type_index<T>();
|
|
if (n > 0) {
|
|
// rebind allocator
|
|
using alloc_type =
|
|
typename std::allocator_traits<allocator_type>::template rebind_alloc<T>;
|
|
alloc_type a(alloc);
|
|
new_ptr = detail::create_buffer(a, n, iter); // may throw
|
|
}
|
|
destroy();
|
|
size = n;
|
|
type = new_type;
|
|
ptr = new_ptr;
|
|
}
|
|
|
|
#if defined(_MSC_VER)
|
|
#pragma warning(pop)
|
|
#endif
|
|
};
|
|
|
|
template <class Buffer>
|
|
class reference_t {
|
|
public:
|
|
reference_t(Buffer* b, std::size_t i) : buffer_(b), idx_(i) {}
|
|
|
|
reference_t(const reference_t&) = default;
|
|
reference_t& operator=(const reference_t&) = delete; // references do not rebind
|
|
reference_t& operator=(reference_t&&) = delete; // references do not rebind
|
|
|
|
// minimal operators for partial ordering
|
|
bool operator<(reference_t rhs) const { return op<detail::less>(rhs); }
|
|
bool operator>(reference_t rhs) const { return op<detail::greater>(rhs); }
|
|
bool operator==(reference_t rhs) const { return op<detail::equal>(rhs); }
|
|
|
|
// adapted copy from boost/operators.hpp for partial ordering
|
|
friend bool operator<=(reference_t x, reference_t y) { return !(y < x); }
|
|
friend bool operator>=(reference_t x, reference_t y) { return !(y > x); }
|
|
friend bool operator!=(reference_t y, reference_t x) { return !(x == y); }
|
|
|
|
template <class U>
|
|
bool operator<(const U& rhs) const {
|
|
return op<detail::less>(rhs);
|
|
}
|
|
|
|
template <class U>
|
|
bool operator>(const U& rhs) const {
|
|
return op<detail::greater>(rhs);
|
|
}
|
|
|
|
template <class U>
|
|
bool operator==(const U& rhs) const {
|
|
return op<detail::equal>(rhs);
|
|
}
|
|
|
|
// adapted copy from boost/operators.hpp
|
|
template <class U>
|
|
friend bool operator<=(reference_t x, const U& y) {
|
|
if (detail::is_unsigned_integral<U>::value) return !(x > y);
|
|
return (x < y) || (x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator>=(reference_t x, const U& y) {
|
|
if (detail::is_unsigned_integral<U>::value) return !(x < y);
|
|
return (x > y) || (x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator>(const U& x, reference_t y) {
|
|
return y < x;
|
|
}
|
|
template <class U>
|
|
friend bool operator<(const U& x, reference_t y) {
|
|
return y > x;
|
|
}
|
|
template <class U>
|
|
friend bool operator<=(const U& x, reference_t y) {
|
|
if (detail::is_unsigned_integral<U>::value) return !(y < x);
|
|
return (y > x) || (y == x);
|
|
}
|
|
template <class U>
|
|
friend bool operator>=(const U& x, reference_t y) {
|
|
if (detail::is_unsigned_integral<U>::value) return !(y > x);
|
|
return (y < x) || (y == x);
|
|
}
|
|
template <class U>
|
|
friend bool operator==(const U& y, reference_t x) {
|
|
return x == y;
|
|
}
|
|
template <class U>
|
|
friend bool operator!=(const U& y, reference_t x) {
|
|
return !(x == y);
|
|
}
|
|
template <class U>
|
|
friend bool operator!=(reference_t y, const U& x) {
|
|
return !(y == x);
|
|
}
|
|
|
|
operator double() const {
|
|
return buffer_->apply(
|
|
[this](const auto* tp) { return static_cast<double>(tp[idx_]); });
|
|
}
|
|
|
|
protected:
|
|
template <class Binary, class U>
|
|
bool op(const reference_t<U>& rhs) const {
|
|
const auto i = idx_;
|
|
const auto j = rhs.idx_;
|
|
return buffer_->apply([i, j, &rhs](const auto* ptr) {
|
|
const auto& pi = ptr[i];
|
|
return rhs.buffer_->apply([&pi, j](const auto* q) { return Binary()(pi, q[j]); });
|
|
});
|
|
}
|
|
|
|
template <class Binary, class U>
|
|
bool op(const U& rhs) const {
|
|
const auto i = idx_;
|
|
return buffer_->apply([i, &rhs](const auto* tp) { return Binary()(tp[i], rhs); });
|
|
}
|
|
|
|
template <class U>
|
|
friend class reference_t;
|
|
|
|
Buffer* buffer_;
|
|
std::size_t idx_;
|
|
};
|
|
|
|
public:
|
|
using const_reference = reference_t<const buffer_type>;
|
|
|
|
class reference : public reference_t<buffer_type> {
|
|
using base_type = reference_t<buffer_type>;
|
|
|
|
public:
|
|
using base_type::base_type;
|
|
|
|
reference operator=(reference t) {
|
|
t.buffer_->apply([this, &t](const auto* otp) { *this = otp[t.idx_]; });
|
|
return *this;
|
|
}
|
|
|
|
reference operator=(const_reference t) {
|
|
t.buffer_->apply([this, &t](const auto* otp) { *this = otp[t.idx_]; });
|
|
return *this;
|
|
}
|
|
|
|
template <class U>
|
|
reference operator=(const U& t) {
|
|
base_type::buffer_->apply([this, &t](auto* tp) {
|
|
tp[this->idx_] = 0;
|
|
adder()(tp, *(this->buffer_), this->idx_, t);
|
|
});
|
|
return *this;
|
|
}
|
|
|
|
template <class U>
|
|
reference operator+=(const U& t) {
|
|
base_type::buffer_->apply(adder(), *base_type::buffer_, base_type::idx_, t);
|
|
return *this;
|
|
}
|
|
|
|
template <class U>
|
|
reference operator*=(const U& t) {
|
|
base_type::buffer_->apply(multiplier(), *base_type::buffer_, base_type::idx_, t);
|
|
return *this;
|
|
}
|
|
|
|
template <class U>
|
|
reference operator-=(const U& t) {
|
|
return operator+=(-t);
|
|
}
|
|
|
|
template <class U>
|
|
reference operator/=(const U& t) {
|
|
return operator*=(1.0 / static_cast<double>(t));
|
|
}
|
|
|
|
reference operator++() {
|
|
base_type::buffer_->apply(incrementor(), *base_type::buffer_, base_type::idx_);
|
|
return *this;
|
|
}
|
|
|
|
// minimal operators for partial ordering
|
|
bool operator<(reference rhs) const { return base_type::operator<(rhs); }
|
|
bool operator>(reference rhs) const { return base_type::operator>(rhs); }
|
|
bool operator==(reference rhs) const { return base_type::operator==(rhs); }
|
|
|
|
// adapted copy from boost/operators.hpp for partial ordering
|
|
friend bool operator<=(reference x, reference y) { return !(y < x); }
|
|
friend bool operator>=(reference x, reference y) { return !(y > x); }
|
|
friend bool operator!=(reference y, reference x) { return !(x == y); }
|
|
};
|
|
|
|
private:
|
|
template <class Value, class Reference, class Buffer>
|
|
class iterator_t
|
|
: public boost::iterator_adaptor<iterator_t<Value, Reference, Buffer>, std::size_t,
|
|
Value, std::random_access_iterator_tag, Reference,
|
|
std::ptrdiff_t> {
|
|
|
|
public:
|
|
iterator_t() = default;
|
|
template <class V, class R, class B>
|
|
iterator_t(const iterator_t<V, R, B>& it)
|
|
: iterator_t::iterator_adaptor_(it.base()), buffer_(it.buffer_) {}
|
|
iterator_t(Buffer* b, std::size_t i) noexcept
|
|
: iterator_t::iterator_adaptor_(i), buffer_(b) {}
|
|
|
|
protected:
|
|
template <class V, class R, class B>
|
|
bool equal(const iterator_t<V, R, B>& rhs) const noexcept {
|
|
return buffer_ == rhs.buffer_ && this->base() == rhs.base();
|
|
}
|
|
Reference dereference() const { return {buffer_, this->base()}; }
|
|
|
|
friend class ::boost::iterator_core_access;
|
|
template <class V, class R, class B>
|
|
friend class iterator_t;
|
|
|
|
private:
|
|
Buffer* buffer_ = nullptr;
|
|
};
|
|
|
|
public:
|
|
using const_iterator = iterator_t<const value_type, const_reference, const buffer_type>;
|
|
using iterator = iterator_t<value_type, reference, buffer_type>;
|
|
|
|
explicit unlimited_storage(allocator_type a = {}) : buffer(std::move(a)) {}
|
|
unlimited_storage(const unlimited_storage&) = default;
|
|
unlimited_storage& operator=(const unlimited_storage&) = default;
|
|
unlimited_storage(unlimited_storage&&) = default;
|
|
unlimited_storage& operator=(unlimited_storage&&) = default;
|
|
|
|
template <class T>
|
|
unlimited_storage(const storage_adaptor<T>& s) {
|
|
using V = detail::remove_cvref_t<decltype(s[0])>;
|
|
constexpr auto ti = type_index<V>();
|
|
detail::static_if_c<(ti < mp11::mp_size<types>::value)>(
|
|
[&](auto) { buffer.template make<V>(s.size(), s.begin()); },
|
|
[&](auto) { buffer.template make<double>(s.size(), s.begin()); }, 0);
|
|
}
|
|
|
|
template <class Iterable, class = detail::requires_iterable<Iterable>>
|
|
unlimited_storage& operator=(const Iterable& s) {
|
|
*this = unlimited_storage(s);
|
|
return *this;
|
|
}
|
|
|
|
allocator_type get_allocator() const { return buffer.alloc; }
|
|
|
|
void reset(std::size_t s) { buffer.template make<uint8_t>(s); }
|
|
|
|
std::size_t size() const noexcept { return buffer.size; }
|
|
|
|
reference operator[](std::size_t i) noexcept { return {&buffer, i}; }
|
|
const_reference operator[](std::size_t i) const noexcept { return {&buffer, i}; }
|
|
|
|
bool operator==(const unlimited_storage& o) const noexcept {
|
|
if (size() != o.size()) return false;
|
|
return buffer.apply([&o](const auto* ptr) {
|
|
return o.buffer.apply([ptr, &o](const auto* optr) {
|
|
return std::equal(ptr, ptr + o.size(), optr, detail::equal{});
|
|
});
|
|
});
|
|
}
|
|
|
|
template <class T>
|
|
bool operator==(const T& o) const {
|
|
if (size() != o.size()) return false;
|
|
return buffer.apply([&o](const auto* ptr) {
|
|
return std::equal(ptr, ptr + o.size(), std::begin(o), detail::equal{});
|
|
});
|
|
}
|
|
|
|
unlimited_storage& operator*=(const double x) {
|
|
buffer.apply(multiplier(), buffer, x);
|
|
return *this;
|
|
}
|
|
|
|
iterator begin() noexcept { return {&buffer, 0}; }
|
|
iterator end() noexcept { return {&buffer, size()}; }
|
|
const_iterator begin() const noexcept { return {&buffer, 0}; }
|
|
const_iterator end() const noexcept { return {&buffer, size()}; }
|
|
|
|
/// @private used by unit tests, not part of generic storage interface
|
|
template <class T>
|
|
unlimited_storage(std::size_t s, const T* p, allocator_type a = {})
|
|
: buffer(std::move(a)) {
|
|
buffer.template make<T>(s, p);
|
|
}
|
|
|
|
template <class Archive>
|
|
void serialize(Archive&, unsigned);
|
|
|
|
private:
|
|
struct incrementor {
|
|
template <class T, class Buffer>
|
|
void operator()(T* tp, Buffer& b, std::size_t i) {
|
|
if (!detail::safe_increment(tp[i])) {
|
|
using U = mp11::mp_at_c<types, (type_index<T>() + 1)>;
|
|
b.template make<U>(b.size, tp);
|
|
++static_cast<U*>(b.ptr)[i];
|
|
}
|
|
}
|
|
|
|
template <class Buffer>
|
|
void operator()(mp_int* tp, Buffer&, std::size_t i) {
|
|
++tp[i];
|
|
}
|
|
|
|
template <class Buffer>
|
|
void operator()(double* tp, Buffer&, std::size_t i) {
|
|
++tp[i];
|
|
}
|
|
};
|
|
|
|
struct adder {
|
|
template <class Buffer, class U>
|
|
void operator()(double* tp, Buffer&, std::size_t i, const U& x) {
|
|
tp[i] += static_cast<double>(x);
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void operator()(T* tp, Buffer& b, std::size_t i, const U& x) {
|
|
U_is_integral(std::is_integral<U>{}, tp, b, i, x);
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void U_is_integral(std::false_type, T* tp, Buffer& b, std::size_t i, const U& x) {
|
|
b.template make<double>(b.size, tp);
|
|
operator()(static_cast<double*>(b.ptr), b, i, x);
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void U_is_integral(std::true_type, T* tp, Buffer& b, std::size_t i, const U& x) {
|
|
U_is_unsigned_integral(std::is_unsigned<U>{}, tp, b, i, x);
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void U_is_unsigned_integral(std::false_type, T* tp, Buffer& b, std::size_t i,
|
|
const U& x) {
|
|
if (x >= 0)
|
|
U_is_unsigned_integral(std::true_type{}, tp, b, i, detail::make_unsigned(x));
|
|
else
|
|
U_is_integral(std::false_type{}, tp, b, i, static_cast<double>(x));
|
|
}
|
|
|
|
template <class Buffer, class U>
|
|
void U_is_unsigned_integral(std::true_type, mp_int* tp, Buffer&, std::size_t i,
|
|
const U& x) {
|
|
tp[i] += x;
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void U_is_unsigned_integral(std::true_type, T* tp, Buffer& b, std::size_t i,
|
|
const U& x) {
|
|
if (detail::safe_radd(tp[i], x)) return;
|
|
using V = mp11::mp_at_c<types, (type_index<T>() + 1)>;
|
|
b.template make<V>(b.size, tp);
|
|
U_is_unsigned_integral(std::true_type{}, static_cast<V*>(b.ptr), b, i, x);
|
|
}
|
|
};
|
|
|
|
struct multiplier {
|
|
template <class T, class Buffer>
|
|
void operator()(T* tp, Buffer& b, const double x) {
|
|
// potential lossy conversion that cannot be avoided
|
|
b.template make<double>(b.size, tp);
|
|
operator()(static_cast<double*>(b.ptr), b, x);
|
|
}
|
|
|
|
template <class Buffer>
|
|
void operator()(double* tp, Buffer& b, const double x) {
|
|
for (auto end = tp + b.size; tp != end; ++tp) *tp *= x;
|
|
}
|
|
|
|
template <class T, class Buffer, class U>
|
|
void operator()(T* tp, Buffer& b, std::size_t i, const U& x) {
|
|
b.template make<double>(b.size, tp);
|
|
operator()(static_cast<double*>(b.ptr), b, i, x);
|
|
}
|
|
|
|
template <class Buffer, class U>
|
|
void operator()(double* tp, Buffer&, std::size_t i, const U& x) {
|
|
tp[i] *= static_cast<double>(x);
|
|
}
|
|
};
|
|
|
|
buffer_type buffer;
|
|
};
|
|
|
|
} // namespace histogram
|
|
} // namespace boost
|
|
|
|
#endif
|