libdpf/include/grotto/exact_steps.hpp

/// @file grotto/exact_steps.hpp
/// @brief Exact Grotto steps: digit lengths, integer logs, and word bits.
/// @details Counts and boolean results are fixed-point integers,
///          `n << fractional_bits`. `ilog16(0)`, `ilog256(0)`, and
///          `logstar` of a non-positive input return `ilog_of_zero`.
///          `dec_width`, `oct_width`, and `b64_width` count digits of
///          `floor(|x|)` in bases 10, 8, and 64; zero has length 1.
///          `bit_width`, `bit_floor`, `bit_ceil`, `countl_one`, and
///          `has_single_bit` use the unsigned `Raw` bit pattern, matching
///          the C++ `<bit>` operations on that width.
///          `deg2rad` and `rad2deg` are the linear scale by `π/180` and
///          `180/π`. They do not wrap the angle.

#ifndef LIBDPF_INCLUDE_GROTTO_EXACT_STEPS_HPP__
#define LIBDPF_INCLUDE_GROTTO_EXACT_STEPS_HPP__

#include "grotto/dyadic_lut.hpp"

#include <cstdint>
#include <stdexcept>
#include <type_traits>
#include <vector>

namespace grotto
{

namespace exact_detail
{

inline int ceil_div(int numerator, int denominator)
{
    if (numerator >= 0)
        return (numerator + denominator - 1) / denominator;
    return -((-numerator) / denominator);
}

template <typename Raw>
int floor_log2_real(std::int64_t raw, unsigned fractional_bits, bool & exact_power)
{
    const detail::u128 mag = detail::magnitude<Raw>(raw);
    exact_power = mag != 0 && (mag & (mag - 1)) == 0;
    return detail::floor_log2_u128(mag) - static_cast<int>(fractional_bits);
}

inline std::int64_t abs_floor(std::int64_t raw, unsigned fractional_bits)
{
    if (fractional_bits >= 63)
        throw std::invalid_argument("exact step: fractional width does not fit");
    const detail::u128 mag = raw < 0
        ? detail::u128(-static_cast<__int128>(raw))
        : detail::u128(raw);
    return static_cast<std::int64_t>(mag >> fractional_bits);
}

inline int positive_length(std::int64_t magnitude, int base)
{
    if (magnitude <= 0)
        return 1;
    if (base == 2 || base == 8 || base == 64)
    {
        int log = 0;
        auto n = static_cast<std::uint64_t>(magnitude);
        while (n > 1)
        {
            n >>= 1;
            ++log;
        }
        const int group = base == 2 ? 1 : (base == 8 ? 3 : 6);
        return log / group + 1;
    }
    int digits = 0;
    while (magnitude > 0)
    {
        magnitude /= base;
        ++digits;
    }
    return digits;
}

}  // namespace exact_detail

/// @brief `ceil(log2(|x|) / 4)`, the integer log base 16.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_ilog16(std::int64_t raw, unsigned fractional_bits)
{
    if (raw == 0)
        return ilog_of_zero;
    bool exact = false;
    const int e = exact_detail::floor_log2_real<Raw>(raw, fractional_bits, exact);
    const int units = exact_detail::ceil_div(exact ? e : e + 1, 4);
    return detail::encode_units(units, fractional_bits);
}

/// @brief `ceil(log2(|x|) / 8)`, the integer log base 256.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_ilog256(std::int64_t raw, unsigned fractional_bits)
{
    if (raw == 0)
        return ilog_of_zero;
    bool exact = false;
    const int e = exact_detail::floor_log2_real<Raw>(raw, fractional_bits, exact);
    const int units = exact_detail::ceil_div(exact ? e : e + 1, 8);
    return detail::encode_units(units, fractional_bits);
}

/// @brief Iterated logarithm. `log*(x) = 0` for `|x| <= 1`, otherwise
///        `1 + log*(log2(|x|))`. Non-positive inputs return `ilog_of_zero`.
/// \complexity Five comparisons of the magnitude against `2^{fractional_bits + s}` for `s` in `{0,1,2,4,16}`. No iterated loop. `Θ(1)`, extra space `Θ(1)`.
/// Non-positive inputs return `ilog_of_zero` (the dyadic sentinel), not a mathematical log-star.
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_logstar(std::int64_t raw, unsigned fractional_bits)
{
    if (raw <= 0)
        return ilog_of_zero;
    if (fractional_bits >= 63)
        throw std::invalid_argument("exact step: fractional width does not fit");
    const detail::u128 mag = detail::magnitude<Raw>(raw);
    const auto below = [&](unsigned shift) {
        if (fractional_bits + shift >= 128)
            return true;
        return mag <= (detail::u128{1} << (fractional_bits + shift));
    };
    int units = 5;
    if (below(0))
        units = 0;
    else if (below(1))
        units = 1;
    else if (below(2))
        units = 2;
    else if (below(4))
        units = 3;
    else if (below(16))
        units = 4;
    return detail::encode_units(units, fractional_bits);
}

/// @brief `std::bit_width` of the unsigned `Raw` pattern.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_bit_width(std::int64_t raw, unsigned fractional_bits)
{
    const auto bits = detail::raw_bits<Raw>(raw);
    int units = 0;
    if (bits != 0)
        units = 64 - __builtin_clzll(bits);
    return detail::encode_units(units, fractional_bits);
}

/// @brief `std::countl_one` of the unsigned `Raw` pattern.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_countl_one(std::int64_t raw, unsigned fractional_bits)
{
    const int units = detail::countl_one_width(
        detail::raw_bits<Raw>(raw), detail::raw_width<Raw>());
    return detail::encode_units(units, fractional_bits);
}

/// @brief `std::has_single_bit` of the unsigned `Raw` pattern.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_has_single_bit(std::int64_t raw, unsigned fractional_bits)
{
    const auto bits = detail::raw_bits<Raw>(raw);
    const bool on = bits != 0 && (bits & (bits - 1)) == 0;
    return detail::encode_units(on ? 1 : 0, fractional_bits);
}

/// @brief `std::bit_floor` of the unsigned `Raw` pattern, as a raw word.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_bit_floor(std::int64_t raw)
{
    const auto bits = detail::raw_bits<Raw>(raw);
    if (bits == 0)
        return 0;
    const int log = 63 - __builtin_clzll(bits);
    return static_cast<std::int64_t>(std::uint64_t{1} << static_cast<unsigned>(log));
}

/// @brief `std::bit_ceil` of the unsigned `Raw` pattern, as a raw word.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
template <typename Raw>
std::int64_t eval_bit_ceil(std::int64_t raw)
{
    const auto bits = detail::raw_bits<Raw>(raw);
    if (bits <= 1)
        return 1;
    if ((bits & (bits - 1)) == 0)
        return static_cast<std::int64_t>(bits);
    const int log = 64 - __builtin_clzll(bits);
    if (log >= 63)
        throw std::overflow_error("exact step: bit_ceil does not fit int64");
    return std::int64_t{1} << log;
}

/// @brief `floor(x)` at the same fractional scale.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_dec_floor(std::int64_t raw, unsigned fractional_bits)
{
    if (fractional_bits >= 63)
        throw std::invalid_argument("exact step: fractional width does not fit");
    const std::int64_t one = std::int64_t{1} << fractional_bits;
    std::int64_t q = raw / one;
    const std::int64_t r = raw - q * one;
    if (r < 0)
        --q;
    return q << fractional_bits;
}

/// @brief `ceil(x)` at the same fractional scale.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_dec_ceil(std::int64_t raw, unsigned fractional_bits)
{
    if (fractional_bits >= 63)
        throw std::invalid_argument("exact step: fractional width does not fit");
    const std::int64_t one = std::int64_t{1} << fractional_bits;
    std::int64_t q = raw / one;
    const std::int64_t r = raw - q * one;
    if (r > 0)
        ++q;
    return q << fractional_bits;
}

/// @brief Digits of `floor(|x|)` in base 10, 8, or 64. Zero has length 1.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_value_length(std::int64_t raw, unsigned fractional_bits, int base)
{
    if (base != 8 && base != 10 && base != 64)
        throw std::invalid_argument("exact step: length base must be 8, 10, or 64");
    const int units = exact_detail::positive_length(
        exact_detail::abs_floor(raw, fractional_bits), base);
    return detail::encode_units(units, fractional_bits);
}
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16

inline std::int64_t eval_dec_width(std::int64_t raw, unsigned fractional_bits)
{
    return eval_value_length(raw, fractional_bits, 10);
}
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16

inline std::int64_t eval_oct_width(std::int64_t raw, unsigned fractional_bits)
{
    return eval_value_length(raw, fractional_bits, 8);
}
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16

inline std::int64_t eval_b64_width(std::int64_t raw, unsigned fractional_bits)
{
    return eval_value_length(raw, fractional_bits, 64);
}

/// @brief True when `floor(|x|)` is a single decimal digit, including 0.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_has_single_digit(std::int64_t raw, unsigned fractional_bits)
{
    const std::int64_t mag = exact_detail::abs_floor(raw, fractional_bits);
    return detail::encode_units(mag <= 9 ? 1 : 0, fractional_bits);
}

namespace exact_detail
{

inline constexpr detail::u128 pi_over_180_64 = 321956420358983237ULL;
inline constexpr detail::u128 rad_to_deg_64 = (detail::u128{57} << 64) | 5456168980075999427ULL;

inline std::int64_t scale_angle(detail::u128 mag64, unsigned fractional_bits)
{
    if (fractional_bits > 64)
        throw std::invalid_argument("exact step: fractional width does not fit");
    const unsigned shift = 64u - fractional_bits;
    detail::u128 mag = mag64;
    if (shift > 0)
    {
        mag += detail::u128{1} << (shift - 1);
        mag >>= shift;
    }
    if (mag > static_cast<detail::u128>(INT64_MAX))
        throw std::overflow_error("exact step: angle scale does not fit");
    return static_cast<std::int64_t>(mag);
}

inline std::int64_t mul_angle(std::int64_t raw, std::int64_t factor, unsigned fractional_bits)
{
    const __int128 prod = static_cast<__int128>(raw) * factor;
    const bool neg = prod < 0;
    auto mag = static_cast<detail::u128>(neg ? -prod : prod);
    if (fractional_bits > 0)
    {
        mag += detail::u128{1} << (fractional_bits - 1);
        mag >>= fractional_bits;
    }
    if (mag > static_cast<detail::u128>(INT64_MAX))
        throw std::overflow_error("exact step: angle does not fit int64");
    const auto out = static_cast<std::int64_t>(mag);
    return neg ? -out : out;
}

}  // namespace exact_detail

/// @brief Degrees to radians, `x * π / 180`.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_deg2rad(std::int64_t raw, unsigned fractional_bits)
{
    return exact_detail::mul_angle(
        raw, exact_detail::scale_angle(exact_detail::pi_over_180_64, fractional_bits),
        fractional_bits);
}

/// @brief Radians to degrees, `x * 180 / π`.
/// \complexity A constant number of shifts and a `countl` / `bit_width` style scan of the `Raw` width (at most 64 iterations of `n >>= 1` in `positive_length` / `floor_log2_real`). `Θ(width)` bit operations, extra space `Θ(1)`.
/// Counts are returned as `units << fractional_bits`.
/// @see grotto::make_exact_constant_lut
/// @see grotto::eval_ilog16
inline std::int64_t eval_rad2deg(std::int64_t raw, unsigned fractional_bits)
{
    return exact_detail::mul_angle(
        raw, exact_detail::scale_angle(exact_detail::rad_to_deg_64, fractional_bits),
        fractional_bits);
}

/// @brief Piecewise-constant LUT of a step whose value is an encoded integer.
/// @details The domain scan is the cut set. `Raw` wider than 16 bits is refused.
/// \complexity Scans every raw value from `numeric_limits<Raw>::min()` to `max()`. The header rejects `width > 16`, so the scan is `Θ(2^w)` with `w ≤ 16`.
/// The result stores one cut per change of `at`. Extra space is that cut vector.
/// @see grotto::eval_dec_floor
/// @see grotto::easy_lut
template <typename Raw, typename At>
easy_lut<Raw> make_exact_step_lut(unsigned fractional_bits, At && at)
{
    (void)fractional_bits;
    constexpr int width = detail::raw_width<Raw>();
    if (width > 16)
        throw std::invalid_argument("exact step lut: raw width must be at most 16");
    using lim = std::numeric_limits<Raw>;
    const auto minv = static_cast<std::int64_t>(lim::min());
    const auto maxv = static_cast<std::int64_t>(lim::max());
    std::vector<std::int64_t> cuts;
    auto previous = at(minv);
    for (std::int64_t raw = minv + 1; raw <= maxv; ++raw)
    {
        const auto value = at(raw);
        if (value != previous)
        {
            cuts.push_back(raw);
            previous = value;
        }
        if (raw == maxv)
            break;
    }
    return detail::steps_from_cuts<Raw>(std::move(cuts), [=](std::int64_t raw) {
        return detail::easy_poly{at(raw), 0, 0, 1};
    });
}

}  // namespace grotto

#endif  // LIBDPF_INCLUDE_GROTTO_EXACT_STEPS_HPP__