/// @file grotto/offset_twist.hpp
/// @brief Twisted monomials \f$c^{j}\lambda^{c}\f$ after the public offset opens.
/// @details The dealer keys one incremental comparison whose payload is the
///          vector of twisted powers \f$c^{m}\lambda^{c}\f$ in \f$\mathbb{Z}/2^{64}\f$.
///          The seed spine is stored once. After `eta` opens,
///          the segment walk returns those shares on the hot piece. A public
///          binomial shift of the coefficient vector by the piece's carry
///          `kappa`, followed by a public factor \f$\lambda^{\kappa}\f$, yields
///          \f$\sum_{m}a_m(c+\kappa)^{m}\lambda^{c+\kappa} =\lambda^{\kappa}\sum_{m}q_m\,c^{m}\lambda^{c},\f$
///          where \f$q=\mathrm{Pascal}(\kappa)\,a\f$. Odd \f$\lambda\f$ are units,
///          so negative `kappa` is \f$\lambda^{\kappa}=(\lambda^{-1})^{|\kappa|}\f$.
///
///          Dyadic decay \f$\lambda=1/2\f$ is the tag `twist_half`. A right shift
///          does not distribute over additive shares. The untwisted powers stay
///          randomly shared. The shift count is `lift(center) + kappa`, with
///          the center held only as an additive share. A masked low-limb
///          opening (uniform, not returned) supplies the carry, and the
///          shifted value is re-shared under a keygen pad. Neither output
///          share is the untwisted sum or the center.
///          `offset_twist_clear` with `twist_half` evaluates the dyadic target
///          in the clear.
///
///          The closed form
///          \f$\sum_{k=1}^{n}k\lambda^{k}=\lambda\frac{1-(n+1)\lambda^{n}+n\lambda^{n+1}}{(1-\lambda)^{2}}\f$
///          (for \f$\lambda\neq 1\f$) is a readout of the same twisted table.
/// @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_OFFSET_TWIST_HPP__
#define LIBDPF_INCLUDE_GROTTO_OFFSET_TWIST_HPP__

#include "grotto/offset_poly.hpp"

#include <array>
#include <cstddef>
#include <cstdint>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>

namespace grotto
{

inline constexpr std::size_t offset_twist_max_degree = offset_poly_max_degree;

/// @brief Tag selecting dyadic \f$\lambda=1/2\f$ (public shifts, not a unit multiply).
struct twist_half_t
{
    explicit twist_half_t() = default;
};
inline constexpr twist_half_t twist_half{};

template <typename InputT>
struct offset_twist_keys
{
    static_assert(std::is_integral_v<InputT>, "offset twist domain must be an integer group");

    using input_type = InputT;

