From fdcb75a0d5025faf4cee898c52eb55568133003f Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Fri, 23 Oct 2026 18:19:00 -0400 Subject: [PATCH 1/2] Accelerate odd-degree polynomial square-root recovery --- .../quotient_ring/polynomial_frobenius.hpp | 85 +++++++++++ .../quotient_ring/polynomial_square_root.hpp | 100 +++++++++++++ .../polynomial_x_norm_reconstruction.hpp | 64 ++++---- libs/math/test/polynomial_frobenius.cpp | 140 ++++++++++++++++++ libs/math/test/polynomial_square_root.cpp | 136 +++++++++++++++++ .../test/polynomial_x_norm_reconstruction.cpp | 118 ++++++++++++++- 6 files changed, 615 insertions(+), 28 deletions(-) diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp index f6616708a..0f016156d 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp @@ -27,6 +27,8 @@ #include #include +#include +#include #include #include @@ -147,6 +149,89 @@ namespace nil::crypto3::math { output = std::move(result); } + /** + * Compute input * F(input) * ... * F^(term_count - 1)(input) modulo B, where F is a power of the coefficient + * field's Frobenius map. The supplied x_to_frobenius_power is the canonical reduced image F(X); it can represent + * X^(Q^s) for any s >= 0, with Q the coefficient-field order. Its algebraic identity is a caller precondition, + * not checked by repeating the exponentiation. B need not be irreducible. + * + * Binary doubling uses O(log(term_count)) compositions and products. Each doubling shares one composition + * precomputation between the partial product and the next Frobenius image. Odd steps reuse the initial image's + * precomputation. Only the initial cache and one doubling cache are live at once. + * + * An empty product is the quotient-ring identity, and a singleton is the reduced input. The result is canonical + * and may alias either input polynomial. Reducing an unreduced input may require more inverse coefficients than + * products of reduced operands; otherwise degree(B) - 1 inverse coefficients suffice. + * + * @throws std::invalid_argument if an input is empty or noncanonical, the supplied image is not reduced, a + * needed composition has a zero cached-power limit, or the divisor inverse has insufficient precision. + */ + template + void frobenius_product(typename Backend::polynomial_type &output, const typename Backend::polynomial_type &input, + std::size_t term_count, const typename Backend::polynomial_type &x_to_frobenius_power, + const polynomial_divisor_context &divisor_context, + polynomial_arithmetic::polynomial_context &arithmetic_context) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename polynomial_type::value_type; + using cache_type = polynomial_composition_precomputation; + 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(x_to_frobenius_power) || + x_to_frobenius_power.size() > std::max(1, degree)) { + throw std::invalid_argument("Frobenius products require canonical inputs and a reduced Frobenius image"); + } + if (degree == 0 || term_count == 0) { + output.assign(1, degree == 0 ? value_type::zero() : value_type::one()); + return; + } + polynomial_type base; + remainder(base, input, divisor_context, arithmetic_context); + if (term_count == 1) { + output = std::move(base); + return; + } + // Every coefficient-field constant is fixed by F, including zero. + if (base.size() == 1) { + powmod(output, base, term_count, divisor_context, arithmetic_context); + return; + } + + 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_image) -> std::pair { + if (count == 1) { + return {base, x_to_frobenius_power}; + } + auto [product, image] = self(self, count / 2, true); + std::optional doubling_cache; + if (count / 2 > 1) { + doubling_cache.emplace(image, degree, divisor_context, arithmetic_context); + } + const auto &cache = doubling_cache ? *doubling_cache : first_cache; + // With P_k = product and H_k = image, P_2k = P_k * P_k(H_k) and H_2k = H_k(H_k). + polynomial_type shifted; + compose_mod(shifted, product, cache, divisor_context, arithmetic_context); + mulmod(product, product, shifted, divisor_context, arithmetic_context); + if (need_image) { + compose_mod(image, image, cache, divisor_context, arithmetic_context); + } + if (count % 2 != 0) { + // P_(2k+1) = base * F(P_2k); use the same initial cache at every odd step. + compose_mod(shifted, product, first_cache, divisor_context, arithmetic_context); + mulmod(product, base, shifted, divisor_context, arithmetic_context); + if (need_image) { + compose_mod(image, image, first_cache, divisor_context, arithmetic_context); + } + } + return {std::move(product), std::move(image)}; + }; + // The final image is not used by the caller, so avoid its compositions. + auto result = chain(chain, term_count, false); + output = std::move(result.first); + } + } // namespace nil::crypto3::math #endif // CRYPTO3_MATH_POLYNOMIAL_FROBENIUS_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp index 2fd072b0b..89c361c88 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp @@ -32,10 +32,13 @@ #include +#include #include #include +#include #include +#include namespace nil::crypto3::math { @@ -248,6 +251,103 @@ namespace nil::crypto3::math { return detail::is_square_mod_impl(input, divisor_context, arithmetic_context, nullptr); } + /** + * Compute a square root of X modulo an irreducible polynomial B of odd degree d over an odd-order field K. + * Write Q = |K|, x = X mod B, and N = Norm(x) = -B(0) / leading_coefficient(B). The element x is square exactly + * when N is square in K. For d = 1 the answer is simply sqrt(N), including zero when B is a multiple of X. + * + * For d > 1, let b = x^((Q + 1) / 2) and m = (d - 1) / 2. The product + * + * D = b^Q * b^(Q^3) * ... * b^(Q^(2m - 1)) + * + * satisfies D^2 = x^(Q + Q^2 + ... + Q^(d - 1)) = N / x. Therefore R = x * D / sqrt(N) satisfies R^2 = x. + * Only the seed exponentiation depends on the bit length of Q. Frobenius products use logarithmically many + * compositions and modular products in d, without constructing a quadratic nonresidue in K[X]/(B). + * The Q-power image of X is recovered as b^2 / x; division by x and final multiplication by x are linear passes. + * + * Return false and set output to zero if x is nonsquare. Otherwise output is a canonical reduced root. + * @throws std::invalid_argument if d is zero or even, K has characteristic two, B(0) is zero with d > 1, + * a needed composition has a zero cached-power limit, or the divisor inverse has insufficient precision. + * @pre B is irreducible. For d > 1, degree(B) - 1 inverse coefficients suffice. K provides is_square() and + * the coefficient-field sqrt_known_square operation. Irreducibility is not checked here. + */ + template + requires algebra::FieldValue && + requires(const typename Backend::polynomial_type::value_type &value) { + { value.is_square() } -> std::convertible_to; + } + bool square_root_mod_x_odd_degree(typename Backend::polynomial_type &output, + const polynomial_divisor_context &divisor_context, + polynomial_arithmetic::polynomial_context &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(); + if (degree % 2 == 0 || algebra::fields::field_characteristic() == 2) { + throw std::invalid_argument("direct square roots of X require an odd degree and odd characteristic"); + } + if (degree > 1 && modulus[0] == value_type::zero()) { + throw std::invalid_argument("an irreducible divisor of degree greater than one has nonzero constant term"); + } + + const value_type inverse_leading = divisor_context.reversed_divisor_inverse()[0]; + const value_type norm = -modulus[0] * inverse_leading; + if (!norm.is_square()) { + output.assign(1, value_type::zero()); + return false; + } + const value_type norm_root = algebra::fields::sqrt_known_square(norm); + if (degree == 1) { + output.assign(1, norm_root); + return true; + } + if (arithmetic_context.options().modular_composition_cached_power_limit == 0) { + throw std::invalid_argument("direct square roots of X require a positive composition cached-power limit"); + } + + const boost::multiprecision::cpp_int exponent = (algebra::fields::field_order() + 1) >> 1; + polynomial_type seed, seed_squared; + 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]. This also handles nonmonic B. + const value_type c = seed_squared[0] * modulus[0].inversed(); + polynomial_type x_to_q(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); + + const std::size_t term_count = (degree - 1) / 2; + polynomial_type product, x_to_q_squared; + { + // Reuse this cache for b^Q and X^(Q^2), then release it before allocating the product-chain caches. + const polynomial_composition_precomputation cache(x_to_q, degree, divisor_context, + arithmetic_context); + compose_mod(product, seed, cache, divisor_context, arithmetic_context); + if (term_count > 1) { + compose_mod(x_to_q_squared, x_to_q, cache, divisor_context, arithmetic_context); + } + } + if (term_count > 1) { + 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]; + } + condense(product); + } + scalar_multiplication(output, product, norm_root.inversed()); + return true; + } + /** * Compute a square root of input in K[X]/(B) with Tonelli-Shanks, using the multiplicative-group decomposition and * quadratic nonresidue cached in square_root_context. The irreducible polynomial B is stored in the divisor context diff --git a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp index 88101ba30..e371abd2d 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp @@ -144,11 +144,15 @@ namespace nil::crypto3::math { * P^2 - X * Q^2 = g, * * with degree(P) at most floor(degree(g) / 2) and degree(Q) at most - * floor((degree(g) - 1) / 2). The coefficient generator supplies field elements used to find a quadratic - * nonresidue in K[X]/(g); it remains owned by the caller. + * floor((degree(g) - 1) / 2). Odd degrees use direct Frobenius root recovery when modular composition is enabled; + * linear inputs always use a coefficient-field root. These paths consume no samples from the coefficient generator. + * Even degrees, and higher odd degrees with modular composition disabled, use Tonelli-Shanks. The coefficient + * generator supplies field elements used to find its quadratic nonresidue in K[X]/(g); it remains owned by the + * caller. * - * Irreducibility of g is a caller precondition and is not tested. The coefficient generator must eventually produce - * coefficients forming a nonsquare canonical representative of degree below degree(g). + * Irreducibility of g is a caller precondition and is not tested. When Tonelli-Shanks is needed, the coefficient + * generator must eventually produce coefficients forming a nonsquare canonical representative of degree below + * degree(g). The selected root, and hence the returned representation, need not match the Tonelli-Shanks result. * * K must be a finite field of odd characteristic. Under these preconditions and the composed arithmetic contracts, * recovery is complete: no value means no exact representation exists. This irreducible-only operation does not @@ -189,31 +193,37 @@ namespace nil::crypto3::math { const std::size_t inverse_precision = std::max(1, degree - 1); polynomial_divisor_context divisor_context(g, inverse_precision, arithmetic_context); - // Work with the canonical representative of X in K[X]/(g). Reduction is necessary when g is linear. - const polynomial_type x = {value_type::zero(), value_type::one()}; - polynomial_type x_mod_g; - remainder(x_mod_g, x, divisor_context, arithmetic_context); - if (!is_square_mod(x_mod_g, divisor_context, arithmetic_context)) { - return std::nullopt; - } - - // Adapt caller-generated field coefficients into canonical representatives of K[X]/(g). No random source or - // polynomial backend is constructed here. - auto quotient_representative_generator = [&]() { - polynomial_type representative(degree); - for (std::size_t index = 0; index < degree; ++index) { - representative[index] = coefficient_generator(); + polynomial_type root; + if (degree % 2 != 0 && + (degree == 1 || arithmetic_context.options().modular_composition_cached_power_limit != 0)) { + // Select this path before constructing a Tonelli-Shanks context: it needs no quotient-field nonresidue. + if (!square_root_mod_x_odd_degree(root, divisor_context, arithmetic_context)) { + return std::nullopt; + } + } else { + const polynomial_type x = {value_type::zero(), value_type::one()}; + polynomial_type x_mod_g; + remainder(x_mod_g, x, divisor_context, arithmetic_context); + if (!is_square_mod(x_mod_g, divisor_context, arithmetic_context)) { + return std::nullopt; } - condense(representative); - return representative; - }; - polynomial_square_root_context square_root_context(divisor_context, arithmetic_context, - quotient_representative_generator); - // Recover R with R^2 = X mod g. A failure here contradicts the successful square test above. - polynomial_type root; - if (!square_root_mod(root, x_mod_g, square_root_context, arithmetic_context)) { - throw std::logic_error("polynomial X-norm recovery failed after X was reported square modulo g"); + // Adapt caller-generated field coefficients into canonical representatives of K[X]/(g). No random + // source or polynomial backend is constructed here. + auto quotient_representative_generator = [&]() { + polynomial_type representative(degree); + for (std::size_t index = 0; index < degree; ++index) { + representative[index] = coefficient_generator(); + } + condense(representative); + return representative; + }; + polynomial_square_root_context square_root_context(divisor_context, arithmetic_context, + quotient_representative_generator); + // A failure here contradicts the successful square test above. + if (!square_root_mod(root, x_mod_g, square_root_context, arithmetic_context)) { + throw std::logic_error("polynomial X-norm recovery failed after X was reported square modulo g"); + } } // Reconstruct P = R * Q mod g with the unique degree bounds required by the norm equation. diff --git a/libs/math/test/polynomial_frobenius.cpp b/libs/math/test/polynomial_frobenius.cpp index 87d1877aa..4f90721a0 100644 --- a/libs/math/test/polynomial_frobenius.cpp +++ b/libs/math/test/polynomial_frobenius.cpp @@ -25,17 +25,20 @@ #define BOOST_TEST_MODULE polynomial_frobenius_test #include +#include #include #include #include #include +#include #include #include #include #include +#include #include namespace { @@ -55,6 +58,21 @@ namespace { } return value; } + + struct counting_product_backend : polynomial_arithmetic::schoolbook_backend { + using base_type = polynomial_arithmetic::schoolbook_backend; + explicit counting_product_backend(std::size_t &operations) : operations(operations) { + } + void square(polynomial_type &output, const polynomial_type &input) const { + ++operations; + base_type::square(output, input); + } + void multiply(polynomial_type &output, const polynomial_type &left, const polynomial_type &right) const { + ++operations; + base_type::multiply(output, left, right); + } + std::size_t &operations; + }; } // namespace BOOST_AUTO_TEST_SUITE(polynomial_frobenius_test_suite) @@ -228,4 +246,126 @@ BOOST_AUTO_TEST_CASE(frobenius_context_rejects_invalid_precomputation_inputs) { std::invalid_argument); } +BOOST_AUTO_TEST_CASE(frobenius_products_match_geometric_exponents_and_support_aliasing) { + using value = fields::babybear::value_type; + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + polynomial_arithmetic::polynomial_context context; + const polynomial modulus {value(5), value(7), value(3), value(2)}; + const math::polynomial_divisor_context divisor(modulus, 4, context); + const polynomial input {value(2), value(5), value(7), value(11), value(13), value(17)}; + const boost::multiprecision::cpp_int q = fields::field_order(); + for (const boost::multiprecision::cpp_int &step : + {boost::multiprecision::cpp_int(q), boost::multiprecision::cpp_int(q * q)}) { + polynomial image; + math::powmod_x(image, step, divisor, context); + for (const std::size_t count : {0, 1, 2, 3, 4, 7, 8, 9, 15, 16, 17, 31, 32, 33}) { + boost::multiprecision::cpp_int exponent = 0, power = 1; + for (std::size_t i = 0; i < count; ++i) { + exponent += power; + power *= step; + } + polynomial expected, actual; + math::powmod(expected, input, exponent, divisor, context); + math::frobenius_product(actual, input, count, image, divisor, context); + BOOST_CHECK(actual == expected); + auto alias = input; + math::frobenius_product(alias, alias, count, image, divisor, context); + BOOST_CHECK(alias == expected); + alias = image; + math::frobenius_product(alias, input, count, alias, divisor, context); + BOOST_CHECK(alias == expected); + } + } +} + +BOOST_AUTO_TEST_CASE(frobenius_products_use_bounded_work_for_maximum_term_count) { + using backend = counting_product_backend; + using polynomial = backend::polynomial_type; + using value = polynomial::value_type; + std::size_t operations = 0; + polynomial_arithmetic::polynomial_context context {backend(operations)}; + const polynomial modulus {-value::one(), value::zero(), value::one()}; + const math::polynomial_divisor_context divisor(modulus, 1, context); + const polynomial x {value::zero(), value::one()}; + operations = 0; + polynomial actual; + // Q is odd and X^2 = 1, so F fixes X and the product is X for every odd count. + math::frobenius_product(actual, x, std::numeric_limits::max(), x, divisor, context); + BOOST_CHECK(actual == x); + BOOST_CHECK_LT(operations, 8 * std::numeric_limits::digits); +} + +BOOST_AUTO_TEST_CASE(frobenius_products_over_fp12_reuse_both_q_and_q_squared_images) { + using backend = polynomial_arithmetic::bn254_small_prime_fft_backend; + using reference_backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + using value = fq12_value_type; + value nonsquare = fq12_value(1); + while (nonsquare.is_square()) { + nonsquare += value::one(); + } + const value scale = fq12_value(31); + const polynomial modulus {-scale * nonsquare, value::zero(), scale}; + const polynomial input {fq12_value(5), fq12_value(17)}; + const polynomial x {value::zero(), value::one()}, minus_x {value::zero(), -value::one()}; + polynomial_arithmetic::polynomial_context context {backend(8)}; + polynomial_arithmetic::polynomial_context reference; + const math::polynomial_divisor_context divisor(modulus, 1, context); + const math::polynomial_divisor_context reference_divisor(modulus, 1, reference); + const value norm = input[0].squared() - nonsquare * input[1].squared(); + // X^2 = nonsquare: Q-power conjugation negates X but fixes all twelve coefficient coordinates. + for (const std::size_t count : {0, 1, 2, 3, 7, 8, 9}) { + polynomial expected, actual; + math::scalar_multiplication(expected, count % 2 == 0 ? polynomial {value::one()} : input, norm.pow(count / 2)); + math::frobenius_product(actual, input, count, minus_x, divisor, context); + BOOST_CHECK(actual == expected); + // The Q^2-power map fixes every element in this quadratic extension. + math::powmod(expected, input, count, reference_divisor, reference); + math::frobenius_product(actual, input, count, x, divisor, context); + BOOST_CHECK(actual == expected); + } +} + +BOOST_AUTO_TEST_CASE(frobenius_products_handle_constants_zero_rings_and_invalid_inputs) { + using value = fields::babybear::value_type; + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 0; + polynomial_arithmetic::polynomial_context context(backend {}, options); + const polynomial modulus {value(5), value(7), value(3), value(2)}; + const math::polynomial_divisor_context divisor(modulus, 4, context); + const polynomial x {value::zero(), value::one()}; + polynomial image; + math::powmod_x(image, fields::field_order(), divisor, context); + for (const auto &constant : {value::zero(), value::one(), value(7)}) { + for (const std::size_t count : {0, 1, 2, 17}) { + polynomial result; + math::frobenius_product(result, polynomial {constant}, count, image, divisor, context); + BOOST_CHECK(result == polynomial({constant.pow(count)})); + } + } + polynomial output {value(19)}; + math::frobenius_product(output, x, 1, image, divisor, context); + BOOST_CHECK(output == x); + BOOST_CHECK_THROW(math::frobenius_product(output, x, 2, image, divisor, context), std::invalid_argument); + BOOST_CHECK(output == x); + polynomial empty; + empty.get_storage().clear(); + polynomial padded(2); + padded[0] = value::one(); + padded[1] = value::zero(); + BOOST_CHECK_THROW(math::frobenius_product(output, empty, 1, image, divisor, context), std::invalid_argument); + BOOST_CHECK_THROW(math::frobenius_product(output, padded, 1, image, divisor, context), std::invalid_argument); + for (const auto &bad_image : {empty, padded, modulus}) { + BOOST_CHECK_THROW(math::frobenius_product(output, x, 1, bad_image, divisor, context), std::invalid_argument); + } + const math::polynomial_divisor_context constant_divisor(polynomial {value(7)}, 1, context); + for (const std::size_t count : {0, 1, 9}) { + math::frobenius_product(output, x, count, polynomial {value::zero()}, constant_divisor, context); + BOOST_CHECK(output == polynomial({value::zero()})); + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/libs/math/test/polynomial_square_root.cpp b/libs/math/test/polynomial_square_root.cpp index 46b75752c..19f96c245 100644 --- a/libs/math/test/polynomial_square_root.cpp +++ b/libs/math/test/polynomial_square_root.cpp @@ -37,6 +37,7 @@ #include #include +#include #include namespace { @@ -344,4 +345,139 @@ BOOST_AUTO_TEST_CASE(tonelli_shanks_recovers_a_square_in_a_quadratic_extension_o BOOST_CHECK(recovered_square == square); } +BOOST_AUTO_TEST_CASE(direct_odd_degree_roots_match_tonelli_shanks_for_sparse_and_dense_divisors) { + polynomial_arithmetic::polynomial_context context; + const polynomial_type x {value_type::zero(), value_type::one()}; + const polynomial_type nonresidue {first_quadratic_non_residue()}; + std::size_t squares = 0, nonsquares = 0; + for (const std::size_t degree : {3, 5, 9, 27}) { + const std::size_t prime_factor = degree == 5 ? 5 : 3; + const auto character_exponent = (field_type::modulus - 1u) / prime_factor; + value_type c(2); + while (c.pow(character_exponent) == value_type::one()) { + c += value_type::one(); + } + c = c.squared(); + // The binomial irreducibility criterion applies: all prime factors of degree divide ord(c), and none + // divide (|K| - 1) / ord(c). Translation and nonzero scaling preserve irreducibility. + for (const std::size_t shift : {0, 1, 2, 3}) { + polynomial_type monic {value_type::one()}; + const polynomial_type linear {value_type(shift), value_type::one()}; + for (std::size_t i = 0; i < degree; ++i) { + context.multiply(monic, monic, linear); + } + monic[0] -= c; + const bool expected_square = (-monic[0]).is_square(); + expected_square ? ++squares : ++nonsquares; + for (const auto &scale : {value_type::one(), value_type(7)}) { + polynomial_type modulus; + math::scalar_multiplication(modulus, monic, scale); + const math::polynomial_divisor_context divisor(modulus, degree - 1, context); + const math::polynomial_square_root_context tonelli(nonresidue, divisor, context); + polynomial_type root, reference; + BOOST_REQUIRE_EQUAL(math::square_root_mod_x_odd_degree(root, divisor, context), expected_square); + BOOST_REQUIRE_EQUAL(math::square_root_mod(reference, x, tonelli, context), expected_square); + if (expected_square) { + polynomial_type squared, negative_reference; + math::squaremod(squared, root, divisor, context); + math::scalar_multiplication(negative_reference, reference, -value_type::one()); + BOOST_CHECK(squared == x); + BOOST_CHECK(root == reference || root == negative_reference); + BOOST_CHECK_LE(root.size(), degree); + BOOST_CHECK(root.size() == 1 || root.back() != value_type::zero()); + } else { + BOOST_CHECK(root == polynomial_type({value_type::zero()})); + } + } + } + } + BOOST_CHECK_GT(squares, 0); + BOOST_CHECK_GT(nonsquares, 0); +} + +BOOST_AUTO_TEST_CASE(direct_linear_roots_handle_zero_nonmonic_divisors_and_disabled_composition) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 0; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + for (const auto &x : {value_type::zero(), value_type::one(), value_type(49), first_quadratic_non_residue()}) { + for (const auto &scale : {value_type::one(), value_type(11)}) { + const polynomial_type modulus {-scale * x, scale}; + const math::polynomial_divisor_context divisor(modulus, 1, context); + polynomial_type root {value_type(19)}; + BOOST_REQUIRE_EQUAL(math::square_root_mod_x_odd_degree(root, divisor, context), x.is_square()); + BOOST_REQUIRE_EQUAL(root.size(), 1); + BOOST_CHECK(root[0].squared() == (x.is_square() ? x : value_type::zero())); + } + } +} + +BOOST_AUTO_TEST_CASE(direct_odd_degree_roots_support_dense_bn254_fp12_coefficients_and_small_prime_products) { + using extension_field = fields::fp12_2over3over2>; + using extension_value = extension_field::value_type; + using backend = polynomial_arithmetic::bn254_small_prime_fft_backend; + using reference_backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + boost::random::mt19937 rng(0x0DD012); + const boost::multiprecision::cpp_int noncube_exponent = (fields::field_order() - 1) / 3; + extension_value c = extension_value::one(); + for (std::size_t i = 0; i < 64 && c.pow(noncube_exponent) == extension_value::one(); ++i) { + c = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE(!c.is_zero()); + BOOST_REQUIRE(c.pow(noncube_exponent) != extension_value::one()); + c = c.squared(); + polynomial_arithmetic::polynomial_context reference; + const polynomial x {extension_value::zero(), extension_value::one()}; + for (const std::size_t degree : {3, 9}) { + extension_value shift = extension_value::one(); + for (std::size_t i = 0; i < 64 && !(c - shift.pow(degree)).is_square(); ++i) { + shift = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE((c - shift.pow(degree)).is_square()); + polynomial modulus {extension_value::one()}; + const polynomial linear {shift, extension_value::one()}; + for (std::size_t i = 0; i < degree; ++i) { + reference.multiply(modulus, modulus, linear); + } + modulus[0] -= c; + const extension_value scale = nil::crypto3::algebra::random_element(rng); + BOOST_REQUIRE(!scale.is_zero()); + math::scalar_multiplication(modulus, modulus, scale); + polynomial_arithmetic::polynomial_context_options options; + // Also exercise the smallest valid composition cache and inverse-based polynomial reduction. + options.modular_composition_cached_power_limit = degree == 3 ? 1 : 3; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context context(backend(32), options); + const math::polynomial_divisor_context divisor(modulus, degree - 1, context); + const math::polynomial_divisor_context reference_divisor(modulus, degree - 1, reference); + polynomial root, squared; + BOOST_REQUIRE(math::square_root_mod_x_odd_degree(root, divisor, context)); + math::squaremod(squared, root, reference_divisor, reference); + BOOST_CHECK(squared == x); + BOOST_CHECK_LE(root.size(), degree); + } +} + +BOOST_AUTO_TEST_CASE(direct_odd_degree_roots_reject_unsupported_contexts) { + polynomial_arithmetic::polynomial_context context; + polynomial_type output {value_type(19)}; + for (const auto &modulus : + {polynomial_type {value_type::one()}, + polynomial_type {-first_quadratic_non_residue(), value_type::zero(), value_type::one()}, + polynomial_type {value_type::zero(), value_type::one(), value_type::zero(), value_type::one()}}) { + const math::polynomial_divisor_context divisor(modulus, 2, context); + BOOST_CHECK_THROW(math::square_root_mod_x_odd_degree(output, divisor, context), std::invalid_argument); + BOOST_CHECK(output == polynomial_type({value_type(19)})); + } + const value_type c = first_cubic_non_residue().squared(); + const polynomial_type cubic {-c, value_type::zero(), value_type::zero(), value_type::one()}; + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 0; + polynomial_arithmetic::polynomial_context disabled(backend_type {}, options); + const math::polynomial_divisor_context divisor(cubic, 2, disabled); + BOOST_CHECK_THROW(math::square_root_mod_x_odd_degree(output, divisor, disabled), std::invalid_argument); + BOOST_CHECK(output == polynomial_type({value_type(19)})); +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/libs/math/test/polynomial_x_norm_reconstruction.cpp b/libs/math/test/polynomial_x_norm_reconstruction.cpp index 7c0ff74b6..7b59c8f23 100644 --- a/libs/math/test/polynomial_x_norm_reconstruction.cpp +++ b/libs/math/test/polynomial_x_norm_reconstruction.cpp @@ -36,6 +36,7 @@ #include #include +#include #include #include #include @@ -183,11 +184,16 @@ BOOST_AUTO_TEST_CASE(recovers_and_normalizes_an_irreducible_quadratic_with_the_s BOOST_REQUIRE(g == polynomial_type({value_type::one(), value_type::zero() - value_type(79), value_type::one()})); boost::random::mt19937 rng(0x584E4F52); - auto coefficient_generator = [&] { return nil::crypto3::algebra::random_element(rng); }; + std::size_t generator_calls = 0; + auto coefficient_generator = [&] { + ++generator_calls; + return nil::crypto3::algebra::random_element(rng); + }; const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); BOOST_REQUIRE(result.has_value()); + BOOST_CHECK_GT(generator_calls, 0); BOOST_CHECK_LE(result->p.size() - 1, 1); BOOST_CHECK_LE(result->q.size() - 1, 0); // Rational reconstruction makes Q monic, initially producing lambda = 1/81. Exact normalization restores Q^2 = 81. @@ -229,6 +235,116 @@ BOOST_AUTO_TEST_CASE(returns_no_value_when_the_scalar_multiple_cannot_be_normali BOOST_CHECK(!result.has_value()); } +BOOST_AUTO_TEST_CASE(odd_degree_norm_recovery_verifies_exact_identities_without_sampling) { + polynomial_arithmetic::polynomial_context context; + value_type c(2); + const auto noncube_exponent = (field_type::modulus - 1u) / 3u; + while (c.pow(noncube_exponent) == value_type::one()) { + c += value_type::one(); + } + c = c.squared(); + std::size_t generator_calls = 0; + auto generator = [&]() -> value_type { + ++generator_calls; + throw std::logic_error("odd-degree recovery must not sample a polynomial nonresidue"); + }; + for (const std::size_t degree : {3, 9, 27}) { + // X^degree - c is irreducible by the binomial criterion; translation preserves irreducibility. + value_type shift = value_type::one(); + while (!(c - shift.pow(degree)).is_square()) { + shift += value_type::one(); + } + polynomial_type monic {value_type::one()}; + const polynomial_type linear {shift, value_type::one()}; + for (std::size_t i = 0; i < degree; ++i) { + context.multiply(monic, monic, linear); + } + monic[0] -= c; + for (const auto &scale : {value_type::one(), value_type(49), first_quadratic_non_residue()}) { + polynomial_type g; + math::scalar_multiplication(g, monic, scale); + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, context, generator); + // -1 is square in BabyBear. The leading coefficient of an odd-degree norm is minus a square. + BOOST_REQUIRE_EQUAL(result.has_value(), scale.is_square()); + if (result) { + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, context) == g); + BOOST_CHECK_LE(result->p.size() - 1, degree / 2); + BOOST_CHECK_LE(result->q.size() - 1, (degree - 1) / 2); + } + } + } + const polynomial_type x {value_type::zero(), value_type::one()}; + const auto zero_root = math::recover_irreducible_polynomial_x_norm_representation(x, context, generator); + BOOST_REQUIRE(zero_root.has_value()); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*zero_root, context) == x); + const polynomial_type nonsquare_linear {-first_quadratic_non_residue(), value_type::one()}; + BOOST_CHECK(!math::recover_irreducible_polynomial_x_norm_representation(nonsquare_linear, context, generator)); + BOOST_CHECK_EQUAL(generator_calls, 0); +} + +BOOST_AUTO_TEST_CASE(odd_degree_recovery_preserves_tonelli_shanks_when_composition_is_disabled) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 0; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + value_type c(2); + const auto noncube_exponent = (field_type::modulus - 1u) / 3u; + while (c.pow(noncube_exponent) == value_type::one()) { + c += value_type::one(); + } + const polynomial_type g {-c.squared(), value_type::zero(), value_type::zero(), value_type::one()}; + const auto nonresidue = first_quadratic_non_residue(); + std::size_t generator_calls = 0; + auto generator = [&] { + // A coefficient-field nonsquare remains nonsquare in an odd-degree extension. + return generator_calls++ % 3 == 0 ? nonresidue : value_type::zero(); + }; + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, context, generator); + BOOST_REQUIRE(result.has_value()); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, context) == g); + BOOST_CHECK_EQUAL(generator_calls, 3); + + const polynomial_type linear {value_type(9), -value_type(4)}; + const auto linear_result = math::recover_irreducible_polynomial_x_norm_representation(linear, context, generator); + BOOST_REQUIRE(linear_result.has_value()); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*linear_result, context) == linear); + BOOST_CHECK_EQUAL(generator_calls, 3); +} + +BOOST_AUTO_TEST_CASE(odd_degree_bn254_fp12_norm_recovery_agrees_across_backends_without_sampling) { + using backend = polynomial_arithmetic::bn254_small_prime_fft_backend; + using polynomial = backend::polynomial_type; + boost::random::mt19937 rng(0xF0120DD); + const boost::multiprecision::cpp_int noncube_exponent = (fields::field_order() - 1) / 3; + fq12_value_type c = fq12_value_type::one(); + for (std::size_t i = 0; i < 64 && c.pow(noncube_exponent) == fq12_value_type::one(); ++i) { + c = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE(!c.is_zero()); + BOOST_REQUIRE(c.pow(noncube_exponent) != fq12_value_type::one()); + c = c.squared(); + fq12_value_type shift = fq12_value_type::one(); + for (std::size_t i = 0; i < 64 && !(c - shift.pow(3)).is_square(); ++i) { + shift = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE((c - shift.pow(3)).is_square()); + polynomial g {shift.pow(3) - c, fq12_scalar(3) * shift.squared(), fq12_scalar(3) * shift, fq12_value_type::one()}; + const auto scale = nil::crypto3::algebra::random_element(rng).squared(); + BOOST_REQUIRE(!scale.is_zero()); + math::scalar_multiplication(g, g, scale); + auto generator = []() -> fq12_value_type { + throw std::logic_error("odd-degree recovery must not sample a polynomial nonresidue"); + }; + polynomial_arithmetic::polynomial_context context {backend(16)}; + polynomial_arithmetic::polynomial_context reference; + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, context, generator); + const auto expected = math::recover_irreducible_polynomial_x_norm_representation(g, reference, generator); + BOOST_REQUIRE(result.has_value()); + BOOST_REQUIRE(expected.has_value()); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, reference) == g); + BOOST_CHECK(result->p == expected->p); + BOOST_CHECK(result->q == expected->q); +} + BOOST_AUTO_TEST_CASE(recovers_a_linear_bn254_fq12_norm_with_the_schoolbook_backend) { polynomial_arithmetic::polynomial_context arithmetic_context; check_bn254_fq12_linear_recovery(arithmetic_context, 0xF0125001); From 50fdee1d18861d51334670b310e9cb111f0b398a Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Fri, 23 Oct 2026 19:02:42 -0400 Subject: [PATCH 2/2] Use characteristic-sized seeds for odd-degree square roots --- .../quotient_ring/polynomial_frobenius.hpp | 152 +++++++++++------- .../quotient_ring/polynomial_square_root.hpp | 81 ++++++++-- libs/math/test/polynomial_frobenius.cpp | 38 +++++ libs/math/test/polynomial_square_root.cpp | 42 +++++ 4 files changed, 239 insertions(+), 74 deletions(-) diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp index 0f016156d..ad3220f6d 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp @@ -149,6 +149,97 @@ namespace nil::crypto3::math { output = std::move(result); } + namespace detail { + + // Shared chain for coefficient-fixed Frobenius products and the characteristic-power square-root seed. + // In the characteristic case the supplied image is X^p, and coefficients must support Frobenius_map. + template + void frobenius_product_impl(typename Backend::polynomial_type &output, + const typename Backend::polynomial_type &input, std::size_t term_count, + const typename Backend::polynomial_type &x_to_frobenius_power, + const polynomial_divisor_context &divisor_context, + polynomial_arithmetic::polynomial_context &arithmetic_context) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename polynomial_type::value_type; + using cache_type = polynomial_composition_precomputation; + 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(x_to_frobenius_power) || + x_to_frobenius_power.size() > std::max(1, degree)) { + throw std::invalid_argument( + "Frobenius products require canonical inputs and a reduced Frobenius image"); + } + if (degree == 0 || term_count == 0) { + output.assign(1, degree == 0 ? value_type::zero() : value_type::one()); + return; + } + polynomial_type base; + remainder(base, input, divisor_context, arithmetic_context); + if (term_count == 1) { + output = std::move(base); + return; + } + if constexpr (!CharacteristicFrobenius) { + // Constants are fixed by Q-Frobenius, but generally not by p-Frobenius in an extension field. + if (base.size() == 1) { + powmod(output, base, term_count, divisor_context, arithmetic_context); + return; + } + } + + const auto apply = [&](polynomial_type &result, const polynomial_type &operand, std::size_t power, + const cache_type &cache) { + if constexpr (CharacteristicFrobenius) { + // Absolute Frobenius is semilinear: f^(p^power) = sigma^power(f)(X^(p^power)). + // Copy before twisting so result may alias operand, including the cached-image input. + polynomial_type twisted = operand; + for (auto &coefficient : twisted) { + coefficient = coefficient.Frobenius_map(power); + } + compose_mod(result, twisted, cache, divisor_context, arithmetic_context); + } else { + compose_mod(result, operand, cache, divisor_context, arithmetic_context); + } + }; + + 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_image) -> std::pair { + if (count == 1) { + return {base, x_to_frobenius_power}; + } + auto [product, image] = self(self, count / 2, true); + std::optional doubling_cache; + if (count / 2 > 1) { + doubling_cache.emplace(image, degree, divisor_context, arithmetic_context); + } + const auto &cache = doubling_cache ? *doubling_cache : first_cache; + // With P_k = product and H_k = image, P_2k = P_k * F^k(P_k) and H_2k = F^k(H_k). + polynomial_type shifted; + apply(shifted, product, count / 2, cache); + mulmod(product, product, shifted, divisor_context, arithmetic_context); + if (need_image) { + apply(image, image, count / 2, cache); + } + if (count % 2 != 0) { + // P_(2k+1) = base * F(P_2k); use the same initial cache at every odd step. + apply(shifted, product, 1, first_cache); + mulmod(product, base, shifted, divisor_context, arithmetic_context); + if (need_image) { + apply(image, image, 1, first_cache); + } + } + return {std::move(product), std::move(image)}; + }; + // The final image is not used by the caller, so avoid its compositions. + auto result = chain(chain, term_count, false); + output = std::move(result.first); + } + + } // namespace detail + /** * Compute input * F(input) * ... * F^(term_count - 1)(input) modulo B, where F is a power of the coefficient * field's Frobenius map. The supplied x_to_frobenius_power is the canonical reduced image F(X); it can represent @@ -171,65 +262,8 @@ namespace nil::crypto3::math { std::size_t term_count, const typename Backend::polynomial_type &x_to_frobenius_power, const polynomial_divisor_context &divisor_context, polynomial_arithmetic::polynomial_context &arithmetic_context) { - using polynomial_type = typename Backend::polynomial_type; - using value_type = typename polynomial_type::value_type; - using cache_type = polynomial_composition_precomputation; - 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(x_to_frobenius_power) || - x_to_frobenius_power.size() > std::max(1, degree)) { - throw std::invalid_argument("Frobenius products require canonical inputs and a reduced Frobenius image"); - } - if (degree == 0 || term_count == 0) { - output.assign(1, degree == 0 ? value_type::zero() : value_type::one()); - return; - } - polynomial_type base; - remainder(base, input, divisor_context, arithmetic_context); - if (term_count == 1) { - output = std::move(base); - return; - } - // Every coefficient-field constant is fixed by F, including zero. - if (base.size() == 1) { - powmod(output, base, term_count, divisor_context, arithmetic_context); - return; - } - - 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_image) -> std::pair { - if (count == 1) { - return {base, x_to_frobenius_power}; - } - auto [product, image] = self(self, count / 2, true); - std::optional doubling_cache; - if (count / 2 > 1) { - doubling_cache.emplace(image, degree, divisor_context, arithmetic_context); - } - const auto &cache = doubling_cache ? *doubling_cache : first_cache; - // With P_k = product and H_k = image, P_2k = P_k * P_k(H_k) and H_2k = H_k(H_k). - polynomial_type shifted; - compose_mod(shifted, product, cache, divisor_context, arithmetic_context); - mulmod(product, product, shifted, divisor_context, arithmetic_context); - if (need_image) { - compose_mod(image, image, cache, divisor_context, arithmetic_context); - } - if (count % 2 != 0) { - // P_(2k+1) = base * F(P_2k); use the same initial cache at every odd step. - compose_mod(shifted, product, first_cache, divisor_context, arithmetic_context); - mulmod(product, base, shifted, divisor_context, arithmetic_context); - if (need_image) { - compose_mod(image, image, first_cache, divisor_context, arithmetic_context); - } - } - return {std::move(product), std::move(image)}; - }; - // The final image is not used by the caller, so avoid its compositions. - auto result = chain(chain, term_count, false); - output = std::move(result.first); + detail::frobenius_product_impl(output, input, term_count, x_to_frobenius_power, divisor_context, + arithmetic_context); } } // namespace nil::crypto3::math diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp index 89c361c88..1daf171c1 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp @@ -251,6 +251,67 @@ namespace nil::crypto3::math { return detail::is_square_mod_impl(input, divisor_context, arithmetic_context, nullptr); } + 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. + template + std::pair + square_root_x_seed(const polynomial_divisor_context &divisor_context, + polynomial_arithmetic::polynomial_context &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; + }) { + 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. + // 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() - 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); + frobenius_product_impl(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() + 1) >> 1; + polynomial_type seed_squared; + 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); + } + return {std::move(seed), std::move(x_to_q)}; + } + + } // namespace detail + /** * Compute a square root of X modulo an irreducible polynomial B of odd degree d over an odd-order field K. * Write Q = |K|, x = X mod B, and N = Norm(x) = -B(0) / leading_coefficient(B). The element x is square exactly @@ -261,9 +322,10 @@ namespace nil::crypto3::math { * D = b^Q * b^(Q^3) * ... * b^(Q^(2m - 1)) * * satisfies D^2 = x^(Q + Q^2 + ... + Q^(d - 1)) = N / x. Therefore R = x * D / sqrt(N) satisfies R^2 = x. - * Only the seed exponentiation depends on the bit length of Q. Frobenius products use logarithmically many - * compositions and modular products in d, without constructing a quadratic nonresidue in K[X]/(B). - * The Q-power image of X is recovered as b^2 / x; division by x and final multiplication by x are linear passes. + * For K = F_(p^e) with coefficient Frobenius support, the seed uses one characteristic-sized exponentiation + * and O(log(e)) semilinear compositions, producing X^Q at the same time. Other coefficient types retain the + * direct Q-sized seed. The remaining Frobenius product uses O(log(d)) compositions and modular products, + * without constructing a quadratic nonresidue in K[X]/(B). Final multiplication by x is a linear pass. * * Return false and set output to zero if x is nonsquare. Otherwise output is a canonical reduced root. * @throws std::invalid_argument if d is zero or even, K has characteristic two, B(0) is zero with d > 1, @@ -306,18 +368,7 @@ namespace nil::crypto3::math { throw std::invalid_argument("direct square roots of X require a positive composition cached-power limit"); } - const boost::multiprecision::cpp_int exponent = (algebra::fields::field_order() + 1) >> 1; - polynomial_type seed, seed_squared; - 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]. This also handles nonmonic B. - const value_type c = seed_squared[0] * modulus[0].inversed(); - polynomial_type x_to_q(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); + auto [seed, x_to_q] = detail::square_root_x_seed(divisor_context, arithmetic_context); const std::size_t term_count = (degree - 1) / 2; polynomial_type product, x_to_q_squared; diff --git a/libs/math/test/polynomial_frobenius.cpp b/libs/math/test/polynomial_frobenius.cpp index 4f90721a0..1677d04ad 100644 --- a/libs/math/test/polynomial_frobenius.cpp +++ b/libs/math/test/polynomial_frobenius.cpp @@ -368,4 +368,42 @@ BOOST_AUTO_TEST_CASE(frobenius_products_handle_constants_zero_rings_and_invalid_ } } +BOOST_AUTO_TEST_CASE(characteristic_product_twists_coefficients_and_matches_direct_powers) { + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + const polynomial modulus {fq12_value(3), fq12_value(19), fq12_value(37), fq12_value(53)}; + const polynomial input {fq12_value(71), fq12_value(89), fq12_value(107)}; + BOOST_REQUIRE(input[0].Frobenius_map(1) != input[0]); + const boost::multiprecision::cpp_int p = fields::field_characteristic(); + for (const std::size_t cache_limit : {1, 3}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = cache_limit; + polynomial_arithmetic::polynomial_context context(backend {}, options); + const math::polynomial_divisor_context divisor(modulus, 2, context); + polynomial image; + math::powmod_x(image, p, divisor, context); + for (const std::size_t count : {0, 1, 2, 3, 6, 12, 13}) { + boost::multiprecision::cpp_int exponent = 0; + for (std::size_t i = 0; i < count; ++i) { + exponent = exponent * p + 1; + } + // Constants also need the coefficient twist: they are not fixed by absolute Frobenius. + for (const auto &base : {input, polynomial {input[0]}, polynomial {fq12_value_type::zero()}}) { + polynomial expected, result; + math::powmod(expected, base, exponent, divisor, context); + math::detail::frobenius_product_impl(result, base, count, image, divisor, context); + BOOST_CHECK(result == expected); + if (count == 12) { + result = base; + math::detail::frobenius_product_impl(result, result, count, image, divisor, context); + BOOST_CHECK(result == expected); + result = image; + math::detail::frobenius_product_impl(result, base, count, result, divisor, context); + BOOST_CHECK(result == expected); + } + } + } + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/libs/math/test/polynomial_square_root.cpp b/libs/math/test/polynomial_square_root.cpp index 19f96c245..05d762cdd 100644 --- a/libs/math/test/polynomial_square_root.cpp +++ b/libs/math/test/polynomial_square_root.cpp @@ -480,4 +480,46 @@ BOOST_AUTO_TEST_CASE(direct_odd_degree_roots_reject_unsupported_contexts) { BOOST_CHECK(output == polynomial_type({value_type(19)})); } +BOOST_AUTO_TEST_CASE(characteristic_seed_and_image_match_independent_full_powers) { + using extension_field = fields::fp12_2over3over2>; + using extension_value = extension_field::value_type; + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + boost::random::mt19937 rng(0x5EED012); + const boost::multiprecision::cpp_int q = fields::field_order(); + const boost::multiprecision::cpp_int exponent = (q + 1) >> 1; + polynomial_arithmetic::polynomial_context context; + for (const std::size_t degree : {3, 5, 9}) { + polynomial modulus(degree + 1, extension_value::zero()); + for (auto &coefficient : modulus) { + coefficient = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE(!modulus[0].is_zero()); + BOOST_REQUIRE(!modulus.back().is_zero()); + // This identity holds even in reducible quotient rings; full coordinates also exercise semilinearity. + BOOST_REQUIRE(modulus[1].Frobenius_map(1) != modulus[1]); + const math::polynomial_divisor_context divisor(modulus, degree - 1, context); + const auto [seed, image] = math::detail::square_root_x_seed(divisor, context); + polynomial expected_seed, expected_image; + math::powmod_x(expected_seed, exponent, divisor, context); + math::powmod_x(expected_image, q, divisor, context); + BOOST_CHECK(seed == expected_seed); + BOOST_CHECK(image == expected_image); + } +} + +BOOST_AUTO_TEST_CASE(prime_field_seed_retains_direct_exponentiation) { + polynomial_arithmetic::polynomial_context context; + const polynomial_type modulus {value_type(5), value_type(7), value_type(11), value_type(13)}; + const math::polynomial_divisor_context divisor(modulus, 2, context); + const auto [seed, image] = math::detail::square_root_x_seed(divisor, context); + const boost::multiprecision::cpp_int q = fields::field_order(); + const boost::multiprecision::cpp_int exponent = (q + 1) >> 1; + polynomial_type expected_seed, expected_image; + math::powmod_x(expected_seed, exponent, divisor, context); + math::powmod_x(expected_image, q, divisor, context); + BOOST_CHECK(seed == expected_seed); + BOOST_CHECK(image == expected_image); +} + BOOST_AUTO_TEST_SUITE_END()