/// @file grotto/range_lut.hpp
/// @brief Full-domain maps built from the principal-domain cubics.
/// @details Each reduced map is one of the elementary range reductions, and
///          the polynomial it evaluates is the matching principal table:
///          `ln` / `lg` / `log10` share the mantissa logarithm;
///          `exp` / `exp2` / `exp10` share the `2^{-13}` exponential;
///          `sin` / `cos` share the quarter-turn sine;
///          `tan` / `cot` share `tanf` and `tang`;
///          `sec` / `csc` share `sec` and `gsec`;
///          `sinh` / `cosh` / `tanh` / `sech` share the hyperbolic addition;
///          `coth` and `csch` use their principal small-argument tables;
///          `sqrt` / `inv` / `rsqrt` / `invsq` are dyadic lifts of `[1/2, 1]`.
///          `expm1` and `log1p` use those same reductions. On `|x| < ln 2`
///          and `|x| <= 1/2` they sum the Taylor series in extra bits so the
///          cancellation in `exp(x)-1` and `ln(1+x)` is not rounded away.
///          Outside that, `expm1` rebuilds `2^n exp(r) - 1` and `log1p`
///          calls the positive logarithm on the exact fixed-point `1+x`.
/// @copyright Copyright (c) 2019-2026 Ryan Henry and [others](@ref authors)
/// @license Released under a GNU General Public v2.0 (GPLv2) license.

#ifndef LIBDPF_INCLUDE_GROTTO_RANGE_LUT_HPP__
#define LIBDPF_INCLUDE_GROTTO_RANGE_LUT_HPP__

#include <cstdint>
#include <stdexcept>

#include "hedley/hedley.h"

#include "grotto/principal_lut.hpp"

