Ryan Henry
0d22946a0e
Ship the TLS mesh, composer, Beaver/Yao/leaf MPC, prep/online paths, apps, and docs so the tree is pushable before elevating share_expr, security_mode, and prep resume. Co-authored-by: Cursor <cursoragent@cursor.com>
533 lines
23 KiB
C++
533 lines
23 KiB
C++
/// @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__
|