libdpf/include/grotto/constant_lut.hpp

/// @file grotto/constant_lut.hpp
/// @brief Exact piecewise-constant lookup tables for the degree-0 gadgets
///        in Appendix D of Storrier, Vadapalli, Lyons, and Henry (ePrint 2023/108).
/// @details Each gadget is a template specialization of `exact_lut`. The
///          specialization holds one canonical program (sign pattern, exponent
///          class, or the shared powers-of-ten table) and projects it onto any
///          signed word width and fractional precision. Piece values are exact
///          integers, including at 0.

#ifndef LIBDPF_INCLUDE_GROTTO_CONSTANT_LUT_HPP__
#define LIBDPF_INCLUDE_GROTTO_CONSTANT_LUT_HPP__

#include "hedley/hedley.h"

#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>

namespace grotto
{

/// @brief Appendix D gadgets whose polynomial degree is 0 and whose max error is 0.
enum class exact_constant
{
    signum,
    positive,
    negative,
    nonneg,
    nonpos,
    zero,
    nonzero,
    ilogb,
    /// @brief `ceil(log2(|x|))`. Exact powers of two agree with `ilogb`; every other
    /// positive magnitude is one larger. Zero uses the same `-64` sentinel.
    ceil_ilogb,
    ilog10,
    clz,
    clrsb
};

template <typename Raw>
struct constant_lut
{
    static_assert(std::is_integral_v<Raw> && std::is_signed_v<Raw>);

    using raw_type = Raw;

    /// @brief Signed piece starts. `bounds.front()` is `numeric_limits<Raw>::min()`,
    /// and the starts are strictly increasing.
    std::vector<Raw> bounds;

    /// @brief `values[i]` is the function on `[bounds[i], next)`, where `next` is
    /// `bounds[i + 1]` or one past `numeric_limits<Raw>::max()` for the last piece.
    std::vector<std::int64_t> values;

    HEDLEY_NO_THROW
    std::size_t linear_parts() const noexcept { return values.size(); }

    /// @brief Pieces after joining the first and last when they carry the same value.
    /// @details Those two meet across the signed wrap, which is how the paper counts
    /// parts for `zero` and `nonzero` (2, not 3).
    /// @return Pieces after joining the first and last when they carry the same value
    HEDLEY_NO_THROW
    std::size_t wrapped_parts() const noexcept
    {
        if (values.size() >= 2 && values.front() == values.back())
            return values.size() - 1;
        return values.size();
    }
    /// \complexity `upper_bound` on `bounds`. `Θ(log P)` comparisons, `P = linear_parts()`. Extra space `Θ(1)`.
    /// @param x raw domain point
    /// @return the piece value

