prrng.h

Go to the documentation of this file.

 
#ifndef PRRNG_H
#define PRRNG_H
 
#define PRRNG_PCG32_INITSTATE 0x853c49e6748fea9bULL
 
#define PRRNG_PCG32_INITSEQ 0xda3e39cb94b95bdbULL
 
#define PRRNG_PCG32_MULT 6364136223846793005ULL
 
#include <array>
#include <xtensor/xarray.hpp>
#include <xtensor/xnoalias.hpp>
#include <xtensor/xtensor.hpp>
 
#ifndef PRRNG_USE_BOOST
#define PRRNG_USE_BOOST 1
#endif
 
#if PRRNG_USE_BOOST
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <xtensor/xvectorize.hpp>
#endif
 
#define PRRNG_QUOTE_HELPER(x) #x
#define PRRNG_QUOTE(x) PRRNG_QUOTE_HELPER(x)
 
#define PRRNG_ASSERT_IMPL(expr, file, line) \
    if (!(expr)) { \
        throw std::runtime_error( \
            std::string(file) + ':' + std::to_string(line) + ": assertion failed (" #expr ") \n\t" \
        ); \
    }
 
#define PRRNG_WARNING_IMPL(message, file, line, function) \
    std::cout << std::string(file) + ":" + std::to_string(line) + " (" + std::string(function) + \
                     ")" + ": " message ") \n\t";
 
#ifndef PRRNG_VERSION
#define PRRNG_VERSION "@PROJECT_VERSION@"
#endif
 
#ifdef PRRNG_ENABLE_ASSERT
#define PRRNG_ASSERT(expr) PRRNG_ASSERT_IMPL(expr, __FILE__, __LINE__)
#else
#define PRRNG_ASSERT(expr)
#endif
 
#ifdef PRRNG_ENABLE_DEBUG
#define PRRNG_DEBUG(expr) PRRNG_ASSERT_IMPL(expr, __FILE__, __LINE__)
#else
#define PRRNG_DEBUG(expr)
#endif
 
#ifdef PRRNG_DISABLE_WARNING
#define PRRNG_WARNING(message)
#else
#define PRRNG_WARNING(message) PRRNG_WARNING_IMPL(message, __FILE__, __LINE__, __FUNCTION__)
#endif
 
#ifdef PRRNG_ENABLE_WARNING_PYTHON
#define PRRNG_WARNING_PYTHON(message) PRRNG_WARNING_IMPL(message, __FILE__, __LINE__, __FUNCTION__)
#else
#define PRRNG_WARNING_PYTHON(message)
#endif
 
namespace prrng {
 
enum distribution {
    random, 
    delta, 
    exponential, 
    power, 
    gamma, 
    pareto, 
    weibull, 
    normal, 
    custom 
};
 
std::vector<double> default_parameters(
    enum distribution distribution,
    const std::vector<double>& parameters = std::vector<double>{}
)
{
    std::vector<double> ret;
    switch (distribution) {
    case distribution::random:
        ret = std::vector<double>{1, 0};
        break;
    case distribution::delta:
        ret = std::vector<double>{1, 0};
        break;
    case distribution::exponential:
        ret = std::vector<double>{1, 0};
        break;
    case distribution::power:
        ret = std::vector<double>{1, 0};
        break;
    case distribution::gamma:
        ret = std::vector<double>{1, 1, 0};
        break;
    case distribution::pareto:
        ret = std::vector<double>{1, 1, 0};
        break;
    case distribution::weibull:
        ret = std::vector<double>{1, 1, 0};
        break;
    case distribution::normal:
        ret = std::vector<double>{1, 0, 0};
        break;
    case distribution::custom:
        std::vector<double>{};
        break;
    }
 
    PRRNG_ASSERT(parameters.size() <= ret.size());
    std::copy(parameters.begin(), parameters.end(), ret.begin());
    return ret;
}
 
namespace detail {
 
bool has_correct_parameters(enum distribution distribution, const std::vector<double>& parameters)
{
    switch (distribution) {
    case distribution::random:
        return parameters.size() == 2;
    case distribution::delta:
        return parameters.size() == 2;
    case distribution::exponential:
        return parameters.size() == 2;
    case distribution::power:
        return parameters.size() == 2;
    case distribution::gamma:
        return parameters.size() == 3;
    case distribution::pareto:
        return parameters.size() == 3;
    case distribution::weibull:
        return parameters.size() == 3;
    case distribution::normal:
        return parameters.size() == 3;
    case distribution::custom:
        return true;
    }
 
    return true;
}
 
inline std::string unquote(const std::string& arg)
{
    std::string ret = arg;
    ret.erase(std::remove(ret.begin(), ret.end(), '\"'), ret.end());
    return ret;
}
 
inline std::string replace(std::string str, const std::string& from, const std::string& to)
{
    size_t start_pos = 0;
    while ((start_pos = str.find(from, start_pos)) != std::string::npos) {
        str.replace(start_pos, from.length(), to);
        start_pos += to.length();
    }
    return str;
}
 
template <class T>
struct is_std_array : std::false_type {};
 
template <class T, size_t N>
struct is_std_array<std::array<T, N>> : std::true_type {};
 
template <size_t N, class T, typename = void>
struct check_fixed_rank {
    constexpr static bool value = false;
};
 
template <size_t N, class T>
struct check_fixed_rank<N, T, typename std::enable_if_t<xt::get_rank<T>::value != SIZE_MAX>> {
    constexpr static bool value = (N == xt::get_rank<T>::value);
};
 
template <size_t N, class T>
struct check_fixed_rank<N, T, typename std::enable_if_t<is_std_array<T>::value>> {
    constexpr static bool value = (N == std::tuple_size<T>::value);
};
 
template <typename R, typename = void>
struct get_value_type {
    using type = typename R::value_type;
};
 
template <typename R>
struct get_value_type<R, typename std::enable_if_t<std::is_arithmetic<R>::value>> {
    using type = R;
};
 
template <typename R, size_t N>
struct return_type_fixed {
    using type = typename xt::xtensor<R, N>;
};
 
template <typename R, class S, typename = void>
struct return_type {
    using type = typename xt::xarray<R>;
};
 
template <typename R, class S>
struct return_type<R, S, typename std::enable_if_t<is_std_array<S>::value>> {
    using type = typename xt::xtensor<R, std::tuple_size<S>::value>;
};
 
template <typename R, class S>
struct return_type<R, S, typename std::enable_if_t<xt::has_fixed_rank_t<S>::value>> {
    using type = typename xt::xtensor<R, xt::get_rank<S>::value>;
};
 
template <typename R, class S, class T, typename = void>
struct composite_return_type {
    using type = typename xt::xarray<R>;
};
 
template <typename R, class S, class T>
struct composite_return_type<
    R,
    S,
    T,
    typename std::enable_if_t<is_std_array<S>::value && is_std_array<T>::value>> {
    constexpr static size_t N = std::tuple_size<S>::value + std::tuple_size<T>::value;
    using type = typename xt::xtensor<R, N>;
};
 
template <typename R, typename = void>
struct allocate_return {
 
    using value_type = typename R::value_type;
    R value;
 
    template <class S>
    allocate_return(const S& shape)
    {
#if (XTENSOR_PYTHON_VERSION_MINOR == 26 && XTENSOR_PYTHON_VERSION_PATCH <= 1) || \
    (XTENSOR_PYTHON_VERSION_MINOR < 26)
        if (shape.size() == 0) {
            std::array<typename R::size_type, 1> shape_ = {1};
            value.resize(shape_);
        }
        else {
            value.resize(shape);
        }
#else
        value.resize(shape);
#endif
    }
 
    template <class I, std::size_t L>
    allocate_return(const I (&shape)[L])
    {
        std::array<I, L> shape_;
        std::copy(shape, shape + L, shape_.begin());
        value.resize(shape_);
    }
 
    typename R::value_type* data()
    {
        return &value.front();
    }
 
    size_t size() const
    {
        return value.size();
    }
};
 
template <typename R>
struct allocate_return<R, typename std::enable_if_t<std::is_arithmetic<R>::value>> {
 
    using value_type = R;
    R value;
 
    template <class S>
    allocate_return(const S& shape)
    {
#ifdef PRRNG_ENABLE_ASSERT
        PRRNG_ASSERT(shape.size() <= 1);
        if (shape.size() == 1) {
            PRRNG_ASSERT(shape[0] == 1);
        }
#else
        (void)(shape);
#endif
    }
 
    R* data()
    {
        return &value;
    }
 
    size_t size() const
    {
        return 1;
    }
};
 
template <class S, class T, typename = void>
struct concatenate {
    static auto two(const S& s, const T& t)
    {
        std::vector<size_t> r;
        r.insert(r.begin(), s.cbegin(), s.cend());
        r.insert(r.end(), t.cbegin(), t.cend());
        return r;
    }
};
 
template <class S, class T>
struct concatenate<
    S,
    T,
    typename std::enable_if_t<is_std_array<S>::value && is_std_array<T>::value>> {
    static auto two(const S& s, const T& t)
    {
        std::array<size_t, std::tuple_size<S>::value + std::tuple_size<T>::value> r;
        std::copy(s.cbegin(), s.cend(), r.begin());
        std::copy(t.cbegin(), t.cend(), r.begin() + std::tuple_size<S>::value);
        return r;
    }
};
 
template <class S>
inline size_t size(const S& shape)
{
    using ST = typename S::value_type;
    return std::accumulate(shape.cbegin(), shape.cend(), ST(1), std::multiplies<ST>());
}
 
template <class I, std::size_t L>
std::array<I, L> to_array(const I (&shape)[L])
{
    std::array<I, L> r;
    std::copy(&shape[0], &shape[0] + L, r.begin());
    return r;
}
} // namespace detail
 
inline std::string version()
{
    return detail::unquote(std::string(PRRNG_QUOTE(PRRNG_VERSION)));
}
 
inline std::vector<std::string> version_dependencies()
{
    std::vector<std::string> ret;
 
    ret.push_back("prrng=" + version());
 
    ret.push_back(
        "xtensor=" + detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_VERSION_MINOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_VERSION_PATCH)))
    );
 
#ifdef XSIMD_VERSION_MAJOR
    ret.push_back(
        "xsimd=" + detail::unquote(std::string(PRRNG_QUOTE(XSIMD_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XSIMD_VERSION_MINOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XSIMD_VERSION_PATCH)))
    );
#endif
 
#ifdef XTL_VERSION_MAJOR
    ret.push_back(
        "xtl=" + detail::unquote(std::string(PRRNG_QUOTE(XTL_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTL_VERSION_MINOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTL_VERSION_PATCH)))
    );
#endif
 
#if defined(XTENSOR_PYTHON_VERSION_MAJOR)
    ret.push_back(
        "xtensor-python=" +
        detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_PYTHON_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_PYTHON_VERSION_MINOR))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(XTENSOR_PYTHON_VERSION_PATCH)))
    );
#endif
 
#ifdef BOOST_VERSION
    ret.push_back(
        "boost=" + detail::unquote(std::to_string(BOOST_VERSION / 100000)) + "." +
        detail::unquote(std::to_string((BOOST_VERSION / 100) % 1000)) + "." +
        detail::unquote(std::to_string(BOOST_VERSION % 100))
    );
#endif
 
    std::sort(ret.begin(), ret.end(), std::greater<std::string>());
 
    return ret;
}
 
inline std::vector<std::string> version_compiler()
{
    std::vector<std::string> ret;
 
#ifdef __DATE__
    std::string date = detail::unquote(std::string(PRRNG_QUOTE(__DATE__)));
    ret.push_back("date=" + detail::replace(detail::replace(date, " ", "-"), "--", "-"));
#endif
 
#ifdef __APPLE__
    ret.push_back("platform=apple");
#endif
 
#ifdef __MINGW32__
    ret.push_back("platform=mingw");
#endif
 
#ifdef __linux__
    ret.push_back("platform=linux");
#endif
 
#ifdef _WIN32
    ret.push_back("platform=windows");
#else
#ifdef WIN32
    ret.push_back("platform=windows");
#endif
#endif
 
#ifdef __clang_version__
    ret.push_back(
        "clang=" + detail::unquote(std::string(PRRNG_QUOTE(__clang_major__))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(__clang_minor__))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(__clang_patchlevel__)))
    );
#endif
 
#ifdef __GNUC__
    ret.push_back(
        "gcc=" + detail::unquote(std::string(PRRNG_QUOTE(__GNUC__))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(__GNUC_MINOR__))) + "." +
        detail::unquote(std::string(PRRNG_QUOTE(__GNUC_PATCHLEVEL__)))
    );
#endif
 
#ifdef _MSC_VER
    ret.push_back("msvc=" + std::to_string(_MSC_VER));
#endif
 
    // c++ version
 
#ifdef __cplusplus
    ret.push_back("c++=" + detail::unquote(std::string(PRRNG_QUOTE(__cplusplus))));
#endif
 
    std::sort(ret.begin(), ret.end(), std::greater<std::string>());
 
    return ret;
}
 
namespace iterator {
 
template <class It, class T, class R = size_t>
inline R lower_bound(const It first, const It last, const T& value, R guess = 0, R proximity = 10)
{
    if (proximity == 0) {
        if (value <= *(first)) {
            return 0;
        }
        return std::lower_bound(first, last, value) - first - 1;
    }
 
    if (*(first + guess) < value && value <= *(first + guess + 1)) {
        return guess;
    }
 
    R l = guess > proximity ? guess - proximity : 0;
    R r = std::min(guess + proximity, static_cast<R>(last - first - 1));
 
    if (*(first + l) < value && *(first + r) >= value) {
        return std::lower_bound(first + l, first + r, value) - first - 1;
    }
    else if (value <= *(first)) {
        return 0;
    }
    else {
        return std::lower_bound(first, last, value) - first - 1;
    }
}
 
} // namespace iterator
 
