From e762534ce8f2669114b932948b57d14f526281f4 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Sun, 25 Oct 2026 18:18:27 -0400 Subject: [PATCH] Accelerate even-degree polynomial square-root recovery --- .../quotient_ring/polynomial_frobenius.hpp | 102 +++++++++ .../quotient_ring/polynomial_square_root.hpp | 144 ++++++++---- .../polynomial_x_norm_reconstruction.hpp | 44 ++-- libs/math/test/polynomial_frobenius.cpp | 157 +++++++++++++ libs/math/test/polynomial_square_root.cpp | 215 +++++++++++++++++- .../test/polynomial_x_norm_reconstruction.cpp | 119 +++++++++- 6 files changed, 720 insertions(+), 61 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 ad3220f6d..90d0bce30 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 @@ -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 + 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 &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(weight) || !canonical(x_to_frobenius_power) || + x_to_frobenius_power.size() > std::max(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 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 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 1daf171c1..aaf066db1 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 @@ -253,7 +253,24 @@ 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 + void square_root_multiply_by_x(typename Backend::polynomial_type &value, + const polynomial_divisor_context &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 std::pair square_root_x_seed(const polynomial_divisor_context &divisor_context, @@ -261,52 +278,28 @@ namespace nil::crypto3::math { 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. + // 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() - 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(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; + const boost::multiprecision::cpp_int exponent = (algebra::fields::field_order() - 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)}; } @@ -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; @@ -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 + requires algebra::FieldValue && + requires(const typename Backend::polynomial_type::value_type &value, Generator &generator) { + { value.is_square() } -> std::convertible_to; + { generator() } -> std::convertible_to; + } + bool square_root_mod_x_even_degree(typename Backend::polynomial_type &output, + const polynomial_divisor_context &divisor_context, + polynomial_arithmetic::polynomial_context &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() == 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; } 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 e371abd2d..9ece608c2 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 @@ -146,13 +146,14 @@ namespace nil::crypto3::math { * with degree(P) at most floor(degree(g) / 2) and degree(Q) at most * 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. + * Even degrees use a weighted Frobenius trace when modular composition is enabled, sampling only if the initial + * trace is zero. With composition disabled, higher degrees use Tonelli-Shanks. The coefficient generator supplies + * field elements for trace retries or the Tonelli-Shanks nonresidue search and remains owned by the caller. * - * 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. + * Irreducibility of g is a caller precondition and is not tested. The coefficient generator must eventually + * produce a representative with nonzero weighted trace when a trace retry is needed, or a nonsquare when + * Tonelli-Shanks is needed. Each sample uses degree(g) coefficients. 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 @@ -201,13 +202,6 @@ namespace nil::crypto3::math { 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; - } - // 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 = [&]() { @@ -218,11 +212,25 @@ namespace nil::crypto3::math { 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"); + if (arithmetic_context.options().modular_composition_cached_power_limit != 0) { + // Even degrees use the trace method before any Tonelli-Shanks context is constructed. + if (!square_root_mod_x_even_degree(root, divisor_context, arithmetic_context, + quotient_representative_generator)) { + 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; + } + 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"); + } } } diff --git a/libs/math/test/polynomial_frobenius.cpp b/libs/math/test/polynomial_frobenius.cpp index 1677d04ad..828ff890d 100644 --- a/libs/math/test/polynomial_frobenius.cpp +++ b/libs/math/test/polynomial_frobenius.cpp @@ -406,4 +406,161 @@ BOOST_AUTO_TEST_CASE(characteristic_product_twists_coefficients_and_matches_dire } } +BOOST_AUTO_TEST_CASE(weighted_sums_match_direct_powers_and_support_all_aliases) { + using value = fields::babybear::value_type; + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_type; + const polynomial modulus {value(5), value(7), value(3), value(2)}; + const polynomial input {value(2), value(5), value(7), value(11), value(13), value(17)}; + const polynomial weight {value(19), value(23), value(29), value(31)}; + const boost::multiprecision::cpp_int q = fields::field_order(); + 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, 4, context); + 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}) { + polynomial expected {value::zero()}, actual; + boost::multiprecision::cpp_int exponent = 0, power = 1; + for (std::size_t i = 0; i < count; ++i) { + polynomial factor, term; + math::powmod(factor, weight, exponent, divisor, context); + math::powmod(term, input, power, divisor, context); + math::mulmod(term, factor, term, divisor, context); + math::addition(expected, expected, term); + exponent += power; + power *= step; + } + math::frobenius_weighted_sum(actual, input, weight, count, image, divisor, context); + BOOST_CHECK(actual == expected); + actual = input; + math::frobenius_weighted_sum(actual, actual, weight, count, image, divisor, context); + BOOST_CHECK(actual == expected); + actual = weight; + math::frobenius_weighted_sum(actual, input, actual, count, image, divisor, context); + BOOST_CHECK(actual == expected); + actual = image; + math::frobenius_weighted_sum(actual, input, weight, count, actual, divisor, context); + BOOST_CHECK(actual == expected); + } + } + } +} + +BOOST_AUTO_TEST_CASE(weighted_sums_use_logarithmic_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()}; + const std::size_t count = std::numeric_limits::max(); + polynomial result; + operations = 0; + // F fixes X and X^2 = 1: the summands alternate X, 1, X, 1, ... . + math::frobenius_weighted_sum(result, x, x, count, x, divisor, context); + // Reduce before constructing this 31-bit field value from a full-width size_t. + const auto characteristic = static_cast(fields::babybear::modulus); + BOOST_CHECK(result == polynomial({value((count / 2) % characteristic), value((count / 2 + 1) % characteristic)})); + BOOST_CHECK_LT(operations, 16 * std::numeric_limits::digits); +} + +BOOST_AUTO_TEST_CASE(weighted_sums_fix_fp12_coefficients_and_match_conjugate_sums) { + using value = fq12_value_type; + using backend = polynomial_arithmetic::bn254_small_prime_fft_backend; + using reference_backend = polynomial_arithmetic::schoolbook_backend; + using polynomial = backend::polynomial_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)}, weight {fq12_value(41), fq12_value(53)}; + 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); + for (const bool conjugate : {false, true}) { + const polynomial image {value::zero(), conjugate ? -value::one() : value::one()}; + for (const std::size_t count : {0, 1, 2, 3, 7, 8, 9}) { + polynomial expected {value::zero()}, product {value::one()}, result; + for (std::size_t i = 0; i < count; ++i) { + auto c = input, w = weight; + if (conjugate && i % 2 != 0) { + c[1] = -c[1]; + w[1] = -w[1]; + } + polynomial term; + math::mulmod(term, product, c, reference_divisor, reference); + math::addition(expected, expected, term); + math::mulmod(product, product, w, reference_divisor, reference); + } + math::frobenius_weighted_sum(result, input, weight, count, image, divisor, context); + BOOST_CHECK(result == expected); + } + } +} + +BOOST_AUTO_TEST_CASE(weighted_sums_handle_degenerate_cases_and_preserve_output_on_errors) { + 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 zero {value::zero()}, x {value::zero(), value::one()}; + polynomial image, output; + math::powmod_x(image, fields::field_order(), divisor, context); + for (const auto &weight : {value::zero(), value::one(), value(7)}) { + for (const std::size_t count : {0, 1, 2, 17}) { + value sum = value::zero(), factor = value::one(); + for (std::size_t i = 0; i < count; ++i) { + sum += factor; + factor *= weight; + } + math::frobenius_weighted_sum(output, polynomial {value(11)}, polynomial {weight}, count, image, divisor, + context); + BOOST_CHECK(output == polynomial({value(11) * sum})); + } + } + math::frobenius_weighted_sum(output, zero, x, 17, image, divisor, context); + BOOST_CHECK(output == zero); + math::frobenius_weighted_sum(output, x, modulus, 17, image, divisor, context); + BOOST_CHECK(output == x); + math::frobenius_weighted_sum(output, x, x, 1, image, divisor, context); + BOOST_CHECK(output == x); + BOOST_CHECK_THROW(math::frobenius_weighted_sum(output, x, 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(); + for (const auto &invalid : {empty, padded}) { + BOOST_CHECK_THROW(math::frobenius_weighted_sum(output, invalid, x, 1, image, divisor, context), + std::invalid_argument); + BOOST_CHECK_THROW(math::frobenius_weighted_sum(output, x, invalid, 1, image, divisor, context), + std::invalid_argument); + BOOST_CHECK_THROW(math::frobenius_weighted_sum(output, x, x, 1, invalid, divisor, context), + std::invalid_argument); + BOOST_CHECK(output == x); + } + BOOST_CHECK_THROW(math::frobenius_weighted_sum(output, x, x, 1, modulus, divisor, context), std::invalid_argument); + BOOST_CHECK(output == x); + const math::polynomial_divisor_context constant_divisor(polynomial {value(7)}, 1, context); + for (const std::size_t count : {0, 1, 9}) { + math::frobenius_weighted_sum(output, x, x, count, zero, constant_divisor, context); + BOOST_CHECK(output == 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 05d762cdd..60204fd2c 100644 --- a/libs/math/test/polynomial_square_root.cpp +++ b/libs/math/test/polynomial_square_root.cpp @@ -36,6 +36,8 @@ #include #include +#include + #include #include #include @@ -487,9 +489,9 @@ BOOST_AUTO_TEST_CASE(characteristic_seed_and_image_match_independent_full_powers 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; + const boost::multiprecision::cpp_int exponent = (q - 1) >> 1; polynomial_arithmetic::polynomial_context context; - for (const std::size_t degree : {3, 5, 9}) { + for (const std::size_t degree : {2, 3, 4, 5, 6, 9}) { polynomial modulus(degree + 1, extension_value::zero()); for (auto &coefficient : modulus) { coefficient = nil::crypto3::algebra::random_element(rng); @@ -514,7 +516,7 @@ BOOST_AUTO_TEST_CASE(prime_field_seed_retains_direct_exponentiation) { 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; + 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); @@ -522,4 +524,211 @@ BOOST_AUTO_TEST_CASE(prime_field_seed_retains_direct_exponentiation) { BOOST_CHECK(image == expected_image); } +BOOST_AUTO_TEST_CASE(trace_even_degree_roots_match_tonelli_shanks_for_powers_of_two_and_mixed_degrees) { + polynomial_arithmetic::polynomial_context context; + const polynomial_type x {value_type::zero(), value_type::one()}; + value_type c(2); + while (c.is_square() || c.pow((field_type::modulus - 1u) / 3u) == value_type::one()) { + c += value_type::one(); + } + boost::random::mt19937 rng(0xE7E002); + std::size_t squares = 0, nonsquares = 0; + for (const std::size_t degree : {2, 4, 6, 8, 12, 16, 32}) { + // The binomial criterion applies: c is neither a square nor a cube and |K| = 1 mod 4. + 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); + std::size_t calls = 0; + auto generator = [&] { + BOOST_REQUIRE_LT(++calls, 64); + polynomial_type result(degree); + for (auto &coefficient : result) { + coefficient = nil::crypto3::algebra::random_element(rng); + } + math::condense(result); + return result; + }; + polynomial_type root {value_type(19)}; + BOOST_REQUIRE_EQUAL(math::square_root_mod_x_even_degree(root, divisor, context, generator), + expected_square); + if (expected_square) { + polynomial_type squared, reference, negative_reference; + const math::polynomial_square_root_context tonelli(divisor, context, generator); + BOOST_REQUIRE(math::square_root_mod(reference, x, tonelli, context)); + 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_EQUAL(calls, 0); + } + } + } + } + BOOST_CHECK_GT(squares, 0); + BOOST_CHECK_GT(nonsquares, 0); +} + +BOOST_AUTO_TEST_CASE(trace_even_degree_roots_retry_zero_traces) { + polynomial_arithmetic::polynomial_context context; + // If r^4 + r + 2 = 0 and x = r^2, B(x) = (x^2 + 2)^2 - x = 0. Newton identities give + // Trace(r) = 0 but Trace(r*x) = -3, forcing the initial c = 1 to fail and c = X to succeed. + const polynomial_type monic {value_type(4), -value_type::one(), value_type(4), value_type::zero(), + value_type::one()}; + const polynomial_type x {value_type::zero(), value_type::one()}; + const polynomial_type expected_root {value_type(2), value_type::zero(), value_type::one()}; + const boost::multiprecision::cpp_int q = fields::field_order(); + 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, 3, context); + // Rabin's criterion independently verifies irreducibility of this quartic fixture. + polynomial_type power, difference, gcd; + const boost::multiprecision::cpp_int q2 = q * q, q4 = q2 * q2; + math::powmod_x(power, q2, divisor, context); + math::subtraction(difference, power, x); + math::gcd(gcd, modulus, difference, context); + BOOST_REQUIRE_EQUAL(gcd.size(), 1); + math::powmod_x(power, q4, divisor, context); + BOOST_REQUIRE(power == x); + std::size_t calls = 0; + auto generator = [&] { + BOOST_REQUIRE_LT(++calls, 4); + return calls == 1 ? polynomial_type {value_type::zero()} : + calls == 2 ? polynomial_type {value_type::one()} : + x; + }; + polynomial_type root, squared, negative_expected; + BOOST_REQUIRE(math::square_root_mod_x_even_degree(root, divisor, context, generator)); + BOOST_CHECK_EQUAL(calls, 3); + math::squaremod(squared, root, divisor, context); + BOOST_CHECK(squared == x); + math::scalar_multiplication(negative_expected, expected_root, -value_type::one()); + BOOST_CHECK(root == expected_root || root == negative_expected); + } +} + +BOOST_AUTO_TEST_CASE(trace_even_degree_roots_support_dense_fp12_coefficients_and_small_prime_products) { + using extension_field = fields::fp12_2over3over2>; + using 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(0xE7E012); + const boost::multiprecision::cpp_int noncube_exponent = (fields::field_order() - 1) / 3; + value c = value::one(); + for (std::size_t i = 0; i < 64 && (c.is_square() || c.pow(noncube_exponent) == value::one()); ++i) { + c = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE(!c.is_zero()); + BOOST_REQUIRE(!c.is_square()); + BOOST_REQUIRE(c.pow(noncube_exponent) != value::one()); + polynomial_arithmetic::polynomial_context reference; + const polynomial x {value::zero(), value::one()}; + for (const std::size_t degree : {2, 4, 6, 8}) { + value shift = value::one(); + for (std::size_t i = 0; i < 64 && !(shift.pow(degree) - c).is_square(); ++i) { + shift = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE((shift.pow(degree) - c).is_square()); + polynomial modulus {value::one()}; + const polynomial linear {shift, value::one()}; + for (std::size_t i = 0; i < degree; ++i) { + reference.multiply(modulus, modulus, linear); + } + modulus[0] -= c; + const value scale = nil::crypto3::algebra::random_element(rng); + BOOST_REQUIRE(!scale.is_zero()); + math::scalar_multiplication(modulus, modulus, scale); + BOOST_REQUIRE(modulus[0].Frobenius_map(1) != modulus[0]); + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = degree == 2 ? 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); + std::size_t calls = 0; + auto generator = [&] { + BOOST_REQUIRE_LT(++calls, 64); + polynomial result(degree); + for (auto &coefficient : result) { + coefficient = nil::crypto3::algebra::random_element(rng); + } + math::condense(result); + return result; + }; + polynomial root, squared; + BOOST_REQUIRE(math::square_root_mod_x_even_degree(root, divisor, context, generator)); + math::squaremod(squared, root, reference_divisor, reference); + BOOST_CHECK(squared == x); + BOOST_CHECK_LE(root.size(), degree); + } +} + +BOOST_AUTO_TEST_CASE(trace_even_degree_roots_reject_unsupported_contexts_without_sampling) { + polynomial_arithmetic::polynomial_context context; + polynomial_type output {value_type(19)}; + auto generator = []() -> polynomial_type { + BOOST_FAIL("invalid contexts must not sample"); + return {}; + }; + for (const auto &modulus : + {polynomial_type {value_type::one()}, polynomial_type {value_type::one(), value_type::one()}, + polynomial_type {value_type(2), value_type::zero(), value_type::zero(), value_type::one()}, + polynomial_type {value_type::zero(), value_type::one(), value_type::one()}}) { + const math::polynomial_divisor_context divisor(modulus, 2, context); + BOOST_CHECK_THROW(math::square_root_mod_x_even_degree(output, divisor, context, generator), + std::invalid_argument); + BOOST_CHECK(output == polynomial_type({value_type(19)})); + } + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 0; + polynomial_arithmetic::polynomial_context disabled(backend_type {}, options); + const polynomial_type modulus {value_type::one(), -value_type(79), value_type::one()}; + const math::polynomial_divisor_context divisor(modulus, 1, disabled); + BOOST_CHECK_THROW(math::square_root_mod_x_even_degree(output, divisor, disabled, generator), std::invalid_argument); + BOOST_CHECK(output == polynomial_type({value_type(19)})); +} + +BOOST_AUTO_TEST_CASE(trace_even_degree_roots_validate_retry_candidates_and_propagate_generator_errors) { + polynomial_arithmetic::polynomial_context context; + const polynomial_type modulus {value_type(4), -value_type::one(), value_type(4), value_type::zero(), + value_type::one()}; + const math::polynomial_divisor_context divisor(modulus, 3, context); + polynomial_type empty, padded(2); + empty.get_storage().clear(); + padded[0] = value_type::one(); + padded[1] = value_type::zero(); + polynomial_type output {value_type(19)}; + for (const auto &invalid : {empty, padded, modulus}) { + std::size_t calls = 0; + auto generator = [&] { + ++calls; + return invalid; + }; + BOOST_CHECK_THROW(math::square_root_mod_x_even_degree(output, divisor, context, generator), + std::invalid_argument); + BOOST_CHECK_EQUAL(calls, 1); + BOOST_CHECK(output == polynomial_type({value_type(19)})); + } + auto failing_generator = []() -> polynomial_type { throw std::runtime_error("generator failure"); }; + BOOST_CHECK_THROW(math::square_root_mod_x_even_degree(output, divisor, context, failing_generator), + std::runtime_error); + 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 7b59c8f23..c39ee5558 100644 --- a/libs/math/test/polynomial_x_norm_reconstruction.cpp +++ b/libs/math/test/polynomial_x_norm_reconstruction.cpp @@ -193,7 +193,7 @@ BOOST_AUTO_TEST_CASE(recovers_and_normalizes_an_irreducible_quadratic_with_the_s coefficient_generator); BOOST_REQUIRE(result.has_value()); - BOOST_CHECK_GT(generator_calls, 0); + BOOST_CHECK_EQUAL(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. @@ -595,4 +595,121 @@ BOOST_AUTO_TEST_CASE(screening_preserves_even_multiplicities_and_rejects_odd_obs } } +BOOST_AUTO_TEST_CASE(even_degree_trace_recovery_normalizes_scalars_and_consumes_only_retry_samples) { + polynomial_arithmetic::polynomial_context context; + // This irreducible quartic's root has trace zero; multiplying it by X gives trace -3. + // The square-root tests independently check irreducibility and the forced retry. + const polynomial_type monic {value_type(4), -value_type::one(), value_type(4), value_type::zero(), + value_type::one()}; + for (const auto &scale : {value_type::one(), value_type(49), first_quadratic_non_residue()}) { + polynomial_type g; + math::scalar_multiplication(g, monic, scale); + std::size_t calls = 0; + auto generator = [&] { + const std::size_t index = calls++; + BOOST_REQUIRE_LT(index, 12); + // Three representatives: zero, one (both have zero weighted trace), then X. + return index == 4 || index == 9 ? value_type::one() : value_type::zero(); + }; + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, context, generator); + BOOST_REQUIRE_EQUAL(result.has_value(), scale.is_square()); + BOOST_CHECK_EQUAL(calls, 12); + if (result) { + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, context) == g); + BOOST_CHECK_LE(result->p.size() - 1, 2); + BOOST_CHECK_LE(result->q.size() - 1, 1); + } + } + std::size_t calls = 0; + auto failing_generator = [&]() -> value_type { + ++calls; + throw std::runtime_error("retry coefficient generation failed"); + }; + BOOST_CHECK_THROW(math::recover_irreducible_polynomial_x_norm_representation(monic, context, failing_generator), + std::runtime_error); + BOOST_CHECK_EQUAL(calls, 1); +} + +BOOST_AUTO_TEST_CASE(even_degree_recovery_preserves_tonelli_shanks_when_composition_is_disabled) { + const polynomial_type g = representable_irreducible_quadratic(); + for (const std::size_t cache_limit : {0, 1, 16}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = cache_limit; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + std::size_t calls = 0; + boost::random::mt19937 rng(0xE7E021); + auto generator = [&] { + BOOST_REQUIRE_LT(++calls, 128); + return nil::crypto3::algebra::random_element(rng); + }; + 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); + if (cache_limit == 0) { + BOOST_CHECK_GT(calls, 0); + BOOST_CHECK_EQUAL(calls % 2, 0); + } else { + BOOST_CHECK_EQUAL(calls, 0); + } + const std::size_t calls_before_rejection = calls; + const polynomial_type nonsquare {-first_quadratic_non_residue(), value_type::zero(), value_type::one()}; + BOOST_CHECK(!math::recover_irreducible_polynomial_x_norm_representation(nonsquare, context, generator)); + BOOST_CHECK_EQUAL(calls, calls_before_rejection); + } +} + +BOOST_AUTO_TEST_CASE(even_degree_bn254_fp12_norm_recovery_agrees_across_backends) { + using value = fq12_value_type; + using polynomial = fq12_schoolbook_backend::polynomial_type; + using fft_backend = polynomial_arithmetic::bn254_small_prime_fft_backend; + boost::random::mt19937 rng(0xE7E022); + value c = value::one(); + for (std::size_t i = 0; i < 64 && c.is_square(); ++i) { + c = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE(!c.is_zero()); + BOOST_REQUIRE(!c.is_square()); + polynomial_arithmetic::polynomial_context reference; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context optimized(fft_backend(16), options); + for (const std::size_t degree : {2, 4}) { + value shift = value::one(); + for (std::size_t i = 0; i < 64 && !(shift.pow(degree) - c).is_square(); ++i) { + shift = nil::crypto3::algebra::random_element(rng); + } + BOOST_REQUIRE((shift.pow(degree) - c).is_square()); + polynomial g {value::one()}; + const polynomial linear {shift, value::one()}; + for (std::size_t i = 0; i < degree; ++i) { + reference.multiply(g, g, linear); + } + g[0] -= c; + const value scale = nil::crypto3::algebra::random_element(rng).squared(); + BOOST_REQUIRE(!scale.is_zero()); + math::scalar_multiplication(g, g, scale); + // Fix independent RNG states for equivalent fallback sampling, if a zero trace occurs. + boost::random::mt19937 rng_reference(0xE7E023), rng_optimized(0xE7E023); + auto reference_generator = [&] { + return nil::crypto3::algebra::random_element(rng_reference); + }; + auto optimized_generator = [&] { + return nil::crypto3::algebra::random_element(rng_optimized); + }; + const auto expected = + math::recover_irreducible_polynomial_x_norm_representation(g, reference, reference_generator); + const auto actual = + math::recover_irreducible_polynomial_x_norm_representation(g, optimized, optimized_generator); + BOOST_REQUIRE(expected.has_value()); + BOOST_REQUIRE(actual.has_value()); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*actual, reference) == g); + BOOST_CHECK(actual->p == expected->p); + BOOST_CHECK(actual->q == expected->q); + BOOST_CHECK_LE(actual->p.size() - 1, degree / 2); + BOOST_CHECK_LE(actual->q.size() - 1, (degree - 1) / 2); + BOOST_CHECK(rng_reference == rng_optimized); + } +} + BOOST_AUTO_TEST_SUITE_END()