    HEDLEY_NO_THROW
    std::int64_t operator()(Raw x) const noexcept
    {
        const auto it = std::upper_bound(bounds.begin(), bounds.end(), x);
        const auto index = static_cast<std::size_t>(it - bounds.begin());
        return values[index - 1];
    }
};

template <exact_constant Which>
struct exact_lut;

namespace detail
{

using u128 = unsigned __int128;

template <typename Raw>
HEDLEY_NO_THROW
constexpr u128 magnitude(Raw raw) noexcept
{
    if (raw >= 0)
        return static_cast<u128>(raw);
    if (raw == std::numeric_limits<Raw>::min())
        return u128{1} << std::numeric_limits<Raw>::digits;
    return static_cast<u128>(-static_cast<__int128>(raw));
}

HEDLEY_NO_THROW
constexpr int floor_log2(u128 mag) noexcept
{
    if (mag <= std::uint64_t(-1))
        return 63 - __builtin_clzll(static_cast<std::uint64_t>(mag));
    return 127 - __builtin_clzll(static_cast<std::uint64_t>(mag >> 64));
}

HEDLEY_NO_THROW
constexpr bool shift_fits(u128 value, unsigned shift) noexcept
{
    return shift < 128 && value <= (~u128{0} >> shift);
}

/// @brief 10^0 .. 10^19. Every ilog10 projection reads this one table.
inline constexpr std::uint64_t pow10[] = {
    1ull,
    10ull,
    100ull,
    1000ull,
    10000ull,
    100000ull,
    1000000ull,
    10000000ull,
    100000000ull,
    1000000000ull,
    10000000000ull,
    100000000000ull,
    1000000000000ull,
    10000000000000ull,
    100000000000000ull,
    1000000000000000ull,
    10000000000000000ull,
    100000000000000000ull,
    1000000000000000000ull,
    10000000000000000000ull,
};

HEDLEY_NO_THROW
inline bool magnitude_ge_pow10(u128 mag, int k, unsigned fractional_bits) noexcept
{
    if (mag == 0 || fractional_bits >= 128)
        return false;
    if (k >= 0)
    {
        if (k >= 20 || !shift_fits(pow10[k], fractional_bits))
            return false;
        return mag >= (u128{pow10[k]} << fractional_bits);
    }
    const int exponent = -k;
    if (exponent >= 20)
        return true;
    const u128 scale = pow10[exponent];
    if (mag > (~u128{0}) / scale)
        return true;
    return mag * scale >= (u128{1} << fractional_bits);
}

/// @brief Smallest positive magnitude whose base-10 log is at least `k`.
/// @param k index or exponent
/// @param fractional_bits the number of fractional bits
/// @return Smallest positive magnitude whose base-10 log is at least `k`
HEDLEY_NO_THROW
inline u128 first_magnitude_at_least_pow10(int k, unsigned fractional_bits) noexcept
{
    if (fractional_bits >= 128)
        return 0;
    if (k >= 0)
    {
        if (k >= 20 || !shift_fits(pow10[k], fractional_bits))
            return 0;
        return u128{pow10[k]} << fractional_bits;
    }
    const int exponent = -k;
    if (exponent >= 20)
        return 1;
    const u128 scale = pow10[exponent];
    const u128 thresh = u128{1} << fractional_bits;
    return (thresh + scale - 1) / scale;
}

template <typename Raw>
void push_negative_magnitude(u128 mag, std::vector<Raw> & cuts)
{
    using lim = std::numeric_limits<Raw>;
    constexpr unsigned digits = lim::digits;
    if (mag == 0)
        return;
    if (mag < (u128{1} << digits))
        cuts.push_back(static_cast<Raw>(-static_cast<__int128>(mag)));
    else
        cuts.push_back(lim::min());
}

template <typename Raw>
void push_both_signs(u128 mag, std::vector<Raw> & cuts)
{
    using lim = std::numeric_limits<Raw>;
    constexpr unsigned digits = lim::digits;
    if (mag == 0)
    {
        cuts.push_back(Raw{0});
        return;
    }
    if (mag <= static_cast<u128>(lim::max()))
        cuts.push_back(static_cast<Raw>(mag));
    if (mag < (u128{1} << digits))
    {
        const __int128 neg = -static_cast<__int128>(mag);
        cuts.push_back(static_cast<Raw>(neg));
        if (neg < static_cast<__int128>(lim::max()))
            cuts.push_back(static_cast<Raw>(neg + 1));
    }
    else
    {
        cuts.push_back(lim::min());
        cuts.push_back(static_cast<Raw>(static_cast<__int128>(lim::min()) + 1));
    }
}

template <typename Raw, typename Fn>
constant_lut<Raw> assemble(std::vector<Raw> cuts, Fn && fn)
{
    using lim = std::numeric_limits<Raw>;
    cuts.push_back(lim::min());
    std::sort(cuts.begin(), cuts.end());
    cuts.erase(std::unique(cuts.begin(), cuts.end()), cuts.end());

    constant_lut<Raw> lut;
    for (std::size_t i = 0; i < cuts.size(); ++i)
    {
        const Raw start = cuts[i];
        const Raw last = (i + 1 < cuts.size())
            ? static_cast<Raw>(cuts[i + 1] - 1)
            : lim::max();
        const std::int64_t at_start = fn(start);
        const std::int64_t at_last = fn(last);
        if (at_start != at_last)
            throw std::logic_error(
                "constant lut endpoints differ start=" + std::to_string(static_cast<long long>(start))
                + " last=" + std::to_string(static_cast<long long>(last))
                + " " + std::to_string(at_start) + " vs " + std::to_string(at_last));
        const __int128 width = static_cast<__int128>(last) - static_cast<__int128>(start);
        if (width > 2)
        {
            const Raw mid = static_cast<Raw>(static_cast<__int128>(start) + width / 2);
            if (fn(mid) != at_start)
                throw std::logic_error(
                    "constant lut midpoint differs start=" + std::to_string(static_cast<long long>(start))
                    + " mid=" + std::to_string(static_cast<long long>(mid)));
        }
        if (!lut.values.empty() && lut.values.back() == at_start)
            continue;
        lut.bounds.push_back(start);
        lut.values.push_back(at_start);
    }
    return lut;
}

template <std::int64_t Neg, std::int64_t Zero, std::int64_t Pos>
struct sign_program
{
    static constexpr std::int64_t canonical[] = { Neg, Zero, Pos };

    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned) noexcept
    {
        if (raw < 0)
            return canonical[0];
        if (raw > 0)
            return canonical[2];
        return canonical[1];
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned)
    {
        using lim = std::numeric_limits<Raw>;
        return assemble<Raw>({ lim::min(), Raw{0}, Raw{1} },
            [](Raw raw) { return eval<Raw>(raw, 0); });
    }
};

}  // namespace detail