namespace inplace {
 
template <class T, class V, class R>
inline void
lower_bound(const T& matrix, const V& value, R& index, typename R::value_type proximity = 10)
{
    PRRNG_ASSERT(value.dimension() == matrix.dimension() - 1);
    PRRNG_ASSERT(value.dimension() == index.dimension());
 
    auto nd = value.dimension();
    auto stride = matrix.shape(nd);
    auto n = value.size();
 
#ifdef PRRNG_ENABLE_ASSERT
    for (decltype(nd) i = 0; i < nd; ++i) {
        PRRNG_ASSERT(matrix.shape(i) == value.shape(i));
        PRRNG_ASSERT(matrix.shape(i) == index.shape(i));
    }
#endif
 
    for (decltype(n) i = 0; i < n; ++i) {
        index.flat(i) = iterator::lower_bound(
            &matrix.flat(i * stride),
            &matrix.flat(i * stride) + stride,
            value.flat(i),
            index.flat(i),
            proximity
        );
    }
}
 
template <class V, class I>
inline void cumsum_chunk(V& cumsum, const V& delta, const I& shift)
{
    PRRNG_ASSERT(cumsum.dimension() >= 1);
    PRRNG_ASSERT(cumsum.dimension() == delta.dimension());
 
    if (delta.size() == 0) {
        return;
    }
 
    size_t dim = cumsum.dimension();
    size_t n = cumsum.shape(dim - 1);
    size_t ndelta = delta.shape(dim - 1);
 
    for (size_t i = 0; i < shift.size(); ++i) {
 
        if (shift.flat(i) == 0) {
            continue;
        }
 
        auto* d = &delta.flat(i * ndelta);
        auto* c = &cumsum.flat(i * n);
 
        if (shift.flat(i) > 0) {
            PRRNG_ASSERT(shift.flat(i) <= static_cast<typename I::value_type>(n));
            PRRNG_ASSERT(static_cast<typename I::value_type>(ndelta) >= shift.flat(i));
            size_t nadd = static_cast<size_t>(shift.flat(i));
            size_t nkeep = n - nadd;
            auto offset = *(c + n - 1);
            std::copy(c + n - nkeep, c + n, c);
            std::copy(d, d + nadd, c + nkeep);
            *(c + nkeep) += offset;
            std::partial_sum(c + nkeep, c + n, c + nkeep);
        }
        else {
            PRRNG_ASSERT(-shift.flat(i) < static_cast<typename I::value_type>(n));
            PRRNG_ASSERT(static_cast<typename I::value_type>(ndelta) > -shift.flat(i));
            size_t nadd = static_cast<size_t>(-shift.flat(i));
            size_t nkeep = n - nadd;
            auto offset = *(c);
            std::copy(c, c + nkeep, c + nadd);
            std::copy(d, d + nadd + 1, c);
            std::partial_sum(c, c + nadd + 1, c);
            offset -= *(c + nadd);
            std::transform(c, c + nadd + 1, c, [&](auto& v) { return v + offset; });
        }
    }
}
 
} // namespace inplace
 
template <class T, class V, class R>
inline R lower_bound(const T& matrix, const V& value, const R& index, size_t proximity = 10)
{
    R ret = index;
    inplace::lower_bound(matrix, value, ret, proximity);
    return ret;
}
 
template <class T, class V, class R>
inline R lower_bound(const T& matrix, const V& value)
{
    R ret = xt::zeros<typename R::value_type>(value.shape());
    inplace::lower_bound(matrix, value, ret, 0);
    return ret;
}
 
template <class V, class I>
inline V cumsum_chunk(const V& cumsum, const V& delta, const I& shift)
{
    V ret = cumsum;
    inplace::cumsum_chunk(ret, delta, shift);
    return ret;
}
 
namespace detail {
 
template <class G, class D, class P>
void chunk_align_at(
    G& generator,
    const D& get_chunk,
    const P& param,
    double* data,
    ptrdiff_t size,
    ptrdiff_t* start,
    ptrdiff_t index
)
{
    ptrdiff_t ichunk = index - *start;
 
    if (ichunk > param.buffer && ichunk < size - param.buffer) {
        return;
    }
 
    using R = decltype(get_chunk(size_t{}));
    ptrdiff_t n = size;
    ptrdiff_t offset = 0;
    ichunk -= param.margin;
 
    if (ichunk < 0 && ichunk > -size) {
        n = -ichunk;
        std::copy(data, data + size + ichunk, data + n);
    }
    else if (ichunk > 0 && ichunk < size) {
        n = ichunk;
        offset = size - ichunk;
        std::copy(data + n, data + size, data);
    }
 
    *start = index - param.margin;
    generator.jump_to(*start + offset);
    R extra = get_chunk(static_cast<size_t>(n));
    generator.drawn(n);
    std::copy(extra.begin(), extra.end(), data + offset);
}
 
template <class G, class D, class S, class P>
void cumsum_align_at(
    G& generator,
    const D& get_chunk,
    const S& get_sum,
    const P& param,
    double* data,
    ptrdiff_t size,
    ptrdiff_t* start,
    ptrdiff_t index
)
{
    ptrdiff_t ichunk = index - *start;
 
    if (ichunk > param.buffer && ichunk < size - param.buffer) {
        return;
    }
 
    using R = decltype(get_chunk(size_t{}));
    ichunk -= param.margin;
 
    if (ichunk > 0 && ichunk < size) {
        ptrdiff_t n = ichunk;
        ptrdiff_t offset = size - ichunk;
        double back = data[size - 1];
        std::copy(data + n, data + size, data);
 
        *start = index - param.margin;
        generator.jump_to(*start + offset);
        R extra = get_chunk(static_cast<size_t>(n));
        generator.drawn(n);
 
        extra.front() += back;
        std::partial_sum(extra.begin(), extra.end(), extra.begin());
        std::copy(extra.begin(), extra.end(), data + offset);
        return;
    }
 
    if (ichunk < 0 && ichunk > -size) {
        ptrdiff_t n = -ichunk;
        double front = data[0];
        std::copy(data, data + size + ichunk, data + n);
 
        *start = index - param.margin;
        generator.jump_to(*start);
        R extra = get_chunk(static_cast<size_t>(n + 1));
        generator.drawn(n + 1);
 
        std::partial_sum(extra.begin(), extra.end(), extra.begin());
        extra -= extra.back() - front;
        std::copy(extra.begin(), extra.end(), data);
        return;
    }
 
    if (ichunk < 0) {
        ptrdiff_t n = *start - (index - param.margin + size) + 1;
        *start = index - param.margin;
        generator.jump_to(*start);
        R extra = get_chunk(static_cast<size_t>(size));
        generator.drawn(size);
 
        double front = data[0] - get_sum(n);
        generator.drawn(n);
 
        std::partial_sum(extra.begin(), extra.end(), extra.begin());
        extra -= extra.back() - front;
        std::copy(extra.begin(), extra.end(), data);
        return;
    }
 
    generator.jump_to(*start + size);
    ptrdiff_t n = index - param.margin - (*start + size);
    *start = index - param.margin;
    double back = get_sum(n) + data[size - 1];
    generator.drawn(n);
 
    R extra = get_chunk(static_cast<size_t>(size));
    generator.drawn(size);
    extra.front() += back;
    std::partial_sum(extra.begin(), extra.end(), extra.begin());
    std::copy(extra.begin(), extra.end(), data);
}
 
template <class G, class D>
void prev(
    G& generator,
    const D& get_chunk,
    ptrdiff_t margin,
    double* data,
    ptrdiff_t size,
    ptrdiff_t* start
)
{
    using R = decltype(get_chunk(size_t{}));
 
    generator.jump_to(*start - size + margin);
 
    double front = data[0];
    ptrdiff_t m = size - margin + 1;
    R extra = get_chunk(static_cast<size_t>(m));
    generator.drawn(m);
    std::partial_sum(extra.begin(), extra.end(), extra.begin());
    extra -= extra.back() - front;
 
    std::copy(data, data + margin, data + size - margin);
    std::copy(extra.begin(), extra.end() - 1, data);
 
    *start -= size - margin;
}
 
template <class G, class D>
void next(
    G& generator,
    const D& get_chunk,
    ptrdiff_t margin,
    double* data,
    ptrdiff_t size,
    ptrdiff_t* start
)
{
    using R = decltype(get_chunk(size_t{}));
    PRRNG_ASSERT(margin < size);
 
    generator.jump_to(*start + size);
 
    double back = data[size - 1];
    ptrdiff_t n = size - margin;
    R extra = get_chunk(static_cast<size_t>(n));
    generator.drawn(n);
    extra.front() += back;
    std::partial_sum(extra.begin(), extra.end(), extra.begin());
    std::copy(data + size - margin, data + size, data);
    std::copy(extra.begin(), extra.end(), data + margin);
    *start += n;
}
 
template <class G, class D, class S, class P>
void align(
    G& generator,
    const D& get_chunk,
    const S& get_sum,
    const P& param,
    double* data,
    ptrdiff_t size,
    ptrdiff_t* start,
    ptrdiff_t* i,
    double target,
    bool recursive = false
)
{
    using R = decltype(get_chunk(size_t{}));
 
    if (target > data[size - 1]) {
        double delta = data[size - 1] - data[0];
        ptrdiff_t n = size;
        generator.jump_to(*start + size);
        double back = data[size - 1];
        double j = (target - data[size - 1]) / delta - (double)(param.margin) / (double)(n);
 
        if (j > 1) {
            ptrdiff_t m = static_cast<ptrdiff_t>((j - 1) * static_cast<double>(n));
            back += get_sum(static_cast<size_t>(m));
            generator.drawn(m);
            *start += m + size;
            R extra = get_chunk(static_cast<size_t>(n));
            generator.drawn(n);
            extra.front() += back;
            std::partial_sum(extra.begin(), extra.end(), data);
            return align(generator, get_chunk, get_sum, param, data, size, start, i, target, true);
        }
 
        next(generator, get_chunk, 1 + param.margin, data, size, start);
        return align(generator, get_chunk, get_sum, param, data, size, start, i, target, true);
    }
 
    if (target < data[0]) {
        prev(generator, get_chunk, 0, data, size, start);
        return align(generator, get_chunk, get_sum, param, data, size, start, i, target, true);
    }
 
    if (recursive || *i >= size) {
        *i = std::lower_bound(data, data + size, target) - data - 1;
    }
    else {
        *i = iterator::lower_bound(data, data + size, target, *i);
    }
 
    if (*i == param.margin) {
        return;
    }
    if (!recursive && param.buffer > 0 && *i >= param.buffer && *i + param.buffer < size) {
        return;
    }
    if (*i < param.margin) {
        if (!param.strict && *i >= param.min_margin) {
            return;
        }
        prev(generator, get_chunk, 0, data, size, start);
        return align(generator, get_chunk, get_sum, param, data, size, start, i, target, true);
    }
 
    generator.jump_to(*start + size);
    ptrdiff_t n = *i - param.margin;
    R extra = get_chunk(static_cast<size_t>(n));
    generator.drawn(n);
    *start += n;
    *i -= n;
    extra.front() += data[size - 1];
    std::partial_sum(extra.begin(), extra.end(), extra.begin());
    std::copy(data + n, data + size, data);
    std::copy(extra.begin(), extra.end(), data + size - n);
}
 
} // namespace detail
 
class exponential_distribution {
public:
    exponential_distribution(double scale = 1)
    {
        m_scale = scale;
    }
 
    template <class T>
    T pdf(const T& x)
    {
        double rate = 1.0 / m_scale;
        return rate * xt::exp(-rate * x);
    }
 
    template <class T>
    T cdf(const T& x)
    {
        double rate = 1.0 / m_scale;
        return 1.0 - xt::exp(-rate * x);
    }
 
    template <class T>
    T quantile(const T& p)
    {
        return -xt::log(1.0 - p) * m_scale;
    }
 
private:
    double m_scale;
};
 
class power_distribution {
public:
    power_distribution(double k = 1)
    {
        m_k = k;
    }
 
    template <class T>
    T pdf(const T& x)
    {
        return m_k * xt::pow(x, m_k - 1.0);
    }
 
    template <class T>
    T cdf(const T& x)
    {
        return xt::pow(x, m_k);
    }
 
    template <class T>
    T quantile(const T& p)
    {
        return xt::pow(1.0 - p, 1.0 / m_k);
    }
 
private:
    double m_k;
};
 
class gamma_distribution {
public:
    gamma_distribution(double k = 1, double scale = 1)
    {
        m_shape = k;
        m_scale = scale;
    }
 
    template <class T>
    T pdf(const T& x)
    {
        using value_type = typename T::value_type;
 
#if PRRNG_USE_BOOST
        auto f = xt::vectorize(boost::math::gamma_p_derivative<value_type, value_type>);
        T ret = f(m_shape, x / m_scale);
        return xt::where(xt::equal(x, 0), 0.0, ret);
#else
        T ret = x;
        ret.fill(std::numeric_limits<value_type>::quiet_NaN());
        return ret;
#endif
    }
 
    template <class T>
    T cdf(const T& x)
    {
        using value_type = typename T::value_type;
 
#if PRRNG_USE_BOOST
        auto f = xt::vectorize(boost::math::gamma_p<value_type, value_type>);
        return f(m_shape, x / m_scale);
#else
        auto ret = x;
        ret.fill(std::numeric_limits<value_type>::quiet_NaN());
        return ret;
#endif
    }
 
    template <class T>
    T quantile(const T& p)
    {
        using value_type = typename T::value_type;
 
#if PRRNG_USE_BOOST
        auto f = xt::vectorize(boost::math::gamma_p_inv<value_type, value_type>);
        return m_scale * f(m_shape, p);
#else
        auto ret = p;
        ret.fill(std::numeric_limits<value_type>::quiet_NaN());
        return ret;
#endif
    }
 
private:
    double m_shape;
    double m_scale;
};
 
class pareto_distribution {
public:
    pareto_distribution(double k = 1, double scale = 1)
    {
        m_k = k;
        m_scale = scale;
        PRRNG_ASSERT(m_k > 0);
        PRRNG_ASSERT(m_scale > 0);
    }
 
    template <class T>
    T pdf(const T& x)
    {
        return m_k * std::pow(m_scale, m_k) * xt::pow(x, -(m_k + 1.0));
    }
 
    template <class T>
    T cdf(const T& x)
    {
        return 1.0 - std::pow(m_scale, m_k) * xt::pow(x, -m_k);
    }
 
    template <class T>
    T quantile(const T& p)
    {
        return m_scale * xt::pow(1.0 - p, -1.0 / m_k);
    }
 
private:
    double m_k;
    double m_scale;
};
 
class weibull_distribution {
public:
    weibull_distribution(double k = 1, double scale = 1)
    {
        m_shape = k;
        m_scale = scale;
    }
 
