Skip to content

Repository navigation

Accelerate even-degree polynomial square-root recovery #185

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
odraccir1 merged 1 commit into from
Oct 26, 2026
Merged

Accelerate even-degree polynomial square-root recovery #185

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 102 additions & 0 deletions libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,108 @@ namespace nil::crypto3::math {
arithmetic_context);
}

/**
* Compute sum_{i=0}^{term_count-1} P_i * F^i(input) modulo B, where P_0 = 1 and
* P_i = weight * F(weight) * ... * F^(i-1)(weight). F fixes the coefficient field and the supplied canonical,
* reduced image is F(X). As with frobenius_product, its algebraic identity is a caller precondition and B need
* not be irreducible. An empty sum is zero; a singleton is the reduced input.
*
* Binary doubling uses O(log(term_count)) compositions and products. The sum, weight product and Frobenius
* image share each doubling cache; only the initial cache and one doubling cache are live at once. The last
* level computes only the sum. The canonical result may alias input, weight or the supplied image.
*
* @throws std::invalid_argument if an input is empty or noncanonical, the image is not reduced, a needed
* composition has a zero cached-power limit, or the divisor inverse has insufficient precision.
*/
template<detail::SupportsDivrem Backend>
void frobenius_weighted_sum(typename Backend::polynomial_type &output,
const typename Backend::polynomial_type &input,
const typename Backend::polynomial_type &weight, std::size_t term_count,
const typename Backend::polynomial_type &x_to_frobenius_power,
const polynomial_divisor_context<Backend> &divisor_context,
polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context) {
using polynomial_type = typename Backend::polynomial_type;
using value_type = typename polynomial_type::value_type;
using cache_type = polynomial_composition_precomputation<Backend>;
const std::size_t degree = divisor_context.degree();
const auto canonical = [](const polynomial_type &value) {
return !value.empty() && (value.size() == 1 || value.back() != value_type::zero());
};
if (!canonical(input) || !canonical(weight) || !canonical(x_to_frobenius_power) ||
x_to_frobenius_power.size() > std::max<std::size_t>(1, degree)) {
throw std::invalid_argument("Frobenius sums require canonical inputs and a reduced Frobenius image");
}
if (degree == 0 || term_count == 0) {
output.assign(1, value_type::zero());
return;
}
polynomial_type base, multiplier;
remainder(base, input, divisor_context, arithmetic_context);
if (term_count == 1 || is_zero(base)) {
output = std::move(base);
return;
}
remainder(multiplier, weight, divisor_context, arithmetic_context);
if (is_zero(multiplier)) {
output = std::move(base);
return;
}
if (base.size() == 1 && multiplier.size() == 1) {
// Coefficient-fixed Frobenius reduces this to a geometric sum, including weight = 1.
value_type sum = value_type::zero(), product = value_type::one();
value_type block_sum = value_type::one(), factor = multiplier[0];
for (std::size_t count = term_count; count != 0; count >>= 1) {
if (count % 2 != 0) {
sum += product * block_sum;
product *= factor;
}
block_sum += factor * block_sum;
factor = factor.squared();
}
output.assign(1, base[0] * sum);
return;
}

struct partial_sum {
polynomial_type sum, product, image;
};
const cache_type first_cache(x_to_frobenius_power, degree, divisor_context, arithmetic_context);
const auto chain = [&](auto &&self, std::size_t count, bool need_product_and_image) -> partial_sum {
if (count == 1) {
return {base, multiplier, x_to_frobenius_power};
}
auto result = self(self, count / 2, true);
std::optional<cache_type> doubling_cache;
if (count / 2 > 1) {
doubling_cache.emplace(result.image, degree, divisor_context, arithmetic_context);
}
const auto &cache = doubling_cache ? *doubling_cache : first_cache;
// W_2k = W_k + P_k * F^k(W_k). Form this before replacing P_k by P_2k.
polynomial_type shifted;
compose_mod(shifted, result.sum, cache, divisor_context, arithmetic_context);
mulmod(shifted, result.product, shifted, divisor_context, arithmetic_context);
addition(result.sum, result.sum, shifted);
if (need_product_and_image) {
compose_mod(shifted, result.product, cache, divisor_context, arithmetic_context);
mulmod(result.product, result.product, shifted, divisor_context, arithmetic_context);
compose_mod(result.image, result.image, cache, divisor_context, arithmetic_context);
}
if (count % 2 != 0) {
// Prepend a term: W_(2k+1) = input + weight * F(W_2k), reusing the initial cache.
compose_mod(shifted, result.sum, first_cache, divisor_context, arithmetic_context);
mulmod(shifted, multiplier, shifted, divisor_context, arithmetic_context);
addition(result.sum, base, shifted);
if (need_product_and_image) {
compose_mod(shifted, result.product, first_cache, divisor_context, arithmetic_context);
mulmod(result.product, multiplier, shifted, divisor_context, arithmetic_context);
compose_mod(result.image, result.image, first_cache, divisor_context, arithmetic_context);
}
}
return result;
};
output = std::move(chain(chain, term_count, false).sum);
}

} // namespace nil::crypto3::math