template <>
struct exact_lut<exact_constant::signum> : detail::sign_program<-1, 0, 1> {};

template <>
struct exact_lut<exact_constant::positive> : detail::sign_program<0, 0, 1> {};

template <>
struct exact_lut<exact_constant::negative> : detail::sign_program<1, 0, 0> {};

template <>
struct exact_lut<exact_constant::nonneg> : detail::sign_program<0, 1, 1> {};

template <>
struct exact_lut<exact_constant::nonpos> : detail::sign_program<1, 1, 0> {};

template <>
struct exact_lut<exact_constant::zero> : detail::sign_program<0, 1, 0> {};

template <>
struct exact_lut<exact_constant::nonzero> : detail::sign_program<1, 0, 1> {};

template <typename Raw>
void push_pow2_cuts(std::vector<Raw> & cuts)
{
    using lim = std::numeric_limits<Raw>;
    cuts.push_back(lim::min());
    cuts.push_back(Raw{0});
    for (unsigned k = 0; k <= static_cast<unsigned>(lim::digits); ++k)
        detail::push_both_signs<Raw>(detail::u128{1} << k, cuts);
}

/// @brief floor(log2(|raw|)) - F, with -64 on the |x| <= 2^{-64} class (including 0).
/// @details One exponent program; fractional precision only shifts the stored exponent.
template <>
struct exact_lut<exact_constant::ilogb>
{
    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned fractional_bits) noexcept
    {
        const detail::u128 mag = detail::magnitude(raw);
        if (mag == 0)
            return -64;
        if (fractional_bits >= 64)
        {
            const unsigned shift = fractional_bits - 64;
            if (shift >= 128 || mag <= (detail::u128{1} << shift))
                return -64;
        }
        return static_cast<std::int64_t>(detail::floor_log2(mag))
            - static_cast<std::int64_t>(fractional_bits);
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned fractional_bits)
    {
        std::vector<Raw> cuts;
        push_pow2_cuts(cuts);
        return detail::assemble<Raw>(std::move(cuts),
            [=](Raw raw) { return eval<Raw>(raw, fractional_bits); });
    }
};

/// @brief ceil(log2(|x|)). Same powers-of-two cuts as `ilogb`; exact powers keep the
/// floor exponent and every other magnitude steps up by one.
template <>
struct exact_lut<exact_constant::ceil_ilogb>
{
    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned fractional_bits) noexcept
    {
        const detail::u128 mag = detail::magnitude(raw);
        const bool sentinel = mag == 0
            || (fractional_bits >= 64
                && mag <= (detail::u128{1} << (fractional_bits - 64)));
        if (sentinel)
            return -64;
        const std::int64_t floor_exp = static_cast<std::int64_t>(detail::floor_log2(mag))
            - static_cast<std::int64_t>(fractional_bits);
        const bool power = (mag & (mag - 1)) == 0;
        return power ? floor_exp : floor_exp + 1;
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned fractional_bits)
    {
        using lim = std::numeric_limits<Raw>;
        std::vector<Raw> cuts;
        push_pow2_cuts(cuts);
        // Exact powers are singletons; the following raw already has ceil + 1.
        for (unsigned k = 0; k < static_cast<unsigned>(lim::digits); ++k)
        {
            const detail::u128 mag = detail::u128{1} << k;
            if (mag >= static_cast<detail::u128>(lim::max()))
                break;
            cuts.push_back(static_cast<Raw>(mag + 1));
        }
        return detail::assemble<Raw>(std::move(cuts),
            [=](Raw raw) { return eval<Raw>(raw, fractional_bits); });
    }
};