    template <class T>
    T pdf(const T& x)
    {
        T ret = xt::exp(-xt::pow(x / m_scale, m_shape));
        ret *= xt::pow(x / m_scale, m_shape - 1.0) * m_shape / m_scale;
 
        if (xt::any(xt::equal(x, 0))) {
            if (m_shape == 1) {
                ret = xt::where(xt::equal(x, 0), 1.0 / m_scale, ret);
            }
            else if (m_shape > 1) {
                ret = xt::where(xt::equal(x, 0), 0.0, ret);
            }
            else {
                throw std::runtime_error("[prrng::weibull_distribution::pdf] Overflow error");
            }
        }
 
        return ret;
    }
 
    template <class T>
    T cdf(const T& x)
    {
        return -xt::expm1(-xt::pow(x / m_scale, m_shape));
    }
 
    template <class T>
    T quantile(const T& p)
    {
        return m_scale * xt::pow(-xt::log1p(-p), 1.0 / m_shape);
    }
 
private:
    double m_shape;
    double m_scale;
};
 
class normal_distribution {
public:
    normal_distribution(double mu = 0, double sigma = 1)
    {
        m_mu = mu;
        m_sigma = sigma;
        m_sigma_sqrt2 = m_sigma * std::sqrt(2.0);
    }
 
    template <class T>
    T pdf(const T& x)
    {
        return xt::exp(-0.5 * xt::pow((x - m_mu) / m_sigma, 2.0)) /
               (m_sigma_sqrt2 * std::sqrt(xt::numeric_constants<double>::PI));
    }
 
    template <class T>
    T cdf(const T& x)
    {
        return 0.5 * (1.0 + xt::erf((x - m_mu) / m_sigma_sqrt2));
    }
 
    template <class T>
    T quantile(const T& p)
    {
        using value_type = typename detail::get_value_type<T>::type;
 
#if PRRNG_USE_BOOST
        auto f = xt::vectorize(boost::math::erf_inv<value_type>);
        return m_mu + m_sigma_sqrt2 * f(2.0 * p - 1.0);
#else
        static_assert(xt::is_xexpression<T>::value, "T must be an xexpression");
        auto ret = p;
        ret.fill(std::numeric_limits<value_type>::quiet_NaN());
        return ret;
#endif
    }
 
private:
    double m_mu;
    double m_sigma;
    double m_sigma_sqrt2;
};
 
template <class Derived>
class GeneratorBase {
public:
    double cumsum_random(size_t n)
    {
        double ret = 0.0;
        for (size_t i = 0; i < n; ++i) {
            ret += static_cast<Derived*>(this)->next_double();
        }
        return ret;
    }
 
    double cumsum_delta(size_t n, double scale = 1)
    {
        return static_cast<double>(n) * scale;
    }
 
    double cumsum_exponential(size_t n, double scale = 1)
    {
        double ret = 0.0;
        for (size_t i = 0; i < n; ++i) {
            ret -= std::log(1.0 - static_cast<Derived*>(this)->next_double());
        }
        return scale * ret;
    }
 
    double cumsum_power(size_t n, double k = 1)
    {
        double ret = 0.0;
        double exponent = 1.0 / k;
        for (size_t i = 0; i < n; ++i) {
            ret += std::pow(1.0 - static_cast<Derived*>(this)->next_double(), exponent);
        }
        return ret;
    }
 
    double cumsum_gamma(size_t n, double k = 1, double scale = 1)
    {
#if PRRNG_USE_BOOST
        double ret = 0.0;
        for (size_t i = 0; i < n; ++i) {
            ret += boost::math::gamma_p_inv<double, double>(
                k, static_cast<Derived*>(this)->next_double()
            );
        }
        return scale * ret;
#else
        return std::numeric_limits<double>::quiet_NaN();
#endif
    }
 
    double cumsum_pareto(size_t n, double k = 1, double scale = 1)
    {
        double ret = 0.0;
        double exponent = -1.0 / k;
        for (size_t i = 0; i < n; ++i) {
            ret += std::pow(1.0 - static_cast<Derived*>(this)->next_double(), exponent);
        }
        return scale * ret;
    }
 
    double cumsum_weibull(size_t n, double k = 1, double scale = 1)
    {
        double ret = 0.0;
        double k_inv = 1.0 / k;
        for (size_t i = 0; i < n; ++i) {
            ret += std::pow(-std::log(1.0 - static_cast<Derived*>(this)->next_double()), k_inv);
        }
        return scale * ret;
    }
 
    double cumsum_normal(size_t n, double mu = 0, double sigma = 1)
    {
#if PRRNG_USE_BOOST
        double ret = 0.0;
        for (size_t i = 0; i < n; ++i) {
            ret += boost::math::erf_inv<double>(
                2.0 * static_cast<Derived*>(this)->next_double() - 1.0
            );
        }
        return static_cast<double>(n) * mu + sigma * std::sqrt(2.0) * ret;
#else
        return std::numeric_limits<double>::quiet_NaN();
#endif
    }
 
    template <typename Iterator>
    void shuffle(Iterator begin, Iterator end)
    {
        for (Iterator it = end - 1; it > begin; --it) {
            std::iter_swap(
                it, begin + static_cast<Derived*>(this)->next_uint32((uint32_t)(it - begin + 1))
            );
        }
    }
 
    template <class P>
    auto decide(const P& p) -> typename detail::return_type<bool, P>::type
    {
        using R = typename detail::return_type<bool, P>::type;
        R ret = xt::empty<bool>(p.shape());
        this->decide(p, ret);
        return ret;
    }
 
    template <class P, class R>
    R decide(const P& p)
    {
        using value_type = typename detail::get_value_type<R>::type;
        R ret = xt::empty<value_type>(p.shape());
        this->decide(p, ret);
        return ret;
    }
 
    template <class P, class R>
    void decide(const P& p, R& ret)
    {
        PRRNG_ASSERT(xt::has_shape(p, ret.shape()));
        using value_type = typename detail::get_value_type<R>::type;
 
        for (size_t i = 0; i < p.size(); ++i) {
            if (static_cast<Derived*>(this)->next_double() < p.flat(i)) {
                ret.flat(i) = static_cast<value_type>(true);
            }
            else {
                ret.flat(i) = static_cast<value_type>(false);
            }
        }
    }
 
    template <class P, class T>
    auto decide_masked(const P& p, const T& mask) -> typename detail::return_type<bool, P>::type
    {
        using R = typename detail::return_type<bool, P>::type;
        R ret = xt::empty<bool>(p.shape());
        this->decide_masked(p, mask, ret);
        return ret;
    }
 
    template <class P, class T, class R>
    R decide_masked(const P& p, const T& mask)
    {
        using value_type = typename detail::get_value_type<R>::type;
        R ret = xt::empty<value_type>(p.shape());
        this->decide_masked(p, mask, ret);
        return ret;
    }
 
    template <class P, class T, class R>
    void decide_masked(const P& p, const T& mask, R& ret)
    {
        PRRNG_ASSERT(xt::has_shape(p, ret.shape()));
        using value_type = typename detail::get_value_type<R>::type;
 
        for (size_t i = 0; i < p.size(); ++i) {
            if (mask.flat(i)) {
                ret.flat(i) = static_cast<value_type>(false);
            }
            else if (static_cast<Derived*>(this)->next_double() < p.flat(i)) {
                ret.flat(i) = static_cast<value_type>(true);
            }
            else {
                ret.flat(i) = static_cast<value_type>(false);
            }
        }
    }
 
    double random()
    {
        return static_cast<Derived*>(this)->next_double();
    }
 