#endif // CRYPTO3_MATH_POLYNOMIAL_FROBENIUS_HPP
144 changes: 105 additions & 39 deletions libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,60 +253,53 @@ namespace nil::crypto3::math {

namespace detail {

// Return {X^((Q+1)/2), X^Q} modulo B. Called after the odd-degree root checks, with d > 1 and B(0) != 0.
// Multiply a canonical reduced representative by X, reducing the single possible leading term.
template<SupportsDivrem Backend>
void square_root_multiply_by_x(typename Backend::polynomial_type &value,
const polynomial_divisor_context<Backend> &divisor_context) {
using value_type = typename Backend::polynomial_type::value_type;
const std::size_t degree = divisor_context.degree();
shift_left(value, value, 1);
if (value.size() > degree) {
const value_type leading = value[degree] * divisor_context.reversed_divisor_inverse()[0];
value.resize(degree);
for (std::size_t i = 0; i < degree; ++i) {
value[i] -= leading * divisor_context.divisor()[i];
}
condense(value);
}
}

// Return {X^((Q-1)/2), X^Q} modulo B. Called after the root checks, with d > 1 and B(0) != 0.
template<SupportsDivrem Backend>
std::pair<typename Backend::polynomial_type, typename Backend::polynomial_type>
square_root_x_seed(const polynomial_divisor_context<Backend> &divisor_context,
polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context) {
using polynomial_type = typename Backend::polynomial_type;
using value_type = typename polynomial_type::value_type;
using field_type = typename value_type::field_type;
const auto &modulus = divisor_context.divisor();
const std::size_t degree = divisor_context.degree();
polynomial_type seed, x_to_q;
if constexpr (field_type::arity > 1 && requires(const value_type &value, std::size_t power) {
{ value.Frobenius_map(power) } -> std::convertible_to<value_type>;
}) {
const auto multiply_by_x = [&](polynomial_type &value) {
shift_left(value, value, 1);
if (value.size() > degree) {
const value_type leading = value[degree] * divisor_context.reversed_divisor_inverse()[0];
value.resize(degree);
for (std::size_t i = 0; i < degree; ++i) {
value[i] -= leading * modulus[i];
}
condense(value);
}
};
// For Q = p^e, a = x^((p-1)/2) gives x^p = x*a^2. Accumulate the e conjugates of a:
// T = a * a^p * ... * a^(p^(e-1)) = x^((Q-1)/2). Then b = x*T and x^Q = x*T^2.
// T = a * a^p * ... * a^(p^(e-1)) = x^((Q-1)/2). Then x^Q = x*T^2.
// The chain shares its composition caches between products and images, avoiding a Q-sized power.
const boost::multiprecision::cpp_int exponent =
(algebra::fields::field_characteristic<field_type>() - 1) >> 1;
polynomial_type characteristic_seed, x_to_p;
powmod_x(characteristic_seed, exponent, divisor_context, arithmetic_context);
squaremod(x_to_p, characteristic_seed, divisor_context, arithmetic_context);
multiply_by_x(x_to_p);
square_root_multiply_by_x(x_to_p, divisor_context);
frobenius_product_impl<true>(seed, characteristic_seed, field_type::arity, x_to_p, divisor_context,
arithmetic_context);
squaremod(x_to_q, seed, divisor_context, arithmetic_context);
multiply_by_x(x_to_q);
multiply_by_x(seed);
} else {
// Prime fields and coefficient types without a Frobenius map retain the direct seed.
const boost::multiprecision::cpp_int exponent = (algebra::fields::field_order<field_type>() + 1) >> 1;
polynomial_type seed_squared;
const boost::multiprecision::cpp_int exponent = (algebra::fields::field_order<field_type>() - 1) >> 1;
powmod_x(seed, exponent, divisor_context, arithmetic_context);
squaremod(seed_squared, seed, divisor_context, arithmetic_context);
// For a reduced y, y / X = (y - c * B) / X with c = y[0] / B[0], also for nonmonic B.
const value_type c = seed_squared[0] * modulus[0].inversed();
x_to_q.assign(degree, value_type::zero());
for (std::size_t i = 0; i < degree; ++i) {
x_to_q[i] =
(i + 1 < seed_squared.size() ? seed_squared[i + 1] : value_type::zero()) - c * modulus[i + 1];
}
condense(x_to_q);
}
squaremod(x_to_q, seed, divisor_context, arithmetic_context);
square_root_multiply_by_x(x_to_q, divisor_context);
return {std::move(seed), std::move(x_to_q)};
}

Expand Down Expand Up @@ -369,6 +362,7 @@ namespace nil::crypto3::math {
}

auto [seed, x_to_q] = detail::square_root_x_seed(divisor_context, arithmetic_context);
detail::square_root_multiply_by_x(seed, divisor_context);

const std::size_t term_count = (degree - 1) / 2;
polynomial_type product, x_to_q_squared;
Expand All @@ -385,17 +379,89 @@ namespace nil::crypto3::math {
frobenius_product(product, product, term_count, x_to_q_squared, divisor_context, arithmetic_context);
}

shift_left(product, product, 1);
if (product.size() > degree) {
// A shifted reduced polynomial exceeds degree(B) by at most one coefficient.
const value_type leading = product[degree] * inverse_leading;
product.resize(degree);
for (std::size_t i = 0; i < degree; ++i) {
product[i] -= leading * modulus[i];
detail::square_root_multiply_by_x(product, divisor_context);
scalar_multiplication(output, product, norm_root.inversed());
return true;
}

/**
* Compute a square root of X modulo an irreducible polynomial B of positive even degree d over an odd-order
* field K. Write Q = |K| and x = X mod B. The norm B(0) / leading_coefficient(B) first decides whether x is square.
*
* With lambda = x^((Q-1)/2), F(u) = u^Q and P_i = lambda * F(lambda) * ... * F^(i-1)(lambda), form
*
* W(c) = sum_{i=0}^{d-1} P_i * F^i(c), with P_0 = 1.
*
* For either root r of x, r * W(c) = Trace(r*c) belongs to K. Hence C = x * W(c)^2 is a coefficient-field square.
* When W(c) is nonzero, x * W(c) / sqrt(C) is a root of x. Try c = 1 first. On a zero trace, ask the caller-owned
* generator for another canonical reduced c and reuse lambda and X^Q. A uniformly sampled c has zero trace with
* probability 1/Q. This avoids nonresidue searches and exponentiations with exponents growing with d.
*
* The seed shares the odd-degree method's characteristic-sized exponentiation where coefficient Frobenius is
* available. Each trace uses O(log(d)) compositions and products; final multiplication by x is a linear pass.
* Return false and set output to zero if x is nonsquare; otherwise return a canonical reduced root.
*
* @pre B is irreducible. K provides is_square() and sqrt_known_square. The generator eventually supplies a c
* with nonzero weighted trace. Irreducibility is not checked; degree(B) - 1 inverse coefficients suffice.
* @throws std::invalid_argument if d is zero or odd, K has characteristic two, B(0) is zero, a generated c is
* empty, noncanonical or unreduced, a needed composition has a zero cached-power limit, or the divisor
* inverse has insufficient precision.
* @throws std::logic_error if a nonzero weighted trace does not produce a nonzero square in K.
*/
template<detail::SupportsDivrem Backend, typename Generator>
requires algebra::FieldValue<typename Backend::polynomial_type::value_type> &&
requires(const typename Backend::polynomial_type::value_type &value, Generator &generator) {
{ value.is_square() } -> std::convertible_to<bool>;
{ generator() } -> std::convertible_to<typename Backend::polynomial_type>;
}
bool square_root_mod_x_even_degree(typename Backend::polynomial_type &output,
const polynomial_divisor_context<Backend> &divisor_context,
polynomial_arithmetic::polynomial_context<Backend> &arithmetic_context,
Generator &generator) {
using polynomial_type = typename Backend::polynomial_type;
using value_type = typename polynomial_type::value_type;
using field_type = typename value_type::field_type;
const auto &modulus = divisor_context.divisor();
const std::size_t degree = divisor_context.degree();
if (degree == 0 || degree % 2 != 0 || algebra::fields::field_characteristic<field_type>() == 2) {
throw std::invalid_argument(
"trace square roots of X require a positive even degree and odd characteristic");
}
if (modulus[0] == value_type::zero()) {
throw std::invalid_argument("an irreducible divisor of degree greater than one has nonzero constant term");
}
const value_type norm = modulus[0] * divisor_context.reversed_divisor_inverse()[0];
if (!norm.is_square()) {
output.assign(1, value_type::zero());
return false;
}
if (arithmetic_context.options().modular_composition_cached_power_limit == 0) {
throw std::invalid_argument("trace square roots of X require a positive composition cached-power limit");
}

const auto [weight, image] = detail::square_root_x_seed(divisor_context, arithmetic_context);
polynomial_type candidate {value_type::one()}, sum;
for (;;) {
frobenius_weighted_sum(sum, candidate, weight, degree, image, divisor_context, arithmetic_context);
if (!is_zero(sum)) {
break;
}
candidate = generator();
if (candidate.empty() || candidate.size() > degree ||
(candidate.size() > 1 && candidate.back() == value_type::zero())) {
throw std::invalid_argument("a trace square-root candidate must be a canonical reduced representative");
}
condense(product);
}
scalar_multiplication(output, product, norm_root.inversed());

polynomial_type constant;
squaremod(constant, sum, divisor_context, arithmetic_context);
detail::square_root_multiply_by_x(constant, divisor_context);
if (constant.size() != 1 || constant[0] == value_type::zero() || !constant[0].is_square()) {
throw std::logic_error("a nonzero weighted trace must yield a nonzero coefficient-field square");
}
const value_type scalar = algebra::fields::sqrt_known_square(constant[0]).inversed();
detail::square_root_multiply_by_x(sum, divisor_context);
scalar_multiplication(output, sum, scalar);
return true;
}

Expand Down
Loading
Loading