/// @brief floor(log10(|x|)), with -19 on |x| <= 10^{-19}. Thresholds come from `pow10`.
template <>
struct exact_lut<exact_constant::ilog10>
{
    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned fractional_bits) noexcept
    {
        const detail::u128 mag = detail::magnitude(raw);
        if (!detail::magnitude_ge_pow10(mag, -19, fractional_bits))
            return -19;
        int lo = -19;
        int hi = 18;
        while (lo < hi)
        {
            const int mid = lo + (hi - lo + 1) / 2;
            if (detail::magnitude_ge_pow10(mag, mid, fractional_bits))
                lo = mid;
            else
                hi = mid - 1;
        }
        return lo;
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned fractional_bits)
    {
        using lim = std::numeric_limits<Raw>;
        std::vector<Raw> cuts{ lim::min(), Raw{0} };
        for (int k = -19; k <= 18; ++k)
            detail::push_both_signs<Raw>(
                detail::first_magnitude_at_least_pow10(k, fractional_bits), cuts);
        return detail::assemble<Raw>(std::move(cuts),
            [=](Raw raw) { return eval<Raw>(raw, fractional_bits); });
    }
};

/// @brief 64-bit leading-zero count of trunc(x). Negatives are 0; a zero integer part is 64.
/// @details Exponent k of the integer part maps to 63-k after a shift of `fractional_bits`.
template <>
struct exact_lut<exact_constant::clz>
{
    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned fractional_bits) noexcept
    {
        if (raw < 0)
            return 0;
        const detail::u128 mag = detail::magnitude(raw);
        if (fractional_bits >= 128 || mag < (detail::u128{1} << fractional_bits))
            return 64;
        return __builtin_clzll(static_cast<std::uint64_t>(mag >> fractional_bits));
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned fractional_bits)
    {
        using lim = std::numeric_limits<Raw>;
        std::vector<Raw> cuts{ lim::min(), Raw{0} };
        for (unsigned k = 0; k < 63; ++k)
        {
            const detail::u128 exponent = detail::u128{1} << k;
            if (!detail::shift_fits(exponent, fractional_bits))
                break;
            const detail::u128 mag = exponent << fractional_bits;
            if (mag > static_cast<detail::u128>(lim::max()))
                break;
            cuts.push_back(static_cast<Raw>(mag));
        }
        return detail::assemble<Raw>(std::move(cuts),
            [=](Raw raw) { return eval<Raw>(raw, fractional_bits); });
    }
};

/// @brief 64-bit redundant sign bits of trunc(x) toward zero.
/// @details Positive q uses 62-floor(log2(q)); negative q uses 62-floor(log2(q-1)).
template <>
struct exact_lut<exact_constant::clrsb>
{
    template <typename Raw>
    HEDLEY_NO_THROW
    static std::int64_t eval(Raw raw, unsigned fractional_bits) noexcept
    {
        const detail::u128 mag = detail::magnitude(raw);
        const detail::u128 quotient = mag >> fractional_bits;
        if (raw >= 0)
        {
            if (quotient == 0)
                return 63;
            return 62 - detail::floor_log2(quotient);
        }
        if (quotient > static_cast<detail::u128>(std::numeric_limits<std::int64_t>::max()))
            return 0;
        if (quotient <= 1)
            return 63;
        return 62 - detail::floor_log2(quotient - 1);
    }

    template <typename Raw>
    static constant_lut<Raw> project(unsigned fractional_bits)
    {
        using lim = std::numeric_limits<Raw>;
        std::vector<Raw> cuts{ lim::min(), Raw{0} };
        cuts.push_back(static_cast<Raw>(static_cast<__int128>(lim::min()) + 1));
        for (unsigned k = 0; k < 63; ++k)
        {
            const detail::u128 exponent = detail::u128{1} << k;
            if (!detail::shift_fits(exponent, fractional_bits))
                break;
            const detail::u128 positive = exponent << fractional_bits;
            if (positive <= static_cast<detail::u128>(lim::max()))
                cuts.push_back(static_cast<Raw>(positive));

            // Negative class q in [2^k+1, 2^{k+1}], most-negative raw first.
            const detail::u128 q_hi = exponent << 1;
            if (detail::shift_fits(q_hi + 1, fractional_bits))
            {
                const detail::u128 mag_hi = ((q_hi + 1) << fractional_bits) - 1;
                detail::push_negative_magnitude<Raw>(mag_hi, cuts);
            }
            else
            {
                cuts.push_back(lim::min());
            }
        }
        if (detail::shift_fits(detail::u128{2}, fractional_bits))
            detail::push_negative_magnitude<Raw>((detail::u128{2} << fractional_bits) - 1, cuts);
        return detail::assemble<Raw>(std::move(cuts),
            [=](Raw raw) { return eval<Raw>(raw, fractional_bits); });
    }
};
/// \complexity One call to `exact_lut<Which>::eval`. Sign maps are a branch. Bit maps scan at most the raw width. `Θ(w)` bit operations, extra space `Θ(1)`.
/// @see grotto::make_exact_constant_lut