    template <class S>
    auto random(const S& shape) -> typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->random_impl<R>(shape);
    }
 
    template <class R, class S>
    R random(const S& shape)
    {
        return this->random_impl<R>(shape);
    }
 
    template <class I, std::size_t L>
    auto random(const I (&shape)[L]) -> typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->random_impl<R>(shape);
    }
 
    template <class R, class I, std::size_t L>
    R random(const I (&shape)[L])
    {
        return this->random_impl<R>(shape);
    }
 
    template <typename T>
    T randint(T high)
    {
        PRRNG_ASSERT(high >= 0);
        PRRNG_ASSERT(static_cast<uint32_t>(high) < std::numeric_limits<uint32_t>::max());
 
        uint32_t ret;
        this->draw_list_uint32(&ret, static_cast<uint32_t>(high), 1);
        return static_cast<T>(ret);
    }
 
    template <class S, typename T>
    auto randint(const S& shape, T high) -> typename detail::return_type<T, S>::type
    {
        using R = typename detail::return_type<T, S>::type;
        return this->randint_impl<R>(shape, high);
    }
 
    template <class R, class S, typename T>
    R randint(const S& shape, T high)
    {
        return this->randint_impl<R>(shape, high);
    }
 
    template <class I, std::size_t L, typename T>
    auto randint(const I (&shape)[L], T high) -> typename detail::return_type_fixed<T, L>::type
    {
        using R = typename detail::return_type_fixed<T, L>::type;
        return this->randint_impl<R>(shape, high);
    }
 
    template <class R, class I, std::size_t L, typename T>
    R randint(const I (&shape)[L], T high)
    {
        return this->randint_impl<R>(shape, high);
    }
 
    template <class S, typename T, typename U>
    auto randint(const S& shape, T low, U high) -> typename detail::return_type<T, S>::type
    {
        using R = typename detail::return_type<T, S>::type;
        return this->randint_impl<R>(shape, low, high);
    }
 
    template <class R, class S, typename T, typename U>
    R randint(const S& shape, T low, U high)
    {
        return this->randint_impl<R>(shape, low, high);
    }
 
    template <class I, std::size_t L, typename T, typename U>
    auto randint(const I (&shape)[L], T low, U high) ->
        typename detail::return_type_fixed<T, L>::type
    {
        using R = typename detail::return_type_fixed<T, L>::type;
        return this->randint_impl<R>(shape, low, high);
    }
 
    template <class R, class I, std::size_t L, typename T, typename U>
    R randint(const I (&shape)[L], T low, U high)
    {
        return this->randint_impl<R>(shape, low, high);
    }
 
    double delta(double scale = 1)
    {
        return scale;
    }
 
    template <class S>
    auto delta(const S& shape, double scale = 1) -> typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->delta_impl<R>(shape, scale);
    }
 
    template <class R, class S>
    R delta(const S& shape, double scale = 1)
    {
        return this->delta_impl<R>(shape, scale);
    }
 
    template <class I, std::size_t L>
    auto delta(const I (&shape)[L], double scale = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->delta_impl<R>(shape, scale);
    }
 
    template <class R, class I, std::size_t L>
    R delta(const I (&shape)[L], double scale = 1)
    {
        return this->delta_impl<R>(shape, scale);
    }
 
    double exponential(double scale = 1)
    {
        return -std::log(1.0 - static_cast<Derived*>(this)->next_double()) * scale;
    }
 
    template <class S>
    auto exponential(const S& shape, double scale = 1) ->
        typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->exponential_impl<R>(shape, scale);
    }
 
    template <class R, class S>
    R exponential(const S& shape, double scale = 1)
    {
        return this->exponential_impl<R>(shape, scale);
    }
 
    template <class I, std::size_t L>
    auto exponential(const I (&shape)[L], double scale = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->exponential_impl<R>(shape, scale);
    }
 
    template <class R, class I, std::size_t L>
    R exponential(const I (&shape)[L], double scale = 1)
    {
        return this->exponential_impl<R>(shape, scale);
    }
 
    double power(double k = 1)
    {
        return std::pow(1.0 - static_cast<Derived*>(this)->next_double(), 1.0 / k);
    }
 
    template <class S>
    auto power(const S& shape, double k = 1) -> typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->power_impl<R>(shape, k);
    }
 
    template <class R, class S>
    R power(const S& shape, double k = 1)
    {
        return this->power_impl<R>(shape, k);
    }
 
    template <class I, std::size_t L>
    auto power(const I (&shape)[L], double k = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->power_impl<R>(shape, k);
    }
 
    template <class R, class I, std::size_t L>
    R power(const I (&shape)[L], double k = 1)
    {
        return this->power_impl<R>(shape, k);
    }
 
    double gamma(double k = 1, double scale = 1)
    {
#if PRRNG_USE_BOOST
        return scale * boost::math::gamma_p_inv<double, double>(
                           k, static_cast<Derived*>(this)->next_double()
                       );
#else
        return std::numeric_limits<double>::quiet_NaN();
#endif
    }
 
    template <class S>
    auto gamma(const S& shape, double k = 1, double scale = 1) ->
        typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->gamma_impl<R>(shape, k, scale);
    }
 
    template <class R, class S>
    R gamma(const S& shape, double k = 1, double scale = 1)
    {
        return this->gamma_impl<R>(shape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto gamma(const I (&shape)[L], double k = 1, double scale = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->gamma_impl<R>(shape, k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R gamma(const I (&shape)[L], double k = 1, double scale = 1)
    {
        return this->gamma_impl<R>(shape, k, scale);
    }
 
    double pareto(double k = 1, double scale = 1)
    {
        return scale * std::pow(1.0 - static_cast<Derived*>(this)->next_double(), -1.0 / k);
    }
 
    template <class S>
    auto pareto(const S& shape, double k = 1, double scale = 1) ->
        typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->pareto_impl<R>(shape, k, scale);
    }
 
    template <class R, class S>
    R pareto(const S& shape, double k = 1, double scale = 1)
    {
        return this->pareto_impl<R>(shape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto pareto(const I (&shape)[L], double k = 1, double scale = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->pareto_impl<R>(shape, k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R pareto(const I (&shape)[L], double k = 1, double scale = 1)
    {
        return this->pareto_impl<R>(shape, k, scale);
    }
 
    double weibull(double k = 1, double scale = 1)
    {
        return scale *
               std::pow(-std::log(1.0 - static_cast<Derived*>(this)->next_double()), 1.0 / k);
    }
 
    template <class S>
    auto weibull(const S& shape, double k = 1, double scale = 1) ->
        typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->weibull_impl<R>(shape, k, scale);
    }
 
    template <class R, class S>
    R weibull(const S& shape, double k = 1, double scale = 1)
    {
        return this->weibull_impl<R>(shape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto weibull(const I (&shape)[L], double k = 1, double scale = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->weibull_impl<R>(shape, k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R weibull(const I (&shape)[L], double k = 1, double scale = 1)
    {
        return this->weibull_impl<R>(shape, k, scale);
    }
 
    double normal(double mu = 0, double sigma = 1)
    {
#if PRRNG_USE_BOOST
        return mu + sigma * std::sqrt(2.0) *
                        boost::math::erf_inv<double>(
                            2.0 * static_cast<Derived*>(this)->next_positive_double() - 1.0
                        );
#else
        return std::numeric_limits<double>::quiet_NaN();
#endif
    }
 
    template <class S>
    auto normal(const S& shape, double mu = 0, double sigma = 1) ->
        typename detail::return_type<double, S>::type
    {
        using R = typename detail::return_type<double, S>::type;
        return this->normal_impl<R>(shape, mu, sigma);
    }
 
    template <class R, class S>
    R normal(const S& shape, double mu = 0, double sigma = 1)
    {
        return this->normal_impl<R>(shape, mu, sigma);
    }
 
    template <class I, std::size_t L>
    auto normal(const I (&shape)[L], double mu = 0, double sigma = 1) ->
        typename detail::return_type_fixed<double, L>::type
    {
        using R = typename detail::return_type_fixed<double, L>::type;
        return this->normal_impl<R>(shape, mu, sigma);
    }
 
    template <class R, class I, std::size_t L>
    R normal(const I (&shape)[L], double mu = 0, double sigma = 1)
    {
        return this->normal_impl<R>(shape, mu, sigma);
    }
 
    double draw(
        enum prrng::distribution distribution,
        std::vector<double> parameters = std::vector<double>{},
        bool append_default = true
    )
    {
        if (append_default) {
            parameters = default_parameters(distribution, parameters);
        }
        else {
            PRRNG_ASSERT(detail::has_correct_parameters(distribution, parameters));
        }
 
        switch (distribution) {
        case prrng::distribution::random:
            return this->random() * parameters[0] + parameters[1];
        case prrng::distribution::delta:
            return this->delta(parameters[0]) + parameters[1];
        case prrng::distribution::exponential:
            return this->exponential(parameters[0]) + parameters[1];
        case prrng::distribution::power:
            return this->power(parameters[0]) + parameters[1];
        case prrng::distribution::pareto:
            return this->pareto(parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::weibull:
            return this->weibull(parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::gamma:
            return this->gamma(parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::normal:
            return this->normal(parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::custom:
            throw std::runtime_error("Unknown distribution");
        }
 
        throw std::runtime_error("Unknown distribution");
    }
 
    template <class R, class S>
    R draw(
        const S& shape,
        enum prrng::distribution distribution,
        std::vector<double> parameters = std::vector<double>{},
        bool append_default = true
    )
    {
        if (append_default) {
            parameters = default_parameters(distribution, parameters);
        }
        else {
            PRRNG_ASSERT(detail::has_correct_parameters(distribution, parameters));
        }
 
        switch (distribution) {
        case prrng::distribution::random:
            return this->random<R>(shape) * parameters[0] + parameters[1];
        case prrng::distribution::delta:
            return this->delta<R>(shape, parameters[0]) + parameters[1];
        case prrng::distribution::exponential:
            return this->exponential<R>(shape, parameters[0]) + parameters[1];
        case prrng::distribution::power:
            return this->power<R>(shape, parameters[0]) + parameters[1];
        case prrng::distribution::pareto:
            return this->pareto<R>(shape, parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::weibull:
            return this->weibull<R>(shape, parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::gamma:
            return this->gamma<R>(shape, parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::normal:
            return this->normal<R>(shape, parameters[0], parameters[1]) + parameters[2];
        case prrng::distribution::custom:
            throw std::runtime_error("Unknown distribution");
        }
 
        throw std::runtime_error("Unknown distribution");
    }
 
    double cumsum(
        size_t n,
        enum prrng::distribution distribution,
        std::vector<double> parameters = std::vector<double>{},
        bool append_default = true
    )
    {
        if (append_default) {
            parameters = default_parameters(distribution, parameters);
        }
        else {
            PRRNG_ASSERT(detail::has_correct_parameters(distribution, parameters));
        }
 
        double m = static_cast<double>(n);
 
        switch (distribution) {
        case prrng::distribution::random:
            return this->cumsum_random(n) * parameters[0] + m * parameters[1];
        case prrng::distribution::delta:
            return this->cumsum_delta(n, parameters[0]) + m * parameters[1];
        case prrng::distribution::exponential:
            return this->cumsum_exponential(n, parameters[0]) + m * parameters[1];
        case prrng::distribution::power:
            return this->cumsum_power(n, parameters[0]) + m * parameters[1];
        case prrng::distribution::pareto:
            return this->cumsum_pareto(n, parameters[0], parameters[1]) + m * parameters[2];
        case prrng::distribution::weibull:
            return this->cumsum_weibull(n, parameters[0], parameters[1]) + m * parameters[2];
        case prrng::distribution::gamma:
            return this->cumsum_gamma(n, parameters[0], parameters[1]) + m * parameters[2];
        case prrng::distribution::normal:
            return this->cumsum_normal(n, parameters[0], parameters[1]) + m * parameters[2];
        case prrng::distribution::custom:
            throw std::runtime_error("Unknown distribution");
        }
 
        throw std::runtime_error("Unknown distribution");
    }
 
private:
    void draw_list_double(double* data, size_t n)
    {
        for (size_t i = 0; i < n; ++i) {
            data[i] = static_cast<Derived*>(this)->next_double();
        }
    }
 
    void draw_list_positive_double(double* data, size_t n)
    {
        for (size_t i = 0; i < n; ++i) {
            data[i] = static_cast<Derived*>(this)->next_positive_double();
        }
    }
 
    void draw_list_uint32(uint32_t* data, uint32_t bound, size_t n)
    {
        for (size_t i = 0; i < n; ++i) {
            data[i] = static_cast<Derived*>(this)->next_uint32(bound);
        }
    }
 
    template <class R, class S>
    R positive_random_impl(const S& shape)
    {
        static_assert(
            std::is_same<typename detail::allocate_return<R>::value_type, double>::value,
            "Return value_type must be double"
        );
 
        detail::allocate_return<R> ret(shape);
        this->draw_list_positive_double(ret.data(), ret.size());
        return std::move(ret.value);
    }
 
    template <class R, class S>
    R random_impl(const S& shape)
    {
        static_assert(
            std::is_same<typename detail::allocate_return<R>::value_type, double>::value,
            "Return value_type must be double"
        );
 
        detail::allocate_return<R> ret(shape);
        this->draw_list_double(ret.data(), ret.size());
        return std::move(ret.value);
    }
 
    template <class R, class S, typename T>
    R randint_impl(const S& shape, T high)
    {
        static_assert(
            std::numeric_limits<typename detail::allocate_return<R>::value_type>::max() >=
                std::numeric_limits<T>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        PRRNG_ASSERT(high >= 0);
        PRRNG_ASSERT(static_cast<uint32_t>(high) < std::numeric_limits<uint32_t>::max());
 
        detail::allocate_return<R> ret(shape);
        std::vector<uint32_t> tmp(ret.size());
        this->draw_list_uint32(&tmp.front(), static_cast<uint32_t>(high), ret.size());
        std::copy(tmp.begin(), tmp.end(), ret.data());
        return std::move(ret.value);
    }
 
    template <class R, class S, typename T, typename U>
    R randint_impl(const S& shape, T low, U high)
    {
        static_assert(
            std::numeric_limits<typename detail::allocate_return<R>::value_type>::min() >=
                std::numeric_limits<T>::min(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename detail::allocate_return<R>::value_type>::max() >=
                std::numeric_limits<T>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename detail::allocate_return<R>::value_type>::min() >=
                std::numeric_limits<U>::min(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename detail::allocate_return<R>::value_type>::max() >=
                std::numeric_limits<U>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        PRRNG_ASSERT(high - low >= 0);
        PRRNG_ASSERT(static_cast<uint32_t>(high - low) < std::numeric_limits<uint32_t>::max());
 
        detail::allocate_return<R> ret(shape);
        std::vector<uint32_t> tmp(ret.size());
        this->draw_list_uint32(&tmp.front(), static_cast<uint32_t>(high - low), ret.size());
        std::copy(tmp.begin(), tmp.end(), ret.data());
        return ret.value + low;
    }
 
    template <class R, class S>
    R delta_impl(const S& shape, double scale)
    {
        detail::allocate_return<R> ret(shape);
        ret.value.fill(scale);
        return std::move(ret.value);
    }
 
    template <class R, class S>
    R exponential_impl(const S& shape, double scale)
    {
        R r = this->random_impl<R>(shape);
        return exponential_distribution(scale).quantile(r);
    }
 
    template <class R, class S>
    R power_impl(const S& shape, double k)
    {
        R r = this->random_impl<R>(shape);
        return power_distribution(k).quantile(r);
    }
 
    template <class R, class S>
    R gamma_impl(const S& shape, double k, double scale)
    {
        R r = this->random_impl<R>(shape);
        return gamma_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R pareto_impl(const S& shape, double k, double scale)
    {
        R r = this->random_impl<R>(shape);
        return pareto_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R weibull_impl(const S& shape, double k, double scale)
    {
        R r = this->random_impl<R>(shape);
        return weibull_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R normal_impl(const S& shape, double mu, double sigma)
    {
        R r = this->positive_random_impl<R>(shape);
        return normal_distribution(mu, sigma).quantile(r);
    }
};
 
class pcg32 : public GeneratorBase<pcg32> {
public:
    template <typename T = uint64_t, typename S = uint64_t>
    pcg32(T initstate = PRRNG_PCG32_INITSTATE, S initseq = PRRNG_PCG32_INITSEQ)
    {
        static_assert(sizeof(uint64_t) >= sizeof(T), "Down-casting not allowed.");
        static_assert(sizeof(uint64_t) >= sizeof(S), "Down-casting not allowed.");
        this->seed(static_cast<uint64_t>(initstate), static_cast<uint64_t>(initseq));
    }
 
    void seed(uint64_t initstate = PRRNG_PCG32_INITSTATE, uint64_t initseq = PRRNG_PCG32_INITSEQ)
    {
        m_initstate = initstate;
        m_initseq = initseq;
 
        m_state = 0U;
        m_inc = (initseq << 1u) | 1u;
        (*this)();
        m_state += initstate;
        (*this)();
    }
 
    uint32_t operator()()
    {
        uint64_t oldstate = m_state;
        m_state = oldstate * PRRNG_PCG32_MULT + m_inc;
        uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
        uint32_t rot = oldstate >> 59u;
        return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
    }
 
    uint32_t next_uint32()
    {
        return (*this)();
    }
 
    uint32_t next_uint32(uint32_t bound)
    {
        // To avoid bias, we need to make the range of the RNG a multiple of
        // bound, which we do by dropping output less than a threshold.
        // A naive scheme to calculate the threshold would be to do
        //
        //     uint32_t threshold = 0x100000000ull % bound;
        //
        // but 64-bit div/mod is slower than 32-bit div/mod (especially on
        // 32-bit platforms).  In essence, we do
        //
        //     uint32_t threshold = (0x100000000ull-bound) % bound;
        //
        // because this version will calculate the same modulus, but the LHS
        // value is less than 2^32.
 
        uint32_t threshold = (~bound + 1u) % bound;
 
        // Uniformity guarantees that this loop will terminate.  In practice, it
        // should usually terminate quickly; on average (assuming all bounds are
        // equally likely), 82.25% of the time, we can expect it to require just
        // one iteration.  In the worst case, someone passes a bound of 2^31 + 1
        // (i.e., 2147483649), which invalidates almost 50% of the range.  In
        // practice, bounds are typically small and only a tiny amount of the range
        // is eliminated.
        for (;;) {
            uint32_t r = next_uint32();
            if (r >= threshold) {
                return r % bound;
            }
        }
    }
 
    float next_float()
    {
        // Trick from MTGP: generate an uniformly distributed
        // single precision number in [1,2) and subtract 1.
        union {
            uint32_t u;
            float f;
        } x;
        x.u = (next_uint32() >> 9) | 0x3f800000u;
        return x.f - 1.0f;
    }
 
    double next_double()
    {
        union {
            uint64_t u;
            double d;
        } x;
 
        x.u = ((uint64_t)next_uint32() << 20) | 0x3ff0000000000000ULL;
 
        return x.d - 1.0;
    }
 
    double next_positive_double()
    {
        union {
            uint64_t u;
            double d;
        } x;
 
        x.u = ((uint64_t)next_uint32() << 20) | 0x3ff0000000000000ULL;
 
        if (x.u == 0x3ff0000000000000ULL) {
            x.u += 1;
        }
 
        return x.d - 1.0;
    }
 
    uint64_t state() const
    {
        return m_state;
    }
 
    template <typename R>
    R state()
    {
        static_assert(
            std::numeric_limits<R>::max() >= std::numeric_limits<decltype(m_state)>::max(),
            "Down-casting not allowed."
        );
 
        static_assert(
            std::numeric_limits<R>::min() <= std::numeric_limits<decltype(m_state)>::min(),
            "Down-casting not allowed."
        );
 
        return static_cast<R>(m_state);
    }
 
    uint64_t initstate() const
    {
        return m_initstate;
    }
 
    template <typename R>
    R initstate() const
    {
        static_assert(
            std::numeric_limits<R>::max() >= std::numeric_limits<decltype(m_initstate)>::max(),
            "Down-casting not allowed."
        );
 
        static_assert(
            std::numeric_limits<R>::min() <= std::numeric_limits<decltype(m_initstate)>::min(),
            "Down-casting not allowed."
        );
 
        return static_cast<R>(m_initstate);
    }
 
    uint64_t initseq() const
    {
        return m_initseq;
    }
 
    template <typename R>
    R initseq() const
    {
        static_assert(
            std::numeric_limits<R>::max() >= std::numeric_limits<decltype(m_initseq)>::max(),
            "Down-casting not allowed."
        );
 
        static_assert(
            std::numeric_limits<R>::min() <= std::numeric_limits<decltype(m_initseq)>::min(),
            "Down-casting not allowed."
        );
 
        return static_cast<R>(m_initseq);
    }
 
    template <typename T>
    void restore(T state)
    {
        static_assert(sizeof(uint64_t) >= sizeof(T), "Down-casting not allowed.");
        m_state = static_cast<uint64_t>(state);
    }
 
    int64_t operator-(const pcg32& other) const
    {
        PRRNG_DEBUG(m_inc == other.m_inc);
 
        uint64_t cur_mult = PRRNG_PCG32_MULT;
        uint64_t cur_plus = m_inc;
        uint64_t cur_state = other.m_state;
        uint64_t the_bit = 1u;
        uint64_t distance = 0u;
 
        while (m_state != cur_state) {
            if ((m_state & the_bit) != (cur_state & the_bit)) {
                cur_state = cur_state * cur_mult + cur_plus;
                distance |= the_bit;
            }
            PRRNG_DEBUG((m_state & the_bit) == (cur_state & the_bit));
            the_bit <<= 1;
            cur_plus = (cur_mult + 1ULL) * cur_plus;
            cur_mult *= cur_mult;
        }
 
        return (int64_t)distance;
    }
 
    template <typename R = int64_t>
    R distance(const pcg32& other) const
    {
        static_assert(sizeof(R) >= sizeof(int64_t), "Down-casting not allowed.");
        int64_t r = this->operator-(other);
 
#ifdef PRRNG_ENABLE_DEBUG
        bool u = std::is_unsigned<R>::value;
        PRRNG_DEBUG((r < 0 && !u) || r >= 0);
#endif
 
        return static_cast<R>(r);
    }
 
    template <
        typename R = int64_t,
        typename T,
        std::enable_if_t<std::is_integral<T>::value, bool> = true>
    R distance(T other_state) const
    {
        static_assert(sizeof(R) >= sizeof(int64_t), "Down-casting not allowed.");
 
        uint64_t cur_mult = PRRNG_PCG32_MULT;
        uint64_t cur_plus = m_inc;
        uint64_t cur_state = other_state;
        uint64_t the_bit = 1u;
        uint64_t distance = 0u;
 
        while (m_state != cur_state) {
            if ((m_state & the_bit) != (cur_state & the_bit)) {
                cur_state = cur_state * cur_mult + cur_plus;
                distance |= the_bit;
            }
            PRRNG_DEBUG((m_state & the_bit) == (cur_state & the_bit));
            the_bit <<= 1;
            cur_plus = (cur_mult + 1ULL) * cur_plus;
            cur_mult *= cur_mult;
        }
 
        int64_t r = (int64_t)distance;
 
#ifdef PRRNG_ENABLE_DEBUG
        bool u = std::is_unsigned<R>::value;
        PRRNG_DEBUG((r < 0 && !u) || r >= 0);
#endif
 
        return static_cast<R>(r);
    }
 
    template <typename T>
    void advance(T distance)
    {
        static_assert(sizeof(int64_t) >= sizeof(T), "Down-casting not allowed.");
 
        int64_t delta_ = static_cast<int64_t>(distance);
 
        uint64_t cur_mult = PRRNG_PCG32_MULT;
        uint64_t cur_plus = m_inc;
        uint64_t acc_mult = 1u;
        uint64_t acc_plus = 0u;
 
        // Even though delta is an unsigned integer, we can pass a signed
        // integer to go backwards, it just goes "the long way round".
        uint64_t delta = (uint64_t)delta_;
 
        while (delta > 0) {
            if (delta & 1) {
                acc_mult *= cur_mult;
                acc_plus = acc_plus * cur_mult + cur_plus;
            }
            cur_plus = (cur_mult + 1) * cur_plus;
            cur_mult *= cur_mult;
            delta /= 2;
        }
        m_state = acc_mult * m_state + acc_plus;
    }
 
    bool operator==(const pcg32& other) const
    {
        return m_state == other.m_state && m_inc == other.m_inc;
    }
 
    bool operator!=(const pcg32& other) const
    {
        return m_state != other.m_state || m_inc != other.m_inc;
    }
 
protected:
    uint64_t m_initstate; 
    uint64_t m_initseq; 
    uint64_t m_state; 
    uint64_t m_inc; 
};
 
class pcg32_index : public pcg32 {
private:
    ptrdiff_t m_index; 
    bool m_delta; 
 
public:
    template <typename T = uint64_t, typename S = uint64_t>
    pcg32_index(
        T initstate = PRRNG_PCG32_INITSTATE,
        S initseq = PRRNG_PCG32_INITSEQ,
        bool delta = false
    )
    {
        static_assert(sizeof(uint64_t) >= sizeof(T), "Down-casting not allowed.");
        static_assert(sizeof(uint64_t) >= sizeof(S), "Down-casting not allowed.");
        this->seed(static_cast<uint64_t>(initstate), static_cast<uint64_t>(initseq));
        m_index = 0;
        m_delta = delta;
    }
 
    uint64_t state_at(ptrdiff_t index)
    {
        if (m_delta) {
            return m_state;
        }
        uint64_t state = m_state;
        this->advance(index - m_index);
        uint64_t ret = m_state;
        this->restore(state);
        return ret;
    }
 
    void jump_to(ptrdiff_t index)
    {
        if (m_delta) {
            return;
        }
        this->advance(index - m_index);
        m_index = index;
    }
 
    void drawn(ptrdiff_t n)
    {
        if (m_delta) {
            return;
        }
        m_index += n;
    }
 
    void set_delta(bool delta)
    {
        m_delta = delta;
    }
 
    ptrdiff_t index() const
    {
        return m_index;
    }
 
    void set_index(ptrdiff_t index)
    {
        m_index = index;
    }
};
 
struct alignment {
    alignment(
        ptrdiff_t buffer = 0,
        ptrdiff_t margin = 0,
        ptrdiff_t min_margin = 0,
        bool strict = false
    )
    {
        this->buffer = buffer;
        this->margin = margin;
        this->min_margin = min_margin;
        this->strict = strict;
    }
 
    ptrdiff_t buffer = 0;
 
    ptrdiff_t margin = 0;
 
    ptrdiff_t min_margin = 0;
 
    bool strict = false;
};
 
template <class Data>
class pcg32_cumsum {
private:
    Data m_data; 
    pcg32_index m_gen; 
    std::function<Data(size_t)> m_draw; 
    std::function<double(size_t)> m_sum; 
    bool m_extendible; 
    alignment m_align; 
    distribution m_distro; 
    std::array<double, 3> m_param; 
    ptrdiff_t m_start; 
    ptrdiff_t m_i; 
 
    void auto_functions()
    {
        m_extendible = true;
 
        switch (m_distro) {
        case random:
            m_draw = [this](size_t n) -> Data {
                return m_gen.random<Data>(std::array<size_t, 1>{n}) * m_param[0] + m_param[1];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_random(n) * m_param[0] + static_cast<double>(n) * m_param[1];
            };
            return;
        case delta:
            m_draw = [this](size_t n) -> Data {
                return m_gen.delta<Data>(std::array<size_t, 1>{n}, m_param[0]) + m_param[1];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_delta(n, m_param[0]) + static_cast<double>(n) * m_param[1];
            };
            return;
        case exponential:
            m_draw = [this](size_t n) -> Data {
                return m_gen.exponential<Data>(std::array<size_t, 1>{n}, m_param[0]) + m_param[1];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_exponential(n, m_param[0]) +
                       static_cast<double>(n) * m_param[1];
            };
            return;
        case power:
            m_draw = [this](size_t n) -> Data {
                return m_gen.power<Data>(std::array<size_t, 1>{n}, m_param[0]) + m_param[1];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_power(n, m_param[0]) + static_cast<double>(n) * m_param[1];
            };
            return;
        case gamma:
            m_draw = [this](size_t n) -> Data {
                return m_gen.gamma<Data>(std::array<size_t, 1>{n}, m_param[0], m_param[1]) +
                       m_param[2];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_gamma(n, m_param[0], m_param[1]) +
                       static_cast<double>(n) * m_param[2];
            };
            return;
        case pareto:
            m_draw = [this](size_t n) -> Data {
                return m_gen.pareto<Data>(std::array<size_t, 1>{n}, m_param[0], m_param[1]) +
                       m_param[2];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_pareto(n, m_param[0], m_param[1]) +
                       static_cast<double>(n) * m_param[2];
            };
            return;
        case weibull:
            m_draw = [this](size_t n) -> Data {
                return m_gen.weibull<Data>(std::array<size_t, 1>{n}, m_param[0], m_param[1]) +
                       m_param[2];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_weibull(n, m_param[0], m_param[1]) +
                       static_cast<double>(n) * m_param[2];
            };
            return;
        case normal:
            m_draw = [this](size_t n) -> Data {
                return m_gen.normal<Data>(std::array<size_t, 1>{n}, m_param[0], m_param[1]) +
                       m_param[2];
            };
            m_sum = [this](size_t n) -> double {
                return m_gen.cumsum_normal(n, m_param[0], m_param[1]) +
                       static_cast<double>(n) * m_param[2];
            };
            return;
        case custom:
            m_extendible = false;
            return;
        }
    }
 
    void copy_from(const pcg32_cumsum& other)
    {
        m_data = other.m_data;
        m_gen = other.m_gen;
        m_align = other.m_align;
        m_distro = other.m_distro;
        m_param = other.m_param;
        m_start = other.m_start;
        m_i = other.m_i;
        this->auto_functions();
    }
 
public:
    template <class R, typename T = uint64_t, typename S = uint64_t>
    pcg32_cumsum(
        const R& shape,
        T initstate = PRRNG_PCG32_INITSTATE,
        S initseq = PRRNG_PCG32_INITSEQ,
        enum distribution distribution = distribution::custom,
        const std::vector<double>& parameters = std::vector<double>{},
        const alignment& align = alignment()
    )
    {
        m_data = xt::empty<typename Data::value_type>(shape);
        m_gen = pcg32_index(initstate, initseq, distribution == distribution::delta);
        m_start = m_gen.index();
        m_i = static_cast<ptrdiff_t>(m_data.size());
        m_align = align;
        m_distro = distribution;
        std::copy(parameters.begin(), parameters.end(), m_param.begin());
        this->auto_functions();
 
        if (!m_extendible) {
            return;
        }
 
        using E = decltype(m_draw(size_t{}));
        E extra = m_draw(m_data.size());
        m_gen.drawn(m_data.size());
        std::partial_sum(extra.begin(), extra.end(), m_data.begin());
    }
 
    void set_functions(
        std::function<Data(size_t)> get_chunk,
        std::function<double(size_t)> get_cumsum,
        bool uses_generator = true
    )
    {
        m_extendible = true;
        m_draw = get_chunk;
        m_sum = get_cumsum;
        m_gen.set_delta(!uses_generator);
 
        using E = decltype(m_draw(size_t{}));
        E extra = m_draw(m_data.size());
        m_gen.drawn(m_data.size());
        std::partial_sum(extra.begin(), extra.end(), m_data.begin());
    }
 
    pcg32_cumsum(const pcg32_cumsum& other)
    {
        this->copy_from(other);
    }
 
    void operator=(const pcg32_cumsum& other)
    {
        this->copy_from(other);
    }
 
    bool is_extendible() const
    {
        return m_extendible;
    }
 
    const pcg32_index& generator() const
    {
        return m_gen;
    }
 
    const auto& shape() const
    {
        return m_data.shape();
    }
 
    auto size() const
    {
        return m_data.size();
    }
 
    template <class T>
    pcg32_cumsum& operator+=(const T& value)
    {
        xt::noalias(m_data) += value;
        return *this;
    }
 
    template <class T>
    pcg32_cumsum& operator-=(const T& value)
    {
        xt::noalias(m_data) -= value;
        return *this;
    }
 
    const Data& data() const
    {
        return m_data;
    }
 
    void set_data(const Data& data)
    {
        PRRNG_ASSERT(xt::has_shape(data, m_data.shape()));
        xt::noalias(m_data) = data;
    }
 
    ptrdiff_t start() const
    {
        return m_start;
    }
 
    void set_start(ptrdiff_t index)
    {
        m_start = index;
    }
 
    ptrdiff_t index_at_align() const
    {
        return m_start + m_i;
    }
 
    ptrdiff_t chunk_index_at_align() const
    {
        return m_i;
    }
 
    double left_of_align() const
    {
        return m_data[m_i];
    }
 
    double right_of_align() const
    {
        return m_data[m_i + 1];
    }
 
    uint64_t state_at(ptrdiff_t index)
    {
        return m_gen.state_at(index);
    }
 
    void restore(uint64_t state, double value, ptrdiff_t index)
    {
        m_gen.set_index(index);
        m_gen.restore(state);
        m_start = index;
 
        using E = decltype(m_draw(size_t{}));
        E extra = m_draw(m_data.size());
        m_gen.drawn(m_data.size());
        extra.front() += value - extra.front();
        std::partial_sum(extra.begin(), extra.end(), m_data.begin());
    }
 
    bool contains(double target) const
    {
        return target >= m_data.front() && target <= m_data.back();
    }
 
    void prev(size_t margin = 0)
    {
        PRRNG_ASSERT(m_extendible);
        m_i = static_cast<ptrdiff_t>(m_data.size());
        detail::prev(m_gen, m_draw, margin, m_data.data(), m_data.size(), &m_start);
    }
 
    void next(size_t margin = 0)
    {
        PRRNG_ASSERT(m_extendible);
        m_i = static_cast<ptrdiff_t>(m_data.size());
        detail::next(m_gen, m_draw, margin, m_data.data(), m_data.size(), &m_start);
    }
 
    void align(double target)
    {
        if (!m_extendible) {
            PRRNG_ASSERT(this->contains(target));
            m_i = iterator::lower_bound(m_data.begin(), m_data.end(), target, m_i);
            return;
        }
 
        detail::align(
            m_gen, m_draw, m_sum, m_align, m_data.data(), m_data.size(), &m_start, &m_i, target
        );
    }
};
 
template <class Derived, class M>
class GeneratorBase_array {
public:
    using shape_type = M; 
    using size_type = typename M::value_type; 
 
    size_type size() const
    {
        return m_size;
    }
 
    const shape_type& strides() const
    {
        return m_strides;
    }
 
    const shape_type& shape() const
    {
        return m_shape;
    }
 
    template <class T>
    auto shape(T axis) const
    {
        return m_shape[axis];
    }
 
    template <class T>
    auto flat_index(const T& index) const
    {
        PRRNG_DEBUG(this->inbounds(index));
        return std::inner_product(index.cbegin(), index.cend(), m_strides.cbegin(), 0);
    }
 
    template <class T>
    bool inbounds(const T& index) const
    {
        if (index.size() != m_strides.size()) {
            return false;
        }
 
        for (size_t i = 0; i < m_strides.size(); ++i) {
            if (index[i] >= m_shape[i]) {
                return false;
            }
        }
 
        return true;
    }
 
    template <class S>
    auto random(const S& ishape) -> typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->random_impl<R>(ishape);
    }
 
    template <class R, class S>
    R random(const S& ishape)
    {
        return this->random_impl<R>(ishape);
    }
 
    template <class I, std::size_t L>
    auto random(const I (&ishape)[L]) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->random_impl<R>(detail::to_array(ishape));
    }
 
    template <class R, class I, std::size_t L>
    R random(const I (&ishape)[L])
    {
        return this->random_impl<R>(detail::to_array(ishape));
    }
 
    template <class S, typename T>
    auto randint(const S& ishape, T high) -> typename detail::composite_return_type<T, M, S>::type
    {
        using R = typename detail::composite_return_type<T, M, S>::type;
        return this->randint_impl<R>(ishape, high);
    }
 
    template <class R, class S, typename T>
    R randint(const S& ishape, T high)
    {
        return this->randint_impl<R>(ishape, high);
    }
 
    template <class I, std::size_t L, typename T>
    auto randint(const I (&ishape)[L], T high) ->
        typename detail::composite_return_type<T, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<T, M, std::array<size_t, L>>::type;
        return this->randint_impl<R>(detail::to_array(ishape), high);
    }
 
    template <class R, class I, std::size_t L, typename T>
    R randint(const I (&ishape)[L], T high)
    {
        return this->randint_impl<R>(detail::to_array(ishape), high);
    }
 
    template <class S, typename T, typename U>
    auto randint(const S& ishape, T low, U high) ->
        typename detail::composite_return_type<T, M, S>::type
    {
        using R = typename detail::composite_return_type<T, M, S>::type;
        return this->randint_impl<R>(ishape, low, high);
    }
 
    template <class R, class S, typename T, typename U>
    R randint(const S& ishape, T low, U high)
    {
        return this->randint_impl<R>(ishape, low, high);
    }
 
    template <class I, std::size_t L, typename T, typename U>
    auto randint(const I (&ishape)[L], T low, U high) ->
        typename detail::composite_return_type<T, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<T, M, std::array<size_t, L>>::type;
        return this->randint_impl<R>(detail::to_array(ishape), low, high);
    }
 
    template <class R, class I, std::size_t L, typename T, typename U>
    R randint(const I (&ishape)[L], T low, U high)
    {
        return this->randint_impl<R>(detail::to_array(ishape), low, high);
    }
 
    template <class S>
    auto delta(const S& ishape, double scale = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->delta_impl<R>(ishape, scale);
    }
 
    template <class R, class S>
    R delta(const S& ishape, double scale = 1)
    {
        return this->delta_impl<R>(ishape, scale);
    }
 
    template <class I, std::size_t L>
    auto delta(const I (&ishape)[L], double scale = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->delta_impl<R>(detail::to_array(ishape), scale);
    }
 
    template <class R, class I, std::size_t L>
    R delta(const I (&ishape)[L], double scale = 1)
    {
        return this->delta_impl<R>(detail::to_array(ishape), scale);
    }
 
    template <class S>
    auto exponential(const S& ishape, double scale = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->exponential_impl<R>(ishape, scale);
    }
 
    template <class R, class S>
    R exponential(const S& ishape, double scale = 1)
    {
        return this->exponential_impl<R>(ishape, scale);
    }
 
    template <class I, std::size_t L>
    auto exponential(const I (&ishape)[L], double scale = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->exponential_impl<R>(detail::to_array(ishape), scale);
    }
 
    template <class R, class I, std::size_t L>
    R exponential(const I (&ishape)[L], double scale = 1)
    {
        return this->exponential_impl<R>(detail::to_array(ishape), scale);
    }
 
    template <class S>
    auto power(const S& ishape, double k = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->power_impl<R>(ishape, k);
    }
 
    template <class R, class S>
    R power(const S& ishape, double k = 1)
    {
        return this->power_impl<R>(ishape, k);
    }
 
    template <class I, std::size_t L>
    auto power(const I (&ishape)[L], double k = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->power_impl<R>(detail::to_array(ishape), k);
    }
 
    template <class R, class I, std::size_t L>
    R power(const I (&ishape)[L], double k = 1)
    {
        return this->power_impl<R>(detail::to_array(ishape), k);
    }
 
    template <class S>
    auto gamma(const S& ishape, double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->gamma_impl<R>(ishape, k, scale);
    }
 
    template <class R, class S>
    R gamma(const S& ishape, double k = 1, double scale = 1)
    {
        return this->gamma_impl<R>(ishape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto gamma(const I (&ishape)[L], double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->gamma_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R gamma(const I (&ishape)[L], double k = 1, double scale = 1)
    {
        return this->gamma_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class S>
    auto pareto(const S& ishape, double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->pareto_impl<R>(ishape, k, scale);
    }
 
    template <class R, class S>
    R pareto(const S& ishape, double k = 1, double scale = 1)
    {
        return this->pareto_impl<R>(ishape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto pareto(const I (&ishape)[L], double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->pareto_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R pareto(const I (&ishape)[L], double k = 1, double scale = 1)
    {
        return this->pareto_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class S>
    auto weibull(const S& ishape, double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->weibull_impl<R>(ishape, k, scale);
    }
 
    template <class R, class S>
    R weibull(const S& ishape, double k = 1, double scale = 1)
    {
        return this->weibull_impl<R>(ishape, k, scale);
    }
 
    template <class I, std::size_t L>
    auto weibull(const I (&ishape)[L], double k = 1, double scale = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->weibull_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class R, class I, std::size_t L>
    R weibull(const I (&ishape)[L], double k = 1, double scale = 1)
    {
        return this->weibull_impl<R>(detail::to_array(ishape), k, scale);
    }
 
    template <class S>
    auto normal(const S& ishape, double mu = 0, double sigma = 1) ->
        typename detail::composite_return_type<double, M, S>::type
    {
        using R = typename detail::composite_return_type<double, M, S>::type;
        return this->normal_impl<R>(ishape, mu, sigma);
    }
 
    template <class R, class S>
    R normal(const S& ishape, double mu = 0, double sigma = 1)
    {
        return this->normal_impl<R>(ishape, mu, sigma);
    }
 
    template <class I, std::size_t L>
    auto normal(const I (&ishape)[L], double mu = 0, double sigma = 1) ->
        typename detail::composite_return_type<double, M, std::array<size_t, L>>::type
    {
        using R = typename detail::composite_return_type<double, M, std::array<size_t, L>>::type;
        return this->normal_impl<R>(detail::to_array(ishape), mu, sigma);
    }
 
    template <class R, class I, std::size_t L>
    R normal(const I (&ishape)[L], double mu = 0, double sigma = 1)
    {
        return this->normal_impl<R>(detail::to_array(ishape), mu, sigma);
    }
 
    template <class T>
    auto cumsum_random(const T& n) -> typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_random_impl(ret.data(), n.data());
        return ret;
    }
 
    template <class R, class T>
    R cumsum_random(const T& n)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_random_impl(ret.data(), n.data());
        return ret;
    }
 
    template <class T>
    auto cumsum_delta(const T& n, double scale = 1) -> typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        for (size_t i = 0; i < ret.size(); ++i) {
            ret.flat(i) = static_cast<double>(n.flat(i)) * scale;
        }
        return ret;
    }
 
    template <class R, class T>
    R cumsum_delta(const T& n, double scale = 1)
    {
        R ret = R::from_shape(m_shape);
        for (size_t i = 0; i < ret.size(); ++i) {
            ret.flat(i) = static_cast<double>(n.flat(i)) * scale;
        }
        return ret;
    }
 
    template <class T>
    auto cumsum_exponential(const T& n, double scale = 1) ->
        typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_exponential_impl(ret.data(), n.data(), scale);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_exponential(const T& n, double scale = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_exponential_impl(ret.data(), n.data(), scale);
        return ret;
    }
 
    template <class T>
    auto cumsum_power(const T& n, double k = 1) -> typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_power_impl(ret.data(), n.data(), k);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_power(const T& n, double k = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_power_impl(ret.data(), n.data(), k);
        return ret;
    }
 
    template <class T>
    auto cumsum_gamma(const T& n, double k = 1, double scale = 1) ->
        typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_gamma_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_gamma(const T& n, double k = 1, double scale = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_gamma_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class T>
    auto cumsum_pareto(const T& n, double k = 1, double scale = 1) ->
        typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_pareto_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_pareto(const T& n, double k = 1, double scale = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_pareto_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class T>
    auto cumsum_weibull(const T& n, double k = 1, double scale = 1) ->
        typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_weibull_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_weibull(const T& n, double k = 1, double scale = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_weibull_impl(ret.data(), n.data(), k, scale);
        return ret;
    }
 
    template <class T>
    auto cumsum_normal(const T& n, double mu = 0, double sigma = 1) ->
        typename detail::return_type<double, M>::type
    {
        using R = typename detail::return_type<double, M>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_normal_impl(ret.data(), n.data(), mu, sigma);
        return ret;
    }
 
    template <class R, class T>
    R cumsum_normal(const T& n, double mu = 0, double sigma = 1)
    {
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->cumsum_normal_impl(ret.data(), n.data(), mu, sigma);
        return ret;
    }
 
    template <class P>
    auto decide(const P& p) -> typename detail::return_type<bool, P>::type
    {
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        using R = typename detail::return_type<bool, P>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->decide_impl(p.data(), ret.data());
        return ret;
    }
 
    template <class P, class R>
    R decide(const P& p)
    {
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->decide_impl(p.data(), ret.data());
        return ret;
    }
 
    template <class P, class R>
    void decide(const P& p, R& ret)
    {
        static_assert(
            std::is_same<typename R::value_type, bool>::value, "Return value_type must be bool"
        );
 
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        PRRNG_ASSERT(xt::has_shape(p, ret.shape()));
        static_cast<Derived*>(this)->decide_impl(p.data(), ret.data());
    }
 
    template <class P, class T>
    auto decide_masked(const P& p, const T& mask) -> typename detail::return_type<bool, P>::type
    {
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        PRRNG_ASSERT(xt::has_shape(p, mask.shape()));
        using R = typename detail::return_type<bool, P>::type;
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->decide_masked_impl(p.data(), mask.data(), ret.data());
        return ret;
    }
 
    template <class P, class T, class R>
    R decide_masked(const P& p, const T& mask)
    {
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        PRRNG_ASSERT(xt::has_shape(p, mask.shape()));
        R ret = R::from_shape(m_shape);
        static_cast<Derived*>(this)->decide_masked_impl(p.data(), mask.data(), ret.data());
        return ret;
    }
 
    template <class P, class T, class R>
    void decide_masked(const P& p, const T& mask, R& ret)
    {
        static_assert(
            std::is_same<typename R::value_type, bool>::value, "Return value_type must be bool"
        );
 
        PRRNG_ASSERT(xt::has_shape(p, m_shape));
        PRRNG_ASSERT(xt::has_shape(p, mask.shape()));
        PRRNG_ASSERT(xt::has_shape(p, ret.shape()));
        static_cast<Derived*>(this)->decide_masked_impl(p.data(), mask.data(), ret.data());
    }
 
private:
    template <class R, class S>
    R positive_random_impl(const S& ishape)
    {
        static_assert(
            std::is_same<typename R::value_type, double>::value, "Return value_type must be double"
        );
 
        auto n = detail::size(ishape);
        R ret = R::from_shape(detail::concatenate<M, S>::two(m_shape, ishape));
        static_cast<Derived*>(this)->draw_list_positive_double(&ret.front(), n);
        return ret;
    }
 
    template <class R, class S>
    R random_impl(const S& ishape)
    {
        static_assert(
            std::is_same<typename R::value_type, double>::value, "Return value_type must be double"
        );
 
        auto n = detail::size(ishape);
        R ret = R::from_shape(detail::concatenate<M, S>::two(m_shape, ishape));
        static_cast<Derived*>(this)->draw_list_double(&ret.front(), n);
        return ret;
    }
 
    template <class R, class S, typename T>
    R randint_impl(const S& ishape, T high)
    {
        static_assert(
            std::numeric_limits<typename R::value_type>::max() >= std::numeric_limits<T>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<T>::max() <= std::numeric_limits<uint32_t>::max(), "Bound too large"
        );
 
        auto n = detail::size(ishape);
        R ret = R::from_shape(detail::concatenate<M, S>::two(m_shape, ishape));
        std::vector<uint32_t> tmp(ret.size());
        static_cast<Derived*>(this)->draw_list_uint32(&tmp.front(), static_cast<uint32_t>(high), n);
        std::copy(tmp.begin(), tmp.end(), ret.begin());
        return ret;
    }
 
    template <class R, class S, typename T, typename U>
    R randint_impl(const S& ishape, T low, U high)
    {
        static_assert(
            std::numeric_limits<typename R::value_type>::max() >= std::numeric_limits<T>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename R::value_type>::min() >= std::numeric_limits<T>::min(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename R::value_type>::max() >= std::numeric_limits<U>::max(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            std::numeric_limits<typename R::value_type>::min() >= std::numeric_limits<U>::min(),
            "Return value_type must must be able to accommodate the bound"
        );
 
        static_assert(
            static_cast<uint32_t>(std::numeric_limits<T>::max()) <
                std::numeric_limits<uint32_t>::max(),
            "Bound too large"
        );
 
        auto n = detail::size(ishape);
        R ret = R::from_shape(detail::concatenate<M, S>::two(m_shape, ishape));
        std::vector<uint32_t> tmp(ret.size());
        static_cast<Derived*>(this)->draw_list_uint32(
            &tmp.front(), static_cast<uint32_t>(high - low), n
        );
        std::copy(tmp.begin(), tmp.end(), ret.begin());
        return ret + low;
    }
 
    template <class R, class S>
    R delta_impl(const S& ishape, double scale)
    {
        R ret = R::from_shape(detail::concatenate<M, S>::two(m_shape, ishape));
        ret.fill(scale);
        return ret;
    }
 
    template <class R, class S>
    R exponential_impl(const S& ishape, double scale)
    {
        R r = this->random_impl<R>(ishape);
        return exponential_distribution(scale).quantile(r);
    }
 
    template <class R, class S>
    R power_impl(const S& ishape, double k)
    {
        R r = this->random_impl<R>(ishape);
        return power_distribution(k).quantile(r);
    }
 
    template <class R, class S>
    R gamma_impl(const S& ishape, double k, double scale)
    {
        R r = this->random_impl<R>(ishape);
        return gamma_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R pareto_impl(const S& ishape, double k, double scale)
    {
        R r = this->random_impl<R>(ishape);
        return pareto_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R weibull_impl(const S& ishape, double k, double scale)
    {
        R r = this->random_impl<R>(ishape);
        return weibull_distribution(k, scale).quantile(r);
    }
 
    template <class R, class S>
    R normal_impl(const S& ishape, double mu, double sigma)
    {
        R r = this->positive_random_impl<R>(ishape);
        return normal_distribution(mu, sigma).quantile(r);
    }
 
protected:
    size_type m_size = 0; 
    shape_type m_shape; 
    shape_type m_strides; 
};
 
template <class Generator, class Shape>
class pcg32_arrayBase : public GeneratorBase_array<pcg32_arrayBase<Generator, Shape>, Shape> {
    friend GeneratorBase_array<pcg32_arrayBase<Generator, Shape>, Shape>;
 
private:
    using derived_type = pcg32_arrayBase<Generator, Shape>;
 
public:
    using size_type = typename Shape::value_type; 
    using shape_type = Shape; 
 
protected:
    template <class T>
    void init(const T& initstate)
    {
        std::copy(initstate.shape().cbegin(), initstate.shape().cend(), m_shape.begin());
        std::copy(initstate.strides().cbegin(), initstate.strides().cend(), m_strides.begin());
        m_size = initstate.size();
        m_gen.reserve(m_size);
 
        for (size_type i = 0; i < m_size; ++i) {
            m_gen.push_back(Generator(initstate.flat(i)));
        }
    }
 
    template <class T, class U>
    void init(const T& initstate, const U& initseq)
    {
        PRRNG_ASSERT(xt::has_shape(initstate, initseq.shape()));
 
        std::copy(initstate.shape().cbegin(), initstate.shape().cend(), m_shape.begin());
        std::copy(initstate.strides().cbegin(), initstate.strides().cend(), m_strides.begin());
        m_size = initstate.size();
        m_gen.reserve(m_size);
 
        for (size_type i = 0; i < m_size; ++i) {
            m_gen.push_back(Generator(initstate.flat(i), initseq.flat(i)));
        }
    }
 
public:
    pcg32_arrayBase() = default;
 
    template <class... Args>
    Generator& operator()(Args... args)
    {
        return m_gen[this->get_item(0, 0, args...)];
    }
 
    template <class... Args>
    const Generator& operator()(Args... args) const
    {
        return m_gen[this->get_item(0, 0, args...)];
    }
 
    Generator& operator[](size_t i)
    {
        PRRNG_DEBUG(i < m_size);
        return m_gen[i];
    }
 
    const Generator& operator[](size_t i) const
    {
        PRRNG_DEBUG(i < m_size);
        return m_gen[i];
    }
 
    Generator& flat(size_t i)
    {
        PRRNG_DEBUG(i < m_size);
        return m_gen[i];
    }
 
    const Generator& flat(size_t i) const
    {
        PRRNG_DEBUG(i < m_size);
        return m_gen[i];
    }
 
    auto state() -> typename detail::return_type<uint64_t, Shape>::type
    {
        using R = typename detail::return_type<uint64_t, Shape>::type;
        return this->state<R>();
    }
 
    template <class R>
    R state()
    {
        using value_type = typename R::value_type;
        R ret = R::from_shape(m_shape);
 
        for (size_type i = 0; i < m_size; ++i) {
            ret.flat(i) = m_gen[i].template state<value_type>();
        }
 
        return ret;
    }
 
    auto initstate() -> typename detail::return_type<uint64_t, Shape>::type
    {
        using R = typename detail::return_type<uint64_t, Shape>::type;
        return this->initstate<R>();
    }
 
    template <class R>
    R initstate()
    {
        using value_type = typename R::value_type;
        R ret = R::from_shape(m_shape);
 
        for (size_type i = 0; i < m_size; ++i) {
            ret.flat(i) = m_gen[i].template initstate<value_type>();
        }
 
        return ret;
    }
 
    auto initseq() -> typename detail::return_type<uint64_t, Shape>::type
    {
        using R = typename detail::return_type<uint64_t, Shape>::type;
        return this->initseq<R>();
    }
 
    template <class R>
    R initseq()
    {
        using value_type = typename R::value_type;
        R ret = R::from_shape(m_shape);
 
        for (size_type i = 0; i < m_size; ++i) {
            ret.flat(i) = m_gen[i].template initseq<value_type>();
        }
 
        return ret;
    }
 
    template <class T>
    auto distance(const T& arg) -> typename detail::return_type<int64_t, Shape>::type
    {
        using R = typename detail::return_type<int64_t, Shape>::type;
        return this->distance<R, T>(arg);
    }
 
    template <class R, class T>
    R distance(const T& arg)
    {
        using value_type = typename R::value_type;
        R ret = R::from_shape(m_shape);
 
        for (size_type i = 0; i < m_size; ++i) {
            ret.flat(i) = m_gen[i].template distance<value_type>(arg.flat(i));
        }
 
        return ret;
    }
 
    template <class T>
    void advance(const T& arg)
    {
        for (size_type i = 0; i < m_size; ++i) {
            m_gen[i].advance(arg.flat(i));
        }
    }
 
    template <class T>
    void restore(const T& arg)
    {
        for (size_type i = 0; i < m_size; ++i) {
            m_gen[i].restore(arg.flat(i));
        }
    }
 
protected:
    void decide_impl(const double* p, bool* ret)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].next_double() < p[i];
        }
    }
 
    void decide_masked_impl(const double* p, const bool* mask, bool* ret)
    {
        for (size_type i = 0; i < m_size; ++i) {
            if (mask[i]) {
                ret[i] = false;
            }
            else {
                ret[i] = m_gen[i].next_double() < p[i];
            }
        }
    }
 
    void cumsum_random_impl(double* ret, const size_t* n)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_random(n[i]);
        }
    }
 
    void cumsum_exponential_impl(double* ret, const size_t* n, double scale)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_exponential(n[i], scale);
        }
    }
 
    void cumsum_power_impl(double* ret, const size_t* n, double k)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_power(n[i], k);
        }
    }
 
    void cumsum_gamma_impl(double* ret, const size_t* n, double k, double scale)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_gamma(n[i], k, scale);
        }
    }
 
    void cumsum_pareto_impl(double* ret, const size_t* n, double k, double scale)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_pareto(n[i], k, scale);
        }
    }
 
    void cumsum_weibull_impl(double* ret, const size_t* n, double k, double scale)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_weibull(n[i], k, scale);
        }
    }
 
    void cumsum_normal_impl(double* ret, const size_t* n, double mu, double sigma)
    {
        for (size_type i = 0; i < m_size; ++i) {
            ret[i] = m_gen[i].cumsum_normal(n[i], mu, sigma);
        }
    }
 
    void draw_list_double(double* data, size_t n)
    {
        for (size_type i = 0; i < m_size; ++i) {
            for (size_t j = 0; j < n; ++j) {
                data[i * n + j] = m_gen[i].next_double();
            }
        }
    }
 
    void draw_list_positive_double(double* data, size_t n)
    {
        for (size_type i = 0; i < m_size; ++i) {
            for (size_t j = 0; j < n; ++j) {
                data[i * n + j] = m_gen[i].next_positive_double();
            }
        }
    }
 
    void draw_list_uint32(uint32_t* data, uint32_t bound, size_t n)
    {
        for (size_type i = 0; i < m_size; ++i) {
            for (size_t j = 0; j < n; ++j) {
                data[i * n + j] = m_gen[i].next_uint32(bound);
            }
        }
    }
 
private:
    template <class T>
    size_t get_item(size_t sum, size_t d, T arg)
    {
        return sum + arg * m_strides[d];
    }
 
    template <class T, class... Args>
    size_t get_item(size_t sum, size_t d, T arg, Args... args)
    {
        return get_item(sum + arg * m_strides[d], d + 1, args...);
    }
 
protected:
    std::vector<Generator> m_gen; 
    using GeneratorBase_array<derived_type, Shape>::m_size;
    using GeneratorBase_array<derived_type, Shape>::m_shape;
    using GeneratorBase_array<derived_type, Shape>::m_strides;
};
 
class pcg32_array : public pcg32_arrayBase<pcg32, std::vector<size_t>> {
private:
    using derived_type = pcg32_arrayBase<pcg32, std::vector<size_t>>;
 
public:
    using size_type = size_t; 
    using shape_type = std::vector<size_t>; 
 
    pcg32_array() = default;
 
    template <class T>
    pcg32_array(const T& initstate)
    {
        m_shape.resize(initstate.dimension());
        m_strides.resize(initstate.dimension());
        this->init(initstate);
    }
 
    template <class T, class U>
    pcg32_array(const T& initstate, const U& initseq)
    {
        m_shape.resize(initstate.dimension());
        m_strides.resize(initstate.dimension());
        this->init(initstate, initseq);
    }
 
protected:
    using pcg32_arrayBase<pcg32, shape_type>::m_gen;
    using GeneratorBase_array<derived_type, shape_type>::m_size;
    using GeneratorBase_array<derived_type, shape_type>::m_shape;
    using GeneratorBase_array<derived_type, shape_type>::m_strides;
};
 
template <size_t N>
class pcg32_tensor : public pcg32_arrayBase<pcg32, std::array<size_t, N>> {
private:
    using derived_type = pcg32_arrayBase<pcg32, std::array<size_t, N>>;
 
public:
    using size_type = size_t; 
    using shape_type = std::array<size_t, N>; 
 
    pcg32_tensor() = default;
 
    template <class T>
    pcg32_tensor(const T& initstate)
    {
        static_assert(detail::check_fixed_rank<N, T>::value, "Ranks to not match");
        this->init(initstate);
    }
 
    template <class T, class U>
    pcg32_tensor(const T& initstate, const U& initseq)
    {
        static_assert(detail::check_fixed_rank<N, T>::value, "Ranks to not match");
        this->init(initstate, initseq);
    }
 
protected:
    using pcg32_arrayBase<pcg32, shape_type>::m_gen;
    using GeneratorBase_array<derived_type, shape_type>::m_size;
    using GeneratorBase_array<derived_type, shape_type>::m_shape;
    using GeneratorBase_array<derived_type, shape_type>::m_strides;
};
 
class pcg32_index_array : public pcg32_arrayBase<pcg32_index, std::vector<size_t>> {
private:
    using derived_type = pcg32_arrayBase<pcg32_index, std::vector<size_t>>;
 
public:
    pcg32_index_array() = default;
 
    template <class T, class U>
    pcg32_index_array(const T& initstate, const U& initseq)
    {
        m_shape.resize(initstate.dimension());
        m_strides.resize(initstate.dimension());
        this->init(initstate, initseq);
    }
 
protected:
    using pcg32_arrayBase<pcg32_index, std::vector<size_t>>::m_gen;
    using GeneratorBase_array<derived_type, std::vector<size_t>>::m_size;
    using GeneratorBase_array<derived_type, std::vector<size_t>>::m_shape;
    using GeneratorBase_array<derived_type, std::vector<size_t>>::m_strides;
};
 
template <size_t N>
class pcg32_index_tensor : public pcg32_arrayBase<pcg32_index, std::array<size_t, N>> {
private:
    using derived_type = pcg32_arrayBase<pcg32_index, std::array<size_t, N>>;
 
public:
    pcg32_index_tensor() = default;
 
    template <class T, class U>
    pcg32_index_tensor(const T& initstate, const U& initseq)
    {
        static_assert(detail::check_fixed_rank<N, T>::value, "Ranks to not match");
        this->init(initstate, initseq);
    }
 
protected:
    using pcg32_arrayBase<pcg32_index, std::array<size_t, N>>::m_gen;
    using GeneratorBase_array<derived_type, std::array<size_t, N>>::m_size;
    using GeneratorBase_array<derived_type, std::array<size_t, N>>::m_shape;
    using GeneratorBase_array<derived_type, std::array<size_t, N>>::m_strides;
};
 
namespace detail {
 
template <class T, typename = void>
struct auto_pcg32 {
    static auto get(const T& initseq)
    {
        return pcg32_array(initseq);
    }
 
    template <class S>
    static auto get(const T& initseq, const S& initstate)
    {
        return pcg32_array(initseq, initstate);
    }
};
 
template <class T>
struct auto_pcg32<T, typename std::enable_if_t<xt::has_fixed_rank_t<T>::value>> {
    static auto get(const T& initseq)
    {
        return pcg32_tensor<xt::get_rank<T>::value>(initseq);
    }
 
    template <class S>
    static auto get(const T& initseq, const S& initstate)
    {
        return pcg32_tensor<xt::get_rank<T>::value>(initseq, initstate);
    }
};
 
template <class T>
struct auto_pcg32<T, typename std::enable_if_t<std::is_integral<T>::value>> {
    static auto get(const T& initseq)
    {
        return pcg32(initseq);
    }
 
    template <class S>
    static auto get(const T& initseq, const S& initstate)
    {
        return pcg32(initseq, initstate);
    }
};
 
} // namespace detail
 
template <class T>
inline auto auto_pcg32(const T& initstate)
{
    return detail::auto_pcg32<T>::get(initstate);
}
 
template <class T, class S>
inline auto auto_pcg32(const T& initstate, const S& initseq)
{
    return detail::auto_pcg32<T>::get(initstate, initseq);
}
 
template <class Generator, class Data, class Index, bool is_cumsum>
class pcg32_arrayBase_chunkBase {
    static_assert(std::is_signed<typename Index::value_type>::value, "Index must be signed");
 
public:
    using size_type = typename Data::size_type; 
 
protected:
    Generator m_gen; 
    Data m_data; 
 
    std::vector<std::function<xt::xtensor<double, 1>(size_t)>> m_draw;
 
    std::vector<std::function<double(size_t)>> m_sum;
 
    bool m_extendible;
 
    alignment m_align; 
    distribution m_distro; 
    std::array<double, 3> m_param; 
    Index m_start; 
    Index m_i; 
    size_t m_n; 
 
protected:
    template <class S, class T, class U>
    void init(
        const S& shape,
        const T& initstate,
        const U& initseq,
        enum distribution distribution,
        const std::vector<double>& parameters,
        const alignment& align = alignment()
    )
    {
        PRRNG_ASSERT(xt::has_shape(initstate, initseq.shape()));
 
        m_align = align;
        m_distro = distribution;
        m_gen = Generator(initstate, initseq);
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            m_gen[i].set_delta(distribution == distribution::delta);
        }
 
        std::vector<size_t> data_shape;
        data_shape.resize(initstate.dimension() + shape.size());
        std::copy(initstate.shape().begin(), initstate.shape().end(), data_shape.begin());
        std::copy(shape.begin(), shape.end(), data_shape.begin() + initstate.dimension());
        m_data = xt::empty<typename Data::value_type>(data_shape);
        m_n = m_data.size() / initstate.size();
 
        m_start = xt::zeros<typename Index::value_type>(m_gen.shape());
        m_i = m_n * xt::ones<typename Index::value_type>(m_gen.shape());
 
        auto par = default_parameters(distribution, parameters);
        std::copy(par.begin(), par.end(), m_param.begin());
 
        this->auto_functions();
 
        using E = decltype(m_draw[size_t{}](size_t{}));
 
        // if possible: draw the first chunk
        if (m_extendible) {
            for (size_t i = 0; i < m_gen.size(); ++i) {
                E extra = m_draw[i](m_n);
                m_gen[i].drawn(m_n);
                if constexpr (is_cumsum) {
                    std::partial_sum(extra.begin(), extra.end(), &m_data.flat(i * m_n));
                }
                else {
                    std::copy(extra.begin(), extra.end(), &m_data.flat(i * m_n));
                }
            }
        }
    }
 
    void auto_functions()
    {
        using R = decltype(m_draw[size_t{}](size_t{}));
        m_extendible = true;
 
        if (m_distro != distribution::custom) {
            m_draw.resize(m_gen.size());
            if constexpr (is_cumsum) {
                m_sum.resize(m_gen.size());
            }
        }
 
        switch (m_distro) {
        case random:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template random<R>(std::array<size_t, 1>{n}) * m_param[0] +
                           m_param[1];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_random(n) * m_param[0] +
                               static_cast<double>(n) * m_param[1];
                    };
                }
            }
            return;
        case delta:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template delta<R>(std::array<size_t, 1>{n}, m_param[0]) +
                           m_param[1];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_delta(n, m_param[0]) +
                               static_cast<double>(n) * m_param[1];
                    };
                }
            }
            return;
        case exponential:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template exponential<R>(std::array<size_t, 1>{n}, m_param[0]) +
                           m_param[1];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_exponential(n, m_param[0]) +
                               static_cast<double>(n) * m_param[1];
                    };
                }
            }
            return;
        case power:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template power<R>(std::array<size_t, 1>{n}, m_param[0]) +
                           m_param[1];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_power(n, m_param[0]) +
                               static_cast<double>(n) * m_param[1];
                    };
                }
            }
            return;
        case gamma:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template gamma<R>(
                               std::array<size_t, 1>{n}, m_param[0], m_param[1]
                           ) +
                           m_param[2];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_gamma(n, m_param[0], m_param[1]) +
                               static_cast<double>(n) * m_param[2];
                    };
                }
            }
            return;
        case pareto:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template pareto<R>(
                               std::array<size_t, 1>{n}, m_param[0], m_param[1]
                           ) +
                           m_param[2];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_pareto(n, m_param[0], m_param[1]) +
                               static_cast<double>(n) * m_param[2];
                    };
                }
            }
            return;
        case weibull:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template weibull<R>(
                               std::array<size_t, 1>{n}, m_param[0], m_param[1]
                           ) +
                           m_param[2];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_weibull(n, m_param[0], m_param[1]) +
                               static_cast<double>(n) * m_param[2];
                    };
                }
            }
            return;
        case normal:
            for (size_t i = 0; i < m_gen.size(); ++i) {
                m_draw[i] = [this, i](size_t n) -> R {
                    return m_gen[i].template normal<R>(
                               std::array<size_t, 1>{n}, m_param[0], m_param[1]
                           ) +
                           m_param[2];
                };
                if constexpr (is_cumsum) {
                    m_sum[i] = [this, i](size_t n) -> double {
                        return m_gen[i].cumsum_normal(n, m_param[0], m_param[1]) +
                               static_cast<double>(n) * m_param[2];
                    };
                }
            }
            return;
        case custom:
            m_extendible = false;
            return;
        }
    }
 
    void copy_from(const pcg32_arrayBase_chunkBase& other)
    {
        m_gen = other.m_gen;
        m_data = other.m_data;
        m_align = other.m_align;
        m_distro = other.m_distro;
        m_param = other.m_param;
        m_start = other.m_start;
        m_i = other.m_i;
        m_n = other.m_n;
        this->auto_functions();
    }
 
