From 013266a26dbe875db8bb5cbdefdfb644fc512b5f Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Wed, 14 Oct 2026 16:35:20 -0400 Subject: [PATCH 01/11] bounds and prime plans --- .../math/detail/small_prime_fft_params.hpp | 178 +++++++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_params.cpp | 238 ++++++++++++++++++ 3 files changed, 417 insertions(+) create mode 100644 libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp create mode 100644 libs/math/test/small_prime_fft_params.cpp diff --git a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp new file mode 100644 index 000000000..08d4e3b03 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp @@ -0,0 +1,178 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#ifndef CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP +#define CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace nil::crypto3::math::detail { + + struct small_prime_fft_prime { + std::uint64_t modulus; + std::uint64_t root; // Exact order 2^two_adicity; shorter roots are obtained by squaring. + unsigned two_adicity; + }; + + // Twelve 50-bit primes with large power-of-two factors in modulus - 1 provide CRT capacity + // for full-width BN254 Fp12 products. The tests certify primality and the exact root orders. + // Prime selection uses a prefix of this table. FMA arithmetic and its range checks are separate + // from these integer reconstruction parameters; no floating-point arithmetic is used here. + inline constexpr std::array small_prime_fft_primes = {{ + {1108307720798209ULL, 194751219211145ULL, 44}, + {659706976665601ULL, 422061946989193ULL, 43}, + {1086317488242689ULL, 813873581740013ULL, 42}, + {910395627798529ULL, 212465101045679ULL, 42}, + {699289395265537ULL, 643876812462936ULL, 42}, + {1022545813831681ULL, 136965991847555ULL, 41}, + {1013749720809473ULL, 549289857321204ULL, 41}, + {868614185943041ULL, 322646709537195ULL, 41}, + {864216139431937ULL, 46940759627233ULL, 41}, + {837827860365313ULL, 128509036413366ULL, 41}, + {815837627809793ULL, 584069832576029ULL, 41}, + {600333348765697ULL, 295748699958414ULL, 41}, + }}; + + struct small_prime_fft_crt_params { + std::size_t prime_count; + boost::multiprecision::cpp_int modulus; + }; + + /** + * Select the shortest nonempty catalogue prefix with M > 2*coefficient_bound, where every + * prime supports transform_size. Then each integer in [-coefficient_bound, coefficient_bound] + * has a unique centered CRT representative. The bound must cover the integer expression being + * reconstructed, before reduction modulo the target field characteristic. + * + * Return nullopt for a negative bound, a non-power-of-two size (including zero), insufficient + * root order, or insufficient CRT capacity. Arbitrary-precision setup arithmetic avoids overflow. + */ + inline std::optional + select_small_prime_fft_crt_params(std::size_t transform_size, + const boost::multiprecision::cpp_int &coefficient_bound) { + if (!std::has_single_bit(transform_size) || coefficient_bound < 0) { + return std::nullopt; + } + const unsigned depth = std::countr_zero(transform_size); + const boost::multiprecision::cpp_int twice_bound = 2 * coefficient_bound; + boost::multiprecision::cpp_int modulus = 1; + for (std::size_t i = 0; i < small_prime_fft_primes.size(); ++i) { + if (depth > small_prime_fft_primes[i].two_adicity) { + return std::nullopt; + } + modulus *= small_prime_fft_primes[i].modulus; + if (modulus > twice_bound) { + return small_prime_fft_crt_params {i + 1, std::move(modulus)}; + } + } + return std::nullopt; + } + + /** + * Bound one coordinate of a sum of term_count BN254 Fp12 products, with input coordinates + * lifted canonically into [0, p-1]. In the tower u^2=-1, v^3=9+u, w^2=v, use the equivalent + * basis w^j*u^k (0<=j<6, 0<=k<2). At output w^j there are j+1 unwrapped pairs of w powers, + * with absolute coordinate weight 2 each, and 5-j wrapped pairs, with weight 2*(9+1) each. + * Thus each output coordinate has weight 2*(j+1)+20*(5-j) <= 102. The triangle inequality + * gives C = 102*term_count*(p-1)^2, including the signs introduced by u^2=-1. + * + * This is a bound on the final integer tower product, not on lazy arithmetic intermediates. + * Cyclic input folding, digit splitting and other input lifts require their own bounds. + */ + inline boost::multiprecision::cpp_int bn254_fp12_convolution_bound(std::size_t term_count) { + const boost::multiprecision::cpp_int largest_coordinate = + algebra::fields::field_characteristic>() - 1; + return 102 * boost::multiprecision::cpp_int(term_count) * largest_coordinate * largest_coordinate; + } + + struct small_prime_fft_params { + std::size_t left_size; + std::size_t right_size; + std::size_t result_size; + std::size_t transform_size; + std::size_t spectrum_bytes; // One spectrum: 12 coordinates * prime_count * transform_size 64-bit cells. + boost::multiprecision::cpp_int coefficient_bound; + small_prime_fft_crt_params crt; + }; + + /** + * Plan a full BN254 Fp12 product, or its first coefficient_count coefficients. Only input prefixes + * that can affect that output are used, but the transform still contains their COMPLETE linear + * product: shortening it to the output count would permit cyclic aliasing into low coefficients. + * Square uses the same plan with equal operand lengths. No transform buffers are allocated here. + * + * Return nullopt when either input length or coefficient_count is zero (the caller handles a zero + * result without a transform), when size arithmetic would overflow, when the catalogue cannot + * supply the roots/CRT capacity, or when one spectrum exceeds max_spectrum_bytes. This byte limit + * covers only one spectrum, not all future operation scratch, twiddle tables or prepared operands. + */ + inline std::optional + make_bn254_small_prime_fft_params(std::size_t left_size, std::size_t right_size, + std::size_t coefficient_count = std::numeric_limits::max(), + std::size_t max_spectrum_bytes = std::numeric_limits::max()) { + left_size = std::min(left_size, coefficient_count); + right_size = std::min(right_size, coefficient_count); + if (left_size == 0 || right_size == 0 || + left_size > std::numeric_limits::max() - (right_size - 1)) { + return std::nullopt; + } + const std::size_t product_size = left_size + (right_size - 1); + // std::bit_ceil requires the result to be representable in size_t. + if (product_size > (std::size_t(1) << (std::numeric_limits::digits - 1))) { + return std::nullopt; + } + const std::size_t transform_size = std::bit_ceil(product_size); + boost::multiprecision::cpp_int bound = bn254_fp12_convolution_bound(std::min(left_size, right_size)); + auto crt = select_small_prime_fft_crt_params(transform_size, bound); + if (!crt) { + return std::nullopt; + } + const std::size_t bytes_per_index = 12 * crt->prime_count * sizeof(std::uint64_t); + if (transform_size > max_spectrum_bytes / bytes_per_index) { + return std::nullopt; + } + return small_prime_fft_params {left_size, + right_size, + std::min(product_size, coefficient_count), + transform_size, + transform_size * bytes_per_index, + std::move(bound), + std::move(*crt)}; + } + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PARAMS_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index 94c805953..043a88b84 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -46,6 +46,7 @@ set(TESTS_NAMES "evaluation_domain" "geometric_sequence_domain" "integer_sqrt" + "small_prime_fft_params" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_params.cpp b/libs/math/test/small_prime_fft_params.cpp new file mode 100644 index 000000000..9c6901341 --- /dev/null +++ b/libs/math/test/small_prime_fft_params.cpp @@ -0,0 +1,238 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_params_test + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace { + namespace detail = nil::crypto3::math::detail; + namespace fields = nil::crypto3::algebra::fields; + using boost::multiprecision::cpp_int; + using fq = fields::alt_bn128_base_field<254>; + using fq12 = fields::fp12_2over3over2; + + // Coordinates are [1,u,v,uv,v^2,uv^2,w,uw,vw,uvw,v^2w,uv^2w]. + std::size_t coordinate_index(unsigned w_power, unsigned u_power) { + return 6 * (w_power % 2) + 2 * (w_power / 2) + u_power; + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_params_test_suite) + +BOOST_AUTO_TEST_CASE(catalogue_has_distinct_certified_primes_and_exact_roots) { + for (std::size_t i = 0; i < detail::small_prime_fft_primes.size(); ++i) { + const auto &entry = detail::small_prime_fft_primes[i]; + BOOST_TEST_CONTEXT("prime index " << i) { + BOOST_CHECK_GE(entry.modulus, std::uint64_t(1) << 49); + BOOST_CHECK_LT(entry.modulus, std::uint64_t(1) << 50); + BOOST_CHECK_EQUAL(std::countr_zero(entry.modulus - 1), entry.two_adicity); + BOOST_CHECK_GT(entry.root, 0); + BOOST_CHECK_LT(entry.root, entry.modulus); + for (std::size_t j = 0; j < i; ++j) { + BOOST_CHECK_NE(entry.modulus, detail::small_prime_fft_primes[j].modulus); + } + + const cpp_int modulus = entry.modulus; + const cpp_int order = cpp_int(1) << entry.two_adicity; + const cpp_int root = entry.root; + const cpp_int half_power = boost::multiprecision::powm(root, cpp_int(order / 2), modulus); + BOOST_REQUIRE_EQUAL(half_power, modulus - 1); + BOOST_REQUIRE_EQUAL(boost::multiprecision::powm(root, order, modulus), 1); + + // Every prime divisor r of modulus must give root order exactly 2^s, hence r >= 2^s+1. + // Since (2^s)^2 > modulus, a composite modulus is impossible. This certifies primality, + // rather than relying on a probabilistic primality test or assuming the table is prime. + BOOST_REQUIRE_GT(order * order, modulus); + for (unsigned depth = 0; depth <= entry.two_adicity; ++depth) { + const cpp_int short_order = cpp_int(1) << depth; + const cpp_int short_root = boost::multiprecision::powm(root, cpp_int(order / short_order), modulus); + BOOST_CHECK_EQUAL(boost::multiprecision::powm(short_root, short_order, modulus), 1); + if (depth != 0) { + BOOST_CHECK_EQUAL(boost::multiprecision::powm(short_root, cpp_int(short_order / 2), modulus), + modulus - 1); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(crt_selection_checks_exact_capacity_boundaries) { + cpp_int product = 1; + for (std::size_t i = 0; i < detail::small_prime_fft_primes.size(); ++i) { + product *= detail::small_prime_fft_primes[i].modulus; + // The odd product M distinguishes all integers from -(M-1)/2 through +(M-1)/2. + const auto fits = detail::select_small_prime_fft_crt_params(1, (product - 1) / 2); + BOOST_REQUIRE(fits); + BOOST_CHECK_EQUAL(fits->prime_count, i + 1); + BOOST_CHECK_EQUAL(fits->modulus, product); + const auto next = detail::select_small_prime_fft_crt_params(1, (product + 1) / 2); + if (i + 1 < detail::small_prime_fft_primes.size()) { + BOOST_REQUIRE(next); + BOOST_CHECK_EQUAL(next->prime_count, i + 2); + } else { + BOOST_CHECK(!next); + } + } + BOOST_CHECK(!detail::select_small_prime_fft_crt_params(0, 1)); + BOOST_CHECK(!detail::select_small_prime_fft_crt_params(3, 1)); + BOOST_CHECK(!detail::select_small_prime_fft_crt_params(1, -1)); + const auto zero_bound = detail::select_small_prime_fft_crt_params(1, 0); + BOOST_REQUIRE(zero_bound); + BOOST_CHECK_EQUAL(zero_bound->prime_count, 1); +} + +BOOST_AUTO_TEST_CASE(tower_bound_matches_all_basis_products_in_crypto3_coordinate_order) { + const cpp_int modulus = fields::field_characteristic(); + std::array weights {}; + for (unsigned a = 0; a < 12; ++a) { + for (unsigned b = 0; b < 12; ++b) { + fq12::value_type left = fq12::value_type::zero(); + fq12::value_type right = fq12::value_type::zero(); + left.coordinate(a) = fq::value_type::one(); + right.coordinate(b) = fq::value_type::one(); + const auto actual = left * right; + + const unsigned w_power = a / 6 + (a % 6 / 2) * 2 + b / 6 + (b % 6 / 2) * 2; + const unsigned u_power = a % 2 + b % 2; + const bool wraps = w_power >= 6; + std::array expected {}; + // Reduce w^6 to 9+u, then reduce each u^2 to -1 using signed integer arithmetic. + for (unsigned term = 0; term < (wraps ? 2U : 1U); ++term) { + const unsigned exponent = u_power + term; + const int magnitude = wraps && term == 0 ? 9 : 1; + expected[coordinate_index(w_power % 6, exponent % 2)] += + (exponent / 2 % 2 == 0 ? magnitude : -magnitude); + } + for (unsigned c = 0; c < 12; ++c) { + cpp_int coordinate(actual.coordinate(c).to_integral().backend().to_cpp_int()); + if (coordinate > modulus / 2) { + coordinate -= modulus; + } + BOOST_CHECK_EQUAL(coordinate, expected[c]); + weights[c] += expected[c] < 0 ? -expected[c] : expected[c]; + } + } + } + const std::array expected_weights = {102, 102, 66, 66, 30, 30, 84, 84, 48, 48, 12, 12}; + BOOST_CHECK_EQUAL_COLLECTIONS(weights.begin(), weights.end(), expected_weights.begin(), expected_weights.end()); + const cpp_int maximum = std::numeric_limits::max(); + const cpp_int expected = + *std::max_element(weights.begin(), weights.end()) * maximum * (modulus - 1) * (modulus - 1); + BOOST_CHECK_EQUAL(detail::bn254_fp12_convolution_bound(std::numeric_limits::max()), expected); + BOOST_CHECK_EQUAL(detail::bn254_fp12_convolution_bound(0), 0); +} + +BOOST_AUTO_TEST_CASE(workload_plan_covers_signed_fp12_products) { + const auto plan = detail::make_bn254_small_prime_fft_params(14939, 14939); + BOOST_REQUIRE(plan); + BOOST_CHECK_EQUAL(plan->left_size, 14939); + BOOST_CHECK_EQUAL(plan->right_size, 14939); + BOOST_CHECK_EQUAL(plan->result_size, 29877); + BOOST_CHECK_EQUAL(plan->transform_size, 32768); + BOOST_CHECK_EQUAL(plan->crt.prime_count, 11); + BOOST_CHECK_EQUAL(plan->spectrum_bytes, 33 * 1024 * 1024); + const cpp_int modulus = fields::field_characteristic(); + const cpp_int twice_bound = cpp_int(204) * 14939 * (modulus - 1) * (modulus - 1); + BOOST_CHECK_EQUAL(2 * plan->coefficient_bound, twice_bound); + BOOST_CHECK_EQUAL(boost::multiprecision::msb(twice_bound) + 1, 529); + BOOST_CHECK_GT(plan->crt.modulus, twice_bound); + BOOST_CHECK_LE(plan->crt.modulus / detail::small_prime_fft_primes[10].modulus, twice_bound); +} + +BOOST_AUTO_TEST_CASE(linear_and_low_shapes_do_not_alias_high_coefficients) { + for (std::size_t left : {1, 2, 3, 8, 9, 16, 17}) { + for (std::size_t right : {1, 2, 7, 8, 9, 16, 33}) { + for (std::size_t count : {1, 2, 7, 8, 9, 16, 17, 64}) { + const auto plan = detail::make_bn254_small_prime_fft_params(left, right, count); + BOOST_REQUIRE(plan); + const std::size_t used_left = std::min(left, count); + const std::size_t used_right = std::min(right, count); + const std::size_t full_size = used_left + used_right - 1; + BOOST_CHECK_EQUAL(plan->left_size, used_left); + BOOST_CHECK_EQUAL(plan->right_size, used_right); + BOOST_CHECK_EQUAL(plan->result_size, std::min(count, full_size)); + BOOST_CHECK(std::has_single_bit(plan->transform_size)); + BOOST_CHECK_GE(plan->transform_size, full_size); + BOOST_CHECK_LT(plan->transform_size / 2, full_size); + BOOST_CHECK_EQUAL(plan->coefficient_bound, + detail::bn254_fp12_convolution_bound(std::min(used_left, used_right))); + } + } + } + const auto low = detail::make_bn254_small_prime_fft_params(100, 100, 8); + BOOST_REQUIRE(low); + BOOST_CHECK_EQUAL(low->result_size, 8); + BOOST_CHECK_EQUAL(low->transform_size, 16); +} + +BOOST_AUTO_TEST_CASE(empty_overflow_and_storage_limits_are_checked_before_allocation) { + constexpr auto maximum = std::numeric_limits::max(); + constexpr auto top_bit = std::size_t(1) << (std::numeric_limits::digits - 1); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(0, 1)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(1, 0)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(1, 1, 0)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(maximum, 2)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(2, maximum)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(top_bit, 2)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(1, maximum)); + // Trimming before size arithmetic makes this small low product valid, even with huge input lengths. + const auto trimmed = detail::make_bn254_small_prime_fft_params(maximum, maximum, 1); + BOOST_REQUIRE(trimmed); + BOOST_CHECK_EQUAL(trimmed->transform_size, 1); + const auto plan = detail::make_bn254_small_prime_fft_params(14939, 14939); + BOOST_REQUIRE(plan); + BOOST_CHECK(detail::make_bn254_small_prime_fft_params(14939, 14939, maximum, plan->spectrum_bytes)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(14939, 14939, maximum, plan->spectrum_bytes - 1)); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(1, 1, 1, 0)); +} + +BOOST_AUTO_TEST_CASE(root_capacity_and_larger_prime_prefix_are_checked_without_buffers) { + // These lengths exercise metadata only; the planner never allocates the associated spectra. + if constexpr (std::numeric_limits::digits > 44) { + const std::size_t largest = std::uint64_t(1) << 41; + const auto plan = detail::make_bn254_small_prime_fft_params(largest / 2, largest / 2); + BOOST_REQUIRE(plan); + BOOST_CHECK_EQUAL(plan->transform_size, largest); + BOOST_CHECK_EQUAL(plan->crt.prime_count, 12); + BOOST_CHECK(!detail::make_bn254_small_prime_fft_params(largest, 2)); + // At depth 42 the first five primes support the root, but cannot cover a single Fp12 product. + BOOST_CHECK(!detail::select_small_prime_fft_crt_params(largest * 2, detail::bn254_fp12_convolution_bound(1))); + BOOST_CHECK(detail::select_small_prime_fft_crt_params(std::uint64_t(1) << 44, 1)); + BOOST_CHECK(!detail::select_small_prime_fft_crt_params(std::uint64_t(1) << 45, 1)); + } +} + +BOOST_AUTO_TEST_SUITE_END() From 89d053a6db64217666ff633c3ca982385f30caa8 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Wed, 14 Oct 2026 17:38:59 -0400 Subject: [PATCH 02/11] scalar residue arithmetic --- .../detail/small_prime_fft_arithmetic.hpp | 182 +++++++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_arithmetic.cpp | 249 ++++++++++++++++++ 3 files changed, 432 insertions(+) create mode 100644 libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp create mode 100644 libs/math/test/small_prime_fft_arithmetic.cpp diff --git a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp new file mode 100644 index 000000000..1b1ca1082 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp @@ -0,0 +1,182 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#ifndef CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP +#define CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace nil::crypto3::math::detail { + + /** + * Scalar residue arithmetic for an odd 50-bit modulus q. Values are integer-valued binary64 + * numbers, not approximate coefficients. Primality is a caller concern (the prime catalogue + * is separately certified); these arithmetic operations also work for odd composite moduli. + * + * Select the kernel once per operation with fma_supported(). The FMA methods require this + * predicate to remain true until the operation finishes, including FE_TONEAREST on the calling + * thread. The integer methods are the exact fallback for other rounding modes/compiler flags. + * Keeping this dispatch outside each scalar operation avoids a rounding-mode query per butterfly. + * No method changes the floating-point environment. + */ + class small_prime_fft_arithmetic { + using wide_type = boost::multiprecision::uint128_modular_t; + + static_assert(std::numeric_limits::radix == 2 && std::numeric_limits::digits == 53 && + std::numeric_limits::max_exponent == 1024, + "small-prime residues require binary64 doubles"); + + public: + explicit small_prime_fft_arithmetic(std::uint64_t modulus) : modulus_(modulus) { + if (modulus < (std::uint64_t(1) << 49) || modulus >= (std::uint64_t(1) << 50) || modulus % 2 == 0) { + throw std::invalid_argument("small_prime_fft_arithmetic: expected an odd 50-bit modulus"); + } + + // Since 2^49 < q < 2^50, the nearest binary64 reciprocal is S/2^102, where + // S = round(2^102/q) fits exactly in 53 bits. An odd q cannot produce a halfway tie. + // Construct S using integers so that setup is independent of the active rounding mode. + const wide_type scale = wide_type(1) << 102; + const wide_type significand = (scale + modulus / 2) / modulus; + reciprocal_ = std::ldexp(static_cast(significand.convert_to()), -102); + + // For h = RN(a*b), l = fma(a,b,-h), and k = nearbyint(h*reciprocal), let + // epsilon = reciprocal - 1/q and |a*b| < 2^e. Then h+l is exact and + // |a*b/q - k| <= |a*b*epsilon| + 2^(e-54)*reciprocal + rounding_error. + // For |a*b| < 2*q^2 the rounding_error is at most 5/8. With + // D = |q*S - 2^102|, the remaining numerator over 2^102 is + // T = 2*q*D + 2^(e-54)*S. Requiring T < 3*2^99 makes the total < 1. + // For |a*b| < 4*q^2, e increases by one and the rounding_error is 3/4: + // the condition becomes 2*T < 3*2^100, giving an output strictly inside (-1.5*q,1.5*q). + // All certificate intermediates fit within 128 bits. + const wide_type scaled_reciprocal = significand * modulus; + const wide_type error = scaled_reciprocal > scale ? scaled_reciprocal - scale : scale - scaled_reciprocal; + const wide_type product_limit = 2 * wide_type(modulus) * modulus; + const unsigned exponent = boost::multiprecision::msb(product_limit) + 1; + const wide_type numerator = 2 * wide_type(modulus) * error + (significand << (exponent - 54)); + reciprocal_satisfies_bounds_ = numerator < (wide_type(3) << 99); + } + + std::uint64_t modulus() const { + return modulus_; + } + + double reciprocal() const { + return reciprocal_; + } + + /** Check once before calling FMA kernels; use the integer kernels when this returns false. */ + bool fma_supported() const { +#if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || defined(__RECIPROCAL_MATH__) || defined(_M_FP_FAST) || \ + (defined(__FINITE_MATH_ONLY__) && __FINITE_MATH_ONLY__ > 0) || FLT_EVAL_METHOD != 0 + return false; +#else + return std::numeric_limits::is_iec559 && reciprocal_satisfies_bounds_ && + std::fegetround() == FE_TONEAREST; +#endif + } + + /** + * Return an integer representative of value modulo q strictly inside (-q,q). + * @pre value is an integer with |value| < 2^53, and fma_supported() is true. + */ + double reduce_fma(double value) const { +#if defined(__clang__) + // Clang does not expose every unsafe floating-point option through predefined macros. +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + assert(valid_integer(value)); + assert(fma_supported()); + // |value/q| < 16, so quotient rounding error is < 1. The FMA remainder is an + // exactly representable integer of magnitude < q, even for a large input representative. + const double quotient = std::nearbyint(value * reciprocal_); + return std::fma(-quotient, static_cast(modulus_), value); + } + + /** + * Return an integer representative of left*right modulo q. The output is strictly inside + * (-q,q) when |left*right| < 2*q^2, and inside (-1.5*q,1.5*q) when |left*right| < 4*q^2. + * @pre both operands are integers of magnitude < 2^53, |left*right| < 4*q^2, and + * fma_supported() is true. These are integer product bounds, not rounded double tests. + */ + double multiply_fma(double left, double right) const { +#if defined(__clang__) + // Preserve the rounded high product even when the caller enables reassociation/contraction. +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + assert(valid_integer(left) && valid_integer(right)); + assert(product_magnitude(left, right) < 4 * wide_type(modulus_) * modulus_); + assert(fma_supported()); + const double high = left * right; + const double low = std::fma(left, right, -high); + const double quotient = std::nearbyint(high * reciprocal_); + // Both the compensated subtraction and the final remainder are small enough to be + // exact integers in binary64. Preserve the split product: reassociation is not allowed. + return std::fma(-quotient, static_cast(modulus_), high) + low; + } + + /** Exact fallback for reduce_fma, with the same input/output bounds and any rounding mode. */ + double reduce_integer(double value) const { + assert(valid_integer(value)); + return static_cast(static_cast(value) % static_cast(modulus_)); + } + + /** + * Exact fallback: both operands must be integers of magnitude < 2^53. Their full product + * fits in 106 bits, so the existing 128-bit integer type suffices. The result is inside (-q,q). + */ + double multiply_integer(double left, double right) const { + assert(valid_integer(left) && valid_integer(right)); + const wide_type residue = product_magnitude(left, right) % modulus_; + const double result = static_cast(residue.convert_to()); + return (left < 0) != (right < 0) ? -result : result; + } + + private: + static bool valid_integer(double value) { + return std::isfinite(value) && std::abs(value) < 0x1p53 && std::trunc(value) == value; + } + + // Call only after checking the integer/range preconditions, so both conversions are exact. + static wide_type product_magnitude(double left, double right) { + return wide_type(static_cast(std::abs(left))) * static_cast(std::abs(right)); + } + + std::uint64_t modulus_; + double reciprocal_; + bool reciprocal_satisfies_bounds_; + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_ARITHMETIC_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index 043a88b84..ebdc4d655 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -47,6 +47,7 @@ set(TESTS_NAMES "geometric_sequence_domain" "integer_sqrt" "small_prime_fft_params" + "small_prime_fft_arithmetic" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_arithmetic.cpp b/libs/math/test/small_prime_fft_arithmetic.cpp new file mode 100644 index 000000000..c3a2a8436 --- /dev/null +++ b/libs/math/test/small_prime_fft_arithmetic.cpp @@ -0,0 +1,249 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_arithmetic_test + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + namespace detail = nil::crypto3::math::detail; + using boost::multiprecision::cpp_int; + using boost::multiprecision::cpp_rational; + using arithmetic_type = detail::small_prime_fft_arithmetic; + + struct restore_rounding { + int saved = std::fegetround(); + ~restore_rounding() { + std::fesetround(saved); + } + }; + + cpp_int positive_residue(cpp_int value, std::uint64_t modulus) { + value %= modulus; + return value < 0 ? cpp_int(value + modulus) : value; + } + + void check_residue(double value, const cpp_int &expected, std::uint64_t modulus, double bound) { + BOOST_REQUIRE(std::isfinite(value)); + BOOST_REQUIRE_EQUAL(value, std::trunc(value)); + BOOST_REQUIRE_LT(std::abs(value), bound); + BOOST_CHECK_EQUAL(positive_residue(static_cast(value), modulus), expected); + } + + void check_product(const arithmetic_type &arithmetic, std::int64_t left, std::int64_t right) { + const auto modulus = arithmetic.modulus(); + const cpp_int product = cpp_int(left) * right; + const cpp_int expected = positive_residue(product, modulus); + const cpp_int magnitude = boost::multiprecision::abs(product); + check_residue(arithmetic.multiply_integer(static_cast(left), static_cast(right)), expected, + modulus, static_cast(modulus)); + if (arithmetic.fma_supported() && magnitude < 4 * cpp_int(modulus) * modulus) { + const double bound = magnitude < 2 * cpp_int(modulus) * modulus ? double(modulus) : 1.5 * modulus; + check_residue(arithmetic.multiply_fma(static_cast(left), static_cast(right)), expected, + modulus, bound); + } + } + + void check_reduction(const arithmetic_type &arithmetic, std::int64_t value) { + const auto modulus = arithmetic.modulus(); + const cpp_int expected = positive_residue(value, modulus); + check_residue(arithmetic.reduce_integer(static_cast(value)), expected, modulus, double(modulus)); + if (arithmetic.fma_supported()) { + check_residue(arithmetic.reduce_fma(static_cast(value)), expected, modulus, double(modulus)); + } + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_arithmetic_test_suite) + +BOOST_AUTO_TEST_CASE(reciprocals_have_exact_error_certificates) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + for (const auto &prime : detail::small_prime_fft_primes) { + BOOST_TEST_CONTEXT("modulus " << prime.modulus) { + const arithmetic_type arithmetic(prime.modulus); + BOOST_CHECK_EQUAL(arithmetic.modulus(), prime.modulus); + // Decode the stored binary64 value and verify its error using independent rational arithmetic. + const auto bits = std::bit_cast(arithmetic.reciprocal()); + BOOST_CHECK_EQUAL(bits >> 63, 0); + const int exponent = int((bits >> 52) & 0x7ff) - 1023 - 52; + BOOST_REQUIRE_EQUAL(exponent, -102); + const cpp_int significand = (bits & ((std::uint64_t(1) << 52) - 1)) | (std::uint64_t(1) << 52); + const cpp_rational reciprocal = cpp_rational(significand) / (cpp_int(1) << -exponent); + const cpp_rational error = boost::multiprecision::abs(reciprocal - cpp_rational(1) / prime.modulus); + BOOST_CHECK_LT(error, cpp_rational(1) / (cpp_int(1) << 103)); + + for (unsigned factor : {2, 4}) { + const cpp_int product_limit = factor * cpp_int(prime.modulus) * prime.modulus; + const unsigned product_exponent = boost::multiprecision::msb(product_limit) + 1; + const cpp_rational rounding_error = factor == 2 ? cpp_rational(5) / 8 : cpp_rational(3) / 4; + const cpp_rational quotient_error = + product_limit * error + (cpp_int(1) << (product_exponent - 54)) * reciprocal + rounding_error; + BOOST_CHECK_LT(quotient_error, factor == 2 ? cpp_rational(1) : cpp_rational(3) / 2); + } +#if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || defined(__RECIPROCAL_MATH__) || defined(_M_FP_FAST) || \ + (defined(__FINITE_MATH_ONLY__) && __FINITE_MATH_ONLY__ > 0) || FLT_EVAL_METHOD != 0 + BOOST_CHECK(!arithmetic.fma_supported()); +#else + BOOST_REQUIRE(arithmetic.fma_supported()); +#endif + } + } +} + +BOOST_AUTO_TEST_CASE(products_cover_signs_zeros_and_lazy_range_boundaries) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + for (const auto &prime : detail::small_prime_fft_primes) { + const arithmetic_type arithmetic(prime.modulus); + const auto q = static_cast(prime.modulus); + std::vector values {0, + 1, + q / 2, + q / 2 + 1, + q - 1, + q, + q + 1, + 3 * q / 2, + 3 * q / 2 + 1, + 2 * q - 1, + 2 * q, + 2 * q + 1, + 3 * q - 1, + 3 * q, + 3 * q + 1, + 4 * q - 1, + (std::int64_t(1) << 53) - 1}; + const auto positive_count = values.size(); + for (std::size_t i = 1; i < positive_count; ++i) { + values.push_back(-values[i]); + } + for (auto left : values) { + for (auto right : values) { + BOOST_TEST_CONTEXT("q=" << q << " left=" << left << " right=" << right) { + // The integer path covers all pairs. The FMA path is checked exactly when its product + // precondition holds, including (2*q-1)*(2*q+1) = 4*q^2-1 and both signs. + check_product(arithmetic, left, right); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(random_products_match_the_integer_oracle) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + std::mt19937_64 random(0xc0ffee); + for (const auto &prime : detail::small_prime_fft_primes) { + const arithmetic_type arithmetic(prime.modulus); + const auto q = static_cast(prime.modulus); + for (unsigned i = 0; i < 2000; ++i) { + const auto left = static_cast(random() % (8 * q - 1)) - (4 * q - 1); + const auto right = static_cast(random() % (2 * q - 1)) - (q - 1); + BOOST_TEST_CONTEXT("q=" << q << " left=" << left << " right=" << right) { + check_product(arithmetic, left, right); + } + } + } +} + +BOOST_AUTO_TEST_CASE(reductions_cover_rounding_boundaries_and_large_integer_representatives) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + std::mt19937_64 random(0x123456); + constexpr std::int64_t limit = std::int64_t(1) << 53; + for (const auto &prime : detail::small_prime_fft_primes) { + const arithmetic_type arithmetic(prime.modulus); + const auto q = static_cast(prime.modulus); + for (std::int64_t multiple = -15; multiple <= 15; ++multiple) { + for (std::int64_t offset : + {-q / 2 - 1, -q / 2, std::int64_t(-1), std::int64_t(0), std::int64_t(1), q / 2, q / 2 + 1}) { + const auto value = multiple * q + offset; + if (value > -limit && value < limit) { + BOOST_TEST_CONTEXT("q=" << q << " value=" << value) { + check_reduction(arithmetic, value); + } + } + } + } + check_reduction(arithmetic, limit - 1); + check_reduction(arithmetic, 1 - limit); + for (unsigned i = 0; i < 2000; ++i) { + check_reduction(arithmetic, static_cast(random() % (2 * limit - 1)) - (limit - 1)); + } + } +} + +BOOST_AUTO_TEST_CASE(integer_fallback_and_setup_work_in_all_rounding_modes) { + restore_rounding restore; + for (const auto &prime : detail::small_prime_fft_primes) { + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + const arithmetic_type original(prime.modulus); + const auto reciprocal_bits = std::bit_cast(original.reciprocal()); + const bool supported_at_nearest = original.fma_supported(); + const auto q = static_cast(prime.modulus); + for (int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + const arithmetic_type current(prime.modulus); + BOOST_CHECK_EQUAL(std::fegetround(), mode); + BOOST_CHECK_EQUAL(std::bit_cast(current.reciprocal()), reciprocal_bits); + BOOST_CHECK_EQUAL(current.fma_supported(), mode == FE_TONEAREST && supported_at_nearest); + BOOST_CHECK_EQUAL(original.fma_supported(), mode == FE_TONEAREST && supported_at_nearest); + check_product(current, -2 * q + 1, 2 * q + 1); + check_product(original, (std::int64_t(1) << 53) - 1, -q); + check_reduction(current, -3 * q - 1); + check_reduction(original, (std::int64_t(1) << 53) - 1); + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } + } +} + +BOOST_AUTO_TEST_CASE(invalid_moduli_are_rejected_and_valid_endpoints_work) { + for (auto modulus : + {std::uint64_t(0), std::uint64_t(1), (std::uint64_t(1) << 49) - 1, std::uint64_t(1) << 49, + std::uint64_t(1) << 50, (std::uint64_t(1) << 50) + 1, std::numeric_limits::max()}) { + BOOST_CHECK_THROW(arithmetic_type {modulus}, std::invalid_argument); + } + for (auto modulus : {(std::uint64_t(1) << 49) + 1, (std::uint64_t(1) << 50) - 1}) { + const arithmetic_type arithmetic(modulus); + const auto q = static_cast(modulus); + check_product(arithmetic, q - 1, -q + 1); + check_reduction(arithmetic, -4 * q + 1); + } +} + +BOOST_AUTO_TEST_SUITE_END() From d152424d2167ee020f2f500013ae4218f800dc4d Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Wed, 14 Oct 2026 18:07:12 -0400 Subject: [PATCH 03/11] coordinate conversion --- .../math/detail/bn254_small_prime_fft.hpp | 188 +++++++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_conversion.cpp | 249 ++++++++++++++++++ 3 files changed, 438 insertions(+) create mode 100644 libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp create mode 100644 libs/math/test/small_prime_fft_conversion.cpp 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 new file mode 100644 index 000000000..f35aea544 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp @@ -0,0 +1,188 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#ifndef CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP +#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace nil::crypto3::math::detail { + + /** + * Convert canonical BN254 Fp12 coordinates to a prefix of the small-prime catalogue. + * Coordinate order is Crypto3's [1,u,v,uv,v^2,uv^2,w,uw,vw,uvw,v^2w,uv^2w]. + * Each output cell is an exact integer in [0,q), congruent to the canonical base-field + * coordinate modulo q. In particular, a base-field -1 is lifted to p-1 before reduction. + * + * Precomputed powers are reusable across operands and sizes. Conversion extracts each + * coordinate once for all primes, using fixed-size scratch and caller-owned output storage. + * Kernel selection occurs once per call; the calling thread's rounding mode must remain + * unchanged until it returns. This conversion has no constant-time guarantee. + */ + class bn254_small_prime_fft_converter { + static constexpr std::size_t chunk_bits = 52; + static constexpr std::size_t chunk_count = 5; + static constexpr std::size_t block_size = 32; + using wide_type = boost::multiprecision::uint128_modular_t; + + struct prime_context { + small_prime_fft_arithmetic arithmetic; + std::array powers; + + explicit prime_context(std::uint64_t modulus) : arithmetic(modulus) { + std::uint64_t power = 1; + for (auto &weight : powers) { + power = ((wide_type(power) << chunk_bits) % modulus).convert_to(); + // Center the weights so that |weight| < q/2. These conversions are exact + // under every rounding mode, and setup products fit in 102 bits. + const auto centered = static_cast(power) - + (power > modulus / 2 ? static_cast(modulus) : 0); + weight = static_cast(centered); + } + } + }; + + public: + using field_type = algebra::fields::fp12_2over3over2>; + using value_type = field_type::value_type; + + explicit bn254_small_prime_fft_converter(std::size_t prime_count) { + if (prime_count == 0 || prime_count > small_prime_fft_primes.size()) { + throw std::invalid_argument("bn254_small_prime_fft_converter: invalid catalogue prefix"); + } + primes_.reserve(prime_count); + for (std::size_t i = 0; i < prime_count; ++i) { + primes_.emplace_back(small_prime_fft_primes[i].modulus); + } + } + + std::size_t prime_count() const { + return primes_.size(); + } + + /** + * Write output[(prime*12 + coordinate)*stride + coefficient], padding each channel + * with zeros from coefficients.size() to stride. Pass a subspan to convert a prefix. + * stride need not be a power of two; the transform planner supplies that restriction. + * + * Require coefficients.size() <= stride and output.size() == prime_count()*12*stride. + * Check dimensions without overflow before writing any output; throw invalid_argument + * on mismatch. Empty input writes zeros, including when reusing a previous output buffer. + */ + void convert(std::span output, std::span coefficients, std::size_t stride) const { + const std::size_t channels = field_type::arity * primes_.size(); + if (channels == 0 || coefficients.size() > stride || stride > output.size() / channels || + output.size() != channels * stride) { + throw std::invalid_argument("bn254_small_prime_fft_converter: invalid buffer dimensions"); + } + if (output.empty()) { + return; + } + const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { + return prime.arithmetic.fma_supported(); + }); + if (use_fma) { + convert_impl(output, coefficients, stride); + } else { + convert_impl(output, coefficients, stride); + } + } + + private: + template + void convert_impl(std::span output, std::span coefficients, + std::size_t stride) const { + for (std::size_t channel = 0; channel < field_type::arity * primes_.size(); ++channel) { + auto padding = output.subspan(channel * stride + coefficients.size(), stride - coefficients.size()); + std::fill(padding.begin(), padding.end(), 0.0); + } + + std::array, block_size> chunks; + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + for (std::size_t offset = 0; offset < coefficients.size();) { + const std::size_t count = std::min(block_size, coefficients.size() - offset); + for (std::size_t i = 0; i < count; ++i) { + std::array words {}; + // Use the canonical value, not the field's internal Montgomery representation. + const auto canonical = coefficients[offset + i].coordinate(coordinate).to_integral(); + boost::multiprecision::export_bits(canonical, words.begin(), chunk_bits, false); + for (std::size_t j = 0; j < chunk_count; ++j) { + chunks[i][j] = static_cast(words[j]); + } + } + 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()); + double *destination = + output.data() + (prime * field_type::arity + coordinate) * stride + offset; + for (std::size_t i = 0; i < count; ++i) { + double sum; + if constexpr (UseFma) { + sum = arithmetic.reduce_fma(chunks[i][0]); + } else { + sum = arithmetic.reduce_integer(chunks[i][0]); + } + // Each chunk is < 2^52 and each weight has magnitude < q/2, so the + // products are < 4*q^2 (q > 2^49). The initial remainder is < q; + // adding four terms of magnitude < 1.5*q keeps every partial sum + // inside (-7*q,7*q), hence below 2^53 and exactly representable. + for (std::size_t j = 1; j < chunk_count; ++j) { + if constexpr (UseFma) { + sum += arithmetic.multiply_fma(chunks[i][j], context.powers[j - 1]); + } else { + sum += arithmetic.multiply_integer(chunks[i][j], context.powers[j - 1]); + } + } + if constexpr (UseFma) { + sum = arithmetic.reduce_fma(sum); + } else { + sum = arithmetic.reduce_integer(sum); + } + destination[i] = sum < 0 ? sum + modulus : sum; + } + } + offset += count; + } + } + } + + std::vector primes_; + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index ebdc4d655..bfab84111 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -48,6 +48,7 @@ set(TESTS_NAMES "integer_sqrt" "small_prime_fft_params" "small_prime_fft_arithmetic" + "small_prime_fft_conversion" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_conversion.cpp b/libs/math/test/small_prime_fft_conversion.cpp new file mode 100644 index 000000000..fe6fcab70 --- /dev/null +++ b/libs/math/test/small_prime_fft_conversion.cpp @@ -0,0 +1,249 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_conversion_test + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + namespace detail = nil::crypto3::math::detail; + using boost::multiprecision::cpp_int; + using fq = nil::crypto3::algebra::fields::alt_bn128_base_field<254>; + using converter_type = detail::bn254_small_prime_fft_converter; + using value_type = converter_type::value_type; + using coordinates_type = std::array; + constexpr double sentinel = -12345.0; + + struct restore_rounding { + int saved = std::fegetround(); + ~restore_rounding() { + std::fesetround(saved); + } + }; + + std::vector make_coefficients(std::span coordinates) { + std::vector result(coordinates.size(), value_type::zero()); + for (std::size_t i = 0; i < result.size(); ++i) { + // Assign the tower components directly, independently of coordinate() in the converter. + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + result[i].data[w].data[v].data[u] = + fq::value_type(fq::integral_type(coordinates[i][6 * w + 2 * v + u].str())); + } + } + } + } + return result; + } + + void check_output(std::span output, std::size_t prime_count, + std::span coordinates, std::size_t stride) { + const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); + BOOST_REQUIRE_EQUAL(output.size(), prime_count * 12 * stride); + for (std::size_t prime = 0; prime < prime_count; ++prime) { + const auto modulus = detail::small_prime_fft_primes[prime].modulus; + for (std::size_t coordinate = 0; coordinate < 12; ++coordinate) { + for (std::size_t i = 0; i < stride; ++i) { + BOOST_TEST_CONTEXT("prime=" << prime << " coordinate=" << coordinate << " coefficient=" << i) { + const double actual = output[(prime * 12 + coordinate) * stride + i]; + BOOST_REQUIRE(std::isfinite(actual)); + BOOST_REQUIRE_EQUAL(actual, std::trunc(actual)); + BOOST_REQUIRE_GE(actual, 0.0); + BOOST_REQUIRE_LT(actual, static_cast(modulus)); + cpp_int expected = 0; + if (i < coordinates.size()) { + expected = (coordinates[i][coordinate] % characteristic) % modulus; + } + BOOST_CHECK_EQUAL(static_cast(actual), expected); + } + } + } + } + } + + std::vector check_conversion(const converter_type &converter, std::span coefficients, + std::span coordinates, std::size_t stride) { + const std::size_t cells = converter.prime_count() * 12 * stride; + std::vector storage(cells + 2, sentinel); + auto output = std::span(storage).subspan(1, cells); + converter.convert(output, coefficients, stride); + BOOST_CHECK_EQUAL(storage.front(), sentinel); + BOOST_CHECK_EQUAL(storage.back(), sentinel); + check_output(output, converter.prime_count(), coordinates, stride); + return storage; + } + + std::vector random_coordinates(std::size_t count) { + std::mt19937_64 random(0x51a11); + std::vector result(count); + for (auto &coefficient : result) { + for (auto &coordinate : coefficient) { + for (unsigned limb = 0; limb < 4; ++limb) { + coordinate <<= 64; + coordinate += random(); + } + coordinate &= (cpp_int(1) << 254) - 1; + } + } + return result; + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_conversion_test_suite) + +BOOST_AUTO_TEST_CASE(basis_and_coordinate_layout_are_preserved_for_every_catalogue_prefix) { + std::vector coordinates(14); + for (std::size_t c = 0; c < 12; ++c) { + coordinates[c][c] = 1; + coordinates[12][c] = 100 + c; + } + const auto coefficients = make_coefficients(coordinates); + for (std::size_t count = 1; count <= detail::small_prime_fft_primes.size(); ++count) { + const converter_type converter(count); + BOOST_CHECK_EQUAL(converter.prime_count(), count); + check_conversion(converter, coefficients, coordinates, 19); + } +} + +BOOST_AUTO_TEST_CASE(canonical_values_and_chunk_boundaries_match_integer_oracle) { + const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); + std::vector values = { + 0, 1, characteristic - 1, characteristic, characteristic + 1, (cpp_int(1) << 254) - 1}; + // Every bit position and its neighbors cover chunk boundaries, carries, and sparse coordinates. + for (unsigned bit = 0; bit < 254; ++bit) { + const cpp_int power = cpp_int(1) << bit; + values.push_back(power - 1); + values.push_back(power); + values.push_back(power + 1); + } + for (const auto &prime : detail::small_prime_fft_primes) { + for (int delta : {-1, 0, 1}) { + values.push_back(cpp_int(prime.modulus) + delta); + values.push_back((characteristic / prime.modulus) * prime.modulus + delta); + } + } + std::vector coordinates((values.size() + 11) / 12); + for (std::size_t i = 0; i < values.size(); ++i) { + coordinates[i / 12][i % 12] = values[i]; + } + auto coefficients = make_coefficients(coordinates); + // Exercise a field negation as well as construction from canonical and noncanonical integers. + coefficients[0].data[0].data[1].data[0] = -fq::value_type::one(); + BOOST_REQUIRE_EQUAL(coordinates[0][2], characteristic - 1); + check_conversion(converter_type(12), coefficients, coordinates, std::bit_ceil(coordinates.size())); +} + +BOOST_AUTO_TEST_CASE(random_full_width_coordinates_match_integer_oracle) { + const auto coordinates = random_coordinates(129); + const auto coefficients = make_coefficients(coordinates); + const converter_type converter(12); + check_conversion(converter, coefficients, coordinates, 129); + check_conversion(converter, coefficients, coordinates, 256); +} + +BOOST_AUTO_TEST_CASE(rounding_modes_preserve_exact_conversion_with_reused_context) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + const auto coordinates = random_coordinates(33); + const auto coefficients = make_coefficients(coordinates); + const converter_type original(12); + const auto baseline = check_conversion(original, coefficients, coordinates, 64); + for (int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + const converter_type current(12); + for (const auto *converter : {&original, ¤t}) { + const auto actual = check_conversion(*converter, coefficients, coordinates, 64); + BOOST_CHECK_EQUAL_COLLECTIONS(actual.begin(), actual.end(), baseline.begin(), baseline.end()); + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } + } +} + +BOOST_AUTO_TEST_CASE(reused_buffers_are_zero_padded_and_empty_inputs_work) { + const auto coordinates = random_coordinates(33); + const auto coefficients = make_coefficients(coordinates); + const converter_type converter(3); + constexpr std::size_t stride = 64; + std::vector output(3 * 12 * stride, sentinel); + for (std::size_t count : {std::size_t(33), std::size_t(1), std::size_t(0)}) { + converter.convert(output, std::span(coefficients).first(count), stride); + check_output(output, 3, std::span(coordinates).first(count), stride); + } + // A zero coefficient also overwrites every active cell of a previous nonzero conversion. + converter.convert(output, coefficients, stride); + const std::array zero {}; + converter.convert(output, make_coefficients(zero), stride); + check_output(output, 3, zero, stride); + converter.convert({}, {}, 0); +} + +BOOST_AUTO_TEST_CASE(invalid_dimensions_and_prime_counts_are_rejected_before_writes) { + for (std::size_t count : {std::size_t(0), std::size_t(13), std::numeric_limits::max()}) { + BOOST_CHECK_THROW(converter_type {count}, std::invalid_argument); + } + const converter_type converter(1); + const std::array coefficients = {value_type::one(), value_type::one()}; + std::vector output(24, sentinel); + BOOST_CHECK_THROW(converter.convert(std::span(output).first(23), coefficients, 2), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert(output, coefficients, 1), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert(output, {}, 1), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert(std::span(output).first(12), coefficients, 1), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert(output, {}, std::numeric_limits::max()), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert(output, {}, 0), std::invalid_argument); + BOOST_CHECK_THROW(converter.convert({}, coefficients, 2), std::invalid_argument); + BOOST_CHECK(std::all_of(output.begin(), output.end(), [](double value) { return value == sentinel; })); +} + +BOOST_AUTO_TEST_CASE(conversion_uses_existing_plan_sizes_and_input_prefixes) { + const auto plan = detail::make_bn254_small_prime_fft_params(35, 17, 19); + BOOST_REQUIRE(plan); + const converter_type converter(plan->crt.prime_count); + const auto coordinates = random_coordinates(35); + const auto coefficients = make_coefficients(coordinates); + std::vector output(plan->spectrum_bytes / sizeof(double), sentinel); + for (std::size_t count : {plan->left_size, plan->right_size}) { + converter.convert(output, std::span(coefficients).first(count), plan->transform_size); + check_output(output, plan->crt.prime_count, std::span(coordinates).first(count), + plan->transform_size); + } +} + +BOOST_AUTO_TEST_SUITE_END() From 0c2023ddafea309bd74940249ce00df96191009b Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Wed, 14 Oct 2026 18:33:28 -0400 Subject: [PATCH 04/11] exact CRT export --- .../math/detail/bn254_small_prime_fft.hpp | 164 +++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_crt.cpp | 314 ++++++++++++++++++ 3 files changed, 479 insertions(+) create mode 100644 libs/math/test/small_prime_fft_crt.cpp 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 f35aea544..34226de40 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 @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -38,6 +39,7 @@ #include #include #include +#include namespace nil::crypto3::math::detail { @@ -183,6 +185,168 @@ namespace nil::crypto3::math::detail { std::vector primes_; }; + /** + * Reconstruct signed integer coordinates from a prefix of the small-prime catalogue, then + * reduce them into BN254 Fp12. The supplied bound C must cover the original integer expression, + * 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. + * This operation has no constant-time guarantee. + */ + class bn254_small_prime_fft_reconstructor { + static constexpr std::size_t block_size = 32; + static constexpr std::size_t correction_count = std::bit_width(small_prime_fft_primes.size() - 1); + // Each CRT term is < M, so their sum is < prime_count*M < 12*2^600 < 2^604. + static constexpr unsigned crt_bits = 50 * small_prime_fft_primes.size() + correction_count; + using integer_type = boost::multiprecision::number>; + + struct prime_context { + small_prime_fft_arithmetic arithmetic; + integer_type cofactor; + double inverse; + + prime_context(std::uint64_t modulus, const integer_type &cofactor) : + arithmetic(modulus), cofactor(cofactor) { + 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()); + } + }; + + public: + using field_type = bn254_small_prime_fft_converter::field_type; + using value_type = field_type::value_type; + + bn254_small_prime_fft_reconstructor(const small_prime_fft_crt_params ¶ms, + const boost::multiprecision::cpp_int &coefficient_bound) { + if (params.prime_count == 0 || params.prime_count > small_prime_fft_primes.size() || + coefficient_bound < 0) { + throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid CRT parameters"); + } + boost::multiprecision::cpp_int product = 1; + for (std::size_t i = 0; i < params.prime_count; ++i) { + product *= small_prime_fft_primes[i].modulus; + } + if (params.modulus != product || 2 * coefficient_bound >= product) { + throw std::invalid_argument( + "bn254_small_prime_fft_reconstructor: inconsistent modulus or insufficient capacity"); + } + + // The checks above establish positivity and capacity before conversion to fixed precision. + using setup_backend = integer_type::backend_type::cpp_int_type; + multiples_[0].backend().from_cpp_int(setup_backend(product.backend())); + bound_.backend().from_cpp_int(setup_backend(coefficient_bound.backend())); + half_modulus_ = multiples_[0] >> 1; + for (std::size_t bit = 1; bit < correction_count; ++bit) { + multiples_[bit] = multiples_[bit - 1] << 1; + } + primes_.reserve(params.prime_count); + 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); + } + } + + std::size_t prime_count() const { + return primes_.size(); + } + + /** + * Read residues[(prime*12 + coordinate)*stride + coefficient] for each requested output. + * Require output.size() <= stride and residues.size() == prime_count()*12*stride. + * Check dimensions without overflow before writing any output; throw invalid_argument on + * mismatch. An empty output does no work; entries after output.size() in each channel are + * not read. Input residues must already include any inverse-transform normalization. + * + * @pre Each consumed residue is an integer-valued double of magnitude < 2^53. The thread's + * rounding mode must remain unchanged during the call; kernel selection occurs once. + * + * Throw overflow_error if a centered coordinate exceeds the declared bound; earlier output + * coordinates may already have been written. This check cannot detect an incorrect original + * 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 { + const std::size_t channels = field_type::arity * primes_.size(); + if (channels == 0 || output.size() > stride || stride > residues.size() / channels || + residues.size() != channels * stride) { + throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid buffer dimensions"); + } + if (output.empty()) { + return; + } + const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { + return prime.arithmetic.fma_supported(); + }); + if (use_fma) { + reconstruct_impl(output, residues, stride); + } else { + reconstruct_impl(output, residues, stride); + } + } + + private: + template + void reconstruct_impl(std::span output, std::span residues, + std::size_t stride) const { + using base_field_type = algebra::fields::alt_bn128_base_field<254>; + 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 {}; + 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; + for (std::size_t i = 0; i < count; ++i) { + double digit; + if constexpr (UseFma) { + digit = arithmetic.multiply_fma(arithmetic.reduce_fma(source[i]), context.inverse); + } else { + digit = + arithmetic.multiply_integer(arithmetic.reduce_integer(source[i]), 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); + sums[i] += context.cofactor * 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]; + } + } + 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 base_field_type::value_type value(magnitude); + output[offset + i].coordinate(coordinate) = negative ? -value : value; + } + offset += count; + } + } + } + + std::vector primes_; + std::array multiples_; + integer_type half_modulus_; + integer_type bound_; + }; + } // namespace nil::crypto3::math::detail #endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index bfab84111..399febec9 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -49,6 +49,7 @@ set(TESTS_NAMES "small_prime_fft_params" "small_prime_fft_arithmetic" "small_prime_fft_conversion" + "small_prime_fft_crt" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_crt.cpp b/libs/math/test/small_prime_fft_crt.cpp new file mode 100644 index 000000000..d2b4f4d09 --- /dev/null +++ b/libs/math/test/small_prime_fft_crt.cpp @@ -0,0 +1,314 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_crt_test + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + namespace detail = nil::crypto3::math::detail; + using boost::multiprecision::cpp_int; + using fq = nil::crypto3::algebra::fields::alt_bn128_base_field<254>; + using reconstructor_type = detail::bn254_small_prime_fft_reconstructor; + using value_type = reconstructor_type::value_type; + using coordinates_type = std::array; + + struct restore_rounding { + int saved = std::fegetround(); + ~restore_rounding() { + std::fesetround(saved); + } + }; + + cpp_int canonical(const cpp_int &value, const cpp_int &modulus) { + cpp_int result = value % modulus; + if (result < 0) { + result += modulus; + } + return result; + } + + detail::small_prime_fft_crt_params make_params(std::size_t prime_count) { + cpp_int product = 1; + for (std::size_t i = 0; i < prime_count; ++i) { + product *= detail::small_prime_fft_primes[i].modulus; + } + return {prime_count, product}; + } + + std::vector make_residues(std::span coordinates, std::size_t prime_count, + std::size_t stride, bool lazy = 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; + 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); + if (lazy) { + // Reach large positive and negative representatives while staying below 2^53. + residue += ((prime + coordinate + i) % 2 == 0 ? 7 : -8) * modulus; + } + result[(prime * 12 + coordinate) * stride + i] = + static_cast(residue.convert_to()); + } + } + } + return result; + } + + void check_output(std::span output, std::span coordinates) { + BOOST_REQUIRE_EQUAL(output.size(), coordinates.size()); + const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); + for (std::size_t i = 0; i < output.size(); ++i) { + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + const std::size_t coordinate = 6 * w + 2 * v + u; + BOOST_TEST_CONTEXT("coefficient=" << i << " coordinate=" << coordinate) { + // Inspect the tower directly, independently of coordinate() in reconstruction. + 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)); + } + } + } + } + } + } + + void check_reconstruction(const reconstructor_type &reconstructor, std::span coordinates, + std::size_t stride, bool lazy = false) { + const auto residues = make_residues(coordinates, reconstructor.prime_count(), stride, lazy); + const auto sentinel = value_type::one(); + std::vector storage(coordinates.size() + 2, sentinel); + auto output = std::span(storage).subspan(1, coordinates.size()); + reconstructor.reconstruct(output, residues, 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) { + std::mt19937_64 random(0xc47); + std::vector result(count); + for (auto &coefficient : result) { + for (auto &coordinate : coefficient) { + for (unsigned limb = 0; limb < 10; ++limb) { + coordinate <<= 64; + coordinate += random(); + } + coordinate = coordinate % (2 * bound + 1) - bound; + } + } + return result; + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_crt_test_suite) + +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) { + BOOST_TEST_CONTEXT("prime count=" << count) { + const auto params = make_params(count); + const cpp_int bound = params.modulus / 2; + const reconstructor_type reconstructor(params, bound); + BOOST_CHECK_EQUAL(reconstructor.prime_count(), count); + std::vector values = {0, 1, -1, bound, bound - 1, -bound, 1 - bound}; + for (cpp_int value : {cpp_int(characteristic - 1), characteristic, cpp_int(characteristic + 1)}) { + if (value <= bound) { + values.push_back(value); + values.push_back(-value); + } + } + for (unsigned bit = 1; bit <= boost::multiprecision::msb(bound); ++bit) { + if (bit % 64 != 0 && bit % 50 != 0) { + continue; + } + for (int delta : {-1, 0, 1}) { + const cpp_int value = (cpp_int(1) << bit) + delta; + if (value <= bound) { + values.push_back(value); + values.push_back(-value); + } + } + } + // All CRT digits q_i-1 produce the largest cofactor sum, just below count*M. + // The associated signed integer is -sum(M/q_i), which exercises every correction bit. + cpp_int largest_sum_value = 0; + for (std::size_t prime = 0; prime < count; ++prime) { + largest_sum_value -= params.modulus / detail::small_prime_fft_primes[prime].modulus; + } + values.push_back(largest_sum_value); + std::vector coordinates(values.size()); + for (std::size_t i = 0; i < coordinates.size(); ++i) { + for (std::size_t c = 0; c < 12; ++c) { + coordinates[i][c] = values[(i + c) % values.size()]; + } + } + check_reconstruction(reconstructor, coordinates, coordinates.size() + 3); + check_reconstruction(reconstructor, coordinates, coordinates.size(), true); + } + } +} + +BOOST_AUTO_TEST_CASE(random_signed_integers_cover_all_prime_counts_and_partial_blocks) { + 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); + const auto coordinates = random_coordinates(65, bound); + check_reconstruction(reconstructor, coordinates, 65); + check_reconstruction(reconstructor, coordinates, 96, true); + } +} + +BOOST_AUTO_TEST_CASE(workload_plan_bound_is_recovered_with_both_signs) { + const auto plan = detail::make_bn254_small_prime_fft_params(14939, 14939); + BOOST_REQUIRE(plan); + const reconstructor_type reconstructor(plan->crt, plan->coefficient_bound); + auto coordinates = random_coordinates(33, plan->coefficient_bound); + coordinates[0][0] = plan->coefficient_bound; + coordinates[0][1] = -plan->coefficient_bound; + check_reconstruction(reconstructor, coordinates, 64, true); +} + +BOOST_AUTO_TEST_CASE(conversion_round_trip_preserves_fp12_coordinates_and_requested_prefixes) { + const auto plan = detail::make_bn254_small_prime_fft_params(65, 33); + BOOST_REQUIRE(plan); + const cpp_int characteristic(fq::modulus.backend().to_cpp_int()); + auto coordinates = random_coordinates(65, characteristic - 1); + std::vector input(coordinates.size(), value_type::zero()); + for (std::size_t i = 0; i < input.size(); ++i) { + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + auto &coordinate = coordinates[i][6 * w + 2 * v + u]; + coordinate = canonical(coordinate, characteristic); + input[i].data[w].data[v].data[u] = fq::value_type(fq::integral_type(coordinate.str())); + } + } + } + } + const detail::bn254_small_prime_fft_converter converter(plan->crt.prime_count); + const reconstructor_type reconstructor(plan->crt, plan->coefficient_bound); + std::vector residues(plan->spectrum_bytes / sizeof(double)); + converter.convert(residues, input, plan->transform_size); + std::vector output(input.size(), value_type::one()); + for (std::size_t count : {std::size_t(65), std::size_t(33), std::size_t(1), std::size_t(0)}) { + reconstructor.reconstruct(std::span(output).first(count), residues, plan->transform_size); + check_output(std::span(output).first(count), + std::span(coordinates).first(count)); + } + converter.convert(residues, {}, plan->transform_size); + reconstructor.reconstruct(output, residues, plan->transform_size); + BOOST_CHECK(std::all_of(output.begin(), output.end(), [](const auto &value) { return value.is_zero(); })); + reconstructor.reconstruct({}, {}, 0); +} + +BOOST_AUTO_TEST_CASE(rounding_modes_preserve_reconstruction_with_reused_and_fresh_contexts) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + const auto params = make_params(12); + const cpp_int bound = params.modulus / 2; + const auto coordinates = random_coordinates(33, bound); + const reconstructor_type original(params, bound); + for (int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + const reconstructor_type current(params, bound); + for (const auto *reconstructor : {&original, ¤t}) { + check_reconstruction(*reconstructor, coordinates, 64, true); + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } + } +} + +BOOST_AUTO_TEST_CASE(inconsistent_parameters_and_insufficient_capacity_are_rejected) { + const auto valid = make_params(1); + for (std::size_t count : {std::size_t(0), std::size_t(13), std::numeric_limits::max()}) { + const detail::small_prime_fft_crt_params invalid {count, valid.modulus}; + BOOST_CHECK_THROW(reconstructor_type(invalid, 0), std::invalid_argument); + } + for (cpp_int modulus : {cpp_int(0), cpp_int(-1), cpp_int(valid.modulus - 1), cpp_int(valid.modulus + 1)}) { + const detail::small_prime_fft_crt_params invalid {1, modulus}; + BOOST_CHECK_THROW(reconstructor_type(invalid, 0), std::invalid_argument); + } + BOOST_CHECK_THROW(reconstructor_type(valid, -1), std::invalid_argument); + BOOST_CHECK_THROW(reconstructor_type(valid, (valid.modulus + 1) / 2), std::invalid_argument); + BOOST_CHECK_THROW(reconstructor_type(valid, cpp_int(1) << 10000), std::invalid_argument); + const reconstructor_type zero_bound(valid, 0); + const std::array zero {}; + check_reconstruction(zero_bound, zero, 1); +} + +BOOST_AUTO_TEST_CASE(centered_coordinates_outside_the_declared_bound_are_rejected) { + const reconstructor_type reconstructor(make_params(1), 10); + std::array coordinates {}; + for (int value : {-10, 10}) { + coordinates[0][0] = value; + check_reconstruction(reconstructor, coordinates, 1); + } + for (int value : {-11, 11}) { + coordinates[0][0] = value; + const auto residues = make_residues(coordinates, 1, 1); + std::array output {value_type::one()}; + BOOST_CHECK_THROW(reconstructor.reconstruct(output, residues, 1), std::overflow_error); + } +} + +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(); })); +} + +BOOST_AUTO_TEST_SUITE_END() From 4e932e0c365a9d2bfc24da71dfc120a87f27fef6 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Wed, 14 Oct 2026 19:16:34 -0400 Subject: [PATCH 05/11] complete scalar transform --- .../math/detail/bn254_small_prime_fft.hpp | 277 +++++++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_transform.cpp | 385 ++++++++++++++++++ 3 files changed, 663 insertions(+) create mode 100644 libs/math/test/small_prime_fft_transform.cpp 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 34226de40..cc9642bc5 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 @@ -28,6 +28,8 @@ #include #include #include +#include +#include #include #include #include @@ -347,6 +349,281 @@ namespace nil::crypto3::math::detail { integer_type bound_; }; + /** + * In-place, power-of-two transforms for one prime from the catalogue. The forward transform + * maps natural-order coefficients to evaluations at omega^k in bit-reversed index order. + * The inverse consumes that same order, so convolution needs no permutation passes. + * + * A reusable table of max(1,N/2) centered twiddles serves both directions. Radix-4 recursion + * fuses pairs of layers, with radix-2 leaves for odd depths and special cases for unit twiddles. + * Transform calls allocate no buffers; recursion uses O(log N) stack space. + * + * @pre Inputs are integer-valued doubles in (-3*q,3*q). The calling thread's rounding mode + * must remain unchanged during a call. Each call selects FMA or integer arithmetic once. + * The forward output is in (-3*q,3*q); the normalized inverse output is in [0,q). + * This operation has no constant-time guarantee. + */ + class small_prime_fft_plan { + public: + small_prime_fft_plan(std::size_t prime_index, std::size_t size) : + arithmetic_(checked_prime(prime_index, size).modulus), size_(size), + twiddles_(std::max(std::size_t(1), size / 2), 1.0) { + using wide_type = boost::multiprecision::uint128_modular_t; + const auto &prime = small_prime_fft_primes[prime_index]; + wide_type root = prime.root; + const unsigned depth = std::countr_zero(size_); + for (unsigned level = prime.two_adicity; level > depth; --level) { + root = root * root % prime.modulus; + } + omega_ = root.convert_to(); + size_inverse_ = centered( + static_cast(boost::multiprecision::inverse_mod(wide_type(size_), wide_type(prime.modulus)) + .convert_to())); + + // At a layer with 2^d blocks, block j uses omega_(2^(d+1))^bitreverse(j,d). + // Each new table segment is the preceding prefix times a primitive root of twice + // the previous order: [1], [i], [w8,i*w8], [w16,i*w16,w8*w16,i*w8*w16], ... . + for (std::size_t half = size_ / 4; half != 0; half >>= 1) { + twiddles_[half] = centered(static_cast(root.convert_to())); + root = root * root % prime.modulus; + } + for (std::size_t half = 1; half < size_ / 2; half <<= 1) { + const double step = twiddles_[half]; + for (std::size_t i = 1; i < half; ++i) { + twiddles_[half + i] = centered(arithmetic_.multiply_integer(twiddles_[i], step)); + } + } + } + + std::size_t size() const { + return size_; + } + + std::uint64_t modulus() const { + return arithmetic_.modulus(); + } + + std::uint64_t omega() const { + return omega_; + } + + /** Exact centered representative of 1/N modulo q, also usable to fuse scaling into pointwise work. */ + double size_inverse() const { + return size_inverse_; + } + + /** Require exactly N entries; invalid dimensions throw before any writes. N=1 is the identity. */ + void fft(std::span values) const { + check_buffer(values); + if (arithmetic_.fma_supported()) { + forward_recursive(values.data(), size_, 0); + } else { + forward_recursive(values.data(), size_, 0); + } + } + + /** + * Require exactly N entries in the forward spectrum order. By default, scale by 1/N and + * return canonical residues. With normalize=false, return N times the inverse in (-2*q,2*q), + * allowing the caller to fuse the scaling elsewhere. Invalid dimensions throw before writes. + */ + void inverse_fft(std::span values, bool normalize = true) const { + check_buffer(values); + if (arithmetic_.fma_supported()) { + inverse_impl(values, normalize); + } else { + inverse_impl(values, normalize); + } + } + + private: + static const small_prime_fft_prime &checked_prime(std::size_t prime_index, std::size_t size) { + if (prime_index >= small_prime_fft_primes.size() || !std::has_single_bit(size) || + std::countr_zero(size) > small_prime_fft_primes[prime_index].two_adicity) { + throw std::invalid_argument("small_prime_fft_plan: invalid prime index or transform size"); + } + return small_prime_fft_primes[prime_index]; + } + + void check_buffer(std::span values) const { + if (values.size() != size_ || twiddles_.empty()) { + throw std::invalid_argument("small_prime_fft_plan: invalid buffer dimensions or moved-from plan"); + } +#ifndef NDEBUG + const double limit = 3 * static_cast(modulus()); + for (double value : values) { + assert(std::isfinite(value) && std::trunc(value) == value && std::abs(value) < limit); + } +#endif + } + + double centered(double value) const { + const double q = static_cast(modulus()); + return value > q / 2 ? value - q : (value < -q / 2 ? value + q : value); + } + + double inverse_twiddle(std::size_t index) const { + if (index == 0) { + return 1.0; + } + // Within [b,2*b), reflected exponents add to half the root order, hence w^-1=-w_reflected. + const std::size_t base = std::bit_floor(index); + return -twiddles_[base + (2 * base - 1 - index)]; + } + + // Keep arithmetic selection outside the butterfly loops, with compile-time dispatch inside them. + template + double reduce(double value) const { + if constexpr (UseFma) { + return arithmetic_.reduce_fma(value); + } else { + return arithmetic_.reduce_integer(value); + } + } + + template + double multiply(double left, double right) const { + if constexpr (UseFma) { + return arithmetic_.multiply_fma(left, right); + } else { + return arithmetic_.multiply_integer(left, right); + } + } + + template + void forward_radix4(double *values, std::size_t quarter, std::size_t node) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + const double w2 = twiddles_[node]; + const double w = twiddles_[2 * node]; + const double iw = twiddles_[2 * node + 1]; + for (std::size_t i = 0; i < quarter; ++i) { + const double a0 = reduce(values[i]); + const double a1 = values[i + quarter]; + const double a2 = UnitTwiddle ? reduce(values[i + 2 * quarter]) : + multiply(values[i + 2 * quarter], w2); + const double a3 = UnitTwiddle ? reduce(values[i + 3 * quarter]) : + multiply(values[i + 3 * quarter], w2); + const double even_sum = a0 + a2; + const double even_difference = a0 - a2; + const double odd_sum = UnitTwiddle ? reduce(a1 + a3) : multiply(a1 + a3, w); + const double odd_difference = multiply(a1 - a3, iw); + // Input magnitudes < 3*q and weights < q/2 imply |a2|,|a3| < q and + // |a1 +/- a3| < 4*q. Every product is < 2*q^2, so the odd terms are < q. + // Thus all outputs are < 3*q, with every integer sum exact below 2^53. + values[i] = even_sum + odd_sum; + values[i + quarter] = even_sum - odd_sum; + values[i + 2 * quarter] = even_difference + odd_difference; + values[i + 3 * quarter] = even_difference - odd_difference; + } + } + + template + void forward_recursive(double *values, std::size_t length, std::size_t node) const { + if (length == 1) { + return; + } + if (length == 2) { + const double left = reduce(values[0]); + const double right = + node == 0 ? reduce(values[1]) : multiply(values[1], twiddles_[node]); + values[0] = left + right; + values[1] = left - right; + return; + } + const std::size_t quarter = length / 4; + if (node == 0) { + forward_radix4(values, quarter, node); + } else { + forward_radix4(values, quarter, node); + } + if (quarter > 1) { + for (std::size_t child = 0; child < 4; ++child) { + forward_recursive(values + child * quarter, quarter, 4 * node + child); + } + } + } + + template + void inverse_radix4(double *values, std::size_t quarter, std::size_t node) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + const double w2 = inverse_twiddle(node); + const double w = inverse_twiddle(2 * node); + const double iw = inverse_twiddle(2 * node + 1); + for (std::size_t i = 0; i < quarter; ++i) { + // Leaf inputs may be forward residues of magnitude < 3*q. Reduce them before + // adding four entries; non-leaf inputs already satisfy the inverse bound < 2*q. + const double b0 = quarter == 1 ? reduce(values[i]) : values[i]; + const double b1 = quarter == 1 ? reduce(values[i + quarter]) : values[i + quarter]; + const double b2 = quarter == 1 ? reduce(values[i + 2 * quarter]) : values[i + 2 * quarter]; + const double b3 = quarter == 1 ? reduce(values[i + 3 * quarter]) : values[i + 3 * quarter]; + const double even_sum = b0 + b1; + const double odd_sum = b2 + b3; + const double left = UnitTwiddle ? reduce(b0 - b1) : multiply(b0 - b1, w); + const double right = multiply(b2 - b3, iw); + // Four inputs < 2*q give sums < 8*q < 2^53. Times a centered weight this is + // < 4*q^2, so the a2 remainder is < 1.5*q. Left and right are each < q; + // their sum is < 2*q and their weighted difference is < q. The invariant closes. + values[i] = reduce(even_sum + odd_sum); + values[i + quarter] = left + right; + values[i + 2 * quarter] = + UnitTwiddle ? reduce(even_sum - odd_sum) : multiply(even_sum - odd_sum, w2); + values[i + 3 * quarter] = + UnitTwiddle ? reduce(left - right) : multiply(left - right, w2); + } + } + + template + void inverse_recursive(double *values, std::size_t length, std::size_t node) const { + if (length == 1) { + values[0] = reduce(values[0]); + return; + } + if (length == 2) { + const double left = reduce(values[0]); + const double right = reduce(values[1]); + values[0] = reduce(left + right); + values[1] = + node == 0 ? reduce(left - right) : multiply(left - right, inverse_twiddle(node)); + return; + } + const std::size_t quarter = length / 4; + if (quarter > 1) { + for (std::size_t child = 0; child < 4; ++child) { + inverse_recursive(values + child * quarter, quarter, 4 * node + child); + } + } + if (node == 0) { + inverse_radix4(values, quarter, node); + } else { + inverse_radix4(values, quarter, node); + } + } + + template + void inverse_impl(std::span values, bool normalize) const { + inverse_recursive(values.data(), size_, 0); + if (normalize) { + const double q = static_cast(modulus()); + for (double &value : values) { + value = multiply(value, size_inverse_); + value = value < 0 ? value + q : value; + } + } + } + + small_prime_fft_arithmetic arithmetic_; + std::size_t size_; + std::uint64_t omega_; + double size_inverse_; + std::vector twiddles_; + }; + } // namespace nil::crypto3::math::detail #endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index 399febec9..ce5ce58fb 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -50,6 +50,7 @@ set(TESTS_NAMES "small_prime_fft_arithmetic" "small_prime_fft_conversion" "small_prime_fft_crt" + "small_prime_fft_transform" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_transform.cpp b/libs/math/test/small_prime_fft_transform.cpp new file mode 100644 index 000000000..bbc540605 --- /dev/null +++ b/libs/math/test/small_prime_fft_transform.cpp @@ -0,0 +1,385 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_transform_test + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + namespace detail = nil::crypto3::math::detail; + using boost::multiprecision::cpp_int; + using plan_type = detail::small_prime_fft_plan; + + struct restore_rounding { + int saved = std::fegetround(); + ~restore_rounding() { + std::fesetround(saved); + } + }; + + cpp_int canonical(const cpp_int &value, const cpp_int &modulus) { + cpp_int result = value % modulus; + return result < 0 ? result + modulus : result; + } + + cpp_int as_integer(double value) { + BOOST_REQUIRE(std::isfinite(value)); + BOOST_REQUIRE(std::abs(value) < 0x1p53); + BOOST_REQUIRE_EQUAL(std::trunc(value), value); + return static_cast(value); + } + + cpp_int root_for(std::size_t prime, std::size_t size) { + const auto &entry = detail::small_prime_fft_primes[prime]; + const cpp_int exponent = cpp_int(1) << (entry.two_adicity - std::countr_zero(size)); + return boost::multiprecision::powm(cpp_int(entry.root), exponent, cpp_int(entry.modulus)); + } + + // Direct O(N^2) integer DFT, independent of the butterfly schedule and floating-point kernels. + std::vector dft(std::span input, std::size_t prime, bool inverse = false) { + const cpp_int q = detail::small_prime_fft_primes[prime].modulus; + cpp_int root = root_for(prime, input.size()); + cpp_int scale = 1; + if (inverse) { + root = boost::multiprecision::powm(root, q - 2, q); + scale = boost::multiprecision::powm(cpp_int(input.size()), q - 2, q); + } + std::vector result(input.size()); + cpp_int step = 1; + for (std::size_t k = 0; k < input.size(); ++k) { + cpp_int sum = 0; + cpp_int power = 1; + for (double value : input) { + sum += static_cast(value) * power; + power = power * step % q; + } + result[k] = canonical(sum * scale, q); + step = step * root % q; + } + return result; + } + + std::vector lazy_input(std::size_t size, std::uint64_t modulus, std::mt19937_64 &random) { + const auto q = static_cast(modulus); + const std::array edges {3 * q - 1, 1 - 3 * q, 0, 1, -1, q - 1, q, + q + 1, -q - 1, -q, -q + 1, 2 * q, -2 * q}; + std::vector result(size); + for (std::size_t i = 0; i < size; ++i) { + const auto value = i % 2 == 0 ? edges[(i / 2) % edges.size()] : + static_cast(random() % (6 * q - 1)) - (3 * q - 1); + result[i] = static_cast(value); + } + return result; + } + + std::vector bit_reversed(std::span input) { + std::vector result(input.size()); + const unsigned depth = std::countr_zero(input.size()); + for (std::size_t i = 0; i < input.size(); ++i) { + result[detail::bitreverse(i, depth)] = input[i]; + } + return result; + } + + void check_residue(double actual, const cpp_int &expected, std::uint64_t modulus, unsigned bound, + bool normalized = false) { + const cpp_int value = as_integer(actual); + BOOST_CHECK(std::abs(actual) < bound * static_cast(modulus)); + if (normalized) { + BOOST_CHECK(actual >= 0); + } + BOOST_CHECK_EQUAL(canonical(value, cpp_int(modulus)), expected); + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_transform_tests) + +BOOST_AUTO_TEST_CASE(forward_matches_integer_dft_for_every_prime) { + std::mt19937_64 random(0x46465401); + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + for (std::size_t size = 1; size <= 64; size *= 2) { + BOOST_TEST_CONTEXT("prime=" << prime << " N=" << size) { + const plan_type plan(prime, size); + const auto q = detail::small_prime_fft_primes[prime].modulus; + BOOST_CHECK_EQUAL(plan.size(), size); + BOOST_CHECK_EQUAL(plan.modulus(), q); + BOOST_CHECK_EQUAL(plan.omega(), root_for(prime, size)); + BOOST_CHECK(std::abs(plan.size_inverse()) < static_cast(q) / 2); + BOOST_CHECK_EQUAL(canonical(as_integer(plan.size_inverse()) * size, cpp_int(q)), 1); + + auto values = lazy_input(size, q, random); + const auto expected = dft(values, prime); + plan.fft(values); + for (std::size_t k = 0; k < size; ++k) { + check_residue(values[detail::bitreverse(k, std::countr_zero(size))], expected[k], q, 3); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(inverse_matches_independent_spectra) { + std::mt19937_64 random(0x46465402); + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + for (std::size_t size = 1; size <= 64; size *= 2) { + BOOST_TEST_CONTEXT("prime=" << prime << " N=" << size) { + const plan_type plan(prime, size); + const auto q = plan.modulus(); + const auto spectrum = lazy_input(size, q, random); + const auto expected = dft(spectrum, prime, true); + auto normalized = bit_reversed(spectrum); + auto unnormalized = normalized; + plan.inverse_fft(normalized); + plan.inverse_fft(unnormalized, false); + for (std::size_t i = 0; i < size; ++i) { + check_residue(normalized[i], expected[i], q, 1, true); + check_residue(unnormalized[i], canonical(expected[i] * size, cpp_int(q)), q, 2); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(basis_vectors_pin_spectrum_order) { + constexpr std::size_t size = 16; + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + const plan_type plan(prime, size); + const cpp_int q = plan.modulus(); + const cpp_int root = root_for(prime, size); + for (std::size_t coefficient = 0; coefficient < size; ++coefficient) { + std::vector values(size, 0.0); + values[coefficient] = 1.0; + plan.fft(values); + for (std::size_t k = 0; k < size; ++k) { + const cpp_int expected = boost::multiprecision::powm(root, cpp_int(coefficient * k), q); + check_residue(values[detail::bitreverse(k, 4)], expected, plan.modulus(), 3); + } + plan.inverse_fft(values); + for (std::size_t i = 0; i < size; ++i) { + BOOST_CHECK_EQUAL(values[i], i == coefficient ? 1.0 : 0.0); + } + } + } +} + +BOOST_AUTO_TEST_CASE(constant_and_alternating_extremes) { + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + const auto q = detail::small_prime_fft_primes[prime].modulus; + const auto signed_q = static_cast(q); + for (std::size_t size = 1; size <= 128; size *= 2) { + const plan_type plan(prime, size); + for (const std::int64_t value : + {std::int64_t(0), std::int64_t(1), signed_q, -signed_q, 3 * signed_q - 1, 1 - 3 * signed_q}) { + for (bool alternating : {false, true}) { + std::vector values(size); + for (std::size_t i = 0; i < size; ++i) { + values[i] = static_cast(alternating && i % 2 ? -value : value); + } + plan.fft(values); + for (std::size_t k = 0; k < size; ++k) { + const std::size_t nonzero_index = alternating ? size / 2 : 0; + const cpp_int expected = k == nonzero_index ? canonical(cpp_int(value) * size, cpp_int(q)) : 0; + check_residue(values[detail::bitreverse(k, std::countr_zero(size))], expected, q, 3); + } + plan.inverse_fft(values); + for (std::size_t i = 0; i < size; ++i) { + const cpp_int expected = canonical(cpp_int(alternating && i % 2 ? -value : value), cpp_int(q)); + check_residue(values[i], expected, q, 1, true); + } + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(round_trips_at_larger_even_and_odd_depths) { + std::mt19937_64 random(0x46465403); + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + for (const std::size_t size : {128, 256, 512, 1024, 16384, 32768}) { + BOOST_TEST_CONTEXT("prime=" << prime << " N=" << size) { + const plan_type plan(prime, size); + const auto input = lazy_input(size, plan.modulus(), random); + auto values = input; + plan.fft(values); + for (double value : values) { + as_integer(value); + BOOST_CHECK(std::abs(value) < 3 * static_cast(plan.modulus())); + } + plan.inverse_fft(values); + for (std::size_t i = 0; i < size; ++i) { + check_residue(values[i], canonical(as_integer(input[i]), cpp_int(plan.modulus())), plan.modulus(), + 1, true); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(convolution_matches_schoolbook_with_both_scaling_paths) { + std::mt19937_64 random(0x46465404); + const std::array, 7> shapes { + {{1, 1}, {1, 17}, {2, 3}, {7, 9}, {17, 33}, {32, 32}, {65, 31}}}; + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + for (const auto [left_size, right_size] : shapes) { + BOOST_TEST_CONTEXT("prime=" << prime << " shape=" << left_size << "x" << right_size) { + const auto q = detail::small_prime_fft_primes[prime].modulus; + const auto left = lazy_input(left_size, q, random); + const auto right = lazy_input(right_size, q, random); + const std::size_t size = std::bit_ceil(left_size + right_size - 1); + std::vector expected(size); + for (std::size_t i = 0; i < left_size; ++i) { + for (std::size_t j = 0; j < right_size; ++j) { + expected[i + j] += as_integer(left[i]) * as_integer(right[j]); + } + } + for (auto &value : expected) { + value = canonical(value, cpp_int(q)); + } + + const plan_type plan(prime, size); + const detail::small_prime_fft_arithmetic arithmetic(q); + std::vector a(size, 0.0), b(size, 0.0); + std::copy(left.begin(), left.end(), a.begin()); + std::copy(right.begin(), right.end(), b.begin()); + plan.fft(a); + plan.fft(b); + for (std::size_t i = 0; i < size; ++i) { + a[i] = arithmetic.multiply_integer(a[i], b[i]); + b[i] = arithmetic.multiply_integer(a[i], plan.size_inverse()); + } + auto unnormalized = a; + plan.inverse_fft(a); + plan.inverse_fft(b, false); + plan.inverse_fft(unnormalized, false); + for (std::size_t i = 0; i < size; ++i) { + check_residue(a[i], expected[i], q, 1, true); + check_residue(b[i], expected[i], q, 2); + check_residue(unnormalized[i], canonical(expected[i] * size, cpp_int(q)), q, 2); + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(rounding_modes_preserve_exactness_and_environment) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + std::mt19937_64 random(0x46465405); + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + for (const std::size_t size : {32, 64}) { + const plan_type original(prime, size); + const auto q = original.modulus(); + const auto input = lazy_input(size, q, random); + const auto expected = dft(input, prime); + for (const int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_TEST_CONTEXT("prime=" << prime << " N=" << size << " rounding=" << mode) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + const plan_type constructed_in_mode(prime, size); + for (const plan_type *plan : {&original, &constructed_in_mode}) { + auto values = input; + plan->fft(values); + for (std::size_t k = 0; k < size; ++k) { + check_residue(values[detail::bitreverse(k, std::countr_zero(size))], expected[k], q, 3); + } + plan->inverse_fft(values); + for (std::size_t i = 0; i < size; ++i) { + check_residue(values[i], canonical(as_integer(input[i]), cpp_int(q)), q, 1, true); + } + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } + } + } + } + } +} + +BOOST_AUTO_TEST_CASE(plan_reuse_and_in_place_buffer_boundaries) { + constexpr std::size_t size = 32; + const plan_type original(0, size); + auto copy = original; + const plan_type moved = std::move(copy); + std::mt19937_64 random(0x46465406); + for (const plan_type *plan : {&original, &moved}) { + for (unsigned trial = 0; trial < 3; ++trial) { + const auto input = lazy_input(size, plan->modulus(), random); + std::vector storage(size + 2, -17.0); + auto values = std::span(storage).subspan(1, size); + std::copy(input.begin(), input.end(), values.begin()); + const double *address = values.data(); + plan->fft(values); + plan->inverse_fft(values); + BOOST_CHECK(values.data() == address); + BOOST_CHECK_EQUAL(storage.front(), -17.0); + BOOST_CHECK_EQUAL(storage.back(), -17.0); + for (std::size_t i = 0; i < size; ++i) { + check_residue(values[i], canonical(as_integer(input[i]), cpp_int(plan->modulus())), plan->modulus(), 1, + true); + } + } + } +} + +BOOST_AUTO_TEST_CASE(invalid_parameters_and_dimensions_throw_before_writes) { + BOOST_CHECK_THROW(plan_type(detail::small_prime_fft_primes.size(), 1), std::invalid_argument); + BOOST_CHECK_THROW(plan_type(std::numeric_limits::max(), 1), std::invalid_argument); + for (const std::size_t size : + {std::size_t(0), std::size_t(3), std::size_t(6), std::numeric_limits::max()}) { + BOOST_CHECK_THROW(plan_type(0, size), std::invalid_argument); + } + for (std::size_t prime = 0; prime < detail::small_prime_fft_primes.size(); ++prime) { + const unsigned depth = detail::small_prime_fft_primes[prime].two_adicity + 1; + if (depth < std::numeric_limits::digits) { + BOOST_CHECK_THROW(plan_type(prime, std::size_t(1) << depth), std::invalid_argument); + } + } + const plan_type plan(0, 8); + for (const std::size_t size : {0, 7, 9}) { + std::vector values(size, 42.0); + const auto before = values; + BOOST_CHECK_THROW(plan.fft(values), std::invalid_argument); + BOOST_CHECK_THROW(plan.inverse_fft(values), std::invalid_argument); + BOOST_CHECK_THROW(plan.inverse_fft(values, false), std::invalid_argument); + BOOST_CHECK_EQUAL_COLLECTIONS(values.begin(), values.end(), before.begin(), before.end()); + } +} + +BOOST_AUTO_TEST_SUITE_END() From 24ccb35cd6417b8649cba5b556ccc6beec3aa7dc Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 13:12:37 -0400 Subject: [PATCH 06/11] complete fp12 backend --- .../math/detail/bn254_small_prime_fft.hpp | 194 +++++++ .../bn254_small_prime_fft_backend.hpp | 208 ++++++++ libs/math/test/CMakeLists.txt | 1 + libs/math/test/small_prime_fft_backend.cpp | 483 ++++++++++++++++++ 4 files changed, 886 insertions(+) create mode 100644 libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp create mode 100644 libs/math/test/small_prime_fft_backend.cpp 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 cc9642bc5..edc36c245 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 @@ -624,6 +624,200 @@ namespace nil::crypto3::math::detail { std::vector twiddles_; }; + /** + * Pointwise products in the BN254 tower u^2=-1, v^3=9+u, w^2=v, over one small modulus q. + * This is a quotient ring; irreducibility of the tower polynomials modulo q is not required. + * Coordinate channels have the same order as the converter, with one prime per call. + * + * Reuse the field library's Karatsuba tower identities and quadratic-extension squaring + * identity with the scalar residue kernels. The existing field elements bind their base + * arithmetic to the target modulus and cannot represent these small-prime intermediates. + * + * @pre Input cells are integer-valued doubles of magnitude < 2^53. scale is an integer with + * |scale| < q/2, normally the FFT plan's centered 1/N. The two buffers are disjoint or identical. + * The thread's rounding mode stays unchanged during the call. Output cells are in (-q,q). + * No buffers are allocated and no constant-time guarantee is provided. + */ + class bn254_small_prime_fft_pointwise { + using quadratic_type = std::array; + using cubic_type = std::array; + using tower_type = std::array; + + template + struct kernel { + const small_prime_fft_arithmetic &arithmetic; + + double reduce(double value) const { + if constexpr (UseFma) { + return arithmetic.reduce_fma(value); + } else { + return arithmetic.reduce_integer(value); + } + } + + double multiply(double a, double b) const { + if constexpr (UseFma) { + return arithmetic.multiply_fma(a, b); + } else { + return arithmetic.multiply_integer(a, b); + } + } + + double add(double a, double b) const { + return reduce(a + b); + } + + double subtract(double a, double b) const { + return reduce(a - b); + } + + template + std::array add(const std::array &a, const std::array &b) const { + std::array result; + for (std::size_t i = 0; i < N; ++i) { + result[i] = add(a[i], b[i]); + } + return result; + } + + template + std::array subtract(const std::array &a, const std::array &b) const { + std::array result; + for (std::size_t i = 0; i < N; ++i) { + result[i] = subtract(a[i], b[i]); + } + return result; + } + + // All tower inputs and returned coordinates are in (-q,q). Pair sums are < 2*q; + // the Karatsuba product is < 4*q^2 and its remainder < 1.5*q. Correcting by the + // two diagonal products gives a sum < 3.5*q, hence exact in binary64 before reduction. + quadratic_type multiply(const quadratic_type &a, const quadratic_type &b) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + const double diagonal0 = multiply(a[0], b[0]); + const double diagonal1 = multiply(a[1], b[1]); + const double cross = multiply(a[0] + a[1], b[0] + b[1]); + return {reduce(diagonal0 - diagonal1), reduce(cross - diagonal0 - diagonal1)}; + } + + quadratic_type multiply_xi(const quadratic_type &a) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + // (9+u)*(a0+a1*u) = (9*a0-a1) + (a0+9*a1)*u. Reduce 8*a first: + // |8*a| < 8*q < 2^53, whereas forming 9*a directly could lose an integer bit. + return {reduce(reduce(8 * a[0]) + a[0] - a[1]), reduce(reduce(8 * a[1]) + a[1] + a[0])}; + } + + cubic_type multiply_v(const cubic_type &a) const { + return {multiply_xi(a[2]), a[0], a[1]}; + } + + cubic_type multiply(const cubic_type &a, const cubic_type &b) const { + const auto diagonal0 = multiply(a[0], b[0]); + const auto diagonal1 = multiply(a[1], b[1]); + const auto diagonal2 = multiply(a[2], b[2]); + const auto cross12 = + subtract(subtract(multiply(add(a[1], a[2]), add(b[1], b[2])), diagonal1), diagonal2); + const auto cross01 = + subtract(subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0), diagonal1); + const auto cross02 = + subtract(subtract(multiply(add(a[0], a[2]), add(b[0], b[2])), diagonal0), diagonal2); + return {add(diagonal0, multiply_xi(cross12)), add(cross01, multiply_xi(diagonal2)), + add(cross02, diagonal1)}; + } + + tower_type multiply(const tower_type &a, const tower_type &b) const { + const auto diagonal0 = multiply(a[0], b[0]); + const auto diagonal1 = multiply(a[1], b[1]); + const auto cross = subtract(subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0), diagonal1); + return {add(diagonal0, multiply_v(diagonal1)), cross}; + } + + tower_type square(const tower_type &a) const { + // The quadratic-extension square needs two cubic products instead of three: + // c0 = (a0+a1)*(a0+v*a1) - a0*a1 - v*a0*a1; c1 = 2*a0*a1. + const auto cross = multiply(a[0], a[1]); + const auto combined = multiply(add(a[0], a[1]), add(a[0], multiply_v(a[1]))); + return {subtract(subtract(combined, cross), multiply_v(cross)), add(cross, cross)}; + } + }; + + public: + explicit bn254_small_prime_fft_pointwise(std::uint64_t modulus) : arithmetic_(modulus) { + } + + /** Overwrite left with scale*left*right. Both buffers must contain exactly 12*stride cells. */ + void multiply(std::span left, std::span right, std::size_t stride, double scale) const { + check_buffers(left, right, stride, scale); + if (arithmetic_.fma_supported()) { + apply(left, right, stride, scale); + } else { + apply(left, right, stride, scale); + } + } + + /** Overwrite values with scale*values^2, loading each tower element once. */ + void square(std::span values, std::size_t stride, double scale) const { + check_buffers(values, values, stride, scale); + if (arithmetic_.fma_supported()) { + apply(values, values, stride, scale); + } else { + apply(values, values, stride, scale); + } + } + + private: + void check_buffers(std::span left, std::span right, std::size_t stride, + double scale) const { + if (left.size() % 12 != 0 || left.size() / 12 != stride || left.size() != right.size()) { + throw std::invalid_argument("bn254_small_prime_fft_pointwise: invalid buffer dimensions"); + } + assert(std::isfinite(scale) && std::trunc(scale) == scale && + std::abs(scale) < static_cast(arithmetic_.modulus()) / 2); + } + + template + void apply(std::span left, std::span right, std::size_t stride, double scale) const { + const kernel operations {arithmetic_}; + for (std::size_t i = 0; i < stride; ++i) { + tower_type a, b; + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + const std::size_t offset = (6 * w + 2 * v + u) * stride + i; + a[w][v][u] = operations.reduce(left[offset]); + if constexpr (!Square) { + b[w][v][u] = operations.reduce(right[offset]); + } + } + } + } + const tower_type result = [&]() { + if constexpr (Square) { + return operations.square(a); + } else { + return operations.multiply(a, b); + } + }(); + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + // Fold inverse-transform normalization into this existing write pass. + left[(6 * w + 2 * v + u) * stride + i] = operations.multiply(result[w][v][u], scale); + } + } + } + } + } + + small_prime_fft_arithmetic arithmetic_; + }; + } // namespace nil::crypto3::math::detail #endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP 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 new file mode 100644 index 000000000..7969a05cb --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp @@ -0,0 +1,208 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#ifndef CRYPTO3_MATH_BN254_SMALL_PRIME_FFT_BACKEND_HPP +#define CRYPTO3_MATH_BN254_SMALL_PRIME_FFT_BACKEND_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace nil::crypto3::math::polynomial_arithmetic { + + /** + * Exact BN254 Fp12 polynomial products through small-prime transforms and signed CRT. + * Configure a power-of-two maximum transform length. All smaller power-of-two plans are + * cached at construction; each operation uses the smallest length containing its full linear + * product. multiply_low first discards input coefficients that cannot affect the output. + * + * One catalogue prefix and coefficient bound cover every supported product, allowing all + * lengths to share the converter, CRT cofactors and pointwise arithmetic. If a+b-1 <= Nmax, + * then min(a,b) <= max(1,Nmax/2), so C=102*max(1,Nmax/2)*(p-1)^2 is sufficient. Setup performs + * the capacity and memory-size checks once; repeated operations require no CRT/FFT setup. + * + * Spectrum buffers retain capacity across operations. Squaring converts/transforms one operand + * and uses the dedicated tower square. Results are assembled separately before swapping into + * output, permitting output to alias either input and preserving it on failure. + * Reuse sequentially; concurrent operations require separate backend instances. No constant-time + * guarantee is provided, and the thread's rounding mode must stay unchanged during an operation. + */ + class bn254_small_prime_fft_backend { + using params_type = math::detail::small_prime_fft_params; + + public: + using field_type = math::detail::bn254_small_prime_fft_converter::field_type; + using value_type = field_type::value_type; + using polynomial_type = math::polynomial; + + /** + * max_spectrum_bytes limits one spectrum at the configured maximum length. A second + * spectrum, cached twiddle tables and coefficient output storage are additional allocations. + * Invalid/unsupported lengths or an insufficient byte limit throw before allocating plans. + */ + explicit bn254_small_prime_fft_backend( + std::size_t max_transform_size, + std::size_t max_spectrum_bytes = std::numeric_limits::max()) : + bn254_small_prime_fft_backend(checked_params(max_transform_size, max_spectrum_bytes)) { + } + + std::size_t max_transform_size() const { + return max_transform_size_; + } + + std::size_t prime_count() const { + return converter_.prime_count(); + } + + void multiply(polynomial_type &output, const polynomial_type &left, const polynomial_type &right) { + multiply_low(output, left, right, std::numeric_limits::max()); + } + + void square(polynomial_type &output, const polynomial_type &input) { + multiply_prefixes(output, input, input, std::numeric_limits::max()); + } + + void multiply_low(polynomial_type &output, const polynomial_type &left, const polynomial_type &right, + std::size_t coefficient_count) { + if (&left == &right) { + multiply_prefixes(output, left, right, coefficient_count); + } else { + multiply_prefixes(output, left, right, coefficient_count); + } + } + + private: + static params_type checked_params(std::size_t max_size, std::size_t max_spectrum_bytes) { + if (!std::has_single_bit(max_size)) { + throw std::invalid_argument("bn254_small_prime_fft_backend: expected a power-of-two maximum size"); + } + // These balanced lengths have product length max_size and attain the maximum term count. + auto params = math::detail::make_bn254_small_prime_fft_params( + std::max(std::size_t(1), max_size / 2), max_size / 2 + 1, std::numeric_limits::max(), + max_spectrum_bytes); + if (!params) { + throw std::invalid_argument( + "bn254_small_prime_fft_backend: insufficient roots, CRT capacity or byte limit"); + } + return std::move(*params); + } + + explicit bn254_small_prime_fft_backend(const params_type ¶ms) : + max_transform_size_(params.transform_size), 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()); + 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); + } + 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); + } + } + } + + template + void multiply_prefixes(polynomial_type &output, const polynomial_type &left, const polynomial_type &right, + std::size_t coefficient_count) { + if (plans_.empty()) { + throw std::invalid_argument("bn254_small_prime_fft_backend: moved-from backend"); + } + const std::size_t left_size = std::min(left.size(), coefficient_count); + const std::size_t right_size = std::min(right.size(), coefficient_count); + if (math::is_zero(left.begin(), left.begin() + left_size) || + math::is_zero(right.begin(), right.begin() + right_size)) { + output.assign(1, value_type::zero()); + return; + } + // Nonzero prefixes are nonempty. Check capacity without overflowing left_size+right_size-1. + if (left_size > max_transform_size_ || right_size > max_transform_size_ - (left_size - 1)) { + throw std::invalid_argument("bn254_small_prime_fft_backend: product exceeds configured transform size"); + } + const std::size_t product_size = left_size + right_size - 1; + const std::size_t size = std::bit_ceil(product_size); + const std::size_t result_size = std::min(product_size, coefficient_count); + const std::size_t channels = field_type::arity * prime_count(); + // Constructor checks guarantee that even the largest spectrum's byte count fits size_t. + left_spectrum_.resize(channels * size); + converter_.convert(left_spectrum_, std::span(left.data(), left_size), size); + if constexpr (!Square) { + right_spectrum_.resize(channels * size); + converter_.convert(right_spectrum_, std::span(right.data(), right_size), size); + } + const std::size_t plan_offset = std::countr_zero(size) * prime_count(); + for (std::size_t prime = 0; prime < prime_count(); ++prime) { + const auto &plan = plans_[plan_offset + prime]; + auto a = std::span(left_spectrum_) + .subspan(prime * field_type::arity * size, field_type::arity * size); + std::span b; + if constexpr (!Square) { + b = std::span(right_spectrum_) + .subspan(prime * field_type::arity * size, field_type::arity * size); + } + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + plan.fft(a.subspan(coordinate * size, size)); + if constexpr (!Square) { + plan.fft(b.subspan(coordinate * size, size)); + } + } + if constexpr (Square) { + pointwise_[prime].square(a, size, plan.size_inverse()); + } else { + pointwise_[prime].multiply(a, b, size, plan.size_inverse()); + } + 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); + math::condense(result_); + output.get_storage().swap(result_.get_storage()); + } + + std::size_t max_transform_size_; + math::detail::bn254_small_prime_fft_converter converter_; + math::detail::bn254_small_prime_fft_reconstructor reconstructor_; + std::vector plans_; + std::vector pointwise_; + std::vector left_spectrum_; + std::vector right_spectrum_; + polynomial_type result_; + }; + +} // namespace nil::crypto3::math::polynomial_arithmetic + +#endif // CRYPTO3_MATH_BN254_SMALL_PRIME_FFT_BACKEND_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index ce5ce58fb..197d3f59a 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -51,6 +51,7 @@ set(TESTS_NAMES "small_prime_fft_conversion" "small_prime_fft_crt" "small_prime_fft_transform" + "small_prime_fft_backend" "matrix" "polynomial_gcd" # "kronecker_substitution" # FIXME: This fails. Disabled for passing CI diff --git a/libs/math/test/small_prime_fft_backend.cpp b/libs/math/test/small_prime_fft_backend.cpp new file mode 100644 index 000000000..864eef514 --- /dev/null +++ b/libs/math/test/small_prime_fft_backend.cpp @@ -0,0 +1,483 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE small_prime_fft_backend_test + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace { + namespace math = nil::crypto3::math; + namespace detail = math::detail; + namespace arithmetic = math::polynomial_arithmetic; + using boost::multiprecision::cpp_int; + using backend_type = arithmetic::bn254_small_prime_fft_backend; + using polynomial_type = backend_type::polynomial_type; + using value_type = backend_type::value_type; + using fq = nil::crypto3::algebra::fields::alt_bn128_base_field<254>; + using context_type = arithmetic::polynomial_context; + using reference_backend = arithmetic::schoolbook_backend; + using coordinates_type = std::array; + + static_assert(arithmetic::PolynomialBackend); + static_assert(math::detail::SupportsDivrem); + + struct restore_rounding { + int saved = std::fegetround(); + ~restore_rounding() { + std::fesetround(saved); + } + }; + + cpp_int canonical(const cpp_int &value, const cpp_int &modulus) { + cpp_int result = value % modulus; + return result < 0 ? result + modulus : result; + } + + // Expand all 144 coordinate products and reduce u^2, then w^2, then v^3. + // This oracle uses neither the Karatsuba identities nor the production field multiplication. + coordinates_type tower_product(const coordinates_type &a, const coordinates_type &b) { + coordinates_type result {}; + for (std::size_t i = 0; i < 12; ++i) { + for (std::size_t j = 0; j < 12; ++j) { + const std::size_t u = i % 2 + j % 2; + const std::size_t w = i / 6 + j / 6; + const std::size_t v = (i % 6) / 2 + (j % 6) / 2 + w / 2; + cpp_int product = a[i] * b[j]; + if (u == 2) { + product = -product; + } + const std::size_t index = 6 * (w % 2) + 2 * (v % 3) + u % 2; + if (v < 3) { + result[index] += product; + } else { + result[index] += 9 * product; + result[index ^ 1] += u % 2 ? -product : product; + } + } + } + return result; + } + + value_type field_value(const coordinates_type &coordinates) { + const cpp_int p(fq::modulus.backend().to_cpp_int()); + value_type result = value_type::zero(); + for (std::size_t c = 0; c < 12; ++c) { + result.data[c / 6].data[(c % 6) / 2].data[c % 2] = + fq::value_type(fq::integral_type(canonical(coordinates[c], p).str())); + } + return result; + } + + polynomial_type random_polynomial(std::size_t size, std::mt19937_64 &random) { + const cpp_int p(fq::modulus.backend().to_cpp_int()); + polynomial_type result(size, value_type::zero()); + for (std::size_t i = 0; i < size; ++i) { + coordinates_type coordinates {}; + for (std::size_t c = 0; c < 12; ++c) { + for (unsigned limb = 0; limb < 4; ++limb) { + coordinates[c] = (coordinates[c] << 64) + random(); + } + switch ((i + c) % 7) { + case 0: + coordinates[c] = p - 1; + break; + case 1: + coordinates[c] = 0; + break; + case 2: + coordinates[c] = 1; + break; + case 3: + coordinates[c] = cpp_int(1) << 253; + break; + case 4: + coordinates[c] = p - 2; + break; + default: + break; + } + } + result[i] = field_value(coordinates); + } + return result; + } + + polynomial_type expected_product(const polynomial_type &left, const polynomial_type &right) { + polynomial_type result; + reference_backend {}.multiply(result, left, right); + return result; + } + + polynomial_type expected_low(polynomial_type product, std::size_t count) { + if (count == 0) { + return {value_type::zero()}; + } + if (product.size() > count) { + product.resize(count); + } + math::condense(product); + return product; + } + + void check_polynomial(const polynomial_type &actual, const polynomial_type &expected) { + BOOST_REQUIRE(!actual.empty()); + BOOST_CHECK(actual == expected); + BOOST_CHECK(actual.size() == 1 || !actual[actual.size() - 1].is_zero()); + } + + void check_pointwise(std::span actual, std::span expected, + std::uint64_t modulus, double scale) { + for (std::size_t i = 0; i < expected.size(); ++i) { + for (std::size_t c = 0; c < 12; ++c) { + const double value = actual[c * expected.size() + i]; + BOOST_REQUIRE(std::isfinite(value)); + BOOST_REQUIRE(std::abs(value) < static_cast(modulus)); + BOOST_REQUIRE_EQUAL(value, std::trunc(value)); + BOOST_CHECK_EQUAL(canonical(cpp_int(static_cast(value)), cpp_int(modulus)), + canonical(expected[i][c] * static_cast(scale), cpp_int(modulus))); + } + } + } +} // namespace + +BOOST_AUTO_TEST_SUITE(small_prime_fft_backend_tests) + +BOOST_AUTO_TEST_CASE(pointwise_products_and_squares_match_expanded_integer_tower) { + constexpr std::size_t stride = 9; + std::mt19937_64 random(0x46465410); + for (const auto &prime : detail::small_prime_fft_primes) { + const auto q = static_cast(prime.modulus); + const detail::bn254_small_prime_fft_pointwise pointwise(prime.modulus); + std::vector left(12 * stride), right(12 * stride); + std::vector products(stride), squares(stride); + for (std::size_t i = 0; i < stride; ++i) { + coordinates_type a {}, b {}; + for (std::size_t c = 0; c < 12; ++c) { + const std::int64_t x = + c % 3 == 0 ? 3 * q - 1 : static_cast(random() % (6 * q - 1)) - (3 * q - 1); + const std::int64_t y = + c % 3 == 1 ? 1 - 3 * q : static_cast(random() % (6 * q - 1)) - (3 * q - 1); + a[c] = x; + b[c] = y; + left[c * stride + i] = static_cast(x); + right[c * stride + i] = static_cast(y); + } + products[i] = tower_product(a, b); + squares[i] = tower_product(a, a); + } + for (double scale : {0.0, 1.0, -1.0, static_cast((q - 1) / 2), -static_cast((q - 1) / 2)}) { + std::vector storage(left.size() + 2, -19.0); + auto actual = std::span(storage).subspan(1, left.size()); + std::copy(left.begin(), left.end(), actual.begin()); + pointwise.multiply(actual, right, stride, scale); + check_pointwise(actual, products, prime.modulus, scale); + std::copy(left.begin(), left.end(), actual.begin()); + pointwise.square(actual, stride, scale); + check_pointwise(actual, squares, prime.modulus, scale); + std::copy(left.begin(), left.end(), actual.begin()); + pointwise.multiply(actual, actual, stride, scale); + check_pointwise(actual, squares, prime.modulus, scale); + BOOST_CHECK_EQUAL(storage.front(), -19.0); + BOOST_CHECK_EQUAL(storage.back(), -19.0); + } + } +} + +BOOST_AUTO_TEST_CASE(pointwise_invalid_dimensions_throw_before_writes) { + const detail::bn254_small_prime_fft_pointwise pointwise(detail::small_prime_fft_primes[0].modulus); + std::vector left(24, 7.0), right(24, 8.0); + const auto before = left; + BOOST_CHECK_THROW(pointwise.multiply(left, right, 1, 1.0), std::invalid_argument); + BOOST_CHECK_THROW(pointwise.multiply(left, std::span(right).first(23), 2, 1.0), + std::invalid_argument); + BOOST_CHECK_THROW(pointwise.square(std::span(left).first(23), 2, 1.0), std::invalid_argument); + BOOST_CHECK_THROW(pointwise.square(left, std::numeric_limits::max(), 1.0), std::invalid_argument); + BOOST_CHECK_EQUAL_COLLECTIONS(left.begin(), left.end(), before.begin(), before.end()); + BOOST_CHECK_NO_THROW(pointwise.multiply({}, {}, 0, 1.0)); + BOOST_CHECK_NO_THROW(pointwise.square({}, 0, 1.0)); +} + +BOOST_AUTO_TEST_CASE(all_tower_basis_products_cross_the_conversion_and_crt_pipeline) { + backend_type backend(1); + polynomial_type actual; + for (std::size_t i = 0; i < 12; ++i) { + coordinates_type a {}; + a[i] = 1; + const polynomial_type left {field_value(a)}; + for (std::size_t j = 0; j < 12; ++j) { + coordinates_type b {}; + b[j] = 1; + const polynomial_type right {field_value(b)}; + backend.multiply(actual, left, right); + check_polynomial(actual, polynomial_type {field_value(tower_product(a, b))}); + } + backend.square(actual, left); + check_polynomial(actual, polynomial_type {field_value(tower_product(a, a))}); + } +} + +BOOST_AUTO_TEST_CASE(full_width_products_match_native_fp12_schoolbook) { + backend_type backend(512); + std::mt19937_64 random(0x46465411); + const std::array, 9> shapes { + {{1, 1}, {1, 17}, {2, 3}, {7, 9}, {8, 9}, {17, 33}, {32, 32}, {65, 31}, {257, 129}}}; + polynomial_type result; + for (const auto [a, b] : shapes) { + BOOST_TEST_CONTEXT("shape=" << a << "x" << b) { + const auto left = random_polynomial(a, random); + const auto right = random_polynomial(b, random); + backend.multiply(result, left, right); + check_polynomial(result, expected_product(left, right)); + } + } + // Canonical p-1 coordinates force large integer lifts and signed CRT reconstruction. + coordinates_type coordinates; + coordinates.fill(cpp_int(fq::modulus.backend().to_cpp_int()) - 1); + const polynomial_type dense(17, field_value(coordinates)); + const polynomial_type separate_copy = dense; + backend.multiply(result, dense, separate_copy); + check_polynomial(result, expected_product(dense, separate_copy)); +} + +BOOST_AUTO_TEST_CASE(dedicated_squares_match_multiplication_and_schoolbook) { + backend_type backend(256); + std::mt19937_64 random(0x46465412); + for (const std::size_t size : {1, 2, 7, 16, 33, 65}) { + const auto input = random_polynomial(size, random); + const polynomial_type separate_copy = input; + const auto expected = expected_product(input, input); + polynomial_type squared, multiplied; + backend.square(squared, input); + backend.multiply(multiplied, input, separate_copy); + check_polynomial(squared, expected); + check_polynomial(multiplied, expected); + squared = input; + backend.square(squared, squared); + check_polynomial(squared, expected); + } +} + +BOOST_AUTO_TEST_CASE(low_products_trim_inputs_and_do_not_wrap_high_coefficients) { + backend_type backend(128); + std::mt19937_64 random(0x46465413); + const auto left = random_polynomial(17, random); + const auto right = random_polynomial(33, random); + const auto full = expected_product(left, right); + polynomial_type actual; + for (const std::size_t count : {0, 1, 2, 7, 8, 9, 16, 17, 32, 49, 64}) { + backend.multiply_low(actual, left, right, count); + check_polynomial(actual, expected_low(full, count)); + backend.multiply_low(actual, left, left, count); + check_polynomial(actual, expected_low(expected_product(left, left), count)); + } + // (1+X^8)^2 mod X^9 = 1+2*X^8. A length-16 cyclic convolution would corrupt its constant term. + polynomial_type sparse(9, value_type::zero()); + sparse[0] = sparse[8] = value_type::one(); + backend.multiply_low(actual, sparse, sparse, 9); + check_polynomial(actual, expected_low(expected_product(sparse, sparse), 9)); + backend_type small_backend(4); + small_backend.multiply_low(actual, left, right, 2); + check_polynomial(actual, expected_low(full, 2)); +} + +BOOST_AUTO_TEST_CASE(context_interface_and_all_output_aliases) { + context_type context {backend_type(32)}; + std::mt19937_64 random(0x46465414); + const auto left = random_polynomial(7, random); + const auto right = random_polynomial(9, random); + const auto expected = expected_product(left, right); + polynomial_type result = left; + math::multiplication(result, result, right, context); + check_polynomial(result, expected); + result = right; + math::multiplication(result, left, result, context); + check_polynomial(result, expected); + result = left; + math::multiplication(result, result, result, context); + check_polynomial(result, expected_product(left, left)); + result = left; + math::multiply_low(result, result, right, 5, context); + check_polynomial(result, expected_low(expected, 5)); + result = right; + math::multiply_low(result, left, result, 5, context); + check_polynomial(result, expected_low(expected, 5)); + result = left; + math::multiply_low(result, result, result, 5, context); + check_polynomial(result, expected_low(expected_product(left, left), 5)); + result = left; + math::square(result, result, context); + check_polynomial(result, expected_product(left, left)); +} + +BOOST_AUTO_TEST_CASE(zeros_constants_and_cancellation_produce_canonical_results) { + backend_type backend(32); + const polynomial_type zero {value_type::zero()}, one {value_type::one()}; + const polynomial_type a {value_type::one(), -value_type::one()}; + const polynomial_type b {value_type::one(), value_type::one()}; + polynomial_type result = a; + backend.multiply(result, zero, b); + check_polynomial(result, zero); + backend.multiply(result, a, zero); + check_polynomial(result, zero); + backend.square(result, zero); + check_polynomial(result, zero); + backend.multiply_low(result, a, b, 0); + check_polynomial(result, zero); + backend.multiply(result, a, one); + check_polynomial(result, a); + backend.multiply_low(result, a, b, 2); + check_polynomial(result, one); + polynomial_type padded = a; + padded.resize(12, value_type::zero()); + backend.multiply(result, padded, b); + check_polynomial(result, expected_product(a, b)); + polynomial_type zero_prefix(17, value_type::zero()); + zero_prefix[16] = value_type::one(); + backend.multiply_low(result, zero_prefix, b, 4); + check_polynomial(result, zero); +} + +BOOST_AUTO_TEST_CASE(repeated_size_changes_and_owned_backend_state) { + backend_type original(256); + std::mt19937_64 random(0x46465415); + polynomial_type actual; + const auto initial = random_polynomial(65, random); + original.square(actual, initial); + auto copy = original; + backend_type moved = std::move(copy); + for (backend_type *backend : {&original, &moved}) { + for (const std::size_t size : {33, 1, 17, 2, 65, 7, 33, 1}) { + const auto left = random_polynomial(size, random); + const auto right = random_polynomial(size + 1, random); + const auto expected = expected_product(left, right); + backend->multiply(actual, left, right); + check_polynomial(actual, expected); + backend->square(actual, right); + check_polynomial(actual, expected_product(right, right)); + backend->multiply_low(actual, left, right, 3); + check_polynomial(actual, expected_low(expected, 3)); + } + } +} + +BOOST_AUTO_TEST_CASE(rounding_modes_preserve_pipeline_results_and_environment) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + backend_type original(16); + std::mt19937_64 random(0x46465416); + const auto left = random_polynomial(3, random); + const auto right = random_polynomial(5, random); + const auto product = expected_product(left, right); + const auto squared = expected_product(right, right); + for (const int mode : {FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO, FE_TONEAREST}) { + BOOST_REQUIRE_EQUAL(std::fesetround(mode), 0); + backend_type constructed_in_mode(16); + for (backend_type *backend : {&original, &constructed_in_mode}) { + polynomial_type actual; + backend->multiply(actual, left, right); + check_polynomial(actual, product); + backend->square(actual, right); + check_polynomial(actual, squared); + backend->multiply_low(actual, left, right, 3); + check_polynomial(actual, expected_low(product, 3)); + BOOST_CHECK_EQUAL(std::fegetround(), mode); + } + } +} + +BOOST_AUTO_TEST_CASE(configured_capacity_and_byte_limits_are_checked) { + for (const std::size_t size : + {std::size_t(0), std::size_t(3), std::size_t(6), std::numeric_limits::max()}) { + BOOST_CHECK_THROW(backend_type {size}, std::invalid_argument); + } + BOOST_CHECK_THROW(backend_type(std::size_t(1) << 20, 0), std::invalid_argument); + const auto params = detail::make_bn254_small_prime_fft_params(16, 17); + BOOST_REQUIRE(params.has_value()); + BOOST_CHECK_THROW(backend_type(32, params->spectrum_bytes - 1), std::invalid_argument); + backend_type backend(32, params->spectrum_bytes); + BOOST_CHECK_EQUAL(backend.max_transform_size(), 32); + BOOST_CHECK_EQUAL(backend.prime_count(), params->crt.prime_count); + std::mt19937_64 random(0x46465417); + const auto large = random_polynomial(17, random); + polynomial_type result {value_type::one()}; + const auto before = result; + BOOST_CHECK_THROW(backend.multiply(result, large, large), std::invalid_argument); + BOOST_CHECK_THROW(backend.square(result, large), std::invalid_argument); + BOOST_CHECK_THROW(backend.multiply_low(result, large, large, 17), std::invalid_argument); + check_polynomial(result, before); + polynomial_type alias = large; + BOOST_CHECK_THROW(backend.square(alias, alias), std::invalid_argument); + check_polynomial(alias, large); + backend.multiply_low(result, large, large, 16); + check_polynomial(result, expected_low(expected_product(large, large), 16)); +} + +BOOST_AUTO_TEST_CASE(existing_modular_multiplication_uses_the_new_backend) { + arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + context_type context {backend_type(64), options}; + arithmetic::polynomial_context reference {reference_backend {}, options}; + std::mt19937_64 random(0x46465418); + auto divisor = random_polynomial(9, random); + divisor[8] = value_type::one(); + math::polynomial_divisor_context fixed(divisor, 8, context); + math::polynomial_divisor_context reference_fixed(divisor, 8, reference); + for (unsigned trial = 0; trial < 3; ++trial) { + const auto left = random_polynomial(8, random); + const auto right = random_polynomial(8, random); + polynomial_type actual, expected; + math::mulmod(expected, left, right, reference_fixed, reference); + math::mulmod(actual, left, right, fixed, context); + check_polynomial(actual, expected); + actual = left; + math::mulmod(actual, actual, right, fixed, context); + check_polynomial(actual, expected); + math::squaremod(expected, left, reference_fixed, reference); + actual = left; + math::squaremod(actual, actual, fixed, context); + check_polynomial(actual, expected); + } +} + +BOOST_AUTO_TEST_SUITE_END() From 58dbc46be916775050f976b0858e7f5365ad8d6f Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 13:42:30 -0400 Subject: [PATCH 07/11] Prepared low-products --- .../bn254_small_prime_fft_backend.hpp | 146 +++++++- libs/math/test/small_prime_fft_backend.cpp | 326 ++++++++++++++++++ 2 files changed, 465 insertions(+), 7 deletions(-) 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 7969a05cb..1b6c3645d 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 @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -66,6 +67,44 @@ namespace nil::crypto3::math::polynomial_arithmetic { using value_type = field_type::value_type; using polynomial_type = math::polynomial; + /** + * Owned, immutable spectrum of a fixed operand's low prefix. The preparation has its own + * allocation and survives changes to the source polynomial or the preparing backend. + * Nonzero preparations identify the transform length and catalogue prefix; compatible + * backend instances can reuse them without retaining references to plans or workspace. + */ + class prepared_low_product_type { + public: + std::size_t coefficient_count() const noexcept { + return coefficient_count_; + } + + // Zero denotes a zero fixed prefix or zero requested precision, requiring no spectrum. + std::size_t transform_size() const noexcept { + return transform_size_; + } + + std::size_t storage_bytes() const noexcept { + return spectrum_.capacity() * sizeof(double); + } + + private: + friend class bn254_small_prime_fft_backend; + + prepared_low_product_type(std::size_t coefficient_count, std::size_t prefix_size, + std::size_t transform_size = 0, std::size_t prime_count = 0, + std::vector spectrum = {}) : + coefficient_count_(coefficient_count), prefix_size_(prefix_size), transform_size_(transform_size), + prime_count_(prime_count), spectrum_(std::move(spectrum)) { + } + + std::size_t coefficient_count_; + std::size_t prefix_size_; + std::size_t transform_size_; + std::size_t prime_count_; + std::vector spectrum_; + }; + /** * max_spectrum_bytes limits one spectrum at the configured maximum length. A second * spectrum, cached twiddle tables and coefficient output storage are additional allocations. @@ -102,6 +141,88 @@ namespace nil::crypto3::math::polynomial_arithmetic { } } + /** + * Prepare the fixed right operand for multiply_low, selecting the same prefixes and smallest + * transform as an ordinary product. Positive precision with no variable coefficients returns + * nullopt. Nonzero preparations also require a live backend and capacity for the complete + * prefix product. Zero precision and zero fixed prefixes need no transform or spectrum allocation. + * The source is a coefficient polynomial; trailing zero coefficients are permitted. + */ + std::optional prepare_low_product(const polynomial_type &fixed, + std::size_t variable_coefficient_count, + std::size_t coefficient_count) { + const std::size_t prefix_size = std::min(fixed.size(), coefficient_count); + if (coefficient_count == 0) { + return prepared_low_product_type(coefficient_count, prefix_size); + } + const std::size_t variable_size = std::min(variable_coefficient_count, coefficient_count); + if (variable_size == 0) { + return std::nullopt; + } + if (math::is_zero(fixed.begin(), fixed.begin() + prefix_size)) { + return prepared_low_product_type(coefficient_count, prefix_size); + } + if (plans_.empty() || variable_size > max_transform_size_ || + prefix_size > max_transform_size_ - (variable_size - 1)) { + return std::nullopt; + } + const std::size_t size = std::bit_ceil(variable_size + prefix_size - 1); + // Allocate just this spectrum, so a small preparation never retains larger scratch capacity. + std::vector spectrum(field_type::arity * prime_count() * size); + converter_.convert(spectrum, std::span(fixed.data(), prefix_size), size); + const std::size_t plan_offset = std::countr_zero(size) * prime_count(); + for (std::size_t prime = 0; prime < prime_count(); ++prime) { + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + plans_[plan_offset + prime].fft( + std::span(spectrum).subspan((prime * field_type::arity + coordinate) * size, size)); + } + } + return prepared_low_product_type(coefficient_count, prefix_size, size, prime_count(), std::move(spectrum)); + } + + /** + * Compute left*fixed mod X^coefficient_count without converting or transforming fixed again. + * Return false, leaving output unchanged, for mismatched precision, a moved-from nonzero + * preparation, or an incompatible backend/shape. Nonzero products require the same selected + * FFT length and catalogue prefix; zero products need neither. Successful output is canonical + * and may alias left. The preparation is never modified or copied into working storage. + * + * Matching primes and a product length fitting the current backend suffice for CRT safety: + * its existing coefficient bound covers every linear product that fits, even when the + * preparation came from a backend with a different configured maximum. + */ + bool try_multiply_low_prepared(polynomial_type &output, const polynomial_type &left, + const prepared_low_product_type &prepared, std::size_t coefficient_count) { + if (coefficient_count != prepared.coefficient_count_) { + return false; + } + if (prepared.transform_size_ == 0) { + output.assign(1, value_type::zero()); + return true; + } + if (prepared.spectrum_.empty()) { + return false; + } + const std::size_t left_size = std::min(left.size(), coefficient_count); + if (math::is_zero(left.begin(), left.begin() + left_size)) { + output.assign(1, value_type::zero()); + return true; + } + if (plans_.empty() || prepared.prime_count_ != prime_count() || left_size > max_transform_size_ || + prepared.prefix_size_ > max_transform_size_ - (left_size - 1)) { + return false; + } + const std::size_t product_size = left_size + prepared.prefix_size_ - 1; + const std::size_t size = std::bit_ceil(product_size); + if (size != prepared.transform_size_) { + return false; + } + left_spectrum_.resize(field_type::arity * prime_count() * size); + converter_.convert(left_spectrum_, std::span(left.data(), left_size), size); + finish_product(output, size, std::min(product_size, coefficient_count), prepared.spectrum_); + return true; + } + private: static params_type checked_params(std::size_t max_size, std::size_t max_spectrum_bytes) { if (!std::has_single_bit(max_size)) { @@ -162,20 +283,31 @@ namespace nil::crypto3::math::polynomial_arithmetic { right_spectrum_.resize(channels * size); converter_.convert(right_spectrum_, std::span(right.data(), right_size), size); } + finish_product(output, size, result_size); + } + + // Keep ordinary multiplication, squaring and prepared reuse on the same inverse/CRT path. + // A prepared right spectrum is read directly and skips all of the right operand's forward FFTs. + template + void finish_product(polynomial_type &output, std::size_t size, std::size_t result_size, + std::span prepared_spectrum = {}) { const std::size_t plan_offset = std::countr_zero(size) * prime_count(); for (std::size_t prime = 0; prime < prime_count(); ++prime) { const auto &plan = plans_[plan_offset + prime]; - auto a = std::span(left_spectrum_) - .subspan(prime * field_type::arity * size, field_type::arity * size); - std::span b; + const std::size_t offset = prime * field_type::arity * size; + auto a = std::span(left_spectrum_).subspan(offset, field_type::arity * size); + std::span b; if constexpr (!Square) { - b = std::span(right_spectrum_) - .subspan(prime * field_type::arity * size, field_type::arity * size); + if constexpr (Prepared) { + b = prepared_spectrum.subspan(offset, field_type::arity * size); + } else { + b = std::span(right_spectrum_).subspan(offset, field_type::arity * size); + } } for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { plan.fft(a.subspan(coordinate * size, size)); - if constexpr (!Square) { - plan.fft(b.subspan(coordinate * size, size)); + if constexpr (!Square && !Prepared) { + plan.fft(std::span(right_spectrum_).subspan(offset + coordinate * size, size)); } } if constexpr (Square) { diff --git a/libs/math/test/small_prime_fft_backend.cpp b/libs/math/test/small_prime_fft_backend.cpp index 864eef514..eea906cbe 100644 --- a/libs/math/test/small_prime_fft_backend.cpp +++ b/libs/math/test/small_prime_fft_backend.cpp @@ -60,6 +60,41 @@ namespace { static_assert(arithmetic::PolynomialBackend); static_assert(math::detail::SupportsDivrem); + static_assert(math::detail::SupportsPreparedLowProduct); + + struct prepared_product_counts { + std::size_t preparations = 0; + std::size_t prepared_products = 0; + std::size_t ordinary_products = 0; + }; + + // Observe dispatch at the backend boundary while keeping all arithmetic in the real backend, + // following the existing division tests' counting-backend convention. + struct counting_prepared_backend : backend_type { + counting_prepared_backend(std::size_t size, prepared_product_counts &counts) : + backend_type(size), counts(counts) { + } + + auto prepare_low_product(const polynomial_type &fixed, std::size_t variable_count, std::size_t count) { + ++counts.preparations; + return backend_type::prepare_low_product(fixed, variable_count, count); + } + + bool try_multiply_low_prepared(polynomial_type &output, const polynomial_type &input, + const prepared_low_product_type &prepared, std::size_t count) { + const bool accepted = backend_type::try_multiply_low_prepared(output, input, prepared, count); + counts.prepared_products += accepted; + return accepted; + } + + void multiply_low(polynomial_type &output, const polynomial_type &left, const polynomial_type &right, + std::size_t count) { + ++counts.ordinary_products; + backend_type::multiply_low(output, left, right, count); + } + + prepared_product_counts &counts; + }; struct restore_rounding { int saved = std::fegetround(); @@ -480,4 +515,295 @@ BOOST_AUTO_TEST_CASE(existing_modular_multiplication_uses_the_new_backend) { } } +BOOST_AUTO_TEST_CASE(prepared_products_match_native_and_ordinary_low_products) { + context_type context {backend_type(256)}; + std::mt19937_64 random(0x46465420); + const std::array, 5> shapes {{{1, 1}, {3, 5}, {7, 9}, {17, 33}, {65, 31}}}; + for (const auto [a, b] : shapes) { + const auto input = random_polynomial(a, random); + const auto fixed = random_polynomial(b, random); + for (const std::size_t count : {0, 1, 2, 7, 8, 9, 16, 17, 32, 49, 96}) { + BOOST_TEST_CONTEXT("shape=" << a << "x" << b << ", count=" << count) { + const auto prepared = context.prepare_low_product(fixed, input.size(), count); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->coefficient_count(), count); + const std::size_t size = count == 0 ? 0 : std::bit_ceil(std::min(a, count) + std::min(b, count) - 1); + BOOST_CHECK_EQUAL(prepared->transform_size(), size); + polynomial_type variable = input; + for (unsigned repeat = 0; repeat < 3; ++repeat) { + const auto expected = expected_low(expected_product(variable, fixed), count); + polynomial_type ordinary, actual; + context.multiply_low(ordinary, variable, fixed, count); + check_polynomial(ordinary, expected); + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *prepared, count)); + check_polynomial(actual, expected); + actual = variable; + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, actual, *prepared, count)); + check_polynomial(actual, expected); + variable[0] += value_type::one(); + } + } + } + } + // Large canonical lifts exercise signed CRT through the prepared path as well. + coordinates_type coordinates; + coordinates.fill(cpp_int(fq::modulus.backend().to_cpp_int()) - 1); + const polynomial_type dense(17, field_value(coordinates)); + const auto prepared = context.prepare_low_product(dense, dense.size(), dense.size()); + BOOST_REQUIRE(prepared.has_value()); + polynomial_type actual; + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, dense, *prepared, dense.size())); + check_polynomial(actual, expected_low(expected_product(dense, dense), dense.size())); +} + +BOOST_AUTO_TEST_CASE(prepared_products_handle_zero_and_truncated_prefixes) { + context_type context {backend_type(4)}; + const polynomial_type zero {value_type::zero()}, one {value_type::one()}; + polynomial_type input(17, value_type::one()), actual; + const auto prepared = context.prepare_low_product(input, input.size(), 2); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 4); + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, input, *prepared, 2)); + check_polynomial(actual, expected_low(expected_product(input, input), 2)); + + context_type tiny {backend_type(1)}; + // Zero variable prefixes need no compatible transform, even in a smaller backend. + BOOST_REQUIRE(tiny.try_multiply_low_prepared(actual, zero, *prepared, 2)); + check_polynomial(actual, zero); + polynomial_type zero_prefix(17, value_type::zero()); + zero_prefix[16] = value_type::one(); + BOOST_REQUIRE(tiny.try_multiply_low_prepared(actual, zero_prefix, *prepared, 2)); + check_polynomial(actual, zero); + for (const auto &fixed : {zero, zero_prefix}) { + const auto empty = context.prepare_low_product(fixed, input.size(), 2); + BOOST_REQUIRE(empty.has_value()); + BOOST_CHECK_EQUAL(empty->transform_size(), 0); + BOOST_CHECK_EQUAL(empty->storage_bytes(), 0); + BOOST_REQUIRE(tiny.try_multiply_low_prepared(actual, input, *empty, 2)); + check_polynomial(actual, zero); + actual = one; + BOOST_CHECK(!tiny.try_multiply_low_prepared(actual, input, *empty, 1)); + check_polynomial(actual, one); + } + const auto empty = context.prepare_low_product(input, 0, 0); + BOOST_REQUIRE(empty.has_value()); + BOOST_CHECK_EQUAL(empty->transform_size(), 0); + BOOST_CHECK_EQUAL(empty->storage_bytes(), 0); + BOOST_REQUIRE(tiny.try_multiply_low_prepared(actual, input, *empty, 0)); + check_polynomial(actual, zero); + BOOST_CHECK(!context.prepare_low_product(zero, 0, 2).has_value()); + + // The retained prefixes can end in zeros, although the original polynomials are canonical. + const polynomial_type fixed {value_type(2), value_type::zero(), value_type(7)}; + const polynomial_type variable {value_type(3), value_type::zero(), value_type::one()}; + const auto trailing = context.prepare_low_product(fixed, variable.size(), 2); + BOOST_REQUIRE(trailing.has_value()); + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *trailing, 2)); + check_polynomial(actual, polynomial_type {value_type(6)}); + + context_type larger {backend_type(32)}; + polynomial_type sparse(9, value_type::zero()); + sparse[0] = sparse[8] = value_type::one(); + const auto no_wrap = larger.prepare_low_product(sparse, sparse.size(), 9); + BOOST_REQUIRE(no_wrap.has_value()); + BOOST_CHECK_EQUAL(no_wrap->transform_size(), 32); + BOOST_REQUIRE(larger.try_multiply_low_prepared(actual, sparse, *no_wrap, 9)); + check_polynomial(actual, expected_low(expected_product(sparse, sparse), 9)); +} + +BOOST_AUTO_TEST_CASE(prepared_products_require_matching_precision_and_transform) { + context_type context {backend_type(32)}; + std::mt19937_64 random(0x46465421); + const auto fixed = random_polynomial(3, random); + const auto prepared = context.prepare_low_product(fixed, 3, 10); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 8); + for (const std::size_t size : {3, 4, 5, 6}) { + const auto variable = random_polynomial(size, random); + polynomial_type actual; + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *prepared, 10)); + check_polynomial(actual, expected_product(variable, fixed)); + } + const auto variable = random_polynomial(3, random); + polynomial_type actual {value_type::one()}; + const auto before = actual; + for (const std::size_t count : {0, 9, 11}) { + BOOST_CHECK(!context.try_multiply_low_prepared(actual, variable, *prepared, count)); + check_polynomial(actual, before); + } + for (const std::size_t size : {1, 2, 7, 33}) { + const auto input = random_polynomial(size, random); + actual = input; + BOOST_CHECK(!context.try_multiply_low_prepared(actual, actual, *prepared, 10)); + check_polynomial(actual, input); + } + context_type too_small {backend_type(4)}; + actual = before; + BOOST_CHECK(!too_small.try_multiply_low_prepared(actual, variable, *prepared, 10)); + check_polynomial(actual, before); + BOOST_CHECK(!too_small.prepare_low_product(fixed, 3, 3).has_value()); + BOOST_CHECK(!context.prepare_low_product(fixed, 0, 3).has_value()); + BOOST_CHECK(!context + .prepare_low_product(fixed, std::numeric_limits::max(), + std::numeric_limits::max()) + .has_value()); +} + +BOOST_AUTO_TEST_CASE(prepared_products_own_spectra_across_contexts_and_moves) { + std::mt19937_64 random(0x46465422); + const auto fixed = random_polynomial(5, random); + const auto variable = random_polynomial(7, random); + const auto prepared = [&] { + context_type temporary {backend_type(64)}; + polynomial_type source = fixed; + auto result = temporary.prepare_low_product(source, variable.size(), 17); + source[0] += value_type::one(); + return result; + }(); + BOOST_REQUIRE(prepared.has_value()); + auto copy = *prepared; + const auto moved = std::move(copy); + const auto expected = expected_product(variable, fixed); + for (const std::size_t size : {16, 128}) { + context_type context {backend_type(size)}; + for (const auto *snapshot : {&*prepared, &moved}) { + polynomial_type actual; + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *snapshot, 17)); + check_polynomial(actual, expected); + auto changed_fixed = fixed; + changed_fixed[0] += value_type::one(); + const auto changed = context.prepare_low_product(changed_fixed, variable.size(), 17); + BOOST_REQUIRE(changed.has_value()); + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *changed, 17)); + check_polynomial(actual, expected_product(variable, changed_fixed)); + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *snapshot, 17)); + check_polynomial(actual, expected); + } + polynomial_type actual = variable; + BOOST_CHECK(!context.try_multiply_low_prepared(actual, actual, copy, 17)); + check_polynomial(actual, variable); + } + backend_type source(32); + backend_type destination = std::move(source); + polynomial_type actual = variable; + BOOST_CHECK(!source.prepare_low_product(fixed, variable.size(), 17).has_value()); + BOOST_CHECK(!source.try_multiply_low_prepared(actual, actual, *prepared, 17)); + check_polynomial(actual, variable); + BOOST_REQUIRE(destination.try_multiply_low_prepared(actual, variable, *prepared, 17)); + check_polynomial(actual, expected); +} + +BOOST_AUTO_TEST_CASE(small_preparations_do_not_retain_large_workspaces) { + backend_type backend(256); + polynomial_type actual, large(65, value_type::one()); + backend.multiply(actual, large, large); + const polynomial_type small(9, value_type::one()); + const auto prepared = backend.prepare_low_product(small, small.size(), small.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 32); + const std::size_t bytes_per_position = 12 * backend.prime_count() * sizeof(double); + BOOST_CHECK_GE(prepared->storage_bytes(), 32 * bytes_per_position); + BOOST_CHECK_LT(prepared->storage_bytes(), 256 * bytes_per_position); + BOOST_REQUIRE(backend.try_multiply_low_prepared(actual, small, *prepared, small.size())); + polynomial_type expected(small.size(), value_type::zero()); + for (std::size_t i = 0; i < expected.size(); ++i) { + expected[i] = value_type(i + 1); + } + check_polynomial(actual, expected); +} + +BOOST_AUTO_TEST_CASE(prepared_products_survive_rounding_mode_changes) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + context_type context {backend_type(32)}; + std::mt19937_64 random(0x46465423); + const auto fixed = random_polynomial(5, random); + const auto variable = random_polynomial(9, random); + const auto expected = expected_low(expected_product(variable, fixed), 9); + for (const int prepare_mode : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO}) { + BOOST_REQUIRE_EQUAL(std::fesetround(prepare_mode), 0); + context_type preparing {backend_type(32)}; + const auto prepared = preparing.prepare_low_product(fixed, variable.size(), 9); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(std::fegetround(), prepare_mode); + for (const int use_mode : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO}) { + BOOST_REQUIRE_EQUAL(std::fesetround(use_mode), 0); + polynomial_type actual; + BOOST_REQUIRE(context.try_multiply_low_prepared(actual, variable, *prepared, 9)); + check_polynomial(actual, expected); + BOOST_CHECK_EQUAL(std::fegetround(), use_mode); + } + } +} + +BOOST_AUTO_TEST_CASE(divisor_context_reuses_prepared_products_and_falls_back) { + prepared_product_counts counts; + arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + using counted_context = arithmetic::polynomial_context; + counted_context preparing {counting_prepared_backend(64, counts), options}; + std::mt19937_64 random(0x46465424); + auto divisor = random_polynomial(9, random); + divisor[8] = value_type::one(); + const math::polynomial_divisor_context fixed(divisor, 7, preparing); + BOOST_CHECK_EQUAL(counts.preparations, 2); + auto copy = fixed; + const auto moved = std::move(copy); + counted_context context {counting_prepared_backend(32, counts), options}; + const auto remainder = random_polynomial(6, random); + for (const auto *snapshot : {&fixed, &moved}) { + for (const std::size_t size : {7, 6, 2, 1, 7}) { + auto quotient = random_polynomial(size, random); + quotient[size - 1] = value_type::one(); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + counts = {}; + polynomial_type actual_quotient, actual_remainder; + math::divrem(actual_quotient, actual_remainder, dividend, *snapshot, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + const std::size_t accepted = size == 7 ? 2 : (size == 1 ? 0 : 1); + BOOST_CHECK_EQUAL(counts.prepared_products, accepted); + BOOST_CHECK_EQUAL(counts.ordinary_products, 2 - accepted); + BOOST_CHECK_EQUAL(counts.preparations, 0); + actual_quotient = dividend; + math::divrem(actual_quotient, actual_remainder, actual_quotient, *snapshot, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + actual_remainder = dividend; + math::divrem(actual_quotient, actual_remainder, actual_remainder, *snapshot, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_divisor_context_handles_sparse_fixed_prefixes) { + prepared_product_counts counts; + arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + arithmetic::polynomial_context context {counting_prepared_backend(32, counts), options}; + std::mt19937_64 random(0x46465425); + const auto quotient = random_polynomial(7, random); + const auto remainder = random_polynomial(3, random); + for (const auto &constant : {value_type::zero(), value_type::one()}) { + polynomial_type divisor(9, value_type::zero()); + divisor[0] = constant; + divisor[8] = value_type::one(); + const math::polynomial_divisor_context fixed(divisor, 7, context); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + counts = {}; + polynomial_type actual_quotient, actual_remainder; + math::divrem(actual_quotient, actual_remainder, dividend, fixed, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + BOOST_CHECK_EQUAL(counts.prepared_products, 2); + BOOST_CHECK_EQUAL(counts.ordinary_products, 0); + BOOST_CHECK_EQUAL(counts.preparations, 0); + } +} + BOOST_AUTO_TEST_SUITE_END() From 46d2d36dea552320c499c309bc037574cbe57b4d Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 15:40:55 -0400 Subject: [PATCH 08/11] Prepared cyclic remainder --- .../bn254_small_prime_fft_backend.hpp | 150 ++++++++- libs/math/test/small_prime_fft_backend.cpp | 294 +++++++++++++++++- 2 files changed, 425 insertions(+), 19 deletions(-) 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 1b6c3645d..848b248bb 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 @@ -45,8 +45,9 @@ namespace nil::crypto3::math::polynomial_arithmetic { /** * Exact BN254 Fp12 polynomial products through small-prime transforms and signed CRT. * Configure a power-of-two maximum transform length. All smaller power-of-two plans are - * cached at construction; each operation uses the smallest length containing its full linear - * product. multiply_low first discards input coefficients that cannot affect the output. + * cached at construction; ordinary products use the smallest length containing their full linear + * product. multiply_low first discards input coefficients that cannot affect the output. Prepared + * cyclic remainder reconstruction can use a shorter transform once the exact quotient is known. * * One catalogue prefix and coefficient bound cover every supported product, allowing all * lengths to share the converter, CRT cofactors and pointwise arithmetic. If a+b-1 <= Nmax, @@ -105,6 +106,37 @@ namespace nil::crypto3::math::polynomial_arithmetic { std::vector spectrum_; }; + // Owns the transformed complete divisor modulo X^N - 1, including its leading coefficient. + // Its distinct type keeps cyclic reconstruction separate from multiplication modulo X^k. + class prepared_cyclic_remainder_type { + public: + std::size_t divisor_degree() const noexcept { + return divisor_degree_; + } + + std::size_t transform_size() const noexcept { + return transform_size_; + } + + std::size_t storage_bytes() const noexcept { + return spectrum_.capacity() * sizeof(double); + } + + private: + friend class bn254_small_prime_fft_backend; + + prepared_cyclic_remainder_type(std::size_t divisor_degree, std::size_t transform_size, + std::size_t prime_count, std::vector spectrum) : + divisor_degree_(divisor_degree), transform_size_(transform_size), prime_count_(prime_count), + spectrum_(std::move(spectrum)) { + } + + std::size_t divisor_degree_; + std::size_t transform_size_; + std::size_t prime_count_; + std::vector spectrum_; + }; + /** * max_spectrum_bytes limits one spectrum at the configured maximum length. A second * spectrum, cached twiddle tables and coefficient output storage are additional allocations. @@ -167,17 +199,9 @@ namespace nil::crypto3::math::polynomial_arithmetic { return std::nullopt; } const std::size_t size = std::bit_ceil(variable_size + prefix_size - 1); - // Allocate just this spectrum, so a small preparation never retains larger scratch capacity. - std::vector spectrum(field_type::arity * prime_count() * size); - converter_.convert(spectrum, std::span(fixed.data(), prefix_size), size); - const std::size_t plan_offset = std::countr_zero(size) * prime_count(); - for (std::size_t prime = 0; prime < prime_count(); ++prime) { - for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { - plans_[plan_offset + prime].fft( - std::span(spectrum).subspan((prime * field_type::arity + coordinate) * size, size)); - } - } - return prepared_low_product_type(coefficient_count, prefix_size, size, prime_count(), std::move(spectrum)); + return prepared_low_product_type( + coefficient_count, prefix_size, size, prime_count(), + prepare_spectrum(std::span(fixed.data(), prefix_size), size)); } /** @@ -223,7 +247,93 @@ namespace nil::crypto3::math::polynomial_arithmetic { return true; } + /** + * Prepare the complete divisor H for cyclic remainder reconstruction in the reduced-product + * range degree(V) <= 2*n-2, n=degree(H). Cap the expected quotient length at n-1. Use the smallest + * N >= n only when it is shorter than the supported ordinary low-product transform. Constant, + * linear and monomial divisors, unsuitable sizes and moved-from backends return nullopt. + * + * Fold H in Fp12 before conversion: if N==n, add its leading coefficient into cell zero. + * Canonical field addition keeps every converted coordinate in [0,p-1], preserving the CRT + * bound used by cyclic_remainder_size. The preparation owns its spectrum and cached metadata. + * + * @pre divisor is a nonempty canonical coefficient polynomial. + */ + std::optional prepare_cyclic_remainder(const polynomial_type &divisor, + std::size_t quotient_coefficient_count) { + if (divisor.size() < 3) { + return std::nullopt; + } + const std::size_t n = divisor.size() - 1; + const std::size_t size = cyclic_remainder_size(n, std::min(quotient_coefficient_count, n - 1)); + if (size == 0 || math::is_zero(divisor.begin(), divisor.begin() + n)) { + return std::nullopt; + } + std::span coefficients(divisor.data(), divisor.size()); + if (size == n) { + result_.get_storage().assign(divisor.begin(), divisor.begin() + n); + result_[0] += divisor[n]; + coefficients = std::span(result_.data(), n); + } + return prepared_cyclic_remainder_type(n, size, prime_count(), prepare_spectrum(coefficients, size)); + } + + /** + * Reconstruct R=fold_N(V)-(Q*H mod (X^N-1)). Since degree(R) n - 1 || + quotient.size() != dividend.size() - n || prepared.prime_count_ != prime_count()) { + return false; + } + const std::size_t size = cyclic_remainder_size(n, quotient.size()); + if (size == 0 || size != prepared.transform_size_) { + return false; + } + // Q has fewer than n coefficients, so it needs only zero padding, never cyclic folding. + left_spectrum_.resize(field_type::arity * prime_count() * size); + converter_.convert(left_spectrum_, std::span(quotient.data(), quotient.size()), size); + finish_product(output, size, size, prepared.spectrum_, + std::span(dividend.data(), dividend.size())); + return true; + } + private: + std::size_t cyclic_remainder_size(std::size_t n, std::size_t quotient_size) const { + if (plans_.empty() || n < 2 || quotient_size == 0 || quotient_size > n - 1 || n > max_transform_size_ || + quotient_size > max_transform_size_ - (n - 1)) { + return 0; + } + // Require the ordinary Q*(H mod X^n) product to fit, as in the existing cyclic backend. + // Then 2*quotient_size <= quotient_size+n-1 <= Nmax. Each cyclic coefficient contains at + // most quotient_size products of canonical Fp12 values, so the existing Nmax/2 CRT bound + // remains sufficient even across backend instances. Do not infer this bound from N alone. + const std::size_t size = std::bit_ceil(n); + return size < std::bit_ceil(quotient_size + n - 1) ? size : 0; + } + + std::vector prepare_spectrum(std::span coefficients, std::size_t size) const { + // Allocate just this spectrum, so a small preparation never retains larger scratch capacity. + std::vector spectrum(field_type::arity * prime_count() * size); + converter_.convert(spectrum, coefficients, size); + const std::size_t plan_offset = std::countr_zero(size) * prime_count(); + for (std::size_t prime = 0; prime < prime_count(); ++prime) { + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + plans_[plan_offset + prime].fft( + std::span(spectrum).subspan((prime * field_type::arity + coordinate) * size, size)); + } + } + return spectrum; + } + static params_type checked_params(std::size_t max_size, std::size_t max_spectrum_bytes) { if (!std::has_single_bit(max_size)) { throw std::invalid_argument("bn254_small_prime_fft_backend: expected a power-of-two maximum size"); @@ -290,7 +400,8 @@ namespace nil::crypto3::math::polynomial_arithmetic { // A prepared right spectrum is read directly and skips all of the right operand's forward FFTs. template void finish_product(polynomial_type &output, std::size_t size, std::size_t result_size, - std::span prepared_spectrum = {}) { + std::span prepared_spectrum = {}, + std::span cyclic_dividend = {}) { const std::size_t plan_offset = std::countr_zero(size) * prime_count(); for (std::size_t prime = 0; prime < prime_count(); ++prime) { const auto &plan = plans_[plan_offset + prime]; @@ -321,6 +432,17 @@ namespace nil::crypto3::math::polynomial_arithmetic { } result_.resize(result_size); reconstructor_.reconstruct(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. + for (std::size_t i = 0; i < size; ++i) { + value_type folded = i < cyclic_dividend.size() ? cyclic_dividend[i] : value_type::zero(); + if (i + size < cyclic_dividend.size()) { + folded += cyclic_dividend[i + size]; + } + result_[i] = folded - result_[i]; + } + } math::condense(result_); output.get_storage().swap(result_.get_storage()); } diff --git a/libs/math/test/small_prime_fft_backend.cpp b/libs/math/test/small_prime_fft_backend.cpp index eea906cbe..43e76d928 100644 --- a/libs/math/test/small_prime_fft_backend.cpp +++ b/libs/math/test/small_prime_fft_backend.cpp @@ -61,11 +61,14 @@ namespace { static_assert(arithmetic::PolynomialBackend); static_assert(math::detail::SupportsDivrem); static_assert(math::detail::SupportsPreparedLowProduct); + static_assert(math::detail::SupportsPreparedCyclicRemainder); struct prepared_product_counts { std::size_t preparations = 0; std::size_t prepared_products = 0; std::size_t ordinary_products = 0; + std::size_t cyclic_preparations = 0; + std::size_t cyclic_remainders = 0; }; // Observe dispatch at the backend boundary while keeping all arithmetic in the real backend, @@ -93,6 +96,18 @@ namespace { backend_type::multiply_low(output, left, right, count); } + auto prepare_cyclic_remainder(const polynomial_type &divisor, std::size_t quotient_count) { + ++counts.cyclic_preparations; + return backend_type::prepare_cyclic_remainder(divisor, quotient_count); + } + + bool try_cyclic_remainder(polynomial_type &output, const polynomial_type ÷nd, + const polynomial_type "ient, const prepared_cyclic_remainder_type &prepared) { + const bool accepted = backend_type::try_cyclic_remainder(output, dividend, quotient, prepared); + counts.cyclic_remainders += accepted; + return accepted; + } + prepared_product_counts &counts; }; @@ -747,7 +762,8 @@ BOOST_AUTO_TEST_CASE(divisor_context_reuses_prepared_products_and_falls_back) { auto divisor = random_polynomial(9, random); divisor[8] = value_type::one(); const math::polynomial_divisor_context fixed(divisor, 7, preparing); - BOOST_CHECK_EQUAL(counts.preparations, 2); + BOOST_CHECK_EQUAL(counts.preparations, 1); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 1); auto copy = fixed; const auto moved = std::move(copy); counted_context context {counting_prepared_backend(32, counts), options}; @@ -763,10 +779,13 @@ BOOST_AUTO_TEST_CASE(divisor_context_reuses_prepared_products_and_falls_back) { math::divrem(actual_quotient, actual_remainder, dividend, *snapshot, context); check_polynomial(actual_quotient, quotient); check_polynomial(actual_remainder, remainder); - const std::size_t accepted = size == 7 ? 2 : (size == 1 ? 0 : 1); - BOOST_CHECK_EQUAL(counts.prepared_products, accepted); - BOOST_CHECK_EQUAL(counts.ordinary_products, 2 - accepted); + const std::size_t inverse_accepted = size == 7 ? 1 : 0; + const std::size_t cyclic_accepted = size > 1 ? 1 : 0; + BOOST_CHECK_EQUAL(counts.prepared_products, inverse_accepted); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, cyclic_accepted); + BOOST_CHECK_EQUAL(counts.ordinary_products, 2 - inverse_accepted - cyclic_accepted); BOOST_CHECK_EQUAL(counts.preparations, 0); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 0); actual_quotient = dividend; math::divrem(actual_quotient, actual_remainder, actual_quotient, *snapshot, context); check_polynomial(actual_quotient, quotient); @@ -800,9 +819,274 @@ BOOST_AUTO_TEST_CASE(prepared_divisor_context_handles_sparse_fixed_prefixes) { math::divrem(actual_quotient, actual_remainder, dividend, fixed, context); check_polynomial(actual_quotient, quotient); check_polynomial(actual_remainder, remainder); - BOOST_CHECK_EQUAL(counts.prepared_products, 2); + BOOST_CHECK_EQUAL(counts.prepared_products, constant.is_zero() ? 2 : 1); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, constant.is_zero() ? 0 : 1); BOOST_CHECK_EQUAL(counts.ordinary_products, 0); BOOST_CHECK_EQUAL(counts.preparations, 0); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 0); + } +} + +BOOST_AUTO_TEST_CASE(cyclic_remainders_match_native_division_and_preserve_aliases) { + context_type context {backend_type(256)}; + std::mt19937_64 random(0x46465430); + for (const std::size_t degree : {4, 6, 7, 8, 12, 15, 16, 24, 31, 32, 48, 63, 64}) { + const auto divisor = random_polynomial(degree + 1, random); + auto quotient = random_polynomial(degree - 1, random); + const auto remainder = random_polynomial(degree, random); + const auto prepared = context.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->divisor_degree(), degree); + BOOST_CHECK_EQUAL(prepared->transform_size(), std::bit_ceil(degree)); + for (unsigned repeat = 0; repeat < 3; ++repeat) { + BOOST_TEST_CONTEXT("degree=" << degree << ", repeat=" << repeat) { + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + polynomial_type expected_quotient, expected_remainder, actual; + math::division(expected_quotient, expected_remainder, dividend, divisor); + check_polynomial(expected_quotient, quotient); + check_polynomial(expected_remainder, remainder); + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, expected_remainder); + // All N cells are reconstructed and condensed, so this also checks cancellation above degree(H). + BOOST_CHECK_LT(actual.size(), divisor.size()); + actual = dividend; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, actual, quotient, *prepared)); + check_polynomial(actual, remainder); + actual = quotient; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, actual, *prepared)); + check_polynomial(actual, remainder); + quotient[0] += value_type::one(); + } + } + } +} + +BOOST_AUTO_TEST_CASE(cyclic_remainders_fold_complete_divisors_with_canonical_coordinates) { + context_type context {backend_type(16)}; + coordinates_type coordinates; + coordinates.fill(cpp_int(fq::modulus.backend().to_cpp_int()) - 1); + const polynomial_type divisor(9, field_value(coordinates)), quotient(7, field_value(coordinates)); + const polynomial_type remainder {value_type(19), field_value(coordinates)}; + const auto prepared = context.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 8); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + polynomial_type actual; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, remainder); + + // H=X^8-1 folds to zero, while V=Q*H+R has nonzero high coefficients that must still be folded. + polynomial_type vanishing(9, value_type::zero()), sparse_quotient(7, value_type::zero()); + vanishing[0] = -value_type::one(); + vanishing[8] = sparse_quotient[6] = value_type::one(); + const auto zero_spectrum = context.prepare_cyclic_remainder(vanishing, sparse_quotient.size()); + BOOST_REQUIRE(zero_spectrum.has_value()); + BOOST_CHECK_EQUAL(zero_spectrum->transform_size(), 8); + dividend = expected_product(sparse_quotient, vanishing); + math::addition(dividend, dividend, remainder); + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, sparse_quotient, *zero_spectrum)); + check_polynomial(actual, remainder); + dividend = expected_product(sparse_quotient, vanishing); + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, sparse_quotient, *zero_spectrum)); + check_polynomial(actual, polynomial_type {value_type::zero()}); +} + +BOOST_AUTO_TEST_CASE(cyclic_remainders_reject_unsupported_shapes_without_writes) { + context_type context {backend_type(16)}; + std::mt19937_64 random(0x46465431); + const auto divisor = random_polynomial(9, random); + const auto prepared = context.prepare_cyclic_remainder(divisor, 7); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(divisor, 0).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(divisor, 1).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(polynomial_type {value_type(7)}, 7).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(polynomial_type {value_type(7), value_type(11)}, 7).has_value()); + polynomial_type monomial(9, value_type::zero()); + monomial[8] = value_type::one(); + BOOST_CHECK(!context.prepare_cyclic_remainder(monomial, 7).has_value()); + // A degree-five divisor and four-coefficient quotient already need only an eight-point low product. + BOOST_CHECK(!context.prepare_cyclic_remainder(random_polynomial(6, random), 4).has_value()); + const auto capped = context.prepare_cyclic_remainder(divisor, std::numeric_limits::max()); + BOOST_REQUIRE(capped.has_value()); + BOOST_CHECK_EQUAL(capped->transform_size(), 8); + + const auto reject = [&](const polynomial_type ÷nd, const polynomial_type "ient) { + polynomial_type actual {value_type(99)}; + BOOST_CHECK(!context.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, polynomial_type {value_type(99)}); + actual = dividend; + BOOST_CHECK(!context.try_cyclic_remainder(actual, actual, quotient, *prepared)); + check_polynomial(actual, dividend); + actual = quotient; + BOOST_CHECK(!context.try_cyclic_remainder(actual, dividend, actual, *prepared)); + check_polynomial(actual, quotient); + }; + for (const std::size_t count : {1, 8, 9}) { + const auto quotient = random_polynomial(count, random); + reject(expected_product(quotient, divisor), quotient); + } + reject(random_polynomial(8, random), polynomial_type {value_type::zero()}); + const auto quotient = random_polynomial(7, random); + const auto dividend = expected_product(quotient, divisor); + reject(dividend, random_polynomial(6, random)); + for (const std::size_t size : {4, 8}) { + context_type too_small {backend_type(size)}; + BOOST_CHECK(!too_small.prepare_cyclic_remainder(divisor, quotient.size()).has_value()); + polynomial_type actual = dividend; + BOOST_CHECK(!too_small.try_cyclic_remainder(actual, actual, quotient, *prepared)); + check_polynomial(actual, dividend); + } + polynomial_type actual; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *capped)); + check_polynomial(actual, polynomial_type {value_type::zero()}); +} + +BOOST_AUTO_TEST_CASE(cyclic_preparations_own_spectra_across_contexts_and_moves) { + std::mt19937_64 random(0x46465432); + const auto divisor = random_polynomial(9, random); + const auto quotient = random_polynomial(7, random); + const auto remainder = random_polynomial(8, random); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + const auto prepared = [&] { + context_type temporary {backend_type(64)}; + auto source = divisor; + auto result = temporary.prepare_cyclic_remainder(source, quotient.size()); + source[0] += value_type::one(); + return result; + }(); + BOOST_REQUIRE(prepared.has_value()); + auto copy = *prepared; + const auto moved = std::move(copy); + for (const std::size_t size : {16, 128}) { + context_type context {backend_type(size)}; + for (const auto *snapshot : {&*prepared, &moved}) { + polynomial_type actual; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *snapshot)); + check_polynomial(actual, remainder); + auto changed_divisor = divisor; + changed_divisor[0] += value_type::one(); + const auto changed = context.prepare_cyclic_remainder(changed_divisor, quotient.size()); + BOOST_REQUIRE(changed.has_value()); + const auto changed_dividend = expected_product(quotient, changed_divisor); + BOOST_REQUIRE(context.try_cyclic_remainder(actual, changed_dividend, quotient, *changed)); + check_polynomial(actual, polynomial_type {value_type::zero()}); + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *snapshot)); + check_polynomial(actual, remainder); + } + polynomial_type actual = dividend; + BOOST_CHECK(!context.try_cyclic_remainder(actual, actual, quotient, copy)); + check_polynomial(actual, dividend); + } + backend_type source(32); + backend_type destination = std::move(source); + polynomial_type actual = dividend; + BOOST_CHECK(!source.prepare_cyclic_remainder(divisor, quotient.size()).has_value()); + BOOST_CHECK(!source.try_cyclic_remainder(actual, actual, quotient, *prepared)); + check_polynomial(actual, dividend); + BOOST_REQUIRE(destination.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, remainder); +} + +BOOST_AUTO_TEST_CASE(small_cyclic_preparations_do_not_retain_large_workspaces) { + backend_type backend(256); + polynomial_type actual, large(65, value_type::one()); + backend.square(actual, large); + const polynomial_type divisor(9, value_type::one()), quotient(7, value_type::one()); + const auto prepared = backend.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 8); + const std::size_t bytes_per_position = 12 * backend.prime_count() * sizeof(double); + BOOST_CHECK_GE(prepared->storage_bytes(), 8 * bytes_per_position); + BOOST_CHECK_LT(prepared->storage_bytes(), 256 * bytes_per_position); + auto dividend = expected_product(quotient, divisor); + dividend[0] += value_type(19); + BOOST_REQUIRE(backend.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, polynomial_type {value_type(19)}); +} + +BOOST_AUTO_TEST_CASE(cyclic_remainders_survive_rounding_mode_changes) { + restore_rounding restore; + BOOST_REQUIRE_EQUAL(std::fesetround(FE_TONEAREST), 0); + context_type context {backend_type(16)}; + std::mt19937_64 random(0x46465433); + const auto divisor = random_polynomial(9, random); + const auto quotient = random_polynomial(7, random); + const auto remainder = random_polynomial(8, random); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + for (const int prepare_mode : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO}) { + BOOST_REQUIRE_EQUAL(std::fesetround(prepare_mode), 0); + context_type preparing {backend_type(16)}; + const auto prepared = preparing.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(std::fegetround(), prepare_mode); + for (const int use_mode : {FE_TONEAREST, FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO}) { + BOOST_REQUIRE_EQUAL(std::fesetround(use_mode), 0); + polynomial_type actual; + BOOST_REQUIRE(context.try_cyclic_remainder(actual, dividend, quotient, *prepared)); + check_polynomial(actual, remainder); + BOOST_CHECK_EQUAL(std::fegetround(), use_mode); + } + } +} + +BOOST_AUTO_TEST_CASE(cyclic_divisor_context_falls_back_outside_reduced_product_range) { + prepared_product_counts counts; + arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + arithmetic::polynomial_context context {counting_prepared_backend(64, counts), options}; + std::mt19937_64 random(0x46465434); + const auto divisor = random_polynomial(9, random); + const math::polynomial_divisor_context fixed(divisor, 9, context); + BOOST_CHECK_EQUAL(counts.preparations, 1); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 1); + const auto remainder = random_polynomial(8, random); + for (const std::size_t count : {7, 8, 9, 1, 7}) { + const auto quotient = random_polynomial(count, random); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + counts = {}; + polynomial_type actual_quotient, actual_remainder; + math::divrem(actual_quotient, actual_remainder, dividend, fixed, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + const std::size_t inverse_accepted = count == 9 ? 1 : 0; + const std::size_t cyclic_accepted = count == 7 ? 1 : 0; + BOOST_CHECK_EQUAL(counts.prepared_products, inverse_accepted); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, cyclic_accepted); + BOOST_CHECK_EQUAL(counts.ordinary_products, 2 - inverse_accepted - cyclic_accepted); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + } +} + +BOOST_AUTO_TEST_CASE(divisor_context_keeps_prepared_low_products_when_cyclic_is_not_shorter) { + prepared_product_counts counts; + arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + arithmetic::polynomial_context context {counting_prepared_backend(16, counts), options}; + std::mt19937_64 random(0x46465435); + const auto divisor = random_polynomial(6, random); + const math::polynomial_divisor_context fixed(divisor, 4, context); + BOOST_CHECK_EQUAL(counts.preparations, 2); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 1); + const auto remainder = random_polynomial(5, random); + for (unsigned repeat = 0; repeat < 3; ++repeat) { + const auto quotient = random_polynomial(4, random); + auto dividend = expected_product(quotient, divisor); + math::addition(dividend, dividend, remainder); + counts = {}; + polynomial_type actual_quotient, actual_remainder; + math::divrem(actual_quotient, actual_remainder, dividend, fixed, context); + check_polynomial(actual_quotient, quotient); + check_polynomial(actual_remainder, remainder); + BOOST_CHECK_EQUAL(counts.prepared_products, 2); + BOOST_CHECK_EQUAL(counts.cyclic_remainders + counts.ordinary_products, 0); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); } } From 073909bfaec368060123205e1682768a1c18bac0 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 16:27:24 -0400 Subject: [PATCH 09/11] benchmarks --- libs/math/test/bench_test/CMakeLists.txt | 12 + .../run_small_prime_fft_benchmark.py | 243 ++++++++++++ .../bench_test/small_prime_fft_benchmark.cpp | 357 ++++++++++++++++++ 3 files changed, 612 insertions(+) create mode 100644 libs/math/test/bench_test/run_small_prime_fft_benchmark.py create mode 100644 libs/math/test/bench_test/small_prime_fft_benchmark.cpp diff --git a/libs/math/test/bench_test/CMakeLists.txt b/libs/math/test/bench_test/CMakeLists.txt index 35f70b7f2..f615c54c8 100644 --- a/libs/math/test/bench_test/CMakeLists.txt +++ b/libs/math/test/bench_test/CMakeLists.txt @@ -40,3 +40,15 @@ target_link_libraries(${MIXED_RADIX_BENCHMARK} ${CMAKE_WORKSPACE_NAME}::algebra ${CMAKE_WORKSPACE_NAME}::multiprecision Boost::random) + +# Compare production BN254 Fp12 backends; correctness checks are also available through --check. +set(SMALL_PRIME_FFT_BENCHMARK ${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME}_small_prime_fft_benchmark) +add_executable(${SMALL_PRIME_FFT_BENCHMARK} EXCLUDE_FROM_ALL small_prime_fft_benchmark.cpp) +set_target_properties(${SMALL_PRIME_FFT_BENCHMARK} PROPERTIES CXX_STANDARD 23) +target_compile_options(${SMALL_PRIME_FFT_BENCHMARK} PRIVATE "-march=native") +target_link_libraries(${SMALL_PRIME_FFT_BENCHMARK} + ${CMAKE_WORKSPACE_NAME}::random + ${CMAKE_WORKSPACE_NAME}::math + ${CMAKE_WORKSPACE_NAME}::algebra + ${CMAKE_WORKSPACE_NAME}::multiprecision + Boost::random) diff --git a/libs/math/test/bench_test/run_small_prime_fft_benchmark.py b/libs/math/test/bench_test/run_small_prime_fft_benchmark.py new file mode 100644 index 000000000..04e02cc4c --- /dev/null +++ b/libs/math/test/bench_test/run_small_prime_fft_benchmark.py @@ -0,0 +1,243 @@ +#!/usr/bin/env python3 +"""Package, verify and measure production BN254 Fp12 backends with reproducible logs.""" + +import argparse +import csv +import datetime +import hashlib +import io +import json +import math +import os +from pathlib import Path +import platform +import shlex +import shutil +import statistics +import subprocess +import sys +import tarfile +import tempfile +import time +import traceback + + +BENCH = Path("libs/math/test/bench_test") +SOURCE = BENCH / "small_prime_fft_benchmark.cpp" +RUNNER = BENCH / "run_small_prime_fft_benchmark.py" +MANIFEST = "benchmark-snapshot.json" + + +def timestamp(): + return datetime.datetime.now(datetime.timezone.utc).isoformat() + + +def save(path, value): + path.write_text(json.dumps(value, indent=2, sort_keys=True) + "\n") + + +def digest(path): + checksum = hashlib.sha256() + with path.open("rb") as source: + for block in iter(lambda: source.read(1024 * 1024), b""): + checksum.update(block) + return checksum.hexdigest() + + +def include_roots(repo): + return sorted([*repo.glob("libs/*/include"), *repo.glob("libs/marshalling/*/include")]) + + +def source_metadata(repo): + snapshot = repo / MANIFEST + if snapshot.is_file(): + metadata = json.loads(snapshot.read_text()) + for name, expected in metadata["files"].items(): + if digest(repo / name) != expected: + raise RuntimeError(f"snapshot source changed: {name}") + return metadata + files = {repo / SOURCE, repo / RUNNER, repo / BENCH / "CMakeLists.txt"} + for root in include_roots(repo): + files.update(path for path in root.rglob("*") if path.is_file()) + + def git(*arguments): + return subprocess.check_output(["git", "-C", str(repo), *arguments], text=True) + + if Path(git("rev-parse", "--show-toplevel").strip()).resolve() != repo: + raise RuntimeError("--repo must be a Crypto3 checkout or an exported benchmark snapshot") + return { + "created_utc": timestamp(), + "revision": git("rev-parse", "HEAD").strip(), + "git_status": git("status", "--short", "--untracked-files=all"), + "git_diff": git("diff", "--no-ext-diff", "--binary", "HEAD"), + "files": {str(path.relative_to(repo)): digest(path) for path in sorted(files)}, + } + + +def package(repo, destination): + metadata = source_metadata(repo) + destination.parent.mkdir(parents=True, exist_ok=True) + # Exclusive creation preserves previous handoffs. Every file is rooted under one directory. + with destination.open("xb") as output, tarfile.open(fileobj=output, mode="w:gz") as archive: + prefix = "crypto3-small-prime-benchmark" + for name in metadata["files"]: + archive.add(repo / name, arcname=f"{prefix}/{name}") + if (repo / "LICENSE").is_file(): + archive.add(repo / "LICENSE", arcname=f"{prefix}/LICENSE") + payload = (json.dumps(metadata, indent=2, sort_keys=True) + "\n").encode() + entry = tarfile.TarInfo(f"{prefix}/{MANIFEST}") + entry.size = len(payload) + archive.addfile(entry, io.BytesIO(payload)) + # Do not produce a seemingly valid snapshot if a source changed during packaging. + for name, expected in metadata["files"].items(): + if digest(repo / name) != expected: + raise RuntimeError(f"source changed during packaging; discard {destination}: {name}") + print(f"PACKAGE={destination}\nSHA256={digest(destination)}") + + +def summarize(samples, trials, warmups): + with samples.open(newline="") as source: + rows = list(csv.DictReader(source)) + results = {} + common = {"multiply", "square", "multiply_low", "prepared_low", "divrem", "mulmod", "squaremod"} + for backend in ("mixed_radix", "small_prime"): + backend_rows = [row for row in rows if row["backend"] == backend] + setup = [row for row in backend_rows if row["phase"] == "setup"] + if {row["operation"] for row in setup} != {"backend", "divisor_context", "low_operand", "cyclic_divisor"}: + raise RuntimeError(f"incomplete setup output for {backend}") + cyclic = next(row for row in setup if row["operation"] == "cyclic_divisor") + operations = common | ({"cyclic_remainder"} if int(cyclic["prepared_transform"]) else set()) + results[backend] = {"setup": setup, "operations": {}} + for operation in sorted(operations): + group = [row for row in backend_rows if row["phase"] == "trial" and row["operation"] == operation] + if len(group) != trials or {int(row["trial"]) for row in group} != set(range(warmups, warmups + trials)): + raise RuntimeError(f"incomplete samples: {backend}/{operation}") + measurements = {} + for metric in ("wall_ms", "cpu_ms"): + values = [float(row[metric]) for row in group] + if any(not math.isfinite(value) or value < 0 for value in values): + raise RuntimeError("invalid timing sample") + measurements[metric] = {"median": statistics.median(values), "min": min(values), "max": max(values)} + results[backend]["operations"][operation] = measurements + speedups = {} + for operation in sorted(results["mixed_radix"]["operations"].keys() & results["small_prime"]["operations"].keys()): + baseline = results["mixed_radix"]["operations"][operation]["wall_ms"]["median"] + candidate = results["small_prime"]["operations"][operation]["wall_ms"]["median"] + speedups[operation] = baseline / candidate if candidate else None + ratio = f"{speedups[operation]:.3f}x" if candidate else "below timer resolution" + print(f"{operation:18} mixed_radix={baseline:.6f} ms small_prime={candidate:.6f} ms ratio={ratio}") + return {"backends": results, "wall_time_speedups": speedups} + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--repo", type=Path, default=Path(__file__).resolve().parents[4]) + parser.add_argument("--package", type=Path, help="export source headers and this runner; do not build or measure") + parser.add_argument("--cxx", default="g++") + parser.add_argument("--boost-include", type=Path) + parser.add_argument("--cpu", type=int, help="Linux CPU affinity; default is the highest permitted CPU") + parser.add_argument("--degree", type=int, default=14938) + parser.add_argument("--trials", type=int, default=5) + parser.add_argument("--warmups", type=int, default=2) + parser.add_argument("--timeout", type=int, default=3600, help="limit for each command, in seconds") + parser.add_argument("--check-only", action="store_true") + args = parser.parse_args() + repo = args.repo.resolve() + if not (repo / SOURCE).is_file(): + parser.error("--repo does not contain the benchmark source") + if not 1 <= args.degree <= 17487 or not 1 <= args.trials <= 1000 or not 0 <= args.warmups <= 1000 or args.timeout < 1: + parser.error("degree must be 1..17487, trials 1..1000, warmups 0..1000, timeout positive") + if args.package: + package(repo, args.package.resolve()) + return 0 + permitted = sorted(os.sched_getaffinity(0)) if hasattr(os, "sched_getaffinity") else [] + cpu = args.cpu if args.cpu is not None else (permitted[-1] if permitted else None) + if cpu is not None and (cpu not in permitted or not shutil.which("taskset")): + parser.error("CPU is outside the permitted affinity or taskset is unavailable") + env = dict(os.environ, OMP_NUM_THREADS="1", OPENBLAS_NUM_THREADS="1", MKL_NUM_THREADS="1") + parent = repo / "build" / "small-prime-fft-benchmarks" + parent.mkdir(parents=True, exist_ok=True) + logs = Path(tempfile.mkdtemp(prefix="run-", dir=parent)) + print(f"LOG_DIRECTORY={logs}\nCPU={cpu if cpu is not None else 'not pinned'}", flush=True) + metadata = { + "started_utc": timestamp(), "exit_status": 1, "platform": platform.platform(), + "invocation": sys.argv, "configuration": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()}, + "cpu": cpu, "permitted_cpus": permitted, + "environment": {key: env.get(key) for key in ("PATH", "CPATH", "CPLUS_INCLUDE_PATH", "OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS")}, + "instrumentation": "none", "commands": [], + } + save(logs / "metadata.json", metadata) + + def command(name, argv): + print(f"{name}: {shlex.join(argv)}", flush=True) + record = {"name": name, "argv": argv, "started_utc": timestamp(), "exit_status": None} + metadata["commands"].append(record) + save(logs / "metadata.json", metadata) + started = time.monotonic() + try: + with (logs / f"{name}.stdout.log").open("w") as out, (logs / f"{name}.stderr.log").open("w") as err: + completed = subprocess.run(argv, cwd=repo, env=env, stdout=out, stderr=err, timeout=args.timeout) + record["exit_status"] = completed.returncode + if completed.returncode: + raise RuntimeError(f"{name} failed; see {logs / (name + '.stderr.log')}") + finally: + record["wall_seconds"] = time.monotonic() - started + save(logs / "metadata.json", metadata) + return (logs / f"{name}.stdout.log").read_text() + + try: + sources = source_metadata(repo) + save(logs / "sources.json", sources) + for source in (SOURCE, RUNNER, BENCH / "CMakeLists.txt"): + shutil.copyfile(repo / source, logs / source.name) + command("machine", ["uname", "-a"]) + if sys.platform.startswith("linux"): + command("cpu", ["lscpu"]) + elif sys.platform == "darwin": + command("cpu", ["sysctl", "-n", "machdep.cpu.brand_string"]) + compiler = shlex.split(args.cxx) + version = command("compiler", compiler + ["--version"]) + binary = logs / "small_prime_fft_benchmark" + compile_command = compiler + ["-O3", "-DNDEBUG", "-std=c++23", "-march=native"] + compile_command += ["-fconstexpr-steps=2147483647" if "clang" in version.lower() else "-fconstexpr-ops-limit=2147483647"] + if args.boost_include: + compile_command += ["-isystem", str(args.boost_include.resolve())] + for root in include_roots(repo): + compile_command += ["-I", str(root)] + compile_command += [str(logs / SOURCE.name), "-o", str(binary)] + command("compile", compile_command) + metadata["binary_sha256"] = digest(binary) + prefix = ["taskset", "-c", str(cpu)] if cpu is not None else [] + command("checks", prefix + [str(binary), "--check", "--input-dir", str(logs / "check-inputs")]) + if not args.check_only: + command("benchmark", prefix + [str(binary), "--degree", str(args.degree), "--trials", str(args.trials), + "--warmups", str(args.warmups), "--input-dir", str(logs / "inputs")]) + shutil.copyfile(logs / "benchmark.stdout.log", logs / "samples.csv") + save(logs / "summary.json", summarize(logs / "samples.csv", args.trials, args.warmups)) + for name, expected in sources["files"].items(): + if digest(repo / name) != expected: + raise RuntimeError(f"source changed during the run: {name}") + metadata["exit_status"] = 0 + print("Correctness checks passed." if args.check_only else "All measured outputs verified.", flush=True) + except BaseException: + (logs / "failure.log").write_text(traceback.format_exc()) + raise + finally: + metadata["input_sha256"] = { + str(path.relative_to(logs)): digest(path) + for directory in (logs / "inputs", logs / "check-inputs") if directory.exists() + for path in sorted(directory.glob("*.txt")) + } + metadata["finished_utc"] = timestamp() + save(logs / "metadata.json", metadata) + archive_path = logs.with_suffix(".tar.gz") + with tarfile.open(archive_path, "w:gz") as archive: + for path in sorted(logs.iterdir()): + if path.name != "small_prime_fft_benchmark": + archive.add(path, arcname=f"{logs.name}/{path.name}") + print(f"EXIT_STATUS={metadata['exit_status']}\nRESULT_ARCHIVE={archive_path}", flush=True) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/libs/math/test/bench_test/small_prime_fft_benchmark.cpp b/libs/math/test/bench_test/small_prime_fft_benchmark.cpp new file mode 100644 index 000000000..5ab14a3f0 --- /dev/null +++ b/libs/math/test/bench_test/small_prime_fft_benchmark.cpp @@ -0,0 +1,357 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + namespace math = nil::crypto3::math; + namespace pa = math::polynomial_arithmetic; + namespace fs = std::filesystem; + using small_backend = pa::bn254_small_prime_fft_backend; + using field_type = small_backend::field_type; + using value_type = small_backend::value_type; + using polynomial = small_backend::polynomial_type; + using root_field = nil::crypto3::algebra::fields::alt_bn128_base_field<254>; + using native_backend = pa::mixed_radix_backend; + constexpr std::size_t native_max_order = 34974; + constexpr unsigned seed = 0xC0FFEE; + + void require(bool condition, const char *message) { + if (!condition) { + throw std::runtime_error(message); + } + } + + struct timing { + double wall_ms, cpu_ms; + }; + + template + timing measure(Function &&function) { + const auto cpu_start = std::clock(); + const auto start = std::chrono::steady_clock::now(); + function(); + const auto finish = std::chrono::steady_clock::now(); + const auto cpu_finish = std::clock(); + require(cpu_start != std::clock_t(-1) && cpu_finish != std::clock_t(-1), "CPU clock unavailable"); + return {std::chrono::duration(finish - start).count(), + 1000.0 * (cpu_finish - cpu_start) / CLOCKS_PER_SEC}; + } + + std::size_t native_order(std::size_t required) { + for (std::size_t order = required; order <= native_max_order; ++order) { + if (native_max_order % order == 0) { + return order; + } + } + throw std::invalid_argument("no supported native transform order"); + } + + struct options { + std::size_t degree = 14938, trials = 5, warmups = 2; + bool check = false; + fs::path input_directory; + }; + + // Export canonical coordinates outside timed regions so the runner can hash the exact inputs. + void write_input(const polynomial &input, const fs::path &path) { + std::ofstream output; + output.exceptions(std::ios::failbit | std::ios::badbit); + output.open(path, std::ios::out | std::ios::noreplace); + output << input.size() << '\n' << std::hex; + for (const auto &coefficient : input) { + for (std::size_t c = 0; c < field_type::arity; ++c) { + output << coefficient.data[c / 6].data[(c % 6) / 2].data[c % 2].data << '\n'; + } + } + } + + struct fixture { + polynomial h, a, b, product, square, low, quotient, remainder, square_remainder; + std::string name; + + fixture(const options &settings, bool sparse) : + name(sparse ? "synthetic_nonmonic_sparse" : "synthetic_monic_dense") { + boost::random::mt19937 rng(seed); + const auto generate = [&](std::size_t size) { + polynomial result(size, value_type::zero()); + for (auto &coefficient : result) { + coefficient = nil::crypto3::algebra::random_element(rng); + } + require(!result.back().is_zero(), "zero leading coefficient in deterministic fixture"); + return result; + }; + h = generate(settings.degree + 1); + a = generate(settings.degree); + b = generate(settings.degree); + if (sparse) { + for (std::size_t i = 1; i + 1 < a.size(); ++i) { + a[i] = value_type::zero(); + } + } else { + h.back() = value_type::one(); + } + if (!settings.input_directory.empty()) { + fs::create_directories(settings.input_directory); + const auto prefix = std::to_string(settings.degree) + "-" + name + "-"; + write_input(h, settings.input_directory / (prefix + "h.txt")); + write_input(a, settings.input_directory / (prefix + "a.txt")); + write_input(b, settings.input_directory / (prefix + "b.txt")); + } + } + + template + void reference(pa::polynomial_context &context, bool classical) { + context.multiply(product, a, b); + context.square(square, a); + polynomial square_quotient; + if (classical) { + math::division(quotient, remainder, product, h); + math::division(square_quotient, square_remainder, square, h); + } else { + math::polynomial_divisor_context divisor(h, h.size() - 2, context); + math::divrem(quotient, remainder, product, divisor, context); + math::divrem(square_quotient, square_remainder, square, divisor, context); + } + low = product; + low.resize(a.size()); + math::condense(low); + polynomial reconstructed; + context.multiply(reconstructed, quotient, h); + math::addition(reconstructed, reconstructed, remainder); + require(reconstructed == product && remainder.size() < h.size(), "product division identity failed"); + context.multiply(reconstructed, square_quotient, h); + math::addition(reconstructed, reconstructed, square_remainder); + require(reconstructed == square && square_remainder.size() < h.size(), "square division identity failed"); + for (unsigned point : {0u, 1u, 7u}) { + const value_type x(point); + const auto av = math::evaluate_polynomial(a, x, a.size()); + const auto bv = math::evaluate_polynomial(b, x, b.size()); + require(math::evaluate_polynomial(product, x, product.size()) == av * bv, + "product evaluation identity failed"); + require(math::evaluate_polynomial(square, x, square.size()) == av * av, + "square evaluation identity failed"); + } + } + }; + + enum class operation { multiply, square, low, prepared_low, cyclic, divrem, mulmod, squaremod }; + constexpr std::array operation_names { + "multiply", "square", "multiply_low", "prepared_low", "cyclic_remainder", "divrem", "mulmod", "squaremod"}; + + template + class benchmark_case { + using context_type = pa::polynomial_context; + using divisor_type = math::polynomial_divisor_context; + const fixture &data_; + std::string_view name_; + std::size_t order_; + std::optional context_; + std::optional divisor_; + std::optional low_; + std::optional cyclic_; + polynomial output_, quotient_; + + void print(std::string_view phase, std::size_t trial, std::string_view operation, timing elapsed, + std::size_t transform = 0, std::size_t bytes = 0) const { + std::cout << name_ << ',' << data_.a.size() << ',' << data_.name << ',' << order_ << ',' << phase << ',' + << trial << ',' << operation << ',' << elapsed.wall_ms << ',' << elapsed.cpu_ms << ',' + << transform << ',' << bytes << '\n'; + } + + public: + benchmark_case(const fixture &data, std::string_view name, std::size_t order) : + data_(data), name_(name), order_(order) { + print("setup", 0, "backend", measure([&] { context_.emplace(Backend(order)); })); + print("setup", 0, "divisor_context", + measure([&] { divisor_.emplace(data.h, std::max(std::size_t(1), data.a.size() - 1), *context_); })); + const auto low_time = + measure([&] { low_ = context_->prepare_low_product(data.b, data.a.size(), data.a.size()); }); + require(low_.has_value(), "fixed low-product preparation failed"); + print("setup", 0, "low_operand", low_time, low_->transform_size(), low_->storage_bytes()); + const auto cyclic_time = + measure([&] { cyclic_ = context_->prepare_cyclic_remainder(data.h, data.a.size() - 1); }); + print("setup", 0, "cyclic_divisor", cyclic_time, cyclic_ ? cyclic_->transform_size() : 0, + cyclic_ ? cyclic_->storage_bytes() : 0); + } + + void sample(operation op, std::string_view phase, std::size_t trial) { + const auto label = operation_names[static_cast(op)]; + if (op == operation::cyclic && !cyclic_) { + if (trial == 0) { + print("unsupported", 0, label, {0, 0}); + } + return; + } + bool accepted = true; + const auto elapsed = measure([&] { + switch (op) { + case operation::multiply: + context_->multiply(output_, data_.a, data_.b); + break; + case operation::square: + context_->square(output_, data_.a); + break; + case operation::low: + context_->multiply_low(output_, data_.a, data_.b, data_.a.size()); + break; + case operation::prepared_low: + accepted = context_->try_multiply_low_prepared(output_, data_.a, *low_, data_.a.size()); + break; + case operation::cyclic: + accepted = context_->try_cyclic_remainder(output_, data_.product, data_.quotient, *cyclic_); + break; + case operation::divrem: + math::divrem(quotient_, output_, data_.product, *divisor_, *context_); + break; + case operation::mulmod: + math::mulmod(output_, data_.a, data_.b, *divisor_, *context_); + break; + case operation::squaremod: + math::squaremod(output_, data_.a, *divisor_, *context_); + break; + } + }); + // Consume and check complete outputs outside the timer on every sample, including warmups. + require(accepted, "prepared operation declined"); + const polynomial &expected = op == operation::multiply ? data_.product : + op == operation::square ? data_.square : + op == operation::low || op == operation::prepared_low ? data_.low : + op == operation::squaremod ? data_.square_remainder : + data_.remainder; + require(output_ == expected && !output_.empty() && (output_.size() == 1 || !output_.back().is_zero()), + "benchmark result mismatch"); + if (op == operation::divrem) { + require(quotient_ == data_.quotient, "benchmark quotient mismatch"); + } + print(phase, trial, label, elapsed, + op == operation::prepared_low ? low_->transform_size() : + op == operation::cyclic ? cyclic_->transform_size() : + 0); + } + }; + + void run(const options &settings, bool sparse) { + fixture data(settings, sparse); + const auto order = native_order(2 * settings.degree - 1); + if (settings.degree <= 67) { + pa::polynomial_context> reference; + data.reference(reference, true); + } else { + pa::polynomial_context reference {native_backend(order)}; + data.reference(reference, false); + } + std::cerr << "Verified reference: degree=" << settings.degree << " fixture=" << data.name + << (settings.degree <= 67 ? " schoolbook/native-division\n" : " mixed-radix/identities\n"); + benchmark_case native(data, "mixed_radix", order); + benchmark_case small(data, "small_prime", std::bit_ceil(2 * settings.degree - 1)); + for (std::size_t op = 0; op < operation_names.size(); ++op) { + for (std::size_t trial = 0; trial < settings.warmups + settings.trials; ++trial) { + const auto phase = trial < settings.warmups ? "warmup" : "trial"; + // Alternate the backend order to reduce systematic timing bias. Operations remain serial. + if ((op + trial) % 2 == 0) { + native.sample(static_cast(op), phase, trial); + small.sample(static_cast(op), phase, trial); + } else { + small.sample(static_cast(op), phase, trial); + native.sample(static_cast(op), phase, trial); + } + } + } + std::cerr << "Verified all samples: degree=" << settings.degree << " fixture=" << data.name << '\n'; + } +} // namespace + +int main(int argc, char **argv) { + try { + options settings; + for (int i = 1; i < argc; ++i) { + const std::string_view argument(argv[i]); + if (argument == "--check") { + settings.check = true; + } else if (argument == "--input-dir") { + require(++i < argc, "missing input directory"); + settings.input_directory = argv[i]; + } else { + auto *target = argument == "--degree" ? &settings.degree : + argument == "--trials" ? &settings.trials : + argument == "--warmups" ? &settings.warmups : + nullptr; + require(target && ++i < argc, + "usage: [--check] [--degree N] [--trials N] [--warmups N] [--input-dir DIR]"); + const std::string_view value(argv[i]); + const auto parsed = std::from_chars(value.data(), value.data() + value.size(), *target); + require(parsed.ec == std::errc {} && parsed.ptr == value.data() + value.size(), "invalid number"); + } + } + require(settings.degree >= 1 && settings.degree <= 17487 && settings.trials >= 1 && settings.trials <= 1000 && + settings.warmups <= 1000, + "require degree 1..17487, trials 1..1000, warmups 0..1000"); + std::cerr << "BN254 Fp12; boost=" << BOOST_LIB_VERSION << "; seed=" << seed + << "; serial; synthetic inputs; setup and checks excluded from operation samples.\n"; + std::cout << std::fixed << std::setprecision(6) << std::unitbuf + << "backend,degree,fixture,max_transform,phase,trial,operation,wall_ms,cpu_ms,prepared_transform," + "prepared_bytes\n"; + if (settings.check) { + settings.trials = 1; + settings.warmups = 0; + for (std::size_t degree : {1, 2, 4, 5, 6, 8, 9, 17, 32, 67}) { + settings.degree = degree; + run(settings, false); + run(settings, true); + } + } else { + run(settings, false); + } + return 0; + } catch (const std::exception &error) { + std::cerr << "ERROR: " << error.what() << '\n'; + return 1; + } +} From a936197ff156aa218f925172c46f616c42ac96d1 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 16:34:05 -0400 Subject: [PATCH 10/11] added explanation for primes choice --- .../math/detail/small_prime_fft_params.hpp | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp index 08d4e3b03..d2f2a70df 100644 --- a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp +++ b/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp @@ -47,10 +47,20 @@ namespace nil::crypto3::math::detail { unsigned two_adicity; }; - // Twelve 50-bit primes with large power-of-two factors in modulus - 1 provide CRT capacity - // for full-width BN254 Fp12 products. The tests certify primality and the exact root orders. - // Prime selection uses a prefix of this table. FMA arithmetic and its range checks are separate - // from these integer reconstruction parameters; no floating-point arithmetic is used here. + // These distinct primes q = k*2^s + 1 (k odd) balance residue arithmetic, FFT length and CRT capacity: + // - 2^49 < q < 2^50 matches the FMA kernel's certified range and leaves three headroom bits in + // binary64's 53-bit significand for exact integer-valued lazy sums with magnitude below 8*q. + // Smaller primes contribute fewer CRT bits, generally requiring more primes and transforms + // to reconstruct the same full-width BN254 Fp12 product. + // - s ranges from 41 to 44, providing roots for every power-of-two transform length up to 2^s. + // Entries are ordered by decreasing s, then decreasing q, favoring root capacity first and + // more CRT bits among primes with the same root capacity. + // - Twelve entries provide spare CRT capacity for larger convolutions; the 32768-point + // full-width BN254 Fp12 workload needs eleven. Selection takes only the shortest prefix whose + // product M exceeds twice the signed coefficient bound, so unused entries incur no transforms. + // Other prime/root choices meeting these conditions can also work. Tests certify this catalogue's + // primality, exact root orders and FMA reciprocal error bounds; arithmetic checks remain separate + // from these integer reconstruction parameters, which use no floating-point arithmetic. inline constexpr std::array small_prime_fft_primes = {{ {1108307720798209ULL, 194751219211145ULL, 44}, {659706976665601ULL, 422061946989193ULL, 43}, From d13640b1a1549c7f2191edffe91e710175e4079a Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Thu, 15 Oct 2026 16:40:07 -0400 Subject: [PATCH 11/11] fix test --- .../crypto3/math/detail/bn254_small_prime_fft.hpp | 12 +++++++----- libs/math/test/small_prime_fft_backend.cpp | 2 ++ 2 files changed, 9 insertions(+), 5 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 edc36c245..c50ac286c 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 @@ -755,9 +755,9 @@ namespace nil::crypto3::math::detail { void multiply(std::span left, std::span right, std::size_t stride, double scale) const { check_buffers(left, right, stride, scale); if (arithmetic_.fma_supported()) { - apply(left, right, stride, scale); + apply(left, right, scale); } else { - apply(left, right, stride, scale); + apply(left, right, scale); } } @@ -765,9 +765,9 @@ namespace nil::crypto3::math::detail { void square(std::span values, std::size_t stride, double scale) const { check_buffers(values, values, stride, scale); if (arithmetic_.fma_supported()) { - apply(values, values, stride, scale); + apply(values, values, scale); } else { - apply(values, values, stride, scale); + apply(values, values, scale); } } @@ -782,7 +782,9 @@ namespace nil::crypto3::math::detail { } template - void apply(std::span left, std::span right, std::size_t stride, double scale) const { + void apply(std::span left, std::span right, double scale) const { + // Derive the indexing bound from the buffer extent checked by both public entry points. + const std::size_t stride = left.size() / 12; const kernel operations {arithmetic_}; for (std::size_t i = 0; i < stride; ++i) { tower_type a, b; diff --git a/libs/math/test/small_prime_fft_backend.cpp b/libs/math/test/small_prime_fft_backend.cpp index 43e76d928..d686fcf7a 100644 --- a/libs/math/test/small_prime_fft_backend.cpp +++ b/libs/math/test/small_prime_fft_backend.cpp @@ -281,6 +281,8 @@ BOOST_AUTO_TEST_CASE(pointwise_invalid_dimensions_throw_before_writes) { BOOST_CHECK_THROW(pointwise.multiply(left, std::span(right).first(23), 2, 1.0), std::invalid_argument); BOOST_CHECK_THROW(pointwise.square(std::span(left).first(23), 2, 1.0), std::invalid_argument); + BOOST_CHECK_THROW(pointwise.multiply(left, right, std::numeric_limits::max(), 1.0), + std::invalid_argument); BOOST_CHECK_THROW(pointwise.square(left, std::numeric_limits::max(), 1.0), std::invalid_argument); BOOST_CHECK_EQUAL_COLLECTIONS(left.begin(), left.end(), before.begin(), before.end()); BOOST_CHECK_NO_THROW(pointwise.multiply({}, {}, 0, 1.0));