From b25f9604c311c282b955746a6405172b711d3e2d Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Tue, 20 Oct 2026 13:00:53 -0400 Subject: [PATCH 1/2] fuse normalization and implmemnt certified direct Montgomery reconstruction --- .../math/detail/bn254_small_prime_fft.hpp | 208 +++++++++++++++--- .../bn254_small_prime_fft_backend.hpp | 18 +- libs/math/test/small_prime_fft_crt.cpp | 155 +++++++++++-- 3 files changed, 327 insertions(+), 54 deletions(-) diff --git a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp index 22ae7501b..7507bd000 100644 --- a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp +++ b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp @@ -34,9 +34,11 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -193,8 +195,10 @@ namespace nil::crypto3::math::detail { * with M > 2*C for the prime product M. Centering must precede reduction into the target field: * the CRT representative M-1 denotes -1, not the field element (M-1) mod p. * - * Cofactors and their inverses are prepared once. Reconstruction uses fixed-size integer scratch - * and caller-owned output, with no arbitrary-precision setup work in the coefficient loop. + * Cofactors, inverses and Montgomery export weights are prepared once. An integer interval + * certificate permits direct export when it proves both centering and the coefficient bound; + * uncertain cases use exact wide reconstruction. Both paths use fixed-size scratch and + * caller-owned output, with no arbitrary-precision setup work in the coefficient loop. * This operation has no constant-time guarantee. */ class bn254_small_prime_fft_reconstructor { @@ -209,18 +213,30 @@ namespace nil::crypto3::math::detail { boost::multiprecision::number>; using export_type = boost::multiprecision::uint512_modular_t; static_assert(50 * small_prime_fft_primes.size() <= 600, "CRT export folding requires M < 2^600"); + using fast_params = bn254_small_prime_fft_converter::field_type::extension_policy::fast_params; + using direct_weight_type = std::array; + static constexpr bool direct_export_supported = base_field_type::integral_type::backend_type::limb_bits == 64; + static constexpr unsigned certificate_bits = 60; + static constexpr std::uint64_t certificate_scale = std::uint64_t(1) << certificate_bits; + static_assert(std::all_of(small_prime_fft_primes.begin(), small_prime_fft_primes.end(), [](const auto &prime) { + return prime.modulus > (std::uint64_t(1) << 49) && prime.modulus < (std::uint64_t(1) << 50); + })); struct prime_context { small_prime_fft_arithmetic arithmetic; integer_type cofactor; double inverse; + std::uint64_t certificate_reciprocal; + direct_weight_type export_weight; - prime_context(std::uint64_t modulus, const integer_type &cofactor) : - arithmetic(modulus), cofactor(cofactor) { + prime_context(std::uint64_t modulus, const integer_type &cofactor, const direct_weight_type &weight) : + arithmetic(modulus), cofactor(cofactor), export_weight(weight) { using wide_type = boost::multiprecision::uint128_modular_t; const auto residue = (cofactor % modulus).convert_to(); inverse = static_cast(boost::multiprecision::inverse_mod(wide_type(residue), wide_type(modulus)) .convert_to()); + certificate_reciprocal = + ((wide_type(1) << (certificate_bits + 50)) / modulus).convert_to(); } }; @@ -256,9 +272,24 @@ namespace nil::crypto3::math::detail { multiples_[bit] = multiples_[bit - 1] << 1; } primes_.reserve(params.prime_count); + const boost::multiprecision::cpp_int characteristic(base_field_type::modulus.backend().to_cpp_int()); + const boost::multiprecision::cpp_int r_squared = + (boost::multiprecision::cpp_int(1) << 512) % characteristic; + scaled_bound_ = ((coefficient_bound << certificate_bits) / product).convert_to(); for (std::size_t i = 0; i < params.prime_count; ++i) { const auto modulus = small_prime_fft_primes[i].modulus; - primes_.emplace_back(modulus, multiples_[0] / modulus); + const boost::multiprecision::cpp_int weight = (product / modulus) * r_squared % characteristic; + direct_weight_type limbs {}; + boost::multiprecision::export_bits(weight, limbs.begin(), 64, false); + primes_.emplace_back(modulus, multiples_[0] / modulus, limbs); + } + // The centering quotient is in [0,prime_count]. Store canonical -k*M*R^2 mod p. + for (std::size_t k = 0; k <= params.prime_count; ++k) { + boost::multiprecision::cpp_int correction = (k * product * r_squared) % characteristic; + if (correction != 0) { + correction = characteristic - correction; + } + boost::multiprecision::export_bits(correction, export_corrections_[k].begin(), 64, false); } } @@ -266,6 +297,23 @@ namespace nil::crypto3::math::detail { return primes_.size(); } + /** + * Combine an inverse-transform scale with (M/q)^-1, returning a centered residue. + * Require an integer-valued scale of magnitude < 2^53; invalid prime indices throw. + */ + double weighted_normalization(std::size_t prime, double scale) const { + const auto &context = primes_.at(prime); + const double q = static_cast(context.arithmetic.modulus()); + double result = + context.arithmetic.multiply_integer(context.arithmetic.reduce_integer(scale), context.inverse); + if (result > q / 2) { + result -= q; + } else if (result < -q / 2) { + result += q; + } + return result; + } + /** * Read residues[(prime*12 + coordinate)*stride + coefficient] for each requested output. * Require output.size() <= stride and residues.size() == prime_count()*12*stride. @@ -281,6 +329,23 @@ namespace nil::crypto3::math::detail { * bound when the integer expression has wrapped by M, so the caller must justify that bound. */ void reconstruct(std::span output, std::span residues, std::size_t stride) const { + reconstruct_dispatch(output, residues, stride); + } + + /** + * Reconstruct residues already multiplied by (M/q)^-1 modulo each q. All dimensions, + * residue bounds and failure semantics are those of reconstruct(). The backend can fuse + * this weight with 1/N in the pointwise pass because the inverse transform is linear. + */ + void reconstruct_weighted(std::span output, std::span residues, + std::size_t stride) const { + reconstruct_dispatch(output, residues, stride); + } + + private: + template + void reconstruct_dispatch(std::span output, std::span residues, + std::size_t stride) const { const std::size_t channels = field_type::arity * primes_.size(); if (channels == 0 || output.size() > stride || stride > residues.size() / channels || residues.size() != channels * stride) { @@ -293,13 +358,12 @@ namespace nil::crypto3::math::detail { return prime.arithmetic.fma_supported(); }); if (use_fma) { - reconstruct_impl(output, residues, stride); + reconstruct_impl(output, residues, stride); } else { - reconstruct_impl(output, residues, stride); + reconstruct_impl(output, residues, stride); } } - private: void add_cofactor_product(integer_type &sum, const integer_type &cofactor, std::uint64_t digit) const { using limb_type = boost::multiprecision::limb_type; using wide_type = boost::multiprecision::double_limb_type; @@ -338,19 +402,90 @@ namespace nil::crypto3::math::detail { return base_field_type::value_type(folded); } - template + base_field_type::value_type export_exact(integer_type magnitude) const { + // Each canonical CRT term is < M. Reduce their sum to [0,M) without division. + for (std::size_t bit = std::bit_width(primes_.size() - 1); bit-- > 0;) { + if (magnitude >= multiples_[bit]) { + magnitude -= multiples_[bit]; + } + } + const bool negative = magnitude > half_modulus_; + if (negative) { + magnitude = multiples_[0] - magnitude; + } + if (magnitude > bound_) { + throw std::overflow_error("bn254_small_prime_fft_reconstructor: coefficient exceeds CRT bound"); + } + const auto value = export_magnitude(magnitude); + return negative ? -value : value; + } + + static void add_export_product(fast_params::limb_array &sum, const direct_weight_type &weight, + std::uint64_t digit) { + using wide_type = algebra::fields::detail::fp12_fast::wide_limb; + wide_type carry = 0; + for (std::size_t limb = 0; limb < weight.size(); ++limb) { + const wide_type value = wide_type(weight[limb]) * digit + sum[limb] + carry; + sum[limb] = static_cast(value); + carry = value >> 64; + } + // At most twelve digits < 2^50 times weights < p, plus one correction < p. + // The total is < 2^54*p < 2^308, so the fifth limb cannot overflow. + sum[weight.size()] += static_cast(carry); + } + + bool certify_export(std::uint64_t lower, std::size_t "ient) const { + // For d < q < 2^50 and a=floor(2^110/q), floor(d*a/2^50) <= 2^60*d/q + // < floor(d*a/2^50)+2. Summing encloses 2^60*sum(d/q) in [lower,upper]. + // With at most twelve primes, all quantities below fit in uint64_t. + const std::uint64_t upper = lower + 2 * primes_.size(); + quotient = (lower + certificate_scale / 2) >> certificate_bits; + const std::uint64_t center = std::uint64_t(quotient) << certificate_bits; + // Containment in k +/- C/M proves the bound AND the unique centering quotient, + // since C/M < 1/2. Using floor(2^60*C/M) makes this test conservative. + return lower + scaled_bound_ >= center && upper <= center + scaled_bound_; + } + + template + base_field_type::value_type export_direct(typename Params::limb_array &sum, std::size_t quotient) const { + add_export_product(sum, export_corrections_[quotient], 1); + // sum = (sum(d_i*M_i)-k*M)*R^2 mod p, and sum < 2^54*p < R*p for R=2^256. + // One existing REDC therefore yields a canonical Montgomery coordinate directly, + // without materializing the wide CRT integer, folding, or a second field conversion. + direct_weight_type reduced; + algebra::fields::detail::fp12_fast::montgomery_reduce(reduced.data(), sum); + base_field_type::value_type result; + std::copy(reduced.begin(), reduced.end(), result.data.backend().base_data().limbs()); + return result; + } + + template void reconstruct_impl(std::span output, std::span residues, std::size_t stride) const { for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { for (std::size_t offset = 0; offset < output.size();) { const std::size_t count = std::min(block_size, output.size() - offset); - std::array sums {}; + using accumulator_type = std::conditional_t; + std::array sums {}; + std::array lower {}; + std::array, small_prime_fft_primes.size()> saved_digits; for (std::size_t prime = 0; prime < primes_.size(); ++prime) { const auto &context = primes_[prime]; const auto &arithmetic = context.arithmetic; const double modulus = static_cast(arithmetic.modulus()); const double *source = residues.data() + (prime * field_type::arity + coordinate) * stride + offset; + const auto accumulate = [&](std::size_t i, std::uint64_t digit) { + if constexpr (Direct) { + using wide_type = algebra::fields::detail::fp12_fast::wide_limb; + saved_digits[prime][i] = digit; + lower[i] += static_cast( + (wide_type(digit) * context.certificate_reciprocal) >> 50); + add_export_product(sums[i], context.export_weight, digit); + } else { + add_cofactor_product(sums[i], context.cofactor, digit); + } + }; std::size_t i = 0; #if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) if constexpr (UseFma) { @@ -359,15 +494,16 @@ namespace nil::crypto3::math::detail { const __m256d inverse = _mm256_set1_pd(context.inverse); const __m256d q = _mm256_set1_pd(modulus); for (; count - i >= 4; i += 4) { - __m256d digits = - packed.multiply_fma(packed.reduce_fma(_mm256_loadu_pd(source + i)), inverse); + __m256d digits = packed.reduce_fma(_mm256_loadu_pd(source + i)); + if constexpr (!Weighted) { + digits = packed.multiply_fma(digits, inverse); + } const __m256d negative = _mm256_cmp_pd(digits, _mm256_setzero_pd(), _CMP_LT_OQ); digits = _mm256_add_pd(digits, _mm256_and_pd(negative, q)); std::array canonical; _mm256_storeu_pd(canonical.data(), digits); for (std::size_t lane = 0; lane < canonical.size(); ++lane) { - add_cofactor_product(sums[i + lane], context.cofactor, - static_cast(canonical[lane])); + accumulate(i + lane, static_cast(canonical[lane])); } } } @@ -376,36 +512,38 @@ namespace nil::crypto3::math::detail { for (; i < count; ++i) { double digit; if constexpr (UseFma) { - digit = arithmetic.multiply_fma(arithmetic.reduce_fma(source[i]), context.inverse); + digit = arithmetic.reduce_fma(source[i]); + if constexpr (!Weighted) { + digit = arithmetic.multiply_fma(digit, context.inverse); + } } else { - digit = - arithmetic.multiply_integer(arithmetic.reduce_integer(source[i]), context.inverse); + digit = arithmetic.reduce_integer(source[i]); + if constexpr (!Weighted) { + digit = arithmetic.multiply_integer(digit, context.inverse); + } } // Both operands have magnitude < q, so the product remainder is in (-q,q). // Canonical digits give 0 <= digit*(M/q) < M and bound the complete sum. const auto canonical = static_cast(digit < 0 ? digit + modulus : digit); - add_cofactor_product(sums[i], context.cofactor, canonical); + accumulate(i, canonical); } } for (std::size_t i = 0; i < count; ++i) { - auto &magnitude = sums[i]; - // The sum is < prime_count*M. Subtract binary multiples of M from high to low - // to obtain [0,M), using at most four comparisons/subtractions and no division. - for (std::size_t bit = std::bit_width(primes_.size() - 1); bit-- > 0;) { - if (magnitude >= multiples_[bit]) { - magnitude -= multiples_[bit]; + auto &destination = output[offset + i].coordinate(coordinate); + if constexpr (Direct) { + std::size_t quotient; + if (certify_export(lower[i], quotient)) { + destination = export_direct(sums[i], quotient); + } else { + integer_type exact_sum = 0; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][i]); + } + destination = export_exact(exact_sum); } + } else { + destination = export_exact(sums[i]); } - const bool negative = magnitude > half_modulus_; - if (negative) { - magnitude = multiples_[0] - magnitude; - } - if (magnitude > bound_) { - throw std::overflow_error( - "bn254_small_prime_fft_reconstructor: coefficient exceeds CRT bound"); - } - const auto value = export_magnitude(magnitude); - output[offset + i].coordinate(coordinate) = negative ? -value : value; } offset += count; } @@ -418,6 +556,8 @@ namespace nil::crypto3::math::detail { integer_type bound_; std::size_t accumulator_limbs_; export_type export_radix_; + std::uint64_t scaled_bound_; + std::array export_corrections_ {}; }; /** diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp index 8843a6a1e..57fa79ad2 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp @@ -288,7 +288,7 @@ namespace nil::crypto3::math::polynomial_arithmetic { auto destination = std::span(left_spectrum_).subspan(offset, count); pointwise_[prime].multiply_sum(destination, spectrum(2 * row), spectrum(4 + column), spectrum(2 * row + 1), spectrum(4 + Columns + column), size, - plan.size_inverse()); + normalizations_[plan_offset + prime]); for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { plan.inverse_fft(destination.subspan(coordinate * size, size), false); } @@ -296,8 +296,8 @@ namespace nil::crypto3::math::polynomial_arithmetic { const std::size_t index = row * Columns + column; auto &result = results[index]; result.resize(result_sizes[index]); - reconstructor_.reconstruct(std::span(result.data(), result.size()), left_spectrum_, - size); + reconstructor_.reconstruct_weighted(std::span(result.data(), result.size()), + left_spectrum_, size); math::condense(result); } } @@ -494,6 +494,7 @@ namespace nil::crypto3::math::polynomial_arithmetic { converter_(params.crt.prime_count), reconstructor_(params.crt, params.coefficient_bound) { const unsigned depth = std::countr_zero(max_transform_size_); plans_.reserve((depth + 1) * prime_count()); + normalizations_.reserve((depth + 1) * prime_count()); pointwise_.reserve(prime_count()); for (std::size_t prime = 0; prime < prime_count(); ++prime) { pointwise_.emplace_back(math::detail::small_prime_fft_primes[prime].modulus); @@ -501,6 +502,9 @@ namespace nil::crypto3::math::polynomial_arithmetic { for (unsigned level = 0; level <= depth; ++level) { for (std::size_t prime = 0; prime < prime_count(); ++prime) { plans_.emplace_back(prime, std::size_t(1) << level); + // Prepared forward spectra stay unweighted; apply this once to the product. + normalizations_.push_back( + reconstructor_.weighted_normalization(prime, plans_.back().size_inverse())); } } } @@ -571,16 +575,17 @@ namespace nil::crypto3::math::polynomial_arithmetic { } } if constexpr (Square) { - pointwise_[prime].square(a, size, plan.size_inverse()); + pointwise_[prime].square(a, size, normalizations_[plan_offset + prime]); } else { - pointwise_[prime].multiply(a, b, size, plan.size_inverse()); + pointwise_[prime].multiply(a, b, size, normalizations_[plan_offset + prime]); } for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { plan.inverse_fft(a.subspan(coordinate * size, size), false); } } result_.resize(result_size); - reconstructor_.reconstruct(std::span(result_.data(), result_size), left_spectrum_, size); + reconstructor_.reconstruct_weighted(std::span(result_.data(), result_size), left_spectrum_, + size); if (!cyclic_dividend.empty()) { // Here result_size==size and degree(V)<=2*n-2<2*size, so at most one high coefficient // folds into each cell. Keep the complete difference: any uncancelled tail stays visible. @@ -601,6 +606,7 @@ namespace nil::crypto3::math::polynomial_arithmetic { math::detail::bn254_small_prime_fft_converter converter_; math::detail::bn254_small_prime_fft_reconstructor reconstructor_; std::vector plans_; + std::vector normalizations_; std::vector pointwise_; std::vector left_spectrum_; std::vector right_spectrum_; diff --git a/libs/math/test/small_prime_fft_crt.cpp b/libs/math/test/small_prime_fft_crt.cpp index 3faf7ccad..6fb091d26 100644 --- a/libs/math/test/small_prime_fft_crt.cpp +++ b/libs/math/test/small_prime_fft_crt.cpp @@ -73,14 +73,20 @@ namespace { } std::vector make_residues(std::span coordinates, std::size_t prime_count, - std::size_t stride, bool lazy = false) { + std::size_t stride, bool lazy = false, bool weighted = false) { // Padding is intentionally not a valid residue: reconstruction must read only the requested prefix. std::vector result(prime_count * 12 * stride, std::numeric_limits::quiet_NaN()); for (std::size_t prime = 0; prime < prime_count; ++prime) { const cpp_int modulus = detail::small_prime_fft_primes[prime].modulus; + cpp_int inverse = 1; + if (weighted) { + // Compute the weight independently of the reconstructor's cached inverses. + const cpp_int cofactor = make_params(prime_count).modulus / modulus; + inverse = boost::multiprecision::powm(cpp_int(cofactor % modulus), cpp_int(modulus - 2), modulus); + } for (std::size_t coordinate = 0; coordinate < 12; ++coordinate) { for (std::size_t i = 0; i < coordinates.size(); ++i) { - cpp_int residue = canonical(coordinates[i][coordinate], modulus); + cpp_int residue = canonical(coordinates[i][coordinate] * inverse, modulus); if (lazy) { // Reach large positive and negative representatives while staying below 2^53. residue += ((prime + coordinate + i) % 2 == 0 ? 7 : -8) * modulus; @@ -106,6 +112,9 @@ namespace { const cpp_int actual( output[i].data[w].data[v].data[u].to_integral().backend().to_cpp_int()); BOOST_CHECK_EQUAL(actual, canonical(coordinates[i][coordinate], characteristic)); + const cpp_int stored( + output[i].data[w].data[v].data[u].data.backend().base_data().to_cpp_int()); + BOOST_CHECK_LT(stored, characteristic); } } } @@ -123,6 +132,13 @@ namespace { BOOST_CHECK(storage.front() == sentinel); BOOST_CHECK(storage.back() == sentinel); check_output(output, coordinates); + + const auto weighted = make_residues(coordinates, reconstructor.prime_count(), stride, lazy, true); + std::fill(output.begin(), output.end(), sentinel); + reconstructor.reconstruct_weighted(output, weighted, stride); + BOOST_CHECK(storage.front() == sentinel); + BOOST_CHECK(storage.back() == sentinel); + check_output(output, coordinates); } std::vector random_coordinates(std::size_t count, const cpp_int &bound) { @@ -143,6 +159,29 @@ namespace { BOOST_AUTO_TEST_SUITE(small_prime_fft_crt_test_suite) +BOOST_AUTO_TEST_CASE(combined_normalization_is_centered_and_exact_in_every_rounding_mode) { + restore_rounding restore; + for (int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + const auto params = make_params(count); + const reconstructor_type reconstructor(params, params.modulus / 2); + for (std::size_t prime = 0; prime < count; ++prime) { + const cpp_int q = detail::small_prime_fft_primes[prime].modulus; + for (std::size_t size : {1, 2, 4, 8, 32, 256}) { + const detail::small_prime_fft_plan plan(prime, size); + const double weight = reconstructor.weighted_normalization(prime, plan.size_inverse()); + BOOST_CHECK_LT(std::abs(weight), q.convert_to() / 2); + const cpp_int exact_weight = static_cast(weight); + BOOST_CHECK_EQUAL(canonical(exact_weight * size * (params.modulus / q), q), 1); + } + } + BOOST_CHECK_THROW(reconstructor.weighted_normalization(count, 1), std::out_of_range); + } + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } +} + BOOST_AUTO_TEST_CASE(every_catalogue_prefix_recovers_signed_boundaries_and_carries) { const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { @@ -200,6 +239,91 @@ BOOST_AUTO_TEST_CASE(random_signed_integers_cover_all_prime_counts_and_partial_b } } +BOOST_AUTO_TEST_CASE(direct_export_covers_every_centering_correction_and_maximum_digits) { + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + const auto params = make_params(count); + const cpp_int bound = params.modulus / 2; + const reconstructor_type reconstructor(params, bound); + for (std::size_t quotient = 0; quotient <= count; ++quotient) { + cpp_int sum = 0; + for (std::size_t prime = 0; prime < count; ++prime) { + const cpp_int q = detail::small_prime_fft_primes[prime].modulus; + const cpp_int digit = (q - 1) * quotient / count; + sum += digit * (params.modulus / q); + } + BOOST_REQUIRE_EQUAL((sum + bound) / params.modulus, quotient); + std::array coordinates; + coordinates[0].fill(sum - quotient * params.modulus); + check_reconstruction(reconstructor, coordinates, 3, true); + } + } +} + +BOOST_AUTO_TEST_CASE(exact_fallback_preserves_tiny_bounds_and_values_adjacent_to_bounds) { + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + const auto params = make_params(count); + // Include bounds smaller than the certificate's resolution and bounds near M/2. + const std::vector bounds { + 0, 1, params.modulus >> 60, (params.modulus >> 60) * 32 + 7, params.modulus / 4, params.modulus / 2}; + for (const auto &bound : bounds) { + const reconstructor_type reconstructor(params, bound); + auto coordinates = random_coordinates(33, bound); + coordinates[0][0] = bound; + coordinates[0][1] = -bound; + check_reconstruction(reconstructor, coordinates, 35, true); + if (bound == params.modulus / 2) { + continue; + } + for (int sign : {-1, 1}) { + coordinates[0][0] = sign * (bound + 1); + for (bool weighted : {false, true}) { + const auto residues = make_residues(coordinates, count, 35, true, weighted); + std::vector output(coordinates.size(), value_type::one()); + if (weighted) { + BOOST_CHECK_THROW(reconstructor.reconstruct_weighted(output, residues, 35), + std::overflow_error); + } else { + BOOST_CHECK_THROW(reconstructor.reconstruct(output, residues, 35), std::overflow_error); + } + BOOST_CHECK(std::all_of(output.begin(), output.end(), + [](const auto &value) { return value == value_type::one(); })); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(bound_failure_preserves_completed_prefixes_across_block_boundaries) { + restore_rounding restore; + const auto plan = detail::make_bn254_small_prime_fft_params(14939, 14939); + BOOST_REQUIRE(plan); + for (const cpp_int &bound : {cpp_int(0), plan->coefficient_bound}) { + const reconstructor_type reconstructor(plan->crt, bound); + for (int mode : {FE_TONEAREST, FE_UPWARD}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + for (std::size_t bad : {0, 1, 31, 32, 33}) { + for (int sign : {-1, 1}) { + std::vector coordinates(35); + coordinates[bad][0] = sign * (bound + 1); + for (bool weighted : {false, true}) { + const auto residues = make_residues(coordinates, plan->crt.prime_count, 37, true, weighted); + std::vector output(coordinates.size(), value_type::one()); + if (weighted) { + BOOST_CHECK_THROW(reconstructor.reconstruct_weighted(output, residues, 37), + std::overflow_error); + } else { + BOOST_CHECK_THROW(reconstructor.reconstruct(output, residues, 37), std::overflow_error); + } + for (std::size_t i = 0; i < output.size(); ++i) { + BOOST_CHECK(output[i] == (i < bad ? value_type::zero() : value_type::one())); + } + } + } + } + } + } +} + BOOST_AUTO_TEST_CASE(wide_export_preserves_split_boundaries_and_field_multiples) { const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); for (std::size_t count = 8; count <= detail::small_prime_fft_primes.size(); ++count) { @@ -408,18 +532,21 @@ BOOST_AUTO_TEST_CASE(invalid_dimensions_are_rejected_before_output_writes) { const reconstructor_type reconstructor(make_params(1), 10); std::vector residues(24, 0); std::vector output(2, value_type::one()); - BOOST_CHECK_THROW(reconstructor.reconstruct(output, std::span(residues).first(23), 2), - std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct(output, residues, 1), std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct({}, residues, 1), std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct(output, std::span(residues).first(12), 1), - std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct({}, residues, std::numeric_limits::max()), - std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct({}, residues, 0), std::invalid_argument); - BOOST_CHECK_THROW(reconstructor.reconstruct(output, {}, 2), std::invalid_argument); - BOOST_CHECK( - std::all_of(output.begin(), output.end(), [](const auto &value) { return value == value_type::one(); })); + for (auto reconstruct : {&reconstructor_type::reconstruct, &reconstructor_type::reconstruct_weighted}) { + BOOST_CHECK_THROW((reconstructor.*reconstruct)(output, std::span(residues).first(23), 2), + std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)(output, residues, 1), std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)({}, residues, 1), std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)(output, std::span(residues).first(12), 1), + std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)({}, residues, std::numeric_limits::max()), + std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)({}, residues, 0), std::invalid_argument); + BOOST_CHECK_THROW((reconstructor.*reconstruct)(output, {}, 2), std::invalid_argument); + BOOST_CHECK( + std::all_of(output.begin(), output.end(), [](const auto &value) { return value == value_type::one(); })); + BOOST_CHECK_NO_THROW((reconstructor.*reconstruct)({}, {}, 0)); + } } BOOST_AUTO_TEST_SUITE_END() From e058e69723c57547e8137b600a67cd50e6858798 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Tue, 20 Oct 2026 13:46:55 -0400 Subject: [PATCH 2/2] fix build --- .../detail/element/fp12_fast/limb_ops_x86.hpp | 32 ++++++++++--------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp index 2573b40b0..1ccbec8a3 100644 --- a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp +++ b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp @@ -254,20 +254,20 @@ #define REDC_8_LIMBS_CANCEL_LOW(I) \ /* m = data[i] * p_dash */ \ "movq " T8(I, 0) ", %%rdx\n" \ - "imulq %[p_dash], %%rdx\n" \ + "imulq " PTR(constants, 4) ", %%rdx\n" \ "xor %[zero], %[zero]\n" \ /* multiply m by each modulus limb and add into the rotating window */ \ - "mulxq %[p0], %[low], %[high]\n" \ + "mulxq " PTR(constants, 0) ", %[low], %[high]\n" \ /* use dual carry chains */ \ "adcx %[low], " T8(I, 0) "\n" \ "adox %[high], " T8(I, 1) "\n" \ - "mulxq %[p1], %[low], %[high]\n" \ + "mulxq " PTR(constants, 1) ", %[low], %[high]\n" \ "adcx %[low], " T8(I, 1) "\n" \ "adox %[high], " T8(I, 2) "\n" \ - "mulxq %[p2], %[low], %[high]\n" \ + "mulxq " PTR(constants, 2) ", %[low], %[high]\n" \ "adcx %[low], " T8(I, 2) "\n" \ "adox %[high], " T8(I, 3) "\n" \ - "mulxq %[p3], %[low], %[high]\n" \ + "mulxq " PTR(constants, 3) ", %[low], %[high]\n" \ "adcx %[low], " T8(I, 3) "\n" \ /* merge OF, prior pending carry, and CF into the next live limb */ \ "adox %[zero], %[high]\n" \ @@ -697,7 +697,13 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { template requires(Params::storage_limb_count == 8) inline void montgomery_reduce_x86(limb *result, const limb *data) { - GET_MODULUS_4_LIMBS(Params::base_field_type); + constexpr auto mod_obj = Params::base_field_type::modulus_params.get_mod_obj(); + // One explicit address operand keeps these constants from competing with scratch + // registers, including default-visibility Mach-O and unoptimized builds. + static constexpr std::array constants = { + limb(mod_obj.get_mod().limbs()[0]), limb(mod_obj.get_mod().limbs()[1]), + limb(mod_obj.get_mod().limbs()[2]), limb(mod_obj.get_mod().limbs()[3]), + limb(mod_obj.get_p_dash())}; limb t0, t1, t2, t3, t4; limb low, high, pending, zero; asm volatile( @@ -723,10 +729,10 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { "movq " T8(3, 4) ", %%rdx\n" // Try q = t - p. - "subq %[p0], %[low]\n" - "sbbq %[p1], %[high]\n" - "sbbq %[p2], %[zero]\n" - "sbbq %[p3], %%rdx\n" + "subq " PTR(constants, 0) ", %[low]\n" + "sbbq " PTR(constants, 1) ", %[high]\n" + "sbbq " PTR(constants, 2) ", %[zero]\n" + "sbbq " PTR(constants, 3) ", %%rdx\n" // If t - p did not borrow, keep q as the canonical result. "cmovnc %[low], " T8(0, 4) "\n" @@ -750,11 +756,7 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { [zero]"=&r"(zero) : [data]"r"(data), [result]"r"(result), - [p0]"m"(p0), - [p1]"m"(p1), - [p2]"m"(p2), - [p3]"m"(p3), - [p_dash]"m"(p_dash) + [constants]"r"(constants.data()) : "rdx", "cc", "memory" ); }