template <exact_constant Which, typename Raw>
std::int64_t evaluate_exact(Raw raw, unsigned fractional_bits)
{
    constexpr unsigned bits = static_cast<unsigned>(std::numeric_limits<Raw>::digits + 1);
    if (fractional_bits > bits)
        throw std::invalid_argument("fractional bits exceed the raw width");
    return exact_lut<Which>::template eval<Raw>(raw, fractional_bits);
}
/// \complexity Dispatches to `exact_lut<Which>::project`. Sign patterns emit a constant number of cuts.
/// `ilogb`, `ceil_ilogb`, `clz`, and `clrsb` push one cut per bit of the raw width (`Θ(w)`). `ilog10` walks powers of ten.
/// Extra space is the cut and value vectors.
/// @see grotto::evaluate_exact
/// @see grotto::dyadic_lut
/// @param which which degree-0 map
/// @param fractional_bits must not exceed the raw width
/// @return the piecewise-constant table
/// @note Following Storrier, Vadapalli, Lyons, and Henry, ePrint 2023/108, Appendix D: the degree-0 exact gadgets.

template <typename Raw>
constant_lut<Raw> make_exact_constant_lut(exact_constant which, unsigned fractional_bits)
{
    constexpr unsigned bits = static_cast<unsigned>(std::numeric_limits<Raw>::digits + 1);
    if (fractional_bits > bits)
        throw std::invalid_argument("fractional bits exceed the raw width");
    switch (which)
    {
    case exact_constant::signum:
        return exact_lut<exact_constant::signum>::template project<Raw>(fractional_bits);
    case exact_constant::positive:
        return exact_lut<exact_constant::positive>::template project<Raw>(fractional_bits);
    case exact_constant::negative:
        return exact_lut<exact_constant::negative>::template project<Raw>(fractional_bits);
    case exact_constant::nonneg:
        return exact_lut<exact_constant::nonneg>::template project<Raw>(fractional_bits);
    case exact_constant::nonpos:
        return exact_lut<exact_constant::nonpos>::template project<Raw>(fractional_bits);
    case exact_constant::zero:
        return exact_lut<exact_constant::zero>::template project<Raw>(fractional_bits);
    case exact_constant::nonzero:
        return exact_lut<exact_constant::nonzero>::template project<Raw>(fractional_bits);
    case exact_constant::ilogb:
        return exact_lut<exact_constant::ilogb>::template project<Raw>(fractional_bits);
    case exact_constant::ceil_ilogb:
        return exact_lut<exact_constant::ceil_ilogb>::template project<Raw>(fractional_bits);
    case exact_constant::ilog10:
        return exact_lut<exact_constant::ilog10>::template project<Raw>(fractional_bits);
    case exact_constant::clz:
        return exact_lut<exact_constant::clz>::template project<Raw>(fractional_bits);
    case exact_constant::clrsb:
        return exact_lut<exact_constant::clrsb>::template project<Raw>(fractional_bits);
    }
    throw std::invalid_argument("unknown exact constant");
}

/// @brief Comparison against a public threshold. Two pieces; the cut sits on `bound`
/// (`lt` / `geq`) or just after it (`leq` / `gt`).
enum class threshold_cmp
{
    lt,
    leq,
    gt,
    geq
};
/// \complexity One comparison. `Θ(1)`.
/// @see grotto::make_threshold_lut

template <typename Raw>
HEDLEY_NO_THROW
std::int64_t evaluate_threshold(Raw raw, Raw bound, threshold_cmp kind) noexcept
{
    switch (kind)
    {
    case threshold_cmp::lt: return raw < bound;
    case threshold_cmp::leq: return raw <= bound;
    case threshold_cmp::gt: return raw > bound;
    case threshold_cmp::geq: return raw >= bound;
    }
    return 0;
}
/// \complexity Two or three cuts (the bound, and `bound + 1` for `leq` / `gt`). `Θ(1)`.
/// @see grotto::evaluate_threshold

template <typename Raw>
constant_lut<Raw> make_threshold_lut(Raw bound, threshold_cmp kind)
{
    using lim = std::numeric_limits<Raw>;
    std::vector<Raw> cuts{ lim::min(), bound };
    if (bound < lim::max()
        && (kind == threshold_cmp::leq || kind == threshold_cmp::gt))
        cuts.push_back(static_cast<Raw>(static_cast<__int128>(bound) + 1));
    return detail::assemble<Raw>(std::move(cuts),
        [=](Raw raw) { return evaluate_threshold(raw, bound, kind); });
}