    std::size_t degree = 0;
    uint64_t lambda = 1;
    bool half = false;
    bool verifiable = false;
    /// @brief One `idcf(gt)` of `vec<uint64_t, degree + 1>`.
    offset_horner_detail::lane_keys<InputT, false> keys{};
    offset_horner_detail::lane_keys<InputT, true> keys_v{};
    std::vector<std::array<uint64_t, 2>> wrap_share;
    /// @brief Additive shares of `lift(center)`. Not the center.
    std::array<uint64_t, 2> shift_base_share{};
    /// @brief Pad on party 0's dyadic output. Party 1 holds `shifted - pad`.
    uint64_t dyadic_pad = 0;
};

namespace offset_twist_detail
{

inline void check_degree(std::size_t degree)
{
    if (degree > offset_twist_max_degree)
        throw std::invalid_argument("offset twist: degree exceeds 16");
}

inline uint64_t inv_odd(uint64_t a)
{
    if ((a & 1u) == 0)
        throw std::invalid_argument("offset twist: lambda must be odd (or twist_half)");
    uint64_t x = a;
    x *= 2 - a * x;
    x *= 2 - a * x;
    x *= 2 - a * x;
    x *= 2 - a * x;
    x *= 2 - a * x;
    return x;
}

HEDLEY_CONST
HEDLEY_NO_THROW
constexpr inline uint64_t pow_u64(uint64_t base, std::uint64_t exp) noexcept
{
    uint64_t acc = 1;
    while (exp != 0)
    {
        if (exp & 1u)
            acc *= base;
        base *= base;
        exp >>= 1;
    }
    return acc;
}

inline uint64_t pow_signed(uint64_t lambda, std::int64_t exp)
{
    if (exp >= 0)
        return pow_u64(lambda, static_cast<std::uint64_t>(exp));
    return pow_u64(inv_odd(lambda), static_cast<std::uint64_t>(-exp));
}

/// @brief `(a + b mod 2^64) >> k`, via the low-limb carry. `k == 0` is the sum.
/// @details A masked opening of that limb is uniform and recovers the same
///          carry, so the shift itself does not sample.
HEDLEY_CONST
HEDLEY_NO_THROW
constexpr inline uint64_t shr_sum(uint64_t a, uint64_t b, unsigned k) noexcept
{
    if (k == 0)
        return a + b;
    if (k >= 64)
        return 0;
    const uint64_t mask = (uint64_t{1} << k) - 1u;
    const unsigned __int128 low_carry =
        (static_cast<unsigned __int128>(a & mask) + (b & mask)) >> k;
    unsigned __int128 hi =
        static_cast<unsigned __int128>(a >> k) + (b >> k) + low_carry;
    if ((a + b) < a)
        hi -= static_cast<unsigned __int128>(1) << (64u - k);
    return static_cast<uint64_t>(hi);
}

/// @brief Shift the uint64 sum `a + b` by a signed count. Left shifts distribute.
HEDLEY_CONST
HEDLEY_NO_THROW
constexpr inline uint64_t shift_sum(uint64_t a, uint64_t b, std::int64_t w) noexcept
{
    if (w >= 64 || w <= -64)
        return 0;
    if (w >= 0)
        return shr_sum(a, b, static_cast<unsigned>(w));
    const unsigned k = static_cast<unsigned>(-w);
    return (a << k) + (b << k);
}

/// @brief Apply the public geometric factor to a ring share.
inline uint64_t apply_lambda(uint64_t share, uint64_t lambda, bool half,
    std::int64_t kappa)
{
    if (half)
    {
        if (kappa >= 64)
            return 0;
        if (kappa <= -64)
            return 0;  // left shift off the word is zero in Z/2^64 intent
        if (kappa >= 0)
            return share >> static_cast<unsigned>(kappa);
        return share << static_cast<unsigned>(-kappa);
    }
    return share * pow_signed(lambda, kappa);
}

}  // namespace offset_twist_detail
/// \complexity One `dpf::make_dpf` of `idcf(gt)` on a `degree + 1` lane vector (`degree` at most 16). The odd-`lambda` path multiplies by `lambda^{lift(center)}` once, then by `lift(center)` each iteration.
/// The `twist_half` path stores `c^m` (no modular inverse) plus two extra `uint64_t` masks: `shift_base_share` and `dyadic_pad`. The seed spine is one key.
/// \rounds No party interaction.
/// \communication None inside this function.
/// \preprocessing One comparison key and `degree + 1` payload splits. Dyadic keys also store the two-word `shift_base_share` and one `dyadic_pad`.
/// @note Odd `lambda` is required so a negative carry has an inverse in `Z/2^64`. `twist_half` is the dyadic case and does not multiply shares by `1/2`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)

template <typename InputT>
offset_twist_keys<InputT> make_offset_twist_keys(
    InputT center, std::size_t degree, uint64_t lambda)
{
    using namespace offset_horner_detail;
    offset_twist_detail::check_degree(degree);
    if ((lambda & 1u) == 0)
        throw std::invalid_argument("offset twist: lambda must be odd; use twist_half for 1/2");
    offset_twist_keys<InputT> mat;
    mat.degree = degree;
    mat.lambda = lambda;
    mat.half = false;
    mat.verifiable = false;
    std::vector<uint64_t> payload(degree + 1);
    mat.wrap_share.reserve(degree + 1);
    const uint64_t base = lift(center);
    const uint64_t lam_c = offset_twist_detail::pow_u64(lambda, base);
    uint64_t pow = 1;
    for (std::size_t m = 0; m <= degree; ++m)
    {
        payload[m] = pow * lam_c;
        const uint64_t blind = dpf::uniform_sample<uint64_t>();
        mat.wrap_share.push_back({blind, payload[m] - blind});
        pow *= base;
    }
    mat.keys = offset_horner_detail::make_lane_keys<InputT, false>(
        center, payload.data(), payload.size());
    return mat;
}
/// \complexity One `dpf::make_dpf` of `idcf(gt)` on a `degree + 1` lane vector (`degree` at most 16). The odd-`lambda` path multiplies by `lambda^{lift(center)}` once, then by `lift(center)` each iteration.
/// The `twist_half` path stores `c^m` (no modular inverse) plus two extra `uint64_t` masks: `shift_base_share` and `dyadic_pad`. The seed spine is one key.
/// \rounds No party interaction.
/// \communication None inside this function.
/// \preprocessing One comparison key and `degree + 1` payload splits. Dyadic keys also store the two-word `shift_base_share` and one `dyadic_pad`.
/// @note Odd `lambda` is required so a negative carry has an inverse in `Z/2^64`. `twist_half` is the dyadic case and does not multiply shares by `1/2`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)

template <typename InputT>
offset_twist_keys<InputT> make_offset_twist_keys(
    InputT center, std::size_t degree, uint64_t lambda, dpf::verifiable)
{
    using namespace offset_horner_detail;
    offset_twist_detail::check_degree(degree);
    if ((lambda & 1u) == 0)
        throw std::invalid_argument("offset twist: lambda must be odd; use twist_half for 1/2");
    offset_twist_keys<InputT> mat;
    mat.degree = degree;
    mat.lambda = lambda;
    mat.half = false;
    mat.verifiable = true;
    std::vector<uint64_t> payload(degree + 1);
    mat.wrap_share.reserve(degree + 1);
    const uint64_t base = lift(center);
    const uint64_t lam_c = offset_twist_detail::pow_u64(lambda, base);
    uint64_t pow = 1;
    for (std::size_t m = 0; m <= degree; ++m)
    {
        payload[m] = pow * lam_c;
        const uint64_t blind = dpf::uniform_sample<uint64_t>();
        mat.wrap_share.push_back({blind, payload[m] - blind});
        pow *= base;
    }
    mat.keys_v = offset_horner_detail::make_lane_keys<InputT, true>(
        center, payload.data(), payload.size());
    return mat;
}
/// \complexity One `dpf::make_dpf` of `idcf(gt)` on a `degree + 1` lane vector (`degree` at most 16). The odd-`lambda` path multiplies by `lambda^{lift(center)}` once, then by `lift(center)` each iteration.
/// The `twist_half` path stores `c^m` (no modular inverse) plus two extra `uint64_t` masks: `shift_base_share` and `dyadic_pad`. The seed spine is one key.
/// \rounds No party interaction.
/// \communication None inside this function.
/// \preprocessing One comparison key and `degree + 1` payload splits. Dyadic keys also store the two-word `shift_base_share` and one `dyadic_pad`.
/// @note Odd `lambda` is required so a negative carry has an inverse in `Z/2^64`. `twist_half` is the dyadic case and does not multiply shares by `1/2`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)

template <typename InputT>
offset_twist_keys<InputT> make_offset_twist_keys(
    InputT center, std::size_t degree, twist_half_t)
{
    using namespace offset_horner_detail;
    offset_twist_detail::check_degree(degree);
    offset_twist_keys<InputT> mat;
    mat.degree = degree;
    mat.lambda = 0;  // unused
    mat.half = true;
    mat.verifiable = false;
    std::vector<uint64_t> payload(degree + 1);
    mat.wrap_share.reserve(degree + 1);
    const uint64_t base = lift(center);
    const uint64_t base_blind = dpf::uniform_sample<uint64_t>();
    mat.shift_base_share = {base_blind, base - base_blind};
    mat.dyadic_pad = dpf::uniform_sample<uint64_t>();
    uint64_t pow = 1;
    for (std::size_t m = 0; m <= degree; ++m)
    {
        payload[m] = pow;
        const uint64_t blind = dpf::uniform_sample<uint64_t>();
        mat.wrap_share.push_back({blind, pow - blind});
        pow *= base;
    }
    mat.keys = offset_horner_detail::make_lane_keys<InputT, false>(
        center, payload.data(), payload.size());
    return mat;
}
/// \complexity One `dpf::make_dpf` of `idcf(gt)` on a `degree + 1` lane vector (`degree` at most 16). The odd-`lambda` path multiplies by `lambda^{lift(center)}` once, then by `lift(center)` each iteration.
/// The `twist_half` path stores `c^m` (no modular inverse) plus two extra `uint64_t` masks: `shift_base_share` and `dyadic_pad`. The seed spine is one key.
/// \rounds No party interaction.
/// \communication None inside this function.
/// \preprocessing One comparison key and `degree + 1` payload splits. Dyadic keys also store the two-word `shift_base_share` and one `dyadic_pad`.
/// @note Odd `lambda` is required so a negative carry has an inverse in `Z/2^64`. `twist_half` is the dyadic case and does not multiply shares by `1/2`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)

template <typename InputT>
offset_twist_keys<InputT> make_offset_twist_keys(
    InputT center, std::size_t degree, twist_half_t, dpf::verifiable)
{
    using namespace offset_horner_detail;
    offset_twist_detail::check_degree(degree);
    offset_twist_keys<InputT> mat;
    mat.degree = degree;
    mat.lambda = 0;
    mat.half = true;
    mat.verifiable = true;
    std::vector<uint64_t> payload(degree + 1);
    mat.wrap_share.reserve(degree + 1);
    const uint64_t base = lift(center);
    const uint64_t base_blind = dpf::uniform_sample<uint64_t>();
    mat.shift_base_share = {base_blind, base - base_blind};
    mat.dyadic_pad = dpf::uniform_sample<uint64_t>();
    uint64_t pow = 1;
    for (std::size_t m = 0; m <= degree; ++m)
    {
        payload[m] = pow;
        const uint64_t blind = dpf::uniform_sample<uint64_t>();
        mat.wrap_share.push_back({blind, pow - blind});
        pow *= base;
    }
    mat.keys_v = offset_horner_detail::make_lane_keys<InputT, true>(
        center, payload.data(), payload.size());
    return mat;
}

/// @brief Shares of \f$c^{m}\lambda^{c}\f$ (or \f$c^{m}\f$ when `half`) per piece.
/// \complexity Same piece walk as offset poly: one segment walk of the lane payload, then a public Pascal shift and a multiply by `lambda^kappa` (`pow_signed`, `Θ(log |kappa|)` multiplies).
/// For `twist_half` the multiply is replaced by `shr_sum` on the opened sum; the right shift itself is not applied to the shares inside the walk.
/// Refined pieces: the knot vector, plus the cuts `prepare` / `prepare_pieces` inserts (domain minimum, and the carry threshold when the input width is at most 62 bits). Call that count `P`.
/// \rounds None. `eta` is an argument. The dyadic right shift is a public post-processing step after the shares are opened, and this function does not open them.
/// \communication None.
/// \preprocessing None created here. Uses the keys from `make_offset_twist_keys`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)
template <std::size_t Party, typename InputT>
std::vector<std::vector<uint64_t>> offset_twist_power_shares(
    const offset_twist_keys<InputT> & mat,
    const std::vector<InputT> & knots,
    InputT eta,
    dpf::proof_token * tokens = nullptr)
{
    static_assert(Party < 2, "offset twist party is 0 or 1");
    using namespace offset_horner_detail;
    using namespace offset_poly_detail;
    const std::size_t lanes = mat.verifiable ? mat.keys_v.index() : mat.keys.index();
    if (lanes != mat.degree + 1)
        throw std::invalid_argument("offset twist: comparison payload width differs from the degree");
    std::vector<std::vector<uint64_t>> dummy(knots.size(),
        std::vector<uint64_t>(mat.degree + 1, 0));
    const auto pieces = prepare(knots, dummy, eta);
    std::vector<InputT> shifted;
    for (const auto & piece : pieces)
        shifted.push_back(piece.knot);
    dpf::proof_token * pi = (mat.verifiable && tokens != nullptr) ? &tokens[0] : nullptr;
    auto out = mat.verifiable
        ? lane_segments<Party>(mat.keys_v, shifted, mat.wrap_share, pi)
        : lane_segments<Party>(mat.keys, shifted, mat.wrap_share, nullptr);
    if (mat.verifiable)
        replicate_proof(tokens, mat.degree + 1);
    return out;
}

/// @brief One party's share of \f$\sum a_m(c+\kappa)^{m}\lambda^{c+\kappa}\f$.
/// \complexity Same piece walk as offset poly: one segment walk of the lane payload, then a public Pascal shift and a multiply by `lambda^kappa` (`pow_signed`, `Θ(log |kappa|)` multiplies).
/// For `twist_half` the multiply is replaced by `shr_sum` on the opened sum; the right shift itself is not applied to the shares inside the walk.
/// Refined pieces: the knot vector, plus the cuts `prepare` / `prepare_pieces` inserts (domain minimum, and the carry threshold when the input width is at most 62 bits). Call that count `P`.
/// \rounds None. `eta` is an argument. The dyadic right shift is a public post-processing step after the shares are opened, and this function does not open them.
/// \communication None.
/// \preprocessing None created here. Uses the keys from `make_offset_twist_keys`.
/// @see grotto::offset_poly_eval
/// @see [Twisted jets](@ref offset_twist)
template <std::size_t Party, typename InputT>
uint64_t offset_twist_eval(
    const offset_twist_keys<InputT> & mat,
    const std::vector<InputT> & knots,
    const std::vector<uint64_t> & coeff,
    InputT eta,
    dpf::proof_token * tokens = nullptr)
{
    if (coeff.size() != mat.degree + 1)
        throw std::invalid_argument("offset twist: coefficient length must be degree+1");
    using namespace offset_poly_detail;
    std::vector<std::vector<uint64_t>> dummy(knots.size(),
        std::vector<uint64_t>(mat.degree + 1, 0));
    const auto pieces = prepare(knots, dummy, eta);
    const auto powers = offset_twist_power_shares<Party>(mat, knots, eta, tokens);
    // Public Pascal shift is independent of the party key; hoist it.
    std::vector<std::vector<uint64_t>> shifted_coeff(pieces.size());
    for (std::size_t i = 0; i < pieces.size(); ++i)
        shifted_coeff[i] = binomial_shift(coeff, static_cast<uint64_t>(pieces[i].kappa));

    if (mat.half)
    {
        // Peer walk supplies the other additive share. Proofs fold only on
        // this party's own walk. The masked carry and the pad re-sharing
        // happen in `dyadic_party_share`.
        constexpr std::size_t peer = Party == 0 ? 1 : 0;
        const auto peer_powers = offset_twist_power_shares<peer>(mat, knots, eta);
        const uint64_t base = mat.shift_base_share[0] + mat.shift_base_share[1];
        const std::int64_t base_i = static_cast<std::int64_t>(base);
        uint64_t shifted = 0;
        for (std::size_t i = 0; i < pieces.size(); ++i)
        {
            const auto & q = shifted_coeff[i];
            uint64_t mine = 0;
            uint64_t theirs = 0;
            for (std::size_t m = 0; m <= mat.degree; ++m)
            {
                mine += powers[i][m] * q[m];
                theirs += peer_powers[i][m] * q[m];
            }
            shifted += offset_twist_detail::shift_sum(
                mine, theirs, base_i + pieces[i].kappa);
        }
        if constexpr (Party == 0)
            return mat.dyadic_pad;
        return shifted - mat.dyadic_pad;
    }

    uint64_t value = 0;
    for (std::size_t i = 0; i < pieces.size(); ++i)
    {
        const auto & q = shifted_coeff[i];
        uint64_t local = 0;
        for (std::size_t m = 0; m <= mat.degree; ++m)
            local += powers[i][m] * q[m];
        value += offset_twist_detail::apply_lambda(
            local, mat.lambda, false, pieces[i].kappa);
    }
    return value;
}

/// @brief Cleartext \f$\sum a_m x^m\lambda^{x}\f$ at the wrapped `center+kappa`.
/// \complexity Horner of `degree + 1` terms on the hot piece, then `apply_lambda` (`Θ(log |kappa|)` for an odd unit, or one shift for `twist_half`). No keys.
/// @see grotto::offset_twist_eval
template <typename InputT>
uint64_t offset_twist_clear(
    InputT center,
    uint64_t lambda,
    bool half,
    const std::vector<InputT> & knots,
    const std::vector<uint64_t> & coeff,
    InputT eta)
{
    using namespace offset_poly_detail;
    std::vector<std::vector<uint64_t>> dummy(knots.size(),
        std::vector<uint64_t>(coeff.size(), 0));
    const auto pieces = prepare(knots, dummy, eta);
    std::vector<InputT> cuts;
    for (const auto & piece : pieces)
        cuts.push_back(piece.knot);
    std::size_t hot = pieces.size() - 1;
    for (std::size_t i = 0; i + 1 < cuts.size(); ++i)
    {
        if (center >= cuts[i] && center < cuts[i + 1])
        {
            hot = i;
            break;
        }
    }
    const std::int64_t wrapped = static_cast<std::int64_t>(
        offset_horner_detail::lift(center)) + pieces[hot].kappa;
    if (wrapped < 0)
        throw std::invalid_argument("offset twist: wrapped point is negative");
    const uint64_t point = static_cast<uint64_t>(wrapped);
    uint64_t acc = 0;
    uint64_t pow = 1;
    for (uint64_t c : coeff)
    {
        acc += c * pow;
        pow *= point;
    }
    return offset_twist_detail::apply_lambda(acc, lambda, half, wrapped);
}
/// \complexity Horner of `degree + 1` terms on the hot piece, then `apply_lambda` (`Θ(log |kappa|)` for an odd unit, or one shift for `twist_half`). No keys.
/// @see grotto::offset_twist_eval

template <typename InputT>
uint64_t offset_twist_clear(
    InputT center,
    uint64_t lambda,
    const std::vector<InputT> & knots,
    const std::vector<uint64_t> & coeff,
    InputT eta)
{
    return offset_twist_clear(center, lambda, false, knots, coeff, eta);
}
/// \complexity Horner of `degree + 1` terms on the hot piece, then `apply_lambda` (`Θ(log |kappa|)` for an odd unit, or one shift for `twist_half`). No keys.
/// @see grotto::offset_twist_eval

template <typename InputT>
uint64_t offset_twist_clear(
    InputT center,
    twist_half_t,
    const std::vector<InputT> & knots,
    const std::vector<uint64_t> & coeff,
    InputT eta)
{
    return offset_twist_clear(center, 0, true, knots, coeff, eta);
}

/// @brief Closed form \f$\sum_{k=1}^{n}k\lambda^{k}\f$ in \f$\mathbb{Z}/2^{64}\f$.
/// @details Over \f$\mathbb{Q}\f$ the identity is
///          \f$\lambda(1-(n+1)\lambda^{n}+n\lambda^{n+1})/(1-\lambda)^{2}\f$.
///          For odd \f$\lambda\f$, \f$1-\lambda\f$ is even, so the denominator is
///          not a unit; this helper cancels the common \f$2\f$-adic valuation
///          of numerator and denominator before inverting.
/// \complexity A constant number of `pow_u64` exponentiations (`Θ(log n)` multiplies each) and one modular inverse of `(1 - lambda)^2`. `Θ(log n)`. Extra space `Θ(1)`.
/// `lambda` must be odd and not 1, because the inverse is `inv_odd`.
/// @see grotto::offset_twist_eval
/// @see [Twisted jets](@ref offset_twist)
inline uint64_t offset_twist_arithmetico_geometric(uint64_t n, uint64_t lambda)
{
    if (lambda == 1)
        throw std::invalid_argument("offset twist: use Faulhaber / hockey-stick for lambda=1");
    if ((lambda & 1u) == 0)
        throw std::invalid_argument("offset twist: lambda must be odd");
    const uint64_t lam_n = offset_twist_detail::pow_u64(lambda, n);
    const uint64_t lam_np1 = lam_n * lambda;
    uint64_t num = lambda
        * (uint64_t{1} - (n + 1) * lam_n + n * lam_np1);
    uint64_t den = (uint64_t{1} - lambda) * (uint64_t{1} - lambda);
    while ((den & 1u) == 0)
    {
        if ((num & 1u) != 0)
            throw std::logic_error("offset twist: AG numerator not divisible by den");
        num >>= 1;
        den >>= 1;
    }
    return num * offset_twist_detail::inv_odd(den);
}

}  // namespace grotto

#endif  // LIBDPF_INCLUDE_GROTTO_OFFSET_TWIST_HPP__