public:
    pcg32_arrayBase_chunkBase() = default;
 
    pcg32_arrayBase_chunkBase(const pcg32_arrayBase_chunkBase& other)
    {
        this->copy_from(other);
    }
 
    void operator=(const pcg32_arrayBase_chunkBase& other)
    {
        this->copy_from(other);
    }
 
    template <class T>
    pcg32_arrayBase_chunkBase& operator+=(const T& values)
    {
        xt::noalias(m_data) += values;
        return *this;
    }
 
    template <class T>
    pcg32_arrayBase_chunkBase& operator-=(const T& values)
    {
        xt::noalias(m_data) -= values;
        return *this;
    }
 
    bool is_extendible() const
    {
        return m_extendible;
    }
 
    size_type chunk_size() const
    {
        return static_cast<size_type>(m_n);
    }
 
    const Generator& generators() const
    {
        return m_gen;
    }
 
    const Data& data() const
    {
        return m_data;
    }
 
    void set_data(const Data& data)
    {
        PRRNG_ASSERT(xt::has_shape(data, m_data.shape()));
        xt::noalias(m_data) = data;
    }
 
    const Index& start() const
    {
        return m_start;
    }
 
    void set_start(const Index& index)
    {
        PRRNG_ASSERT(xt::has_shape(index, m_gen.shape()));
        xt::noalias(m_start) = index;
    }
 
    void align_at(const Index& index)
    {
        PRRNG_ASSERT(xt::has_shape(index, m_gen.shape()));
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            if constexpr (!is_cumsum) {
                detail::chunk_align_at(
                    m_gen[i],
                    m_draw[i],
                    m_align,
                    &m_data.flat(i * m_n),
                    m_n,
                    &m_start.flat(i),
                    index.flat(i)
                );
            }
            else {
                detail::cumsum_align_at(
                    m_gen[i],
                    m_draw[i],
                    m_sum[i],
                    m_align,
                    &m_data.flat(i * m_n),
                    m_n,
                    &m_start.flat(i),
                    index.flat(i)
                );
            }
        }
 
        xt::noalias(m_i) = index - m_start;
    }
 
    Index index_at_align() const
    {
        return m_start + m_i;
    }
 
    const Index& chunk_index_at_align() const
    {
        return m_i;
    }
 
    template <class R>
    void left_of_align(R& ret) const
    {
        PRRNG_ASSERT(xt::has_shape(ret, m_gen.shape()));
        using value_type = typename R::value_type;
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            ret.flat(i) = static_cast<value_type>(m_data.flat(i * m_n + m_i.flat(i)));
        }
    }
 
    template <class R>
    void right_of_align(R& ret) const
    {
        PRRNG_ASSERT(xt::has_shape(ret, m_gen.shape()));
        using value_type = typename R::value_type;
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            ret.flat(i) = static_cast<value_type>(m_data.flat(i * m_n + m_i.flat(i) + 1));
        }
    }
 
    template <class R>
    R left_of_align() const
    {
        R ret = R::from_shape(m_gen.shape());
        this->left_of_align(ret);
        return ret;
    }
 
    template <class R>
    R right_of_align() const
    {
        R ret = R::from_shape(m_gen.shape());
        this->right_of_align(ret);
        return ret;
    }
 
    template <class R, class T>
    R state_at(const T& index)
    {
        PRRNG_ASSERT(xt::has_shape(index, m_gen.shape()));
 
        using value_type = typename R::value_type;
        R ret = R::from_shape(m_gen.shape());
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            ret.flat(i) = static_cast<value_type>(m_gen[i].state_at(index.flat(i)));
        }
 
        return ret;
    }
};
 
template <class Generator, class Data, class Index>
class pcg32_arrayBase_chunk : public pcg32_arrayBase_chunkBase<Generator, Data, Index, false> {
protected:
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_data;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_draw;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_gen;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_n;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_i;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_start;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, false>::m_align;
 
public:
    using size_type = typename Data::size_type; 
    using value_type = typename Data::value_type; 
 
    pcg32_arrayBase_chunk() = default;
 
    template <class S, class T>
    void restore(const S& state, const T& index)
    {
        PRRNG_ASSERT(xt::has_shape(state, m_gen.shape()));
        PRRNG_ASSERT(xt::has_shape(index, m_gen.shape()));
        xt::noalias(m_start) = index;
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            m_gen[i].set_index(index.flat(i));
            m_gen[i].restore(state.flat(i));
 
            using E = decltype(m_draw[i](size_t{}));
            E extra = m_draw[i](m_n);
            m_gen[i].drawn(m_n);
            std::copy(extra.begin(), extra.end(), &m_data.flat(i * m_n));
        }
    }
};
 
template <class Generator, class Data, class Index>
class pcg32_arrayBase_cumsum : public pcg32_arrayBase_chunkBase<Generator, Data, Index, true> {
protected:
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_align;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_data;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_draw;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_extendible;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_gen;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_i;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_n;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_start;
    using pcg32_arrayBase_chunkBase<Generator, Data, Index, true>::m_sum;
 
public:
    using size_type = typename Data::size_type; 
 
    pcg32_arrayBase_cumsum() = default;
 
    // TODO: rename align_at_value ?
    template <class T>
    void align(const T& target)
    {
        PRRNG_ASSERT(xt::has_shape(target, m_gen.shape()));
 
        if (!m_extendible) {
            PRRNG_ASSERT(this->contains(target));
            inplace::lower_bound(m_data, target, m_i);
            return;
        }
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            detail::align(
                m_gen[i],
                m_draw[i],
                m_sum[i],
                m_align,
                &m_data.flat(i * m_n),
                m_n,
                &m_start.flat(i),
                &m_i.flat(i),
                target.flat(i)
            );
        }
    }
 
    // todo: overload with array index
 
    void align(size_t i, double target)
    {
        if (!m_extendible) {
            PRRNG_ASSERT(
                target >= m_data.flat(i * m_n) && target <= m_data.flat((i + 1) * m_n - 1)
            );
            m_i.flat(i) = iterator::lower_bound(
                &m_data.flat(i * m_n), &m_data.flat(i * m_n) + m_n, target, m_i.flat(i)
            );
            return;
        }
 
        detail::align(
            m_gen[i],
            m_draw[i],
            m_sum[i],
            m_align,
            &m_data.flat(i * m_n),
            m_n,
            &m_start.flat(i),
            &m_i.flat(i),
            target
        );
    }
 
    template <class S, class V, class T>
    void restore(const S& state, const V& value, const T& index)
    {
        PRRNG_ASSERT(xt::has_shape(state, m_gen.shape()));
        PRRNG_ASSERT(xt::has_shape(value, m_gen.shape()));
        PRRNG_ASSERT(xt::has_shape(index, m_gen.shape()));
        xt::noalias(m_start) = index;
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            m_gen[i].set_index(index.flat(i));
            m_gen[i].restore(state.flat(i));
 
            using E = decltype(m_draw[i](size_t{}));
            E extra = m_draw[i](m_n);
            m_gen[i].drawn(m_n);
            extra.front() += value.flat(i) - extra.front();
            std::partial_sum(extra.begin(), extra.end(), &m_data.flat(i * m_n));
        }
    }
 
    template <class T>
    bool contains(const T& target) const
    {
        PRRNG_ASSERT(xt::has_shape(target, m_gen.shape()));
 
        for (size_t i = 0; i < m_gen.size(); ++i) {
            if (target.flat(i) < m_data.flat(i * m_n) ||
                target.flat(i) > m_data.flat((i + 1) * m_n - 1)) {
                return false;
            }
        }
 
        return true;
    }
};
 