/// @brief `1` on the inclusive clip window `[low, high]`, `0` outside it.
/// @tparam Raw underlying representation
/// @param low the lower endpoint
/// @param high the upper endpoint
/// @return `1` on the inclusive clip window `[low, high]`, `0` outside it
/// @throws std::invalid_argument if `low > high`
/// \complexity Two or three cuts. `Θ(1)`.
/// @see grotto::make_threshold_lut
template <typename Raw>
constant_lut<Raw> make_interval_lut(Raw low, Raw high)
{
    if (low > high)
        throw std::invalid_argument("interval lut: low > high");
    using lim = std::numeric_limits<Raw>;
    std::vector<Raw> cuts{ lim::min(), low };
    if (high < lim::max())
        cuts.push_back(static_cast<Raw>(static_cast<__int128>(high) + 1));
    return detail::assemble<Raw>(std::move(cuts),
        [=](Raw raw) { return raw >= low && raw <= high ? std::int64_t{1} : std::int64_t{0}; });
}

/// @brief `floor(min(max(raw, low), high) / modulus)`, division toward -infinity.
///
/// `modulus`, `low`, and `high` are in the same raw units as the domain, so
/// one program covers every fractional precision: a mathematical step `M`
/// with `F` fractional bits is the raw modulus `M << F`. The paper's
/// `quot(M, T1, T2)` is this function. Piece count is about `(high-low)/modulus`;
/// the build rejects windows that would need more than 2^16 pieces.
/// @tparam Raw underlying representation
/// @param raw the underlying integer
/// @param modulus the public modulus
/// @param low the lower endpoint
/// @param high the upper endpoint
/// @return `floor(min(max(raw, low), high) / modulus)`, division toward -infinity
/// @throws std::invalid_argument if `modulus must be positive`
/// \complexity One clamp and one division. `Θ(1)`. Division is toward -infinity, as the comment on the declaration says.
/// @see grotto::make_clipped_quotient_lut
template <typename Raw>
std::int64_t evaluate_clipped_quotient(Raw raw, Raw modulus, Raw low, Raw high)
{
    if (modulus <= 0)
        throw std::invalid_argument("clipped quotient: modulus must be positive");
    if (low > high)
        throw std::invalid_argument("clipped quotient: low > high");
    const Raw clipped = raw < low ? low : (raw > high ? high : raw);
    const __int128 n = clipped;
    const __int128 d = modulus;
    if (n >= 0)
        return static_cast<std::int64_t>(n / d);
    const __int128 neg = -n;
    return static_cast<std::int64_t>(-((neg + d - 1) / d));
}
/// \complexity The cut loop steps by `modulus` from `low` through `high`. The header rejects spans with more than `2^16` pieces.
/// Time and extra space `Θ((high - low) / modulus)` in raw units.
/// @see grotto::evaluate_clipped_quotient

template <typename Raw>
constant_lut<Raw> make_clipped_quotient_lut(Raw modulus, Raw low, Raw high)
{
    if (modulus <= 0)
        throw std::invalid_argument("clipped quotient: modulus must be positive");
    if (low > high)
        throw std::invalid_argument("clipped quotient: low > high");
    const __int128 span = static_cast<__int128>(high) - static_cast<__int128>(low);
    if (span / modulus > (1 << 16))
        throw std::invalid_argument("clipped quotient: modulus is too small for the window");

    using lim = std::numeric_limits<Raw>;
    std::vector<Raw> cuts{ lim::min(), low };
    const __int128 d = modulus;
    __int128 q = low >= 0
        ? static_cast<__int128>(low) / d
        : -((-static_cast<__int128>(low) + d - 1) / d);
    // First multiple of `modulus` strictly above `low`, through `high`.
    for (__int128 boundary = (q + 1) * d; boundary <= high; boundary += d)
    {
        if (boundary > static_cast<__int128>(lim::max()))
            break;
        if (boundary >= static_cast<__int128>(lim::min()))
            cuts.push_back(static_cast<Raw>(boundary));
    }
    return detail::assemble<Raw>(std::move(cuts),
        [=](Raw raw) { return evaluate_clipped_quotient(raw, modulus, low, high); });
}

}  // namespace grotto

#endif  // LIBDPF_INCLUDE_GROTTO_CONSTANT_LUT_HPP__