namespace grotto
{

enum class reduced : unsigned
{
    ln = 0,
    lg,
    log10,
    exp,
    exp2,
    exp10,
    sin,
    cos,
    tan,
    cot,
    sec,
    csc,
    sinh,
    cosh,
    tanh,
    coth,
    sech,
    csch,
    sqrt,
    inv,
    rsqrt,
    invsq,
    expm1,
    log1p,
};
/// \complexity The `switch` does a constant amount of range reduction and a constant number of `eval_principal` cubics (`Θ(log P)` each).
/// On the small interval, `expm1_series` loops `n = 1 .. 24` and `log1p_series` loops `n = 1 .. 80`, and both stop when the running power is 0.
/// Extra space `Θ(1)`.
/// @see grotto::eval_principal
/// @see grotto::eval_window
/// @see grotto::eval_closed
/// @param which the reduced map
/// @param fractional_bits one of 8, 12, ..., 32
/// @param raw fixed-point argument, value `raw / 2^{fractional_bits}`
/// @return fixed-point result at the same scale

HEDLEY_WARN_UNUSED_RESULT
inline std::int64_t eval_reduced(reduced which, unsigned fractional_bits, std::int64_t raw);

namespace range_detail
{

using u128 = unsigned __int128;

HEDLEY_CONST
HEDLEY_NO_THROW
constexpr u128 words(std::uint64_t hi, std::uint64_t lo) noexcept
{
    return (u128{hi} << 64) | lo;
}

/// @brief `value * 2^64`, rounded half away from zero. Values above `2^64` keep the
/// high limb so the constant is not truncated.
inline constexpr u128 ln2_64 = words(0, 12786308645202655660ULL);
inline constexpr u128 inv_ln2_64 = words(1, 8166282121979093367ULL);
inline constexpr u128 log10_2_64 = words(0, 5553023288523357132ULL);
inline constexpr u128 ln10_64 = words(2, 5581709770980765788ULL);
inline constexpr u128 inv_ln10_64 = words(0, 8011319160293570763ULL);
inline constexpr u128 sqrt2_64 = words(1, 7640891576956012809ULL);
inline constexpr u128 rsqrt2_64 = words(0, 13043817825332782212ULL);
inline constexpr u128 two_over_pi_64 = words(0, 11743562013128004906ULL);
inline constexpr u128 four_over_pi_64 = words(1, 5040379952546458196ULL);
inline constexpr u128 pi_over_4_64 = words(0, 14488038916154245685ULL);

/// @brief `exp(2^{i-13}) * 2^64`.
inline constexpr u128 exp_chunk_64[13] = {
    words(1, 2251937258231296ULL),
    words(1, 4504149427926357ULL),
    words(1, 9009398635954180ULL),
    words(1, 18023197466514910ULL),
    words(1, 36064004308734226ULL),
    words(1, 72198514957318099ULL),
    words(1, 144679606912572172ULL),
    words(1, 290493950045950331ULL),
    words(1, 585562514163419534ULL),
    words(1, 1189712777830127574ULL),
    words(1, 2456155437534072733ULL),
    words(1, 5239344172067481206ULL),
    words(1, 11966795255776918679ULL),
};

inline std::int64_t round_mag(u128 mag, unsigned shift, bool neg)
{
    if (shift >= 128)
        return 0;
    if (shift > 0)
    {
        mag += u128{1} << (shift - 1);
        mag >>= shift;
    }
    if (mag > static_cast<u128>(INT64_MAX))
        throw std::overflow_error("range lut: value does not fit int64");
    const auto out = static_cast<std::int64_t>(mag);
    return neg ? -out : out;
}

inline std::int64_t round_i128(__int128 value, unsigned shift)
{
    const bool neg = value < 0;
    const auto mag = static_cast<u128>(neg ? -value : value);
    return round_mag(mag, shift, neg);
}

inline std::int64_t scale_unit(u128 mag64, unsigned fractional_bits)
{
    return round_mag(mag64, 64u - fractional_bits, false);
}

inline std::int64_t mul_raw(std::int64_t lhs, std::int64_t rhs, unsigned fractional_bits)
{
    return round_i128(static_cast<__int128>(lhs) * rhs, fractional_bits);
}

inline std::int64_t div_raw(std::int64_t num, std::int64_t den, unsigned fractional_bits)
{
    if (den == 0)
        throw std::domain_error("range lut: division by zero");
    const bool neg = (num < 0) != (den < 0);
    auto n = static_cast<u128>(num < 0 ? -static_cast<__int128>(num) : num);
    auto d = static_cast<u128>(den < 0 ? -static_cast<__int128>(den) : den);
    n <<= fractional_bits;
    const u128 quot = (n + d / 2) / d;
    return round_mag(quot, 0, neg);
}

inline std::int64_t shift_pow2(std::int64_t value, int places)
{
    if (places == 0 || value == 0)
        return value;
    if (places > 0)
    {
        if (places >= 62)
            throw std::overflow_error("range lut: exponent overflow");
        const __int128 wide = static_cast<__int128>(value) << places;
        if (wide > INT64_MAX || wide < INT64_MIN)
            throw std::overflow_error("range lut: exponent overflow");
        return static_cast<std::int64_t>(wide);
    }
    return round_i128(value, static_cast<unsigned>(-places));
}

HEDLEY_CONST
HEDLEY_NO_THROW
constexpr std::int64_t one_raw(unsigned fractional_bits) noexcept
{
    return std::int64_t{1} << fractional_bits;
}

HEDLEY_CONST
HEDLEY_NO_THROW
constexpr u128 magnitude_of(std::int64_t raw) noexcept
{
    if (raw >= 0)
        return static_cast<u128>(raw);
    return static_cast<u128>(-static_cast<__int128>(raw));
}

inline std::int64_t abs_raw(std::int64_t raw)
{
    const u128 mag = magnitude_of(raw);
    if (mag > static_cast<u128>(INT64_MAX))
        throw std::overflow_error("range lut: magnitude does not fit int64");
    return static_cast<std::int64_t>(mag);
}

struct dyadic
{
    std::int64_t mantissa_raw;
    int power;
};

inline dyadic split_positive(std::int64_t raw, unsigned fractional_bits)
{
    if (raw <= 0)
        throw std::domain_error("range lut: reduction requires a positive input");
    const auto mag = static_cast<unsigned long long>(raw);
    const int floor_log = 63 - __builtin_clzll(mag);
    const int shift = static_cast<int>(fractional_bits) - floor_log - 1;
    std::int64_t mantissa = shift >= 0
        ? raw << shift
        : round_i128(raw, static_cast<unsigned>(-shift));
    int power = floor_log + 1 - static_cast<int>(fractional_bits);
    const std::int64_t one = one_raw(fractional_bits);
    const std::int64_t half = one >> 1;
    if (mantissa >= one)
    {
        mantissa >>= 1;
        ++power;
    }
    if (mantissa < half)
        mantissa = half;
    return dyadic{mantissa, power};
}

inline std::int64_t ln2_raw(unsigned fractional_bits)
{
    return scale_unit(ln2_64, fractional_bits);
}

inline std::int64_t eval_ln_positive(unsigned fractional_bits, std::int64_t raw);
inline std::int64_t eval_log10_positive(unsigned fractional_bits, std::int64_t raw);
inline u128 exp_scale64(unsigned fractional_bits, std::int64_t raw, std::int64_t & n_bin);
inline std::int64_t finish_wide(u128 wide, int right_shift);

inline std::int64_t eval_exp_at_scale(unsigned fractional_bits, std::int64_t raw)
{
    if (fractional_bits < 13)
    {
        const int lift = static_cast<int>(16u - fractional_bits);
        const __int128 lifted_arg = static_cast<__int128>(raw) << lift;
        if (lifted_arg > INT64_MAX || lifted_arg < INT64_MIN)
            throw std::overflow_error("range lut: exponent overflow");
        const std::int64_t lifted = eval_exp_at_scale(16, static_cast<std::int64_t>(lifted_arg));
        return round_i128(lifted, static_cast<unsigned>(lift));
    }
    std::int64_t n_bin = 0;
    const u128 wide = exp_scale64(fractional_bits, raw, n_bin);
    const int shift = static_cast<int>(64u - fractional_bits) - static_cast<int>(n_bin);
    return finish_wide(wide, shift);
}

HEDLEY_NO_THROW
constexpr std::int64_t fractional_raw(std::int64_t raw, unsigned fractional_bits,
    std::int64_t & whole) noexcept
{
    const std::int64_t one = one_raw(fractional_bits);
    std::int64_t q = raw / one;
    std::int64_t f = raw - q * one;
    if (f < 0)
    {
        f += one;
        --q;
    }
    whole = q;
    return f;
}

inline std::int64_t pow10_raw(int exponent, unsigned fractional_bits)
{
    const std::int64_t one = one_raw(fractional_bits);
    if (exponent == 0)
        return one;
    if (exponent < 0)
        return div_raw(one, pow10_raw(-exponent, fractional_bits), fractional_bits);
    u128 acc = static_cast<u128>(one);
    for (int i = 0; i < exponent; ++i)
    {
        if (acc > static_cast<u128>(INT64_MAX) / 10)
            throw std::overflow_error("range lut: exponent overflow");
        acc *= 10;
    }
    return static_cast<std::int64_t>(acc);
}

struct angle
{
    unsigned index;
    std::int64_t frac_raw;
};

/// @brief `{ |x| * multiplier }` at this precision, with the integer part reduced
/// only as far as the low bits the quadrant logic reads.
/// @param fractional_bits the number of fractional bits
/// @param raw the underlying integer
/// @param multiplier_64 multiplier already scaled by `2^64`
/// @return `{ |x| * multiplier }` at this precision, with the integer part reduced only as far as
/// the low bits the quadrant logic reads
inline angle reduce_positive(unsigned fractional_bits, std::int64_t raw, u128 multiplier_64)
{
    const u128 scaled = magnitude_of(raw) * multiplier_64;
    const u128 rounded = (scaled + (u128{1} << 63)) >> 64;
    const u128 one = u128{1} << fractional_bits;
    return angle{
        static_cast<unsigned>(rounded >> fractional_bits),
        static_cast<std::int64_t>(rounded & (one - 1)),
    };
}

inline std::int64_t principal_sin_fraction(unsigned fractional_bits, std::int64_t fraction_raw, bool complement)
{
    const std::int64_t one = one_raw(fractional_bits);
    std::int64_t argument = complement ? one - fraction_raw : fraction_raw;
    if (argument < 0)
        argument = 0;
    if (argument > one)
        argument = one;
    return eval_principal(principal::sin, fractional_bits, argument);
}

inline std::int64_t sin_from_angle(unsigned fractional_bits, const angle & turned, int sign)
{
    const unsigned which = turned.index & 3u;
    const bool complement = which == 1 || which == 3;
    const int quadrant_sign = (which == 2 || which == 3) ? -1 : 1;
    const std::int64_t magnitude = principal_sin_fraction(
        fractional_bits, turned.frac_raw, complement);
    return magnitude * quadrant_sign * sign;
}

inline std::int64_t cos_from_angle(unsigned fractional_bits, const angle & turned)
{
    angle shifted = turned;
    shifted.index += 1;
    return sin_from_angle(fractional_bits, shifted, 1);
}

inline std::int64_t pi_over_4_raw(unsigned fractional_bits)
{
    return scale_unit(pi_over_4_64, fractional_bits);
}

inline std::int64_t tan_positive(unsigned fractional_bits, std::int64_t magnitude, int quarter_shift)
{
    const angle turned = reduce_positive(fractional_bits, magnitude, four_over_pi_64);
    const unsigned q = (turned.index + static_cast<unsigned>(quarter_shift)) & 3u;
    const std::int64_t one = one_raw(fractional_bits);
    std::int64_t t = (q == 0 || q == 2) ? turned.frac_raw : one - turned.frac_raw;
    if (t < 0)
        t = 0;
    if (t > one)
        t = one;
    const std::int64_t z = mul_raw(t, pi_over_4_raw(fractional_bits), fractional_bits);
    if (q == 0 || q == 3)
    {
        const std::int64_t tanf = eval_principal(principal::tanf, fractional_bits, t);
        const std::int64_t y = mul_raw(z, tanf, fractional_bits);
        return q == 3 ? -y : y;
    }
    if (z == 0)
        throw std::domain_error("range lut: tan pole");
    const std::int64_t tang = eval_principal(principal::tang, fractional_bits, t);
    const std::int64_t y = div_raw(one, z, fractional_bits) + tang;
    return q == 2 ? -y : y;
}

inline std::int64_t sec_positive(unsigned fractional_bits, std::int64_t magnitude, int octant_shift)
{
    const angle turned = reduce_positive(fractional_bits, magnitude, four_over_pi_64);
    const unsigned q8 = (turned.index + static_cast<unsigned>(octant_shift)) & 7u;
    const unsigned q = q8 & 3u;
    const int sigma = (q8 & 4u) == 0 ? 1 : -1;
    const std::int64_t one = one_raw(fractional_bits);
    std::int64_t t = (q == 0 || q == 2) ? turned.frac_raw : one - turned.frac_raw;
    if (t < 0)
        t = 0;
    if (t > one)
        t = one;
    if (q == 0 || q == 3)
    {
        const std::int64_t sec = eval_principal(principal::sec, fractional_bits, t);
        const int sign = (q == 3 ? -1 : 1) * sigma;
        return sec * sign;
    }
    const std::int64_t z = mul_raw(t, pi_over_4_raw(fractional_bits), fractional_bits);
    if (z == 0)
        throw std::domain_error("range lut: sec pole");
    const std::int64_t gsec = eval_principal(principal::gsec, fractional_bits, t);
    std::int64_t y = div_raw(one, z, fractional_bits) + gsec;
    if (q == 2)
        y = -y;
    return y * sigma;
}

inline void quotient_2_13(unsigned fractional_bits, std::int64_t magnitude,
    std::int64_t & quotient, std::int64_t & remainder)
{
    if (fractional_bits >= 13)
    {
        const unsigned shift = fractional_bits - 13;
        quotient = magnitude >> shift;
        const std::int64_t mask = shift >= 63 ? INT64_MAX : (std::int64_t{1} << shift) - 1;
        remainder = shift == 0 ? 0 : magnitude & mask;
        return;
    }
    const int lift = static_cast<int>(13u - fractional_bits);
    const __int128 wide = static_cast<__int128>(magnitude) << lift;
    if (wide > INT64_MAX)
        throw std::overflow_error("range lut: exponent overflow");
    quotient = static_cast<std::int64_t>(wide);
    remainder = 0;
}

inline std::int64_t exp_of_quotient(unsigned fractional_bits, std::int64_t quotient, std::int64_t magnitude)
{
    if (quotient == 0)
        return one_raw(fractional_bits);
    __int128 argument;
    if (fractional_bits >= 13)
        argument = static_cast<__int128>(quotient) << (fractional_bits - 13);
    else
        argument = magnitude;
    if (argument > INT64_MAX)
        throw std::overflow_error("range lut: exponent overflow");
    return eval_exp_at_scale(fractional_bits, static_cast<std::int64_t>(argument));
}

struct hyp
{
    std::int64_t sinh_raw;
    std::int64_t cosh_raw;
};

inline hyp sinh_cosh(unsigned fractional_bits, std::int64_t raw)
{
    const bool neg = raw < 0;
    const auto mag_wide = magnitude_of(raw);
    if (mag_wide > static_cast<u128>(INT64_MAX))
        throw std::overflow_error("range lut: exponent overflow");
    const std::int64_t mag = static_cast<std::int64_t>(mag_wide);
    std::int64_t quotient = 0;
    std::int64_t remainder = 0;
    quotient_2_13(fractional_bits, mag, quotient, remainder);

    const std::int64_t one = one_raw(fractional_bits);
    std::int64_t table = 0;
    if (fractional_bits >= 13 && remainder != 0)
    {
        const __int128 lifted = static_cast<__int128>(remainder) << 13;
        table = lifted > one ? one : static_cast<std::int64_t>(lifted);
    }
    const std::int64_t sr = eval_principal(principal::sinh, fractional_bits, table);
    const std::int64_t cr = eval_principal(principal::cosh, fractional_bits, table);
    std::int64_t sh = sr;
    std::int64_t ch = cr;
    if (quotient != 0)
    {
        const std::int64_t grown = exp_of_quotient(fractional_bits, quotient, mag);
        std::int64_t inv = 0;
        if (grown != 0)
            inv = div_raw(one, grown, fractional_bits);
        const std::int64_t sq = round_i128(static_cast<__int128>(grown) - inv, 1);
        const std::int64_t cq = round_i128(static_cast<__int128>(grown) + inv, 1);
        const __int128 sinh_sum = static_cast<__int128>(sq) * cr + static_cast<__int128>(cq) * sr;
        const __int128 cosh_sum = static_cast<__int128>(cq) * cr + static_cast<__int128>(sq) * sr;
        sh = round_i128(sinh_sum, fractional_bits);
        ch = round_i128(cosh_sum, fractional_bits);
    }
    if (neg)
        sh = -sh;
    return hyp{sh, ch};
}

/// @brief `ln(2^{k+1} ± 1) / 2`, the saturation threshold used by `tanh` and `coth`.
/// @param fractional_bits the number of fractional bits
/// @param plus true for the plus saturation threshold, false for the minus threshold
/// @return `ln(2^{k+1} ± 1) / 2`, the saturation threshold used by `tanh` and `coth`
inline std::int64_t beta_raw(unsigned fractional_bits, bool plus)
{
    u128 ln = u128{fractional_bits + 1} * ln2_64;
    const u128 eps = u128{1} << (63u - fractional_bits);
    if (plus)
        ln += eps;
    else
        ln -= eps;
    return round_mag(ln, 65u - fractional_bits, false);
}

HEDLEY_CONST
HEDLEY_NO_THROW
constexpr int half_pow_of(int power) noexcept
{
    return (power & 1) != 0 ? (power - 1) / 2 : power / 2;
}

struct u256
{
    u128 lo;
    u128 hi;
};

inline u256 mul_u128(u128 a, u128 b)
{
    const auto a0 = static_cast<std::uint64_t>(a);
    const auto a1 = static_cast<std::uint64_t>(a >> 64);
    const auto b0 = static_cast<std::uint64_t>(b);
    const auto b1 = static_cast<std::uint64_t>(b >> 64);
    const u128 p00 = u128{a0} * b0;
    const u128 p01 = u128{a0} * b1;
    const u128 p10 = u128{a1} * b0;
    const u128 p11 = u128{a1} * b1;
    const u128 col = (p00 >> 64) + static_cast<std::uint64_t>(p01) + static_cast<std::uint64_t>(p10);
    u256 out;
    out.lo = static_cast<std::uint64_t>(p00) | (col << 64);
    out.hi = p11 + (p01 >> 64) + (p10 >> 64) + (col >> 64);
    return out;
}

inline u128 round_u256(u256 value, unsigned shift)
{
    if (shift == 0)
        return value.lo;
    if (shift >= 256)
        return 0;
    u256 bumped = value;
    const unsigned bit = shift - 1;
    if (bit < 128)
    {
        const u128 before = bumped.lo;
        bumped.lo += u128{1} << bit;
        if (bumped.lo < before)
            ++bumped.hi;
    }
    else
        bumped.hi += u128{1} << (bit - 128);
    if (shift < 128)
    {
        if (shift == 0)
            return bumped.lo;
        return (bumped.lo >> shift) | (bumped.hi << (128 - shift));
    }
    return bumped.hi >> (shift - 128);
}

/// @brief `ln(m) * 2^64` for `m` in `[1/2, 1]`, via `2 artanh((m-1)/(m+1))`.
/// @details `|z| <= 1/3`, so forty odd powers sit well below `2^{-64}`.
inline u128 ln_mantissa_scale64(unsigned fractional_bits, std::int64_t mantissa_raw)
{
    const u128 one = u128{1} << 64;
    const u128 m = static_cast<u128>(mantissa_raw) << (64u - fractional_bits);
    if (m >= one)
        return 0;
    const u128 num = one - m;
    const u128 den = one + m;
    const u128 z = ((num << 64) + den / 2) / den;
    const u128 z2 = round_u256(mul_u128(z, z), 64);
    u128 acc = z;
    u128 power = z;
    for (int n = 1; n <= 40; ++n)
    {
        power = round_u256(mul_u128(power, z2), 64);
        const unsigned denom = static_cast<unsigned>(2 * n + 1);
        const u128 term = (power + denom / 2) / denom;
        if (term == 0)
            break;
        acc += term;
    }
    return acc << 1;
}

inline std::int64_t round_scale64_to_k(u128 mag, bool neg, unsigned fractional_bits)
{
    return round_mag(mag, 64u - fractional_bits, neg);
}

inline void ln_magnitude_scale64(unsigned fractional_bits, std::int64_t raw, u128 & mag, bool & neg)
{
    const dyadic part = split_positive(raw, fractional_bits);
    // `ln(m) <= 0` on `[1/2, 1]`, so `ln(m * 2^e) = e·ln 2 - |ln m|`.
    const u128 ln_m = ln_mantissa_scale64(fractional_bits, part.mantissa_raw);
    if (part.power >= 0)
    {
        const u128 lift = ln2_64 * static_cast<u128>(part.power);
        if (lift >= ln_m)
        {
            mag = lift - ln_m;
            neg = false;
        }
        else
        {
            mag = ln_m - lift;
            neg = true;
        }
    }
    else
    {
        mag = ln2_64 * static_cast<u128>(-part.power) + ln_m;
        neg = true;
    }
}

inline std::int64_t eval_ln_positive(unsigned fractional_bits, std::int64_t raw)
{
    u128 mag = 0;
    bool neg = false;
    ln_magnitude_scale64(fractional_bits, raw, mag, neg);
    return round_scale64_to_k(mag, neg, fractional_bits);
}

/// @brief `(rem << 64) / den`, rounded. `rem < den` and `den < 2^96`.
inline u128 div_rem_lshift64(u128 rem, u128 den)
{
    const u128 hi = (rem << 32) / den;
    const u128 mid = (rem << 32) % den;
    const u128 lo = (mid << 32) / den;
    const u128 leftover = (mid << 32) % den;
    u128 out = (hi << 32) + lo;
    if (leftover >= den / 2)
        ++out;
    return out;
}

inline std::int64_t eval_log10_positive(unsigned fractional_bits, std::int64_t raw)
{
    u128 ln_mag = 0;
    bool neg = false;
    ln_magnitude_scale64(fractional_bits, raw, ln_mag, neg);
    const u128 quot = ln_mag / ln10_64;
    const u128 rem = ln_mag % ln10_64;
    const u128 log_mag = (quot << 64) + div_rem_lshift64(rem, ln10_64);
    return round_scale64_to_k(log_mag, neg, fractional_bits);
}

/// @brief `exp(x) * 2^64`. The `ln 2` split and the `2^{-13}` chunks stay at
///        scale 64 and are rounded once into the caller's precision.
inline u128 exp_scale64(unsigned fractional_bits, std::int64_t raw, std::int64_t & n_bin)
{
    const u128 x64 = static_cast<u128>(raw < 0 ? -static_cast<__int128>(raw) : raw)
        << (64u - fractional_bits);
    const bool neg = raw < 0;
    u128 mag = x64;
    n_bin = 0;
    if (mag >= ln2_64)
    {
        n_bin = static_cast<std::int64_t>(mag / ln2_64);
        mag -= ln2_64 * static_cast<u128>(n_bin);
    }
    if (neg)
    {
        if (mag == 0)
            n_bin = -n_bin;
        else
        {
            n_bin = -n_bin - 1;
            mag = ln2_64 - mag;
        }
    }
    const u128 step = u128{1} << 51;
    const u128 chunks = mag / step;
    u128 tiny = mag % step;
    u128 acc = u128{1} << 64;
    u128 power = tiny;
    for (int n = 1; n <= 16; ++n)
    {
        const u128 term = (power + static_cast<u128>(n) / 2) / static_cast<u128>(n);
        if (term == 0)
            break;
        acc += term;
        power = round_u256(mul_u128(term, tiny), 64);
    }
    for (unsigned bit = 0; bit < 13; ++bit)
    {
        if (((chunks >> bit) & 1u) == 0)
            continue;
        acc = round_u256(mul_u128(acc, exp_chunk_64[bit]), 64);
    }
    return acc;
}

/// @brief Two Newton steps at scale `2k`, then the exact power-of-two lift.
/// @details The principal cubic is half an ulp at scale `k`. Shifting that
///          rounded word left multiplies the error. Refining before the shift
///          leaves an absolute error below one output ulp across the domain.
inline u128 newton_inv(unsigned fractional_bits, std::int64_t mantissa_raw, std::int64_t seed_raw)
{
    const unsigned K = fractional_bits * 2u;
    u128 m = static_cast<u128>(mantissa_raw) << fractional_bits;
    u128 y = static_cast<u128>(seed_raw) << fractional_bits;
    const u128 two = u128{2} << K;
    for (int step = 0; step < 2; ++step)
    {
        const u128 my = round_u256(mul_u128(m, y), K);
        if (my >= two)
            break;
        y = round_u256(mul_u128(y, two - my), K);
    }
    return y;
}

inline u128 newton_rsqrt(unsigned fractional_bits, std::int64_t mantissa_raw, std::int64_t seed_raw)
{
    const unsigned K = fractional_bits * 2u;
    u128 m = static_cast<u128>(mantissa_raw) << fractional_bits;
    u128 y = static_cast<u128>(seed_raw) << fractional_bits;
    const u128 three = u128{3} << K;
    for (int step = 0; step < 2; ++step)
    {
        const u128 yy = round_u256(mul_u128(y, y), K);
        const u128 myy = round_u256(mul_u128(m, yy), K);
        if (myy >= three)
            break;
        const u128 corr = round_u256(mul_u128(y, three - myy), K);
        y = (corr + 1) >> 1;
    }
    return y;
}

inline std::int64_t finish_wide(u128 wide, int right_shift)
{
    if (right_shift >= 256)
        return 0;
    if (right_shift >= 0)
    {
        u256 value{wide, 0};
        const u128 rounded = round_u256(value, static_cast<unsigned>(right_shift));
        if (rounded > static_cast<u128>(INT64_MAX))
            throw std::overflow_error("range lut: reciprocal does not fit int64");
        return static_cast<std::int64_t>(rounded);
    }
    const int left = -right_shift;
    if (left >= 127)
        throw std::overflow_error("range lut: reciprocal does not fit int64");
    const u128 shifted = wide << static_cast<unsigned>(left);
    if (shifted > static_cast<u128>(INT64_MAX))
        throw std::overflow_error("range lut: reciprocal does not fit int64");
    return static_cast<std::int64_t>(shifted);
}

inline __int128 div_round_i128(__int128 num, int den)
{
    const bool neg = num < 0;
    const auto mag = static_cast<u128>(neg ? -num : num);
    const auto d = static_cast<u128>(den);
    const u128 quot = (mag + d / 2) / d;
    return neg ? -static_cast<__int128>(quot) : static_cast<__int128>(quot);
}

inline __int128 shr_round_i128(__int128 num, unsigned shift)
{
    if (shift == 0)
        return num;
    const bool neg = num < 0;
    auto mag = static_cast<u128>(neg ? -num : num);
    mag = (mag + (u128{1} << (shift - 1))) >> shift;
    return neg ? -static_cast<__int128>(mag) : static_cast<__int128>(mag);
}

/// @brief `expm1` on `|x| < ln 2`, summed at `k+48` fractional bits.
/// @param fractional_bits the number of fractional bits
/// @param raw the underlying integer
/// @return `expm1` on `|x| < ln 2`, summed at `k+48` fractional bits
inline std::int64_t expm1_series(unsigned fractional_bits, std::int64_t raw)
{
    constexpr unsigned extra = 48;
    __int128 power = static_cast<__int128>(raw) << extra;
    __int128 acc = 0;
    for (int n = 1; n <= 24; ++n)
    {
        const __int128 term = div_round_i128(power, n);
        acc += term;
        power = shr_round_i128(term * static_cast<__int128>(raw), fractional_bits);
        if (power == 0)
            break;
    }
    return round_i128(acc, extra);
}

inline std::int64_t eval_expm1(unsigned fractional_bits, std::int64_t raw)
{
    if (raw == 0)
        return 0;
    // Match `exp`: precisions below the 2^{-13} reduction evaluate one
    // scale up and round once, so the power-of-two lift is not rounded early.
    if (fractional_bits < 13)
    {
        const int lift = static_cast<int>(16u - fractional_bits);
        const __int128 lifted_arg = static_cast<__int128>(raw) << lift;
        if (lifted_arg > INT64_MAX || lifted_arg < INT64_MIN)
            throw std::overflow_error("range lut: exponent overflow");
        const std::int64_t lifted = eval_expm1(16, static_cast<std::int64_t>(lifted_arg));
        return round_i128(lifted, static_cast<unsigned>(lift));
    }
    const std::int64_t ln2 = ln2_raw(fractional_bits);
    if (raw > -ln2 && raw < ln2)
        return expm1_series(fractional_bits, raw);

    std::int64_t n_bin = 0;
    const u128 wide = exp_scale64(fractional_bits, raw, n_bin);
    if (n_bin >= 63)
        throw std::overflow_error("range lut: exponent overflow");
    u128 exp64 = wide;
    if (n_bin >= 0)
        exp64 <<= static_cast<unsigned>(n_bin);
    else if (-n_bin >= 128)
        exp64 = 0;
    else
        exp64 >>= static_cast<unsigned>(-n_bin);
    const u128 unit = u128{1} << 64;
    const bool below = exp64 < unit;
    const u128 diff = below ? unit - exp64 : exp64 - unit;
    return round_mag(diff, 64u - fractional_bits, below);
}

/// @brief `log1p` on `|x| <= 1/2`. Every term of a negative argument is negative.
/// @param fractional_bits the number of fractional bits
/// @param raw the underlying integer
/// @return `log1p` on `|x| <= 1/2`
inline std::int64_t log1p_series(unsigned fractional_bits, std::int64_t raw)
{
    constexpr unsigned extra = 48;
    const bool xneg = raw < 0;
    const std::int64_t mag_raw = xneg ? -raw : raw;
    __int128 power = static_cast<__int128>(mag_raw) << extra;
    __int128 acc = 0;
    for (int n = 1; n <= 80; ++n)
    {
        const __int128 term = div_round_i128(power, n);
        const bool neg = xneg || (n % 2 == 0);
        acc += neg ? -term : term;
        power = shr_round_i128(power * static_cast<__int128>(mag_raw), fractional_bits);
        if (power == 0)
            break;
    }
    return round_i128(acc, extra);
}

inline std::int64_t eval_log1p(unsigned fractional_bits, std::int64_t raw)
{
    const std::int64_t one = one_raw(fractional_bits);
    if (raw == 0)
        return 0;
    if (raw <= -one)
        throw std::domain_error("range lut: log1p argument is <= -1");
    const std::int64_t half = one >> 1;
    if (raw >= -half && raw <= half)
        return log1p_series(fractional_bits, raw);
    if (raw > INT64_MAX - one)
    {
        const std::int64_t ln_x = eval_ln_positive(fractional_bits, raw);
        const u128 num = u128{1} << (2u * fractional_bits);
        const auto corr = static_cast<std::int64_t>((num + static_cast<u128>(raw) / 2) / static_cast<u128>(raw));
        return ln_x + corr;
    }
    return eval_ln_positive(fractional_bits, one + raw);
}

}  // namespace range_detail
/// \complexity The `switch` does a constant amount of range reduction and a constant number of `eval_principal` cubics (`Θ(log P)` each).
/// On the small interval, `expm1_series` loops `n = 1 .. 24` and `log1p_series` loops `n = 1 .. 80`, and both stop when the running power is 0.
/// Extra space `Θ(1)`.
/// @see grotto::eval_principal
/// @see grotto::eval_window
/// @see grotto::eval_closed
/// @param which the reduced map
/// @param fractional_bits one of 8, 12, ..., 32
/// @param raw fixed-point argument, value `raw / 2^{fractional_bits}`
/// @return fixed-point result at the same scale

HEDLEY_WARN_UNUSED_RESULT
inline std::int64_t eval_reduced(reduced which, unsigned fractional_bits, std::int64_t raw)
{
    using namespace range_detail;
    if (!principal_precision(fractional_bits))
        throw std::invalid_argument("range lut: precision must be 8, 12, ..., 32");
    const std::int64_t one = one_raw(fractional_bits);
    switch (which)
    {
    case reduced::ln:
        return eval_ln_positive(fractional_bits, raw);
    case reduced::lg:
    {
        const dyadic part = split_positive(raw, fractional_bits);
        const std::int64_t ln_m = eval_principal(principal::ln, fractional_bits, part.mantissa_raw);
        const std::int64_t lg_m = mul_raw(
            ln_m, scale_unit(inv_ln2_64, fractional_bits), fractional_bits);
        return lg_m + (static_cast<std::int64_t>(part.power) << fractional_bits);
    }
    case reduced::log10:
        return eval_log10_positive(fractional_bits, raw);
    case reduced::exp:
        return eval_exp_at_scale(fractional_bits, raw);
    case reduced::exp2:
    {
        std::int64_t whole = 0;
        const std::int64_t frac = fractional_raw(raw, fractional_bits, whole);
        const std::int64_t natural = mul_raw(frac, ln2_raw(fractional_bits), fractional_bits);
        return shift_pow2(eval_exp_at_scale(fractional_bits, natural), static_cast<int>(whole));
    }
    case reduced::exp10:
    {
        std::int64_t whole = 0;
        const std::int64_t frac = fractional_raw(raw, fractional_bits, whole);
        const std::int64_t natural = mul_raw(
            frac, scale_unit(ln10_64, fractional_bits), fractional_bits);
        if (whole > 18 || whole < -18)
            throw std::overflow_error("range lut: exponent overflow");
        return mul_raw(
            eval_exp_at_scale(fractional_bits, natural),
            pow10_raw(static_cast<int>(whole), fractional_bits),
            fractional_bits);
    }
    case reduced::sin:
        return sin_from_angle(
            fractional_bits,
            reduce_positive(fractional_bits, raw, two_over_pi_64),
            raw < 0 ? -1 : 1);
    case reduced::cos:
        return cos_from_angle(
            fractional_bits, reduce_positive(fractional_bits, raw, two_over_pi_64));
    case reduced::tan:
    {
        const std::int64_t y = tan_positive(fractional_bits, abs_raw(raw), 0);
        return raw < 0 ? -y : y;
    }
    case reduced::cot:
    {
        if (raw == 0)
            throw std::domain_error("range lut: cot pole");
        const std::int64_t y = -tan_positive(fractional_bits, abs_raw(raw), 2);
        return raw < 0 ? -y : y;
    }
    case reduced::sec:
        return sec_positive(fractional_bits, abs_raw(raw), 0);
    case reduced::csc:
    {
        if (raw == 0)
            throw std::domain_error("range lut: csc pole");
        const std::int64_t y = sec_positive(fractional_bits, abs_raw(raw), -2);
        return raw < 0 ? -y : y;
    }
    case reduced::sinh:
        return sinh_cosh(fractional_bits, raw).sinh_raw;
    case reduced::cosh:
        return sinh_cosh(fractional_bits, raw).cosh_raw;
    case reduced::tanh:
    {
        if (raw == 0)
            return 0;
        const std::int64_t limit = beta_raw(fractional_bits, false);
        const std::int64_t mag = abs_raw(raw);
        if (mag >= limit)
            return raw < 0 ? -one : one;
        const hyp pair = sinh_cosh(fractional_bits, raw);
        return div_raw(pair.sinh_raw, pair.cosh_raw, fractional_bits);
    }
    case reduced::coth:
    {
        if (raw == 0)
            throw std::domain_error("range lut: coth pole");
        const std::int64_t limit = beta_raw(fractional_bits, true);
        const std::int64_t mag = abs_raw(raw);
        std::int64_t y;
        if (mag >= limit)
            y = one;
        else
        {
            const std::int64_t t = div_raw(mag, limit, fractional_bits);
            const std::int64_t argument = t > one ? one : t;
            const std::int64_t removed = eval_principal(principal::coth, fractional_bits, argument);
            y = removed + div_raw(one, mag, fractional_bits);
        }
        return raw < 0 ? -y : y;
    }
    case reduced::sech:
    {
        const std::int64_t ch = sinh_cosh(fractional_bits, raw).cosh_raw;
        return div_raw(one, ch, fractional_bits);
    }
    case reduced::csch:
    {
        if (raw == 0)
            throw std::domain_error("range lut: csch pole");
        const bool neg = raw < 0;
        const std::int64_t mag = abs_raw(raw);
        std::int64_t y;
        if (mag <= one)
        {
            const std::int64_t removed = eval_principal(principal::csch, fractional_bits, mag);
            y = removed + div_raw(one, mag, fractional_bits);
        }
        else
        {
            y = div_raw(one, sinh_cosh(fractional_bits, mag).sinh_raw, fractional_bits);
        }
        return neg ? -y : y;
    }
    case reduced::sqrt:
    {
        if (raw == 0)
            return 0;
        const dyadic part = split_positive(raw, fractional_bits);
        std::int64_t root = eval_principal(principal::sqrt, fractional_bits, part.mantissa_raw);
        if ((part.power & 1) != 0)
            root = mul_raw(root, scale_unit(sqrt2_64, fractional_bits), fractional_bits);
        return shift_pow2(root, half_pow_of(part.power));
    }
    case reduced::inv:
    {
        const dyadic part = split_positive(raw, fractional_bits);
        const std::int64_t seed = eval_principal(
            principal::inv, fractional_bits, part.mantissa_raw);
        const u128 wide = newton_inv(fractional_bits, part.mantissa_raw, seed);
        return finish_wide(wide, static_cast<int>(fractional_bits) + part.power);
    }
    case reduced::rsqrt:
    {
        const dyadic part = split_positive(raw, fractional_bits);
        const std::int64_t seed = eval_principal(
            principal::rsqrt, fractional_bits, part.mantissa_raw);
        u128 wide = newton_rsqrt(fractional_bits, part.mantissa_raw, seed);
        if ((part.power & 1) != 0)
            wide = round_u256(mul_u128(wide, rsqrt2_64), 64);
        return finish_wide(wide,
            static_cast<int>(fractional_bits) + half_pow_of(part.power));
    }
    case reduced::invsq:
    {
        const dyadic part = split_positive(raw, fractional_bits);
        const std::int64_t seed = eval_principal(
            principal::inv, fractional_bits, part.mantissa_raw);
        const u128 inv = newton_inv(fractional_bits, part.mantissa_raw, seed);
        const unsigned K = fractional_bits * 2u;
        const u128 wide = round_u256(mul_u128(inv, inv), K);
        return finish_wide(wide, static_cast<int>(K) - static_cast<int>(fractional_bits)
            + 2 * part.power);
    }
    case reduced::expm1:
        return eval_expm1(fractional_bits, raw);
    case reduced::log1p:
        return eval_log1p(fractional_bits, raw);
    }
    throw std::invalid_argument("range lut: unknown map");
}

}  // namespace grotto

#endif  // LIBDPF_INCLUDE_GROTTO_RANGE_LUT_HPP__