static_assert(std::is_unsigned_v<typename R::result_type>, "random engine result_type must be unsigned");

static_assert(R::max() > R::min(), "random engine must have non-zero range");

// why the fuck is std::mersenne_twister_engine::operator() noexcept(false)??

// static_assert(noexcept(std::declval<R&>()()), "random engine operator() must be noexcept");

static_assert(std::is_unsigned_v<U>, "unsigned integer type required");

using result_type = U;

constexpr explicit uniform_uint_max_distribution(U max = std::numeric_limits<U>::max()) noexcept

: m_max(max)

{}

constexpr U min() const noexcept { return 0; }

constexpr U max() const noexcept { return m_max; }

template <typename R>

constexpr static U draw(U max, R& rng) {

impl::do_rng_checks(rng);

using r_t = typename R::result_type;

ITLIB_RAND_DIST_CONSTEXPR r_t rng_range = R::max() - R::min();

if ITLIB_RAND_DIST_CONSTEXPR(rng_range < std::numeric_limits<U>::max()) {

// desired max might be bigger than rng's range

if (max <= U(rng_range)) {

// it fits

return draw_in_rng_range(max, rng);

}

else {

// we must draw multiple times

// we represent max as a number in base (rng_range + 1)

// we draw digit by digit, keeping track if we are within the limit

// if we go outside the limit, we reject and stat over

ITLIB_RAND_DIST_CONSTEXPR U base = U(rng_range) + 1;

// collect base_(rng_range+1) digits of max

// array size enough even if rng_range is 2

// (we *could* be fancy and compute this more precisely, but stack is cheap)

r_t max_digits[sizeof(U) * 8];

int max_digit_count = 0;

{

U temp = max;

while (temp > 0) {

max_digits[max_digit_count++] = r_t(temp % base);

temp /= base;

}

}

// rejection loop

while (true) {

// start with most significant digit

U result = 0;

bool tight = true;

for (int i = max_digit_count - 1; ; --i) {

const auto random_digit = rng_range_draw(rng); // unconstrained draw

if (tight) {

if (random_digit > max_digits[i]) {

break;

}

if (random_digit < max_digits[i]) {

tight = false; // we are now free to choose any digit

}

}

result = result * base + U(random_digit);

if (i == 0) {

// done

return result;

}

}

}

}

}

else {

// rng range is guaranteed to cover desired max

return draw_in_rng_range(max, rng);

}

if (umax == 0) return 0;

using r_t = typename R::result_type;

ITLIB_RAND_DIST_CONSTEXPR r_t rng_range = R::max() - R::min();

const auto max = r_t(umax);

if (rng_range == max) {

return U(rng_range_draw(rng));

}

// after the check above, now we have a gurantee that max < rng_range

// this means that result_range_size will not overflow and reject_count will not underflow

const auto result_range_size = r_t(max + 1);

// mathematically we should use just rng_range + 1 here, but if rng_range is limits<r_t>::max(),

// it will overflow and wrap to 0 leading to reject_count == 0 and no rejections ever.

const auto reject_count = r_t(rng_range - result_range_size + 1) % result_range_size;

const auto accept_max = r_t(rng_range - reject_count);

while (true) {

if (auto v = rng_range_draw(rng); v <= accept_max) {

return U(v % result_range_size);

}

// food for thought:

// is it worth checking for broken (malicious?) RNG that might lead to an infinite loop here?

}

static_assert(std::is_integral_v<I>, "integral type required");

using result_type = I;

using U = std::make_unsigned_t<I>;

constexpr uniform_int_distribution(I min, I max) noexcept

: m_min(U(min))

, m_range(U(max) - U(min)) // as per the standard: UB if max < min

{}

constexpr explicit uniform_int_distribution(I max = std::numeric_limits<I>::max()) noexcept

: uniform_int_distribution(0, max)

{}

constexpr I min() const noexcept { return I(m_min); }

constexpr I max() const noexcept { return I(U(m_min + m_range.max())); }

template <typename R>

constexpr static I draw(I min, I max, R& rng) {

// multiple seemingly redundant casts to U to handle 8-bit types which auto-promote to int in expressions

const auto umin = U(min);

const auto urange = U(U(max) - umin); // as per the standard: UB if max < min

return I(U(umin + uniform_uint_max_distribution<U>::draw(urange, rng)));

}

template <typename R>

constexpr I operator()(R& rng) const {

return I(U(m_min + m_range(rng)));

}

const U m_min;

const uniform_uint_max_distribution<U> m_range;

static_assert(std::is_floating_point_v<F>, "floating point type required");

static_assert(std::numeric_limits<F>::radix == 2, "only binary floating point types are supported");

static_assert(std::numeric_limits<F>::digits <= 64, "max floating point precision must fit uint64_t");

using result_type = F;

constexpr fast_uniform_real_distribution(F min = F(0), F max = F(1)) noexcept

: m_min(min)

, m_scale(max - min)

{}

constexpr F min() const noexcept { return m_min; }

constexpr F scale() const noexcept { return m_scale; }

constexpr F max() const noexcept { return m_min + m_scale; }

// draws from [0, 1)

template <typename R>

constexpr static F draw_01(R& rng) {

impl::do_rng_checks(rng);

// (1 << d) - 1 is more readable, but might overflow

constexpr uint64_t max_int = ~uint64_t(0) >> (64 - std::numeric_limits<F>::digits);

using r_t = typename R::result_type;

ITLIB_RAND_DIST_CONSTEXPR uint64_t rng_range = R::max() - R::min();

if ITLIB_RAND_DIST_CONSTEXPR(rng_range >= max_int) {

// rng_range is enough to saturate our float precision

// we slice off the needed bits (and hope that rng is uniform over all bits)

const auto random_value = r_t(rng() - R::min()) & r_t(max_int);

// ideally we would use this here:

// return std::ldexp(F(random_value), -std::numeric_limits<F>::digits);

// ... but in reality no compiler would would optimize it well enough

// every single one does a `call ldexp[f]`!

// at best they would do value * exp2[f](-digits) which is literally the same as below

// well, not literally, as exp2 only gets an overload with C++23, but still...

return F(random_value) / F(max_int + 1);

}

else {

// "stretching" to rng_range

// some F values are unreachable

return F(rng() - R::min()) / (F(rng_range) + 1);

}