template <class Data, class Index>
class pcg32_array_chunk : public pcg32_arrayBase_chunk<pcg32_index_array, Data, Index> {
public:
    pcg32_array_chunk() = default;
 
    template <class S, class T, class U>
    pcg32_array_chunk(
        const S& shape,
        const T& initstate,
        const U& initseq,
        enum distribution distribution,
        const std::vector<double>& parameters,
        const alignment& align = alignment()
    )
    {
        this->init(shape, initstate, initseq, distribution, parameters, align);
    }
};
 
template <class Data, class Index, size_t N>
class pcg32_tensor_chunk : public pcg32_arrayBase_chunk<pcg32_index_tensor<N>, Data, Index> {
public:
    pcg32_tensor_chunk() = default;
 
    template <class S, class T, class U>
    pcg32_tensor_chunk(
        const S& shape,
        const T& initstate,
        const U& initseq,
        enum distribution distribution,
        const std::vector<double>& parameters,
        const alignment& align = alignment()
    )
    {
        this->init(shape, initstate, initseq, distribution, parameters, align);
    }
};
 
template <class Data, class Index>
class pcg32_array_cumsum : public pcg32_arrayBase_cumsum<pcg32_index_array, Data, Index> {
public:
    pcg32_array_cumsum() = default;
 
    template <class S, class T, class U>
    pcg32_array_cumsum(
        const S& shape,
        const T& initstate,
        const U& initseq,
        enum distribution distribution,
        const std::vector<double>& parameters,
        const alignment& align = alignment()
    )
    {
        this->init(shape, initstate, initseq, distribution, parameters, align);
    }
};
 
template <class Data, class Index, size_t N>
class pcg32_tensor_cumsum : public pcg32_arrayBase_cumsum<pcg32_index_tensor<N>, Data, Index> {
public:
    pcg32_tensor_cumsum() = default;
 
    template <class S, class T, class U>
    pcg32_tensor_cumsum(
        const S& shape,
        const T& initstate,
        const U& initseq,
        enum distribution distribution,
        const std::vector<double>& parameters,
        const alignment& align = alignment()
    )
    {
        this->init(shape, initstate, initseq, distribution, parameters, align);
    }
};
 
} // namespace prrng
 
#endif