From 271a937d27b3210521a782b8f3a68f62b0cd7f64 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Fri, 16 Oct 2026 18:02:25 -0400 Subject: [PATCH] optimize coefficent combinations --- cmake/CheckAVX.cmake | 28 +- cmake/CheckHostCPUFeature.cmake | 4 + libs/algebra/CMakeLists.txt | 8 + .../detail/element/fp12_fast/fp12_fast.hpp | 103 ++++ .../detail/element/fp12_fast/limb_ops.hpp | 241 +++++++++- .../detail/element/fp12_fast/limb_ops_x86.hpp | 156 +++++++ libs/algebra/test/fp12.cpp | 355 ++++++++++++++ .../backends/polynomial_backend.hpp | 9 + .../quotient_ring/polynomial_composition.hpp | 442 ++++++++++++++++-- libs/math/test/polynomial_composition.cpp | 422 +++++++++++++++++ 10 files changed, 1737 insertions(+), 31 deletions(-) diff --git a/cmake/CheckAVX.cmake b/cmake/CheckAVX.cmake index e7861c64b..f39da4c2f 100644 --- a/cmake/CheckAVX.cmake +++ b/cmake/CheckAVX.cmake @@ -37,6 +37,19 @@ set(AVX512_CODE " } ") +set(AVX512IFMA_CODE " + #include + int main() { + __m256i a = _mm256_set1_epi64x(1); + a = _mm256_madd52lo_epu64(a, a, a); + a = _mm256_madd52hi_epu64(a, a, a); + __m512i b = _mm512_set1_epi64(1); + b = _mm512_madd52lo_epu64(b, b, b); + b = _mm512_madd52hi_epu64(b, b, b); + return _mm256_extract_epi64(a, 0) + _mm_cvtsi128_si32(_mm512_castsi512_si128(b)); + } +") + macro(check_avx_lang lang type flags) set(__FLAG_I 1) set(CMAKE_REQUIRED_FLAGS_SAVE ${CMAKE_REQUIRED_FLAGS}) @@ -45,7 +58,7 @@ macro(check_avx_lang lang type flags) if(NOT ${lang}_${type}_FOUND) set(CMAKE_REQUIRED_QUIET TRUE) set(CMAKE_REQUIRED_FLAGS ${__FLAG}) - if(lang STREQUAL "CXX") + if("${lang}" STREQUAL "CXX") check_cxx_source_compiles("${${type}_CODE}" ${lang}_HAS_${type}_${__FLAG_I}) else() check_c_source_compiles("${${type}_CODE}" ${lang}_HAS_${type}_${__FLAG_I}) @@ -69,6 +82,19 @@ macro(check_avx_lang lang type flags) endmacro() +# Both widths use the same host contract; VL is required by the four-lane kernel. +macro(check_avx512ifma) + check_host_cpu_feature(AVX512F "avx512f") + check_host_cpu_feature(AVX512IFMA "avx512ifma") + check_host_cpu_feature(AVX512VL "avx512vl") + if(HOST_CPU_AVX512F_FOUND AND HOST_CPU_AVX512IFMA_FOUND AND HOST_CPU_AVX512VL_FOUND) + check_avx_lang(CXX "AVX512IFMA" " ;-mavx512f -mavx512vl -mavx512ifma;/arch:AVX512") + else() + set(CXX_AVX512IFMA_FOUND FALSE CACHE BOOL "CXX AVX512IFMA support" FORCE) + set(CXX_AVX512IFMA_FLAGS "" CACHE STRING "CXX AVX512IFMA flags" FORCE) + endif() +endmacro() + macro(check_avx) check_host_cpu_feature(AVX "avx") check_host_cpu_feature(AVX2 "avx2") diff --git a/cmake/CheckHostCPUFeature.cmake b/cmake/CheckHostCPUFeature.cmake index 71f6ceed9..27e644aa3 100644 --- a/cmake/CheckHostCPUFeature.cmake +++ b/cmake/CheckHostCPUFeature.cmake @@ -30,6 +30,10 @@ function(check_host_cpu_feature feature builtin_name) set(msvc_test "max_leaf >= 7 && (leaf1[2] & (1 << 27)) && (leaf1[2] & (1 << 28)) && ((_xgetbv(0) & 0x6) == 0x6) && (leaf7[1] & (1 << 5))") elseif(feature_upper STREQUAL "AVX512F") set(msvc_test "max_leaf >= 7 && (leaf1[2] & (1 << 27)) && ((_xgetbv(0) & 0xE6) == 0xE6) && (leaf7[1] & (1 << 16))") + elseif(feature_upper STREQUAL "AVX512IFMA") + set(msvc_test "max_leaf >= 7 && (leaf1[2] & (1 << 27)) && ((_xgetbv(0) & 0xE6) == 0xE6) && (leaf7[1] & (1 << 16)) && (leaf7[1] & (1 << 21))") + elseif(feature_upper STREQUAL "AVX512VL") + set(msvc_test "max_leaf >= 7 && (leaf1[2] & (1 << 27)) && ((_xgetbv(0) & 0xE6) == 0xE6) && (leaf7[1] & (1 << 16)) && (leaf7[1] & (1u << 31))") elseif(feature_upper STREQUAL "BMI2") set(msvc_test "max_leaf >= 7 && (leaf7[1] & (1 << 8))") elseif(feature_upper STREQUAL "ADX") diff --git a/libs/algebra/CMakeLists.txt b/libs/algebra/CMakeLists.txt index 7340e3812..b781294b5 100644 --- a/libs/algebra/CMakeLists.txt +++ b/libs/algebra/CMakeLists.txt @@ -14,6 +14,7 @@ cm_project(algebra WORKSPACE_NAME ${CMAKE_WORKSPACE_NAME} LANGUAGES ASM C CXX) include(CMDeploy) include(CheckADX) include(CheckBMI2) +include(CheckAVX) # Blueprint components are using point {0, 0} as a point in infinity, while the standard is {0, 1} # for most curves. We will use compatibility mode with blueprint as default. Setting this flag to @@ -46,6 +47,7 @@ cm_deploy(TARGETS ${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} if(CMAKE_TARGET_ARCHITECTURE STREQUAL "x86_64" OR CMAKE_TARGET_ARCHITECTURE STREQUAL "x86") check_bmi2() check_adx() + check_avx512ifma() if(CXX_BMI2_FOUND) target_compile_definitions(${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} INTERFACE @@ -57,6 +59,12 @@ if(CMAKE_TARGET_ARCHITECTURE STREQUAL "x86_64" OR CMAKE_TARGET_ARCHITECTURE STRE "${CMAKE_UPPER_WORKSPACE_NAME}_HAS_ADX") target_compile_options(${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} INTERFACE ${CXX_ADX_FLAGS}) endif() + if(CXX_AVX512IFMA_FOUND) + target_compile_definitions(${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} INTERFACE + "${CMAKE_UPPER_WORKSPACE_NAME}_HAS_AVX512IFMA") + separate_arguments(CXX_AVX512IFMA_OPTIONS NATIVE_COMMAND "${CXX_AVX512IFMA_FLAGS}") + target_compile_options(${CMAKE_WORKSPACE_NAME}_${CURRENT_PROJECT_NAME} INTERFACE ${CXX_AVX512IFMA_OPTIONS}) + endif() endif() include(CMTest) diff --git a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/fp12_fast.hpp b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/fp12_fast.hpp index 506e9bc4c..6a92fe87f 100644 --- a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/fp12_fast.hpp +++ b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/fp12_fast.hpp @@ -2,7 +2,12 @@ #include #include +#include #include +#include +#include + +#include #include #include @@ -12,6 +17,15 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { template struct fp12_fast { using limb_array = Params::limb_array; + using base_value_type = typename Params::base_field_type::value_type; + using base_backend_type = + std::remove_cvref_t().data.backend().base_data())>; + + // Three Fp6 products, six Fp2 products each, and three base-field products per Fp2 product. + // Every entry is a canonical Montgomery residue. These are reusable linear forms of one operand; + // multiplying corresponding entries and applying recombine_products gives the Fp12 product. + static constexpr std::size_t product_coordinate_count = 54; + using prepared_product_operand = std::array; // Packed input-side Fp2 value: two canonical Fp coefficients. struct fp2_base { @@ -63,6 +77,15 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { fp2_base::add_mod(*this, *this, other); return *this; } + + void to_non_residue(Params::non_residue_type &result) const { + constexpr std::size_t N = Params::base_value_limb_count; + for (std::size_t coefficient = 0; coefficient < 2; ++coefficient) { + for (std::size_t i = 0; i < N; ++i) { + result.data[coefficient].data.backend().base_data().limbs()[i] = data[coefficient * N + i]; + } + } + } }; struct fp2_dbl { @@ -271,6 +294,86 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { } }; + template + static prepared_product_operand prepare_product_operand(const Fp12Value &x) { + using fp2_value_type = typename Params::non_residue_type; + using fp6_value_type = typename Params::underlying_type; + const std::array outer = {x.data[0], x.data[1], x.data[0] + x.data[1]}; + prepared_product_operand result; + for (std::size_t w = 0; w < outer.size(); ++w) { + const auto &a = outer[w].data; + const std::array forms = {a[0], a[1], a[2], a[0] + a[1], a[0] + a[2], a[1] + a[2]}; + for (std::size_t v = 0; v < forms.size(); ++v) { + const auto &b = forms[v].data; + const base_value_type sum = b[0] + b[1]; + const std::size_t offset = 18 * w + 3 * v; + result[offset] = b[0].data.backend().base_data(); + result[offset + 1] = b[1].data.backend().base_data(); + result[offset + 2] = sum.data.backend().base_data(); + } + } + return result; + } + + // This map is linear over Fp, so it also recombines sums of the 54 coordinate products. All corrections + // take place after the dot products, with canonical field arithmetic rather than unbounded lazy values. + template + static Fp12Value recombine_products(const std::array &products) { + using fp2_value_type = typename Params::non_residue_type; + using fp6_value_type = typename Params::underlying_type; + const auto mul_xi = [](const fp2_value_type &x) { + fp2_base packed; + fp2_base::mul_xi(packed, fp2_base(x)); + fp2_value_type result; + packed.to_non_residue(result); + return result; + }; + std::array outer; + for (std::size_t w = 0; w < outer.size(); ++w) { + std::array p; + for (std::size_t v = 0; v < p.size(); ++v) { + const std::size_t offset = 18 * w + 3 * v; + const auto &ac = products[offset]; + const auto &bd = products[offset + 1]; + auto real = ac - bd; + if constexpr (Params::u_squared == -5) { + real -= bd.doubled().doubled(); + } else { + static_assert(Params::u_squared == -1, "unsupported quadratic non-residue"); + } + p[v] = fp2_value_type(real, products[offset + 2] - ac - bd); + } + // p[0..2] are the diagonal Fp2 products; p[3..5] are the pairwise sum products. + outer[w] = fp6_value_type(p[0] + mul_xi(p[5] - p[1] - p[2]), p[3] - p[0] - p[1] + mul_xi(p[2]), + p[4] - p[0] - p[2] + p[1]); + } + const auto &ac = outer[0]; + const auto &bd = outer[1]; + const fp6_value_type v_bd(mul_xi(bd.data[2]), bd.data[0], bd.data[1]); + return Fp12Value(ac + v_bd, outer[2] - ac - bd); + } + + // Preconditions: both ranges contain prepared canonical operands; the right range is at least as long + // as the left. Multipass iterators are required because each coordinate traverses the same operands. + template + static Fp12Value inner_product(LeftIterator first, LeftIterator last, RightIterator right) { + if (first == last) { + return Fp12Value::zero(); + } + std::array products; + for (std::size_t coordinate = 0; coordinate < product_coordinate_count; ++coordinate) { + const auto select = [coordinate](const prepared_product_operand &operand) -> const base_backend_type & { + return operand[coordinate]; + }; + // Reuse the existing 2*N+1-limb accumulator: canonical inputs are < p < R, so any size_t-length + // sum fits. Its final canonicalization handles sums larger than p*R, including 123-term dots. + montgomery_inner_product( + products[coordinate].data.backend().base_data(), boost::make_transform_iterator(first, select), + boost::make_transform_iterator(last, select), boost::make_transform_iterator(right, select)); + } + return recombine_products(products); + } + template static Fp12Value square(const Fp12Value &x) { // Write x = a + b*w in Fp6[w]/(w^2 - v). The quadratic-extension diff --git a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops.hpp b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops.hpp index 9a6bc150e..82ab40e01 100644 --- a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops.hpp +++ b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops.hpp @@ -1,8 +1,10 @@ #pragma once +#include + #include -#if defined(CRYPTO3_HAS_BMI2) && defined(CRYPTO3_HAS_ADX) +#if (defined(CRYPTO3_HAS_BMI2) && defined(CRYPTO3_HAS_ADX)) || defined(CRYPTO3_HAS_AVX512IFMA) #include #endif @@ -22,6 +24,37 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { return load_limbs(Params::base_field_type::modulus_params.get_mod_obj().get_mod()); } + // Sum canonical Montgomery residues, retaining the multiprecision implementation for other limb counts + // and targets. The four-limb kernel changes only product accumulation; reduction still uses the existing + // padded backend and canonicalization, so no dot-length-dependent narrowing is introduced. + template + inline void montgomery_inner_product(Backend &result, LeftIterator first, LeftIterator last, RightIterator right) { + const auto &modulus = Params::base_field_type::modulus_params.get_mod_obj(); +#if defined(CRYPTO3_HAS_BMI2) && defined(CRYPTO3_HAS_ADX) + if constexpr (Params::base_value_limb_count == 4) { + using policy = typename std::remove_cvref_t::policy_type; + using accumulator_type = typename policy::Backend_doubled_padded_limbs; + static_assert(Backend::internal_limb_count == 4 && Backend::limb_bits == 64); + static_assert(sizeof(std::size_t) <= sizeof(limb)); + accumulator_type accumulator; + for (; first != last; ++first, ++right) { + BOOST_ASSERT(eval_lt(*first, modulus.get_mod()) && eval_lt(*right, modulus.get_mod())); + multiply_add_4x4_x86(accumulator.limbs(), (*first).limbs(), (*right).limbs()); + } + // With B=2^64 and R=B^4, at most B-1 canonical products sum to less than (B-1)*p^2. + // REDC adds less than p*R; because p inline void add_limbs_portable(limb *z, const limb *x, const limb *y) { limb carry = 0u; @@ -416,4 +449,210 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { fp2_mul_pre(z, x, y); } + // Five radix-2^52 digits, with Montgomery factor S=2^260 rather than R=2^256. + // Restrict this representation to the four-limb fields whose modulus is below 2^254. + template + requires(Params::base_value_limb_count == 4 && Params::base_field_type::modulus_bits <= 254) + struct montgomery_radix52 { + static constexpr limb mask = (limb(1) << 52) - 1; + static constexpr std::size_t digits = 5; + static constexpr std::size_t max_terms = 256; + using value_type = std::array; + + static value_type pack(const limb *x) { + return {x[0] & mask, ((x[0] >> 52) | (x[1] << 12)) & mask, ((x[1] >> 40) | (x[2] << 24)) & mask, + ((x[2] >> 28) | (x[3] << 36)) & mask, x[3] >> 16}; + } + + static const value_type &modulus() { + static const value_type p = [] { + const auto words = modulus_limbs(); + return pack(words.data()); + }(); + return p; + } + + static limb p_dash() { + return Params::base_field_type::modulus_params.get_mod_obj().get_p_dash() & mask; + } + + template + static value_type from_montgomery64(const Backend &x) { + auto words = load_limbs(x); + // S/R=16: four canonical doublings convert the Montgomery scaling exactly. + for (unsigned i = 0; i < 4; ++i) { + add_low_limbs_mod_portable(words.data(), words.data(), words.data()); + } + return pack(words.data()); + } + + template + static void to_montgomery64(Backend &result, value_type x) { + // Choose m in [0,15] so x+m*p is divisible by 16. Since x> 52); + } + BOOST_ASSERT(carry == 0); + for (std::size_t i = 0; i < digits - 1; ++i) { + x[i] = (x[i] >> 4) | ((x[i + 1] & 15) << 48); + } + x[4] >>= 4; + result = Backend(); + result.limbs()[0] = x[0] | (x[1] << 52); + result.limbs()[1] = (x[1] >> 12) | (x[2] << 40); + result.limbs()[2] = (x[2] >> 24) | (x[3] << 28); + result.limbs()[3] = (x[3] >> 36) | (x[4] << 16); + } + + static void subtract_modulus(value_type &x) { + value_type difference; + limb borrow = 0; + for (std::size_t i = 0; i < digits; ++i) { + const limb subtrahend = modulus()[i] + borrow; + difference[i] = (x[i] - subtrahend) & mask; + borrow = x[i] < subtrahend; + } + if (!borrow) { + x = difference; + } + } + + // Scalar reference: full 104-bit products in 128-bit columns, independently of IFMA's split products. + // At most 256 terms contribute five products to a column, requiring fewer than 115 bits. + template + static value_type inner_product(LeftIterator first, LeftIterator last, RightIterator right) { + value_type result {}; + while (first != last) { + std::array t {}; + std::size_t count = 0; + for (; first != last && count < max_terms; ++first, ++right, ++count) { + for (std::size_t i = 0; i < digits; ++i) { + for (std::size_t j = 0; j < digits; ++j) { + t[i + j] += wide_limb((*first)[i]) * (*right)[j]; + } + } + } + for (std::size_t i = 0; i < digits; ++i) { + const limb m = (limb(t[i]) * p_dash()) & mask; + for (std::size_t j = 0; j < digits; ++j) { + t[i + j] += wide_limb(m) * modulus()[j]; + } + t[i + 1] += t[i] >> 52; + } + value_type reduced; + for (std::size_t i = 5; i < 9; ++i) { + reduced[i - 5] = limb(t[i]) & mask; + t[i + 1] += t[i] >> 52; + } + reduced[4] = limb(t[9]); + // T<256*p^2 and p<2^254 imply REDC(T)<5*p> 52; + } + BOOST_ASSERT(carry == 0); + subtract_modulus(result); + } + return result; + } + }; + + // Left is count consecutive five-digit scalars. Right is [term][digit][lane]; indices selects + // nonzero outer terms without repacking the cached panels. All digits are canonical S-Montgomery inputs. + // Result has [digit][lane] layout. Padded coefficient lanes are zero and never read past a panel. + template + inline void montgomery_inner_product52(limb *result, const limb *left, const limb *right, + const std::size_t *indices, std::size_t count) { + using radix = montgomery_radix52; + using vector = typename Operations::value_type; + constexpr std::size_t lanes = Operations::lanes; + const auto zero = Operations::broadcast(0); + const auto mask = Operations::broadcast(radix::mask); + const auto one = Operations::broadcast(1); + const auto inverse = Operations::broadcast(radix::p_dash()); + std::array modulus, total; + for (std::size_t i = 0; i < 5; ++i) { + modulus[i] = Operations::broadcast(radix::modulus()[i]); + total[i] = zero; + } + const auto subtract = [&](std::array &x) { + std::array difference; + auto borrow = zero; + for (std::size_t i = 0; i < 5; ++i) { + const auto subtrahend = Operations::add(modulus[i], borrow); + difference[i] = Operations::bit_and(Operations::subtract(x[i], subtrahend), mask); + borrow = Operations::select(Operations::less(x[i], subtrahend), one, zero); + } + const auto keep = Operations::less(zero, borrow); + for (std::size_t i = 0; i < 5; ++i) { + x[i] = Operations::select(keep, x[i], difference[i]); + } + }; + for (std::size_t first = 0; first < count;) { + const std::size_t terms = std::min(radix::max_terms, count - first); + std::array t; + t.fill(zero); + // Each column receives at most 9*256 split products, each <2^52: below 2^64. + // Normalize between chunks, never permitting an unbounded IFMA accumulator. + for (std::size_t term = first; term < first + terms; ++term) { + const limb *y = right + indices[term] * 5 * lanes; + std::array b; + for (std::size_t j = 0; j < 5; ++j) { + b[j] = Operations::load(y + j * lanes); + } + for (std::size_t i = 0; i < 5; ++i) { + const auto a = Operations::broadcast(left[5 * term + i]); + for (std::size_t j = 0; j < 5; ++j) { + t[i + j] = Operations::multiply_add_low(t[i + j], a, b[j]); + t[i + j + 1] = Operations::multiply_add_high(t[i + j + 1], a, b[j]); + } + } + } + for (std::size_t i = 0; i < 9; ++i) { + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + t[i] = Operations::bit_and(t[i], mask); + } + // Five cancellation steps divide by S=2^260. IFMA ignores upper operand bits, which is + // exactly the reduction modulo 2^52 needed for m; the carry of t[i] goes to t[i+1]. + for (std::size_t i = 0; i < 5; ++i) { + const auto m = Operations::multiply_add_low(zero, t[i], inverse); + for (std::size_t j = 0; j < 5; ++j) { + t[i + j] = Operations::multiply_add_low(t[i + j], m, modulus[j]); + t[i + j + 1] = Operations::multiply_add_high(t[i + j + 1], m, modulus[j]); + } + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + std::array reduced; + for (std::size_t i = 5; i < 9; ++i) { + reduced[i - 5] = Operations::bit_and(t[i], mask); + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + reduced[4] = t[9]; + for (unsigned i = 0; i < 4; ++i) { + subtract(reduced); + } + auto carry = zero; + for (std::size_t i = 0; i < 5; ++i) { + const auto sum = Operations::add(Operations::add(total[i], reduced[i]), carry); + total[i] = Operations::bit_and(sum, mask); + carry = Operations::shift52(sum); + } + subtract(total); + first += terms; + } + for (std::size_t i = 0; i < 5; ++i) { + Operations::store(result + i * lanes, total[i]); + } + } + } // namespace nil::crypto3::algebra::fields::detail::fp12_fast diff --git a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp index b26f43b20..2573b40b0 100644 --- a/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp +++ b/libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops_x86.hpp @@ -3,6 +3,10 @@ #include #include +#if defined(CRYPTO3_HAS_AVX512IFMA) +#include +#endif + // clang-format off #define STR_IMPL(X) #X @@ -51,6 +55,29 @@ "mov " D4(2, 4) ", " PTR2(Z, Z_BASE, 6) "\n" \ "mov " D4(3, 4) ", " PTR2(Z, Z_BASE, 7) "\n" +// Finish one product limb before adding it to the nine-limb accumulator. The product's CF/OF chains +// and the accumulator's carry are independent: save the latter in pending between schoolbook rounds. +#define MULTIPLY_ADD_4x4_ROUND(ROUND) \ + "xor %[zero], %[zero]\n" \ + "mov " PTR(y, ROUND) ", %%rdx\n" \ + "mulx " PTR(x, 0) ", %[low], %[high]\n" \ + "adox %[low], " D4(0, ROUND) "\n" \ + "mov " D4(0, ROUND) ", %[saved]\n" \ + "adcx %[high], " D4(1, ROUND) "\n" \ + "mulx " PTR(x, 1) ", %[low], %[high]\n" \ + "adox %[low], " D4(1, ROUND) "\n" \ + "adcx %[high], " D4(2, ROUND) "\n" \ + "mulx " PTR(x, 2) ", %[low], %[high]\n" \ + "adox %[low], " D4(2, ROUND) "\n" \ + "adcx %[high], " D4(3, ROUND) "\n" \ + "mulx " PTR(x, 3) ", %[low], " D4(4, ROUND) "\n" \ + "adox %[low], " D4(3, ROUND) "\n" \ + "adox %[zero], " D4(4, ROUND) "\n" \ + "adcx %[zero], " D4(4, ROUND) "\n" \ + "btq $0, %[pending]\n" \ + "adcq %[saved], " PTR(z, ROUND) "\n" \ + "setc %b[pending]\n" + // Get the (i + j) mod 6 d register in the 6-live-limb schoolbook window. #define D6(I, J) "%[d" STR(BOOST_PP_MOD(BOOST_PP_ADD(I, J), 6)) "]" @@ -628,6 +655,45 @@ "adcq %[" #T5 "], " PTR2(Z, Z_BASE, 11) "\n" namespace nil::crypto3::algebra::fields::detail::fp12_fast { + // z += x*y, with four unsigned 64-bit input limbs and nine accumulator limbs. z must be disjoint + // from both inputs, and the exact sum must fit in nine limbs. x and y may alias each other. + template + inline void multiply_add_4x4_x86(Limb *z, const Limb *x, const Limb *y) { + static_assert(std::is_unsigned_v && sizeof(Limb) == sizeof(limb)); + limb d0, d1, d2, d3, low, high, saved, pending, zero; + asm volatile( + "mov " PTR(y, 0) ", %%rdx\n" + "mulx " PTR(x, 0) ", %[d0], %[d1]\n" + "mov %[d0], %[saved]\n" + "mulx " PTR(x, 1) ", %[low], %[d2]\n" + "add %[low], %[d1]\n" + "mulx " PTR(x, 2) ", %[low], %[d3]\n" + "adc %[low], %[d2]\n" + "mulx " PTR(x, 3) ", %[low], %[d0]\n" + "adc %[low], %[d3]\n" + "adc $0, %[d0]\n" + "addq %[saved], " PTR(z, 0) "\n" + "setc %b[pending]\n" + + MULTIPLY_ADD_4x4_ROUND(1) + MULTIPLY_ADD_4x4_ROUND(2) + MULTIPLY_ADD_4x4_ROUND(3) + + // The last round leaves the accumulator carry in CF. Propagate it through the high product + // limbs and into the extra limb, without truncating sums that exceed 512 bits. + "adcq %[d0], " PTR(z, 4) "\n" + "adcq %[d1], " PTR(z, 5) "\n" + "adcq %[d2], " PTR(z, 6) "\n" + "adcq %[d3], " PTR(z, 7) "\n" + "adcq $0, " PTR(z, 8) "\n" + : [d0]"=&r"(d0), [d1]"=&r"(d1), [d2]"=&r"(d2), [d3]"=&r"(d3), + [low]"=&r"(low), [high]"=&r"(high), [saved]"=&r"(saved), + [pending]"=&r"(pending), [zero]"=&r"(zero) + : [z]"r"(z), [x]"r"(x), [y]"r"(y) + : "rdx", "cc", "memory" + ); + } + template requires(Params::storage_limb_count == 8) inline void montgomery_reduce_x86(limb *result, const limb *data) { @@ -1297,6 +1363,7 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { #undef D6 #undef SCHOOLBOOK_4x4_ROUND #undef SCHOOLBOOK_4x4 +#undef MULTIPLY_ADD_4x4_ROUND #undef SCHOOLBOOK_6x6_ROUND #undef SCHOOLBOOK_6x6 #undef DOUBLE_MOD_P @@ -1325,3 +1392,92 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast { #undef SUB_12_LIMBS_MOD // clang-format on + +#if defined(CRYPTO3_HAS_AVX512IFMA) +namespace nil::crypto3::algebra::fields::detail::fp12_fast { + template + struct ifma52_operations; + + template<> + struct ifma52_operations<4> { + static constexpr std::size_t lanes = 4; + struct value_type { + __m256i data; + }; + static value_type broadcast(limb x) { + return {_mm256_set1_epi64x(x)}; + } + static value_type load(const limb *x) { + return {_mm256_loadu_si256(reinterpret_cast(x))}; + } + static void store(limb *x, value_type a) { + _mm256_storeu_si256(reinterpret_cast<__m256i *>(x), a.data); + } + static value_type add(value_type a, value_type b) { + return {_mm256_add_epi64(a.data, b.data)}; + } + static value_type subtract(value_type a, value_type b) { + return {_mm256_sub_epi64(a.data, b.data)}; + } + static value_type bit_and(value_type a, value_type b) { + return {_mm256_and_si256(a.data, b.data)}; + } + static value_type shift52(value_type a) { + return {_mm256_srli_epi64(a.data, 52)}; + } + static __mmask8 less(value_type a, value_type b) { + return _mm256_cmp_epu64_mask(a.data, b.data, _MM_CMPINT_LT); + } + static value_type select(__mmask8 mask, value_type yes, value_type no) { + return {_mm256_mask_blend_epi64(mask, no.data, yes.data)}; + } + static value_type multiply_add_low(value_type sum, value_type a, value_type b) { + return {_mm256_madd52lo_epu64(sum.data, a.data, b.data)}; + } + static value_type multiply_add_high(value_type sum, value_type a, value_type b) { + return {_mm256_madd52hi_epu64(sum.data, a.data, b.data)}; + } + }; + + template<> + struct ifma52_operations<8> { + static constexpr std::size_t lanes = 8; + struct value_type { + __m512i data; + }; + static value_type broadcast(limb x) { + return {_mm512_set1_epi64(x)}; + } + static value_type load(const limb *x) { + return {_mm512_loadu_si512(x)}; + } + static void store(limb *x, value_type a) { + _mm512_storeu_si512(x, a.data); + } + static value_type add(value_type a, value_type b) { + return {_mm512_add_epi64(a.data, b.data)}; + } + static value_type subtract(value_type a, value_type b) { + return {_mm512_sub_epi64(a.data, b.data)}; + } + static value_type bit_and(value_type a, value_type b) { + return {_mm512_and_si512(a.data, b.data)}; + } + static value_type shift52(value_type a) { + return {_mm512_srli_epi64(a.data, 52)}; + } + static __mmask8 less(value_type a, value_type b) { + return _mm512_cmp_epu64_mask(a.data, b.data, _MM_CMPINT_LT); + } + static value_type select(__mmask8 mask, value_type yes, value_type no) { + return {_mm512_mask_blend_epi64(mask, no.data, yes.data)}; + } + static value_type multiply_add_low(value_type sum, value_type a, value_type b) { + return {_mm512_madd52lo_epu64(sum.data, a.data, b.data)}; + } + static value_type multiply_add_high(value_type sum, value_type a, value_type b) { + return {_mm512_madd52hi_epu64(sum.data, a.data, b.data)}; + } + }; +} // namespace nil::crypto3::algebra::fields::detail::fp12_fast +#endif diff --git a/libs/algebra/test/fp12.cpp b/libs/algebra/test/fp12.cpp index 99948fbd3..342ef634b 100644 --- a/libs/algebra/test/fp12.cpp +++ b/libs/algebra/test/fp12.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -149,6 +150,148 @@ namespace { return value; } + // Execute the vector algorithm on scalar lanes as well, independently of the 128-bit-column reference. + // These operations model IFMA's ignored upper operand bits and reject any accumulator overflow. + template + struct scalar_ifma52_operations { + using limb = fields::detail::fp12_fast::limb; + using wide = fields::detail::fp12_fast::wide_limb; + static constexpr std::size_t lanes = Lanes; + static constexpr limb mask = (limb(1) << 52) - 1; + using value_type = std::array; + using mask_type = std::array; + static value_type broadcast(limb x) { + value_type r; + r.fill(x); + return r; + } + static value_type load(const limb *x) { + value_type r; + std::copy_n(x, lanes, r.begin()); + return r; + } + static void store(limb *x, const value_type &a) { + std::copy(a.begin(), a.end(), x); + } + static value_type add(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) { + if (a[i] > std::numeric_limits::max() - b[i]) { + throw std::overflow_error("IFMA lane overflow"); + } + a[i] += b[i]; + } + return a; + } + static value_type subtract(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) + a[i] -= b[i]; + return a; + } + static value_type bit_and(value_type a, const value_type &b) { + for (std::size_t i = 0; i < lanes; ++i) + a[i] &= b[i]; + return a; + } + static value_type shift52(value_type a) { + for (auto &x : a) + x >>= 52; + return a; + } + static mask_type less(const value_type &a, const value_type &b) { + mask_type r; + for (std::size_t i = 0; i < lanes; ++i) + r[i] = a[i] < b[i]; + return r; + } + static value_type select(const mask_type &mask, value_type yes, const value_type &no) { + for (std::size_t i = 0; i < lanes; ++i) + if (!mask[i]) + yes[i] = no[i]; + return yes; + } + static value_type multiply_add_low(value_type sum, const value_type &a, const value_type &b) { + value_type product; + for (std::size_t i = 0; i < lanes; ++i) + product[i] = limb(wide(a[i] & mask) * (b[i] & mask)) & mask; + return add(sum, product); + } + static value_type multiply_add_high(value_type sum, const value_type &a, const value_type &b) { + value_type product; + for (std::size_t i = 0; i < lanes; ++i) + product[i] = limb((wide(a[i] & mask) * (b[i] & mask)) >> 52); + return add(sum, product); + } + }; + + template + void check_radix52_dots() { + namespace fast = fields::detail::fp12_fast; + using params = alt_bn128_254_fp12::extension_policy::fast_params; + using kernel = fast::fp12_fast; + using radix = fast::montgomery_radix52; + using backend = kernel::base_backend_type; + constexpr std::size_t lanes = Operations::lanes, terms = 1025; + const auto &mod = params::base_field_type::modulus_params.get_mod_obj(); + auto maximum = radix::modulus(); + --maximum[0]; + boost::random::mt19937 rng(0x52F00D); + for (bool maximal : {true, false}) { + std::vector left(terms); + std::array, lanes> right; + std::vector packed_left(terms * 5), packed_right(terms * 5 * lanes); + std::vector indices(terms); + for (auto &lane : right) + lane.resize(terms); + for (std::size_t j = 0; j < terms; ++j) { + left[j] = + maximal ? + maximum : + radix::from_montgomery64( + random_fp12(rng).data[0].data[0].data[0].data.backend().base_data()); + indices[j] = terms - j - 1; + for (std::size_t lane = 0; lane < lanes; ++lane) { + right[lane][j] = lane == 0 ? radix::value_type {} : + maximal ? maximum : + radix::from_montgomery64(random_fp12(rng) + .data[0] + .data[0] + .data[0] + .data.backend() + .base_data()); + for (std::size_t i = 0; i < 5; ++i) { + packed_left[j * 5 + i] = left[j][i]; + packed_right[(indices[j] * 5 + i) * lanes + lane] = right[lane][j][i]; + } + } + } + for (std::size_t count : {0, 1, 2, 7, 8, 9, 27, 28, 123, 255, 256, 257, 511, 512, 513, 1025}) { + std::array actual; + actual.fill(0xAABBCCDDEEFF0011); + fast::montgomery_inner_product52(actual.data() + 1, packed_left.data(), + packed_right.data(), indices.data(), count); + BOOST_CHECK_EQUAL(actual.front(), 0xAABBCCDDEEFF0011); + BOOST_CHECK_EQUAL(actual.back(), 0xAABBCCDDEEFF0011); + for (std::size_t lane = 0; lane < lanes; ++lane) { + radix::value_type result; + for (std::size_t i = 0; i < 5; ++i) + result[i] = actual[1 + i * lanes + lane]; + const auto expected = radix::inner_product(left.begin(), left.begin() + count, right[lane].begin()); + BOOST_CHECK(result == expected); + std::vector l(count), r(count); + for (std::size_t j = 0; j < count; ++j) { + radix::to_montgomery64(l[j], left[j]); + radix::to_montgomery64(r[j], right[lane][j]); + } + backend original, converted; + mod.montgomery_inner_product(original, l.begin(), l.end(), r.begin()); + radix::to_montgomery64(converted, result); + BOOST_CHECK(eval_eq(original, converted)); + BOOST_CHECK(eval_lt(converted, mod.get_mod())); + } + } + } + } + } // namespace BOOST_AUTO_TEST_SUITE(fp12_tests) @@ -334,6 +477,218 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(multiplication_matches_generic_implementation, Fp1 } } +BOOST_AUTO_TEST_CASE_TEMPLATE(prepared_product_forms_recombine_exactly, Fp12Field, fp12_field_types) { + using kernel = fields::detail::fp12_fast::fp12_fast; + using base_value_type = typename kernel::base_value_type; + boost::random::mt19937 rng(0x5412); + const auto check = [&](const auto &x, const auto &y) { + const auto left = kernel::prepare_product_operand(x); + const auto right = kernel::prepare_product_operand(y); + std::array products; + for (std::size_t i = 0; i < products.size(); ++i) { + const auto &modulus = base_value_type::modulus_params.get_mod_obj().get_mod(); + BOOST_CHECK(eval_lt(left[i], modulus)); + BOOST_CHECK(eval_lt(right[i], modulus)); + base_value_type a, b; + a.data.backend().base_data() = left[i]; + b.data.backend().base_data() = right[i]; + products[i] = a * b; + } + BOOST_CHECK_EQUAL(kernel::template recombine_products(products), + generic_fp12_mul(x, y)); + }; + for (std::size_t i = 0; i < 6; ++i) { + for (std::size_t j = 0; j < 6; ++j) { + check(boundary_fp12(i), boundary_fp12(j)); + } + } + for (std::size_t i = 0; i < random_samples; ++i) { + check(random_fp12(rng), random_fp12(rng)); + } + check(Fp12Field::value_type::zero(), Fp12Field::value_type::one()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(prepared_inner_products_match_generic_arithmetic, Fp12Field, fp12_field_types) { + using value_type = typename Fp12Field::value_type; + using kernel = fields::detail::fp12_fast::fp12_fast; + using base_value_type = typename kernel::base_value_type; + using prepared_type = typename kernel::prepared_product_operand; + boost::random::mt19937 rng(0x12354); + std::vector left, right; + std::vector prepared_left, prepared_right; + value_type maximum_residues; + auto maximum = base_value_type::modulus_params.get_mod_obj().get_mod(); + typename kernel::base_backend_type one; + one.limbs()[0] = 1; + eval_subtract(maximum, one); + for (std::size_t i = 0; i < Fp12Field::arity; ++i) { + maximum_residues.coordinate(i).data.backend().base_data() = maximum; + } + for (std::size_t i = 0; i < 513; ++i) { + left.push_back(i % 3 == 0 ? maximum_residues : random_fp12(rng)); + right.push_back(i % 3 == 0 ? maximum_residues : boundary_fp12(i % 6)); + prepared_left.push_back(kernel::prepare_product_operand(left.back())); + prepared_right.push_back(kernel::prepare_product_operand(right.back())); + } + for (std::size_t count : {0, 1, 2, 3, 4, 5, 7, 8, 9, 16, 31, 32, 33, 123, 124, 256, 257, 513}) { + value_type expected = value_type::zero(); + for (std::size_t i = 0; i < count; ++i) { + expected += generic_fp12_mul(left[i], right[i]); + } + const auto actual = kernel::template inner_product( + prepared_left.begin(), prepared_left.begin() + count, prepared_right.begin()); + BOOST_CHECK_EQUAL(actual, expected); + } + + // Cancellation across a long dot must remain canonical, including the empty result. + for (std::size_t i = 0; i < 246; i += 2) { + prepared_left[i + 1] = kernel::prepare_product_operand(-left[i]); + prepared_right[i + 1] = prepared_right[i]; + } + BOOST_CHECK_EQUAL(kernel::template inner_product(prepared_left.begin(), prepared_left.begin() + 246, + prepared_right.begin()), + value_type::zero()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(prepared_base_dots_match_portable_with_wide_carries, Fp12Field, fp12_field_types) { + using params = typename Fp12Field::extension_policy::fast_params; + using kernel = fields::detail::fp12_fast::fp12_fast; + using backend = typename kernel::base_backend_type; + const auto &modulus = params::base_field_type::modulus_params.get_mod_obj(); + backend maximum = modulus.get_mod(), one; + one.limbs()[0] = 1; + eval_subtract(maximum, one); + std::vector left(1025, maximum), right(1025, maximum); + // For BN254, 28 products of (p-1)^2 first carry beyond the eight ordinary product limbs. + // Also exercise several REDC canonicalization subtractions and the six-limb portable fallback. + for (std::size_t count : {0, 1, 2, 3, 27, 28, 29, 31, 32, 33, 122, 123, 124, 255, 256, 257, 513, 1025}) { + backend actual, expected; + modulus.montgomery_inner_product(expected, left.begin(), left.begin() + count, right.begin()); + fields::detail::fp12_fast::montgomery_inner_product(actual, left.begin(), left.begin() + count, + right.begin()); + BOOST_CHECK(eval_eq(actual, expected)); + BOOST_CHECK(eval_lt(actual, modulus.get_mod())); + // Writing the result must wait until all input terms have been consumed. + fields::detail::fp12_fast::montgomery_inner_product(left[0], left.begin(), left.begin() + count, + right.begin()); + BOOST_CHECK(eval_eq(left[0], expected)); + left[0] = maximum; + } +} + +BOOST_AUTO_TEST_CASE(radix52_conversion_preserves_canonical_montgomery_values) { + namespace fast = fields::detail::fp12_fast; + using params = alt_bn128_254_fp12::extension_policy::fast_params; + using radix = fast::montgomery_radix52; + using backend = fast::fp12_fast::base_backend_type; + using boost::multiprecision::cpp_int; + const auto integer = [](const auto &x, unsigned bits) { + cpp_int result = 0; + for (auto i = x.rbegin(); i != x.rend(); ++i) { + result <<= bits; + result += *i; + } + return result; + }; + const auto &mod = params::base_field_type::modulus_params.get_mod_obj(); + const cpp_int p = integer(fast::load_limbs(mod.get_mod()), 64); + boost::random::mt19937 rng(0x5264); + for (std::size_t i = 0; i < 300; ++i) { + backend original = i < 16 ? + mod.get_mod() : + random_fp12(rng).data[0].data[0].data[0].data.backend().base_data(); + if (i < 16) { + backend delta; + delta.limbs()[0] = i + 1; + eval_subtract(original, delta); + } + if (i == 16) + original = backend(); + const auto packed = radix::from_montgomery64(original); + BOOST_CHECK(integer(packed, 52) == (integer(fast::load_limbs(original), 64) * 16) % p); + for (auto digit : packed) + BOOST_CHECK_LE(digit, radix::mask); + backend restored; + radix::to_montgomery64(restored, packed); + BOOST_CHECK(eval_eq(restored, original)); + } +} + +BOOST_AUTO_TEST_CASE(radix52_four_and_eight_lane_arithmetic_matches_reference) { + check_radix52_dots>(); + check_radix52_dots>(); +} + +#if defined(CRYPTO3_HAS_AVX512IFMA) +BOOST_AUTO_TEST_CASE(radix52_ifma_matches_reference_in_every_lane) { + check_radix52_dots>(); + check_radix52_dots>(); +} +#endif + +#if defined(CRYPTO3_HAS_BMI2) && defined(CRYPTO3_HAS_ADX) +BOOST_AUTO_TEST_CASE(four_limb_multiply_accumulate_preserves_all_nine_limbs) { + using limb = fields::detail::fp12_fast::limb; + using boost::multiprecision::cpp_int; + constexpr limb maximum = std::numeric_limits::max(); + const auto integer = [](const limb *data, std::size_t size) { + cpp_int result = 0; + while (size > 0) { + result <<= 64; + result += data[--size]; + } + return result; + }; + const auto check = [&](const std::array &left, const std::array &right, + const std::array &initial, bool square = false) { + std::array storage; + storage.front() = 0xAABBCCDDEEFF0011; + storage.back() = 0x1122334455667788; + std::copy(initial.begin(), initial.end(), storage.begin() + 1); + const auto *right_data = square ? left.data() : right.data(); + const cpp_int expected = integer(initial.data(), initial.size()) + + integer(left.data(), left.size()) * integer(right_data, right.size()); + fields::detail::fp12_fast::multiply_add_4x4_x86(storage.data() + 1, left.data(), right_data); + BOOST_CHECK_EQUAL(integer(storage.data() + 1, 9), expected); + BOOST_CHECK_EQUAL(storage.front(), 0xAABBCCDDEEFF0011); + BOOST_CHECK_EQUAL(storage.back(), 0x1122334455667788); + }; + const std::array full {maximum, maximum, maximum, maximum}, zero {}; + std::array initial; + initial.fill(maximum); + check(zero, full, initial); + initial.back() = maximum - 1; // Carry into a nearly full ninth limb without overflowing it. + check(full, full, initial); + check(full, full, initial, true); + initial.back() = 0; + check(full, full, initial); + check(full, full, {}); + for (std::size_t i = 0; i < 4; ++i) { + for (std::size_t j = 0; j < 4; ++j) { + for (unsigned bit : {0, 1, 31, 32, 63}) { + std::array left {}, right {}; + left[i] = limb(1) << bit; + right[j] = maximum; + check(left, right, initial); + } + } + } + boost::random::mt19937 rng(0x4ADC); + const auto random_limb = [&]() { return (limb(rng()) << 32) | rng(); }; + for (std::size_t sample = 0; sample < 1024; ++sample) { + std::array left, right; + for (auto &value : left) + value = random_limb(); + for (auto &value : right) + value = random_limb(); + for (auto &value : initial) + value = random_limb(); + initial.back() >>= 1; // Leave room for the product's carry. + check(left, right, initial, sample % 2 == 0); + } +} +#endif + BOOST_AUTO_TEST_CASE_TEMPLATE(square_matches_generic_implementation, Fp12Field, fp12_field_types) { using fp12_value_type = typename Fp12Field::value_type; boost::random::mt19937 rng(0x2545); diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp index d57a83a5e..f628ca09d 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp @@ -57,6 +57,15 @@ namespace nil::crypto3::math::polynomial_arithmetic { // Zero keeps the unscreened path. The limit is capped at the degree-block size so the computed Frobenius // powers can be reused as the start of that table if factorization continues. std::size_t distinct_degree_screening_limit = 0; + + // Optional prepared Fp12 coefficient combinations. This limits the additional owned coefficient payload, + // excluding the ordinary baby steps and per-call scratch. Zero keeps the scalar path; enable a budget + // explicitly until benchmarks establish the crossover and acceptable memory use for the workload. + std::size_t modular_composition_prepared_byte_limit = 0; + + // Prepared BN254 combinations may use four or eight IFMA lanes when the build supports them. + // Zero retains the scalar prepared kernel. The byte budget above remains the opt-in switch. + std::size_t modular_composition_simd_lanes = 8; }; /** diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp index 7da858649..a2ae2080e 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp @@ -26,16 +26,314 @@ #define CRYPTO3_MATH_POLYNOMIAL_COMPOSITION_HPP #include +#include +#include #include +#include #include +#include #include #include +#include #include #include namespace nil::crypto3::math { + namespace detail { + template + class prepared_composition_coefficients { + public: + static constexpr bool supported = false; + }; + + template + class prepared_composition_coefficients< + Value, std::void_t> { + using params = typename Value::field_type::extension_policy::fast_params; + using kernel = algebra::fields::detail::fp12_fast::fp12_fast; + using prepared_type = typename kernel::prepared_product_operand; + using limb = algebra::fields::detail::fp12_fast::limb; + + public: + static constexpr bool supported = true; + static constexpr std::size_t block_tile_size = 4; + static constexpr std::size_t coefficient_tile_size = 4; + + // Scratch belongs to a composition call, so concurrent readers never mutate the shared cache. + struct workspace { + std::array, block_tile_size> scalars; + std::array, block_tile_size> exponents; + std::array, block_tile_size> simd_scalars; + }; + + prepared_composition_coefficients() = default; + + template + prepared_composition_coefficients(const std::vector &powers, std::size_t byte_limit, + std::size_t simd_lanes = 8) { + if (simd_lanes != 0 && simd_lanes != 4 && simd_lanes != 8) { + throw std::invalid_argument("composition SIMD lanes must be 0, 4 or 8"); + } + if (byte_limit == 0 || powers.size() < 2) { + return; + } + const std::size_t terms = powers.size() - 1; // The constant power is handled by direct addition. + std::size_t width = 0; + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + width = std::max(width, powers[exponent].size()); + } +#if defined(CRYPTO3_HAS_AVX512IFMA) + if constexpr (params::base_value_limb_count == 4 && params::base_field_type::modulus_bits <= 254) { + // Cache dimensions rule out short dots and narrow polynomials. use_for_block also checks + // the actual nonzero outer terms before selecting either prepared kernel. + if (simd_lanes != 0 && width >= 4 && terms >= minimum_prepared_terms && + prepare_simd(powers, width, terms, byte_limit, simd_lanes)) { + return; + } + } +#endif + const std::size_t capacity = std::min(byte_limit / sizeof(prepared_type), coefficients_.max_size()); + // Division before multiplication checks both the budget and size_t overflow without allocating. + if (width == 0 || width > capacity / terms) { + return; + } + coefficients_.resize(width * terms); + for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + if (coefficient < powers[exponent].size()) { + coefficients_[coefficient * terms + exponent - 1] = + kernel::prepare_product_operand(powers[exponent][coefficient]); + } + // Default-constructed backends are zero, padding shorter powers without borrowed views. + } + } + terms_ = terms; + width_ = width; + } + + bool available() const { + return !coefficients_.empty() || !simd_coefficients_.empty(); + } + + std::size_t bytes() const { + return coefficients_.size() * sizeof(prepared_type) + simd_coefficients_.size() * sizeof(limb); + } + + std::size_t simd_lanes() const { + return simd_lanes_; + } + + // Short or sparse sums reuse ordinary Fp12 multiplication instead of paying for operand preparation + // and recombination. The constant coefficient needs only an addition and does not count as a term. + template + bool use_for_block(const Polynomial &outer, std::size_t block_index) const { + if (!available()) { + return false; + } + const std::size_t block_size = terms_ + 1; + const std::size_t first = block_index * block_size; + assert(first < outer.size()); + const std::size_t count = std::min(block_size, outer.size() - first); + if (count <= minimum_prepared_terms) { + return false; + } + std::size_t active_terms = 0; + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero() && ++active_terms == minimum_prepared_terms) { + return true; + } + } + return false; + } + + // Evaluate up to four adjacent outer blocks. Each scalar is prepared once and reused across all + // output coefficients. The coefficient-major cache keeps each dot contiguous; tiles reuse the same + // cached columns across adjacent blocks without storing a matrix of every block's output. + template + void evaluate_blocks(std::array &output, const Polynomial &outer, + std::size_t first_block, std::size_t block_count, workspace &scratch) const { + assert(available() && block_count > 0 && block_count <= block_tile_size); +#if defined(CRYPTO3_HAS_AVX512IFMA) + if constexpr (params::base_value_limb_count == 4 && params::base_field_type::modulus_bits <= 254) { + if (simd_lanes_ == 4) { + evaluate_simd_blocks<4>(output, outer, first_block, block_count, scratch); + return; + } + if (simd_lanes_ == 8) { + evaluate_simd_blocks<8>(output, outer, first_block, block_count, scratch); + return; + } + } +#endif + const std::size_t block_size = terms_ + 1; + for (std::size_t block = 0; block < block_count; ++block) { + const std::size_t first = (first_block + block) * block_size; + assert(first < outer.size()); + const std::size_t count = std::min(block_size, outer.size() - first); + auto &scalars = scratch.scalars[block]; + auto &exponents = scratch.exponents[block]; + scalars.clear(); + exponents.clear(); + scalars.reserve(count - 1); + exponents.reserve(count - 1); + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero()) { + scalars.push_back(kernel::prepare_product_operand(outer[first + exponent])); + exponents.push_back(exponent - 1); + } + } + output[block].assign(scalars.empty() ? 1 : width_, Value::zero()); + output[block][0] = outer[first]; + } + for (std::size_t first_coefficient = 0; first_coefficient < width_;) { + const std::size_t end = + first_coefficient + std::min(coefficient_tile_size, width_ - first_coefficient); + for (std::size_t coefficient = first_coefficient; coefficient < end; ++coefficient) { + const auto select_power = [this, coefficient](std::size_t exponent) -> const prepared_type & { + return coefficients_[coefficient * terms_ + exponent]; + }; + for (std::size_t block = 0; block < block_count; ++block) { + const auto &scalars = scratch.scalars[block]; + if (!scalars.empty()) { + const auto powers = scratch.exponents[block] | std::views::transform(select_power); + output[block][coefficient] += kernel::template inner_product( + scalars.begin(), scalars.end(), powers.begin()); + } + } + } + first_coefficient = end; + } + for (std::size_t block = 0; block < block_count; ++block) { + condense(output[block]); + } + } + + private: + // Conservative dispatch heuristic shared with the SIMD cache cutoff, not an arithmetic bound. + static constexpr std::size_t minimum_prepared_terms = 8; +#if defined(CRYPTO3_HAS_AVX512IFMA) + template + bool prepare_simd(const std::vector &powers, std::size_t width, std::size_t terms, + std::size_t byte_limit, std::size_t lanes) { + using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; + constexpr std::size_t coordinates = kernel::product_coordinate_count; + const std::size_t panel_count = 1 + (width - 1) / lanes; + const std::size_t stride = coordinates * 5 * lanes; + const std::size_t capacity = std::min(byte_limit / sizeof(limb), simd_coefficients_.max_size()); + if (terms > capacity / stride || panel_count > capacity / stride / terms) { + return false; + } + simd_coefficients_.resize(panel_count * stride * terms); + // Owned [panel][coordinate][term][digit][lane] storage. Only this representation is retained; + // 40 bytes per form replace the scalar backend object, and padding is included in the budget. + for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { + const std::size_t panel = coefficient / lanes, lane = coefficient % lanes; + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + if (coefficient >= powers[exponent].size()) { + continue; + } + const auto prepared = kernel::prepare_product_operand(powers[exponent][coefficient]); + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const auto packed = radix::from_montgomery64(prepared[coordinate]); + const std::size_t offset = ((panel * coordinates + coordinate) * terms + exponent - 1) * 5; + for (std::size_t digit = 0; digit < 5; ++digit) { + simd_coefficients_[(offset + digit) * lanes + lane] = packed[digit]; + } + } + } + } + terms_ = terms; + width_ = width; + simd_lanes_ = lanes; + return true; + } + + template + void evaluate_simd_blocks(std::array &output, const Polynomial &outer, + std::size_t first_block, std::size_t block_count, workspace &scratch) const { + namespace fast = algebra::fields::detail::fp12_fast; + using radix = fast::montgomery_radix52; + constexpr std::size_t coordinates = kernel::product_coordinate_count; + bool any_active = false; + for (std::size_t block = 0; block < block_count; ++block) { + const std::size_t first = (first_block + block) * (terms_ + 1); + assert(first < outer.size()); + const std::size_t count = std::min(terms_ + 1, outer.size() - first); + auto &indices = scratch.exponents[block]; + indices.clear(); + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero()) { + indices.push_back(exponent - 1); + } + } + auto &scalars = scratch.simd_scalars[block]; + if (indices.size() > scalars.max_size() / (coordinates * 5)) { + throw std::length_error("composition SIMD scratch exceeds maximum size"); + } + scalars.resize(indices.size() * coordinates * 5); + for (std::size_t term = 0; term < indices.size(); ++term) { + const auto prepared = kernel::prepare_product_operand(outer[first + indices[term] + 1]); + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const auto packed = radix::from_montgomery64(prepared[coordinate]); + std::copy(packed.begin(), packed.end(), + scalars.begin() + (coordinate * indices.size() + term) * 5); + } + } + any_active = any_active || !indices.empty(); + output[block].assign(indices.empty() ? 1 : width_, Value::zero()); + output[block][0] = outer[first]; + } + for (std::size_t first = 0; any_active && first < width_;) { + const std::size_t active_lanes = std::min(Lanes, width_ - first); + std::array, Lanes>, + block_tile_size> + products; + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const limb *right = simd_coefficients_.data() + + ((first / Lanes * coordinates + coordinate) * terms_) * 5 * Lanes; + // Reuse this coordinate panel across neighboring blocks before moving to the next one. + for (std::size_t block = 0; block < block_count; ++block) { + const auto &indices = scratch.exponents[block]; + if (indices.empty()) { + continue; + } + std::array reduced; + fast::montgomery_inner_product52>( + reduced.data(), scratch.simd_scalars[block].data() + coordinate * indices.size() * 5, + right, indices.data(), indices.size()); + for (std::size_t lane = 0; lane < active_lanes; ++lane) { + typename radix::value_type packed; + for (std::size_t digit = 0; digit < 5; ++digit) { + packed[digit] = reduced[digit * Lanes + lane]; + } + radix::to_montgomery64(products[block][lane][coordinate].data.backend().base_data(), + packed); + } + } + } + for (std::size_t block = 0; block < block_count; ++block) { + if (!scratch.exponents[block].empty()) { + for (std::size_t lane = 0; lane < active_lanes; ++lane) { + output[block][first + lane] += + kernel::template recombine_products(products[block][lane]); + } + } + } + first += active_lanes; + } + for (std::size_t block = 0; block < block_count; ++block) { + condense(output[block]); + } + } +#endif + std::size_t terms_ = 0, width_ = 0, simd_lanes_ = 0; + std::vector coefficients_; + std::vector simd_coefficients_; + }; + } // namespace detail + /** * Reusable baby-step precomputation for Brent-Kung composition with a fixed inner polynomial and divisor. For a * block size k, the object stores the reduced powers @@ -46,6 +344,11 @@ namespace nil::crypto3::math { * maximum_outer_coefficient_count coefficients. A divisor context used for later compositions must represent the * same polynomial B and contain sufficient inverse precision. * + * Supported Fp12 coefficients may additionally cache their 54 prepared base-field forms. The context's + * modular_composition_prepared_byte_limit bounds this optional payload; zero, insufficient space, or an + * unsupported coefficient type keeps the scalar path. Prepared data is owned by this object and read-only + * during composition. Copies own independent data; scratch belongs to each composition call. + * * @throws std::invalid_argument if maximum_outer_coefficient_count or the configured cached-power limit is zero, * or if reducing inner requires more inverse coefficients than divisor_context stores. * @pre inner is a nonempty canonical coefficient polynomial. @@ -56,6 +359,7 @@ namespace nil::crypto3::math { using backend_type = Backend; using polynomial_type = typename backend_type::polynomial_type; using value_type = typename polynomial_type::value_type; + using coefficient_cache_type = detail::prepared_composition_coefficients; polynomial_composition_precomputation( const polynomial_type &inner, std::size_t maximum_outer_coefficient_count, @@ -110,6 +414,11 @@ namespace nil::crypto3::math { mulmod(giant_step_, baby_steps_[block_size_ / 2], baby_steps_[block_size_ / 2 + 1], divisor_context, arithmetic_context); } + if constexpr (coefficient_cache_type::supported) { + prepared_coefficients_ = coefficient_cache_type( + baby_steps_, arithmetic_context.options().modular_composition_prepared_byte_limit, + arithmetic_context.options().modular_composition_simd_lanes); + } } std::size_t maximum_outer_coefficient_count() const { @@ -129,13 +438,71 @@ namespace nil::crypto3::math { return giant_step_; } + std::size_t prepared_coefficient_bytes() const { + if constexpr (coefficient_cache_type::supported) { + return prepared_coefficients_.bytes(); + } else { + return 0; + } + } + + const coefficient_cache_type &prepared_coefficients() const { + return prepared_coefficients_; + } + private: std::size_t maximum_outer_coefficient_count_; std::size_t block_size_ = 0; std::vector baby_steps_; polynomial_type giant_step_; + [[no_unique_address]] coefficient_cache_type prepared_coefficients_; }; + namespace detail { + // Add one evaluated outer block to an existing accumulator. The cache must have a nonconstant divisor; + // outer must fit it, block_index must exist, and result must not alias outer or the cached powers. + template + void add_composition_block(typename Backend::polynomial_type &result, + const typename Backend::polynomial_type &outer, + const polynomial_composition_precomputation &precomputation, + std::size_t block_index) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename polynomial_type::value_type; + using cache_type = typename polynomial_composition_precomputation::coefficient_cache_type; + if constexpr (cache_type::supported) { + if (precomputation.prepared_coefficients().use_for_block(outer, block_index)) { + typename cache_type::workspace scratch; + std::array blocks; + precomputation.prepared_coefficients().evaluate_blocks(blocks, outer, block_index, 1, scratch); + addition(result, result, blocks[0]); + return; + } + } + + const std::size_t block_size = precomputation.block_size(); + const std::size_t first = block_index * block_size; + const std::size_t count = std::min(block_size, outer.size() - first); + for (std::size_t exponent = 0; exponent < count; ++exponent) { + const value_type &scalar = outer[first + exponent]; + if (scalar == value_type {}) { + continue; + } + if (exponent == 0) { + result[0] = result[0] + scalar; + continue; + } + const polynomial_type &power = precomputation.baby_step(exponent); + if (result.size() < power.size()) { + result.resize(power.size(), value_type {}); + } + for (std::size_t coefficient = 0; coefficient < power.size(); ++coefficient) { + result[coefficient] = result[coefficient] + power[coefficient] * scalar; + } + } + condense(result); + } + } // namespace detail + /** * Compute output = outer(inner(X)) mod B by Horner's rule, where B is the nonzero polynomial stored in * divisor_context. This Horner algorithm performs one modular multiplication for each coefficient of outer except @@ -206,6 +573,11 @@ namespace nil::crypto3::math { * quadratic work can dominate when L and degree(B) are both large. The cached-power limit may select a smaller k, * storing at most k * degree(B) coefficients at the cost of more giant-step multiplications. * + * When the precomputation contains prepared coefficient forms, up to four adjacent blocks with sufficiently + * many nonzero nonconstant terms are evaluated together using fused base-field dot products. Short or sparse + * blocks use scalar multiplication. Both paths preserve the descending Horner order. At most four block + * polynomials are buffered; each selected outer coefficient is prepared once, inside the block kernel. + * * @throws std::invalid_argument if outer exceeds the maximum coefficient count selected when the precomputation * was constructed. * @pre outer is a nonempty canonical coefficient polynomial. divisor_context represents the same divisor used to @@ -232,40 +604,52 @@ namespace nil::crypto3::math { } const std::size_t block_size = precomputation.block_size(); - // Evaluate F_block(inner) as a linear combination of the baby steps and add it directly to the Horner - // accumulator. Streaming one block at a time avoids storing the full matrix of evaluated block polynomials. - // TODO: This coefficientwise phase costs O(outer.size() * degree(B)), which is quadratic when both dimensions - // grow together. A faster matrix-product or backend-provided linear-combination kernel could replace it. - auto add_block = [&](polynomial_type &result, std::size_t block_index) { - const std::size_t first_coefficient = block_index * block_size; - const std::size_t block_coefficient_count = std::min(block_size, outer.size() - first_coefficient); - for (std::size_t exponent = 0; exponent < block_coefficient_count; ++exponent) { - const value_type &scalar = outer[first_coefficient + exponent]; - if (scalar == value_type {}) { - continue; - } - if (exponent == 0) { - result[0] = result[0] + scalar; - continue; - } - - const polynomial_type &power = precomputation.baby_step(exponent); - if (result.size() < power.size()) { - result.resize(power.size(), value_type {}); - } - for (std::size_t coefficient = 0; coefficient < power.size(); ++coefficient) { - result[coefficient] = result[coefficient] + power[coefficient] * scalar; + const std::size_t block_count = 1 + (outer.size() - 1) / block_size; + using cache_type = typename polynomial_composition_precomputation::coefficient_cache_type; + if constexpr (cache_type::supported) { + if (precomputation.prepared_coefficients().available()) { + typename cache_type::workspace scratch; + std::array blocks; + polynomial_type result; + bool leading = true; + for (std::size_t end = block_count; end > 0;) { + std::size_t first = end; + while (first > 0 && end - first < cache_type::block_tile_size && + precomputation.prepared_coefficients().use_for_block(outer, first - 1)) { + --first; + } + if (first == end) { + if (!leading) { + mulmod(result, result, precomputation.giant_step(), divisor_context, arithmetic_context); + } + detail::add_composition_block(result, outer, precomputation, end - 1); + leading = false; + --end; + continue; + } + const std::size_t count = end - first; + precomputation.prepared_coefficients().evaluate_blocks(blocks, outer, first, count, scratch); + // Evaluate blocks together, then consume them in the original descending Horner order. + for (std::size_t block = count; block-- > 0;) { + if (leading) { + result = std::move(blocks[block]); + leading = false; + } else { + mulmod(result, result, precomputation.giant_step(), divisor_context, arithmetic_context); + addition(result, result, blocks[block]); + } + } + end = first; } + output = std::move(result); + return; } - condense(result); - }; - - const std::size_t block_count = 1 + (outer.size() - 1) / block_size; + } polynomial_type result; - add_block(result, block_count - 1); + detail::add_composition_block(result, outer, precomputation, block_count - 1); for (std::size_t block = block_count - 1; block-- > 0;) { mulmod(result, result, precomputation.giant_step(), divisor_context, arithmetic_context); - add_block(result, block); + detail::add_composition_block(result, outer, precomputation, block); } output = std::move(result); diff --git a/libs/math/test/polynomial_composition.cpp b/libs/math/test/polynomial_composition.cpp index 35d91af72..468106878 100644 --- a/libs/math/test/polynomial_composition.cpp +++ b/libs/math/test/polynomial_composition.cpp @@ -25,6 +25,8 @@ #define BOOST_TEST_MODULE polynomial_composition_test #include +#include +#include #include @@ -35,6 +37,7 @@ #include #include #include +#include namespace { namespace math = nil::crypto3::math; @@ -318,4 +321,423 @@ BOOST_AUTO_TEST_CASE(composition_matches_between_schoolbook_and_mixed_radix_back BOOST_CHECK(result == expected); } +BOOST_AUTO_TEST_CASE(prepared_composition_cache_obeys_its_budget_and_pads_short_powers) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + using precomputation_type = math::polynomial_composition_precomputation; + using kernel = fields::detail::fp12_fast::fp12_fast; + // k=5, inner=X: the nonconstant baby steps have lengths 2,3,4,5. Padding is included in the budget. + const std::size_t required = 4 * 5 * sizeof(kernel::prepared_product_operand); + const polynomial_type inner = {fq12_value_type::zero(), fq12_value_type::one()}; + polynomial_type divisor(8, fq12_value_type::zero()); + divisor[0] = fq12_value(3); + divisor.back() = fq12_value(17); // Nonmonic divisor. + polynomial_type outer(21); + for (std::size_t i = 0; i < outer.size(); ++i) { + outer[i] = fq12_value(23 + 12 * i); + } + for (std::size_t limit : + {std::size_t(0), required - 1, required, required + 1, std::numeric_limits::max()}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 5; + options.modular_composition_prepared_byte_limit = limit; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + math::polynomial_divisor_context division(divisor, 6, context); + const precomputation_type cache(inner, outer.size(), division, context); + BOOST_CHECK_EQUAL(cache.block_size(), 5); + BOOST_CHECK_EQUAL(cache.prepared_coefficient_bytes(), limit >= required ? required : 0); + polynomial_type expected, actual; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(actual, outer, cache, division, context); + BOOST_CHECK(actual == expected); + + // Check the shared single-block entry point too, including an existing accumulator and cancellation. + for (std::size_t block = 0; block < 5; ++block) { + const std::size_t first = block * 5; + const std::size_t end = std::min(first + 5, outer.size()); + polynomial_type slice(outer.begin() + first, outer.begin() + end), evaluated; + math::compose_mod_reference(evaluated, slice, inner, division, context); + actual = {fq12_value(77), fq12_value(91)}; + math::addition(expected, actual, evaluated); + math::detail::add_composition_block(actual, outer, cache, block); + BOOST_CHECK(actual == expected); + math::subtraction(actual, polynomial_type {}, evaluated); + math::detail::add_composition_block(actual, outer, cache, block); + BOOST_CHECK(actual == polynomial_type {fq12_value_type::zero()}); + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_composition_cache_owns_its_data_and_survives_copy_and_move) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + using precomputation_type = math::polynomial_composition_precomputation; + const polynomial_type inner = {fq12_value(1), fq12_value(13), fq12_value(25)}; + const polynomial_type divisor = {fq12_value(37), fq12_value(49), fq12_value(61), fq12_value(73), fq12_value(85)}; + std::optional saved; + { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_prepared_byte_limit = 1024 * 1024; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + math::polynomial_divisor_context division(divisor, 3, context); + polynomial_type temporary_inner = inner; + saved.emplace(temporary_inner, 25, division, context); + temporary_inner.assign(1, fq12_value(999)); + } + BOOST_REQUIRE_GT(saved->prepared_coefficient_bytes(), 0); + auto copy = *saved; + const std::size_t bytes = copy.prepared_coefficient_bytes(); + saved.reset(); + auto moved = std::move(copy); + BOOST_CHECK_EQUAL(moved.prepared_coefficient_bytes(), bytes); + polynomial_arithmetic::polynomial_context context; + math::polynomial_divisor_context division(divisor, 3, context); + for (std::size_t count : {1, 4, 5, 6, 19, 20, 21, 25}) { + polynomial_type outer(count); + for (std::size_t i = 0; i < count; ++i) { + outer[i] = i % 3 == 0 ? fq12_value_type::zero() : fq12_value(i * 12 + count); + } + outer.back() = fq12_value(123); + polynomial_type expected, actual; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(actual, outer, moved, division, context); + BOOST_CHECK(actual == expected); + math::compose_mod(outer, outer, moved, division, context); + BOOST_CHECK(outer == expected); + } + polynomial_type oversized(26, fq12_value(1)), unchanged {fq12_value(91)}; + BOOST_CHECK_THROW(math::compose_mod(unchanged, oversized, moved, division, context), std::invalid_argument); + BOOST_CHECK(unchanged == polynomial_type {fq12_value(91)}); +} + +BOOST_AUTO_TEST_CASE(prepared_composition_dispatch_counts_nonzero_nonconstant_terms) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + constexpr std::size_t block_size = 17; + polynomial_type divisor(10), inner(9); + for (std::size_t i = 0; i < divisor.size(); ++i) + divisor[i] = fq12_value(12 * i + 1); + for (std::size_t i = 0; i < inner.size(); ++i) + inner[i] = fq12_value(12 * i + 121); + + for (std::size_t lanes : {0, 4, 8}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = block_size; + options.modular_composition_prepared_byte_limit = 4 * 1024 * 1024; + options.modular_composition_simd_lanes = lanes; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + math::polynomial_divisor_context division(divisor, 8, context); + const math::polynomial_composition_precomputation cache(inner, block_size * block_size, division, + context); + BOOST_REQUIRE_EQUAL(cache.block_size(), block_size); + BOOST_REQUIRE(cache.prepared_coefficients().available()); +#if defined(CRYPTO3_HAS_AVX512IFMA) + BOOST_CHECK_EQUAL(cache.prepared_coefficients().simd_lanes(), lanes); +#else + BOOST_CHECK_EQUAL(cache.prepared_coefficients().simd_lanes(), 0); +#endif + for (std::size_t active_terms : {0, 1, 7, 8, 9, 16}) { + polynomial_type outer(2 * block_size, fq12_value_type::zero()); + outer[0] = fq12_value(13); // A nonzero constant must not move a seven-term sum above the cutoff. + outer.back() = fq12_value(25); + // Place active terms at the end of the block, including a lone high-exponent term. + for (std::size_t j = 1; j <= active_terms; ++j) + outer[block_size - j] = fq12_value(12 * j + 37); + BOOST_CHECK_EQUAL(cache.prepared_coefficients().use_for_block(outer, 0), active_terms >= 8); + BOOST_CHECK(!cache.prepared_coefficients().use_for_block(outer, 1)); + + polynomial_type slice(outer.begin(), outer.begin() + block_size), evaluated; + math::condense(slice); + math::compose_mod_reference(evaluated, slice, inner, division, context); + polynomial_type actual(11, fq12_value(61)), expected; + math::addition(expected, actual, evaluated); + math::detail::add_composition_block(actual, outer, cache, 0); + BOOST_CHECK(actual == expected); + math::subtraction(actual, polynomial_type {}, evaluated); + math::detail::add_composition_block(actual, outer, cache, 0); + BOOST_CHECK(actual == polynomial_type {fq12_value_type::zero()}); + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_composition_mixes_scalar_blocks_with_dense_batches_and_partial_tails) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + constexpr std::size_t block_size = 17; + // Five adjacent dense blocks cross a four-block batch boundary. Subsequent blocks alternate between paths, + // including an all-zero block, a constant block, and seven-term sums immediately below the cutoff. + const std::array active_terms = {8, 16, 8, 9, 16, 7, 8, 0, 1, 8, 16, 7, 0, 8}; + polynomial_type divisor(10), inner(9); + for (std::size_t i = 0; i < divisor.size(); ++i) + divisor[i] = fq12_value(12 * i + 1); + for (std::size_t i = 0; i < inner.size(); ++i) + inner[i] = fq12_value(12 * i + 121); + + for (std::size_t lanes : {0, 4, 8}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = block_size; + options.modular_composition_prepared_byte_limit = 4 * 1024 * 1024; + options.modular_composition_simd_lanes = lanes; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + math::polynomial_divisor_context division(divisor, 8, context); + const math::polynomial_composition_precomputation cache(inner, block_size * block_size, division, + context); + BOOST_REQUIRE_EQUAL(cache.block_size(), block_size); + BOOST_REQUIRE(cache.prepared_coefficients().available()); + for (std::size_t tail_size : {1, 2, 8, 9, 17}) { + polynomial_type outer(active_terms.size() * block_size + tail_size, fq12_value_type::zero()); + for (std::size_t block = 0; block < active_terms.size(); ++block) { + const std::size_t first = block * block_size; + if (block != 7) + outer[first] = fq12_value(12 * first + 217); + for (std::size_t j = 1; j <= active_terms[block]; ++j) + outer[first + block_size - j] = fq12_value(12 * (first + j) + 229); + BOOST_CHECK_EQUAL(cache.prepared_coefficients().use_for_block(outer, block), active_terms[block] >= 8); + } + for (std::size_t i = active_terms.size() * block_size; i < outer.size(); ++i) + outer[i] = fq12_value(12 * i + 241); + BOOST_CHECK_EQUAL(cache.prepared_coefficients().use_for_block(outer, active_terms.size()), tail_size >= 9); + polynomial_type expected, actual; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(actual, outer, cache, division, context); + BOOST_CHECK(actual == expected); + math::compose_mod(outer, outer, cache, division, context); + BOOST_CHECK(outer == expected); + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_composition_tiles_match_horner_at_block_and_coefficient_boundaries) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + // Dense full blocks contain eight nonconstant terms, so they still exercise the prepared kernels. + constexpr std::size_t block_size = 9; + for (std::size_t degree : {1, 3, 4, 5, 9}) { + polynomial_type divisor(degree + 1), inner(degree); + for (std::size_t i = 0; i <= degree; ++i) { + divisor[i] = fq12_value(i * 12 + 1); + if (i < degree) { + inner[i] = fq12_value(i * 12 + 73); + } + } + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = block_size; + options.modular_composition_prepared_byte_limit = 1024 * 1024; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + math::polynomial_divisor_context division(divisor, std::max(std::size_t(1), degree - 1), context); + const math::polynomial_composition_precomputation cache(inner, 82, division, context); + BOOST_REQUIRE_EQUAL(cache.block_size(), block_size); + BOOST_REQUIRE_GT(cache.prepared_coefficient_bytes(), 0); + for (std::size_t count : {1, 8, 9, 10, 35, 36, 37, 44, 45, 46, 71, 72, 73, 81, 82}) { + for (bool sparse : {false, true}) { + polynomial_type outer(count, fq12_value_type::zero()); + for (std::size_t i = 0; i < count; ++i) { + if (!sparse || i % block_size == 0 || i + 1 == count) { + outer[i] = fq12_value(1 + 12 * i); + } + } + polynomial_type expected, actual; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(actual, outer, cache, division, context); + BOOST_CHECK(actual == expected); + } + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_composition_supports_123_terms_and_the_small_prime_backend) { + using backend_type = polynomial_arithmetic::bn254_small_prime_fft_backend; + using polynomial_type = typename backend_type::polynomial_type; + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 123; + options.modular_composition_prepared_byte_limit = 4 * 1024 * 1024; + polynomial_arithmetic::polynomial_context context(backend_type(32), options); + const polynomial_type divisor = {fq12_value(1), fq12_value(13), fq12_value(25), + fq12_value(37), fq12_value(49), fq12_value(61)}; + const polynomial_type inner = {fq12_value(73), fq12_value(85), fq12_value(97), fq12_value(109), fq12_value(121)}; + math::polynomial_divisor_context division(divisor, 4, context); + const math::polynomial_composition_precomputation cache(inner, 122 * 122 + 1, division, context); + BOOST_REQUIRE_EQUAL(cache.block_size(), 123); + BOOST_REQUIRE_GT(cache.prepared_coefficient_bytes(), 0); + using reference_backend = polynomial_arithmetic::schoolbook_backend; + polynomial_arithmetic::polynomial_context reference_context; + math::polynomial_divisor_context reference_division(divisor, 4, reference_context); + for (std::size_t count : {122, 123, 124, 491, 492, 493, 615}) { + polynomial_type outer(count); + for (std::size_t i = 0; i < count; ++i) { + outer[i] = fq12_value(i * 12 + 1); + } + polynomial_type expected, actual; + math::compose_mod_reference(expected, outer, inner, reference_division, reference_context); + math::compose_mod(actual, outer, cache, division, context); + BOOST_CHECK(actual == expected); + } +} + +BOOST_AUTO_TEST_CASE(prepared_composition_fallbacks_handle_constant_moduli_and_unsupported_fields) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_prepared_byte_limit = 1024 * 1024; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + const polynomial_type inner = {fq12_value(1), fq12_value(13)}; + const polynomial_type outer = {fq12_value(25), fq12_value(37), fq12_value(49)}; + math::polynomial_divisor_context constant_division({fq12_value(7)}, 1, context); + const math::polynomial_composition_precomputation constant_cache(inner, 3, constant_division, + context); + BOOST_CHECK_EQUAL(constant_cache.prepared_coefficient_bytes(), 0); + polynomial_type result = outer; + math::compose_mod(result, result, constant_cache, constant_division, context); + BOOST_CHECK(result == polynomial_type {fq12_value_type::zero()}); + + options.modular_composition_cached_power_limit = 1; + polynomial_arithmetic::polynomial_context scalar_context(backend_type {}, options); + math::polynomial_divisor_context division({fq12_value(1), fq12_value(13), fq12_value(25)}, 1, + scalar_context); + const math::polynomial_composition_precomputation scalar_cache(inner, 3, division, scalar_context); + BOOST_CHECK_EQUAL(scalar_cache.prepared_coefficient_bytes(), 0); + polynomial_type expected; + math::compose_mod_reference(expected, outer, inner, division, scalar_context); + math::compose_mod(result, outer, scalar_cache, division, scalar_context); + BOOST_CHECK(result == expected); + + using base_backend = polynomial_arithmetic::schoolbook_backend; + options.modular_composition_cached_power_limit = 3; + polynomial_arithmetic::polynomial_context base_context(base_backend {}, options); + math::polynomial_divisor_context base_division({fq_value_type(1), fq_value_type(1)}, 1, base_context); + const math::polynomial_composition_precomputation base_cache({fq_value_type(2)}, 3, base_division, + base_context); + BOOST_CHECK_EQUAL(base_cache.prepared_coefficient_bytes(), 0); + base_backend::polynomial_type base_result; + math::compose_mod(base_result, {fq_value_type(1), fq_value_type(2), fq_value_type(3)}, base_cache, base_division, + base_context); + BOOST_CHECK(base_result == base_backend::polynomial_type {fq_value_type(17)}); +} + +BOOST_AUTO_TEST_CASE(prepared_composition_handles_degenerate_inputs_and_one_off_aliases) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_prepared_byte_limit = 1024 * 1024; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + const polynomial_type divisor = {fq12_value(1), fq12_value(13), fq12_value(25), fq12_value(37), fq12_value(49)}; + math::polynomial_divisor_context division(divisor, 3, context); + for (const polynomial_type inner : {polynomial_type {fq12_value_type::zero()}, polynomial_type {fq12_value(61)}, + polynomial_type {fq12_value(73), fq12_value(85), fq12_value(97)}}) { + const math::polynomial_composition_precomputation cache(inner, 25, division, context); + BOOST_REQUIRE_GT(cache.prepared_coefficient_bytes(), 0); + for (const polynomial_type outer : + {polynomial_type {fq12_value_type::zero()}, polynomial_type {fq12_value(109)}, + polynomial_type {-inner[0], fq12_value_type::one()}, + polynomial_type {fq12_value(121), fq12_value_type::zero(), fq12_value(133)}}) { + polynomial_type expected, result; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(result, outer, cache, division, context); + BOOST_CHECK(result == expected); + auto outer_alias = outer; + math::compose_mod(outer_alias, outer_alias, inner, division, context); + BOOST_CHECK(outer_alias == expected); + auto inner_alias = inner; + math::compose_mod(inner_alias, outer, inner_alias, division, context); + BOOST_CHECK(inner_alias == expected); + } + } +} + +BOOST_AUTO_TEST_CASE(prepared_simd_panels_match_scalar_at_lane_and_term_boundaries) { + using polynomial_type = math::polynomial; + using cache_type = math::detail::prepared_composition_coefficients; + constexpr std::size_t budget = 16 * 1024 * 1024; + for (std::size_t width : {3, 4, 5, 7, 8, 9, 15, 17}) { + std::vector powers(18); + powers[0] = {fq12_value_type::one()}; + for (std::size_t j = 1; j < powers.size(); ++j) { + powers[j].resize(j % 3 == 0 ? 1 : width); + for (std::size_t c = 0; c < powers[j].size(); ++c) + powers[j][c] = fq12_value(12 * (j * width + c) + 1); + } + cache_type scalar(powers, budget, 0); + for (std::size_t lanes : {4, 8}) { + cache_type candidate(powers, budget, lanes); + BOOST_REQUIRE(candidate.available()); +#if defined(CRYPTO3_HAS_AVX512IFMA) + BOOST_CHECK_EQUAL(candidate.simd_lanes(), width >= 4 ? lanes : 0); + if (width >= 4) { + const std::size_t expected = + ((width + lanes - 1) / lanes) * lanes * 54 * 17 * 5 * sizeof(std::uint64_t); + BOOST_CHECK_EQUAL(candidate.bytes(), expected); + } +#else + BOOST_CHECK_EQUAL(candidate.simd_lanes(), 0); +#endif + BOOST_CHECK_LE(candidate.bytes(), budget); + const auto bytes = candidate.bytes(); + cache_type exact(powers, bytes, lanes), small(powers, bytes - 1, lanes); + BOOST_CHECK(exact.available()); + BOOST_CHECK_LE(small.bytes(), bytes - 1); + auto copied = candidate; + auto moved = std::move(copied); + typename cache_type::workspace scratch, reference_scratch; + std::array actual, expected; + for (bool sparse : {false, true}) { + polynomial_type outer(5 * powers.size() - 7); + for (std::size_t i = 0; i < outer.size(); ++i) { + outer[i] = sparse && i % 7 ? fq12_value_type::zero() : fq12_value(12 * i + 1); + } + // Include a wholly constant block amid active blocks, then a short final block. + for (std::size_t j = 1; j < powers.size(); ++j) + outer[powers.size() + j] = fq12_value_type::zero(); + for (std::size_t first : {0, 1, 4}) { + const std::size_t count = std::min(std::size_t(4), std::size_t(5) - first); + moved.evaluate_blocks(actual, outer, first, count, scratch); + scalar.evaluate_blocks(expected, outer, first, count, reference_scratch); + for (std::size_t i = 0; i < count; ++i) + BOOST_CHECK(actual[i] == expected[i]); + } + } + } + BOOST_CHECK_THROW(cache_type(powers, budget, 3), std::invalid_argument); + } +} + +BOOST_AUTO_TEST_CASE(prepared_simd_complete_composition_matches_horner_and_aliases) { + using backend = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = backend::polynomial_type; + polynomial_type divisor(10), inner(8), outer(83); + for (std::size_t i = 0; i < divisor.size(); ++i) + divisor[i] = fq12_value(12 * i + 1); + for (std::size_t i = 0; i < inner.size(); ++i) + inner[i] = fq12_value(12 * i + 121); + for (std::size_t i = 0; i < outer.size(); ++i) + outer[i] = i % 3 ? fq12_value(12 * i + 217) : fq12_value_type::zero(); + for (std::size_t lanes : {0, 4, 8}) { + polynomial_arithmetic::polynomial_context_options options; + options.modular_composition_cached_power_limit = 17; + options.modular_composition_prepared_byte_limit = 4 * 1024 * 1024; + options.modular_composition_simd_lanes = lanes; + polynomial_arithmetic::polynomial_context context(backend {}, options); + math::polynomial_divisor_context division(divisor, 8, context); + const math::polynomial_composition_precomputation cache(inner, 257, division, context); + polynomial_type expected, actual = outer; + math::compose_mod_reference(expected, outer, inner, division, context); + math::compose_mod(actual, actual, cache, division, context); + BOOST_CHECK(actual == expected); + } +} + +BOOST_AUTO_TEST_CASE(prepared_simd_configuration_keeps_six_limb_fields_on_scalar_kernel) { + using value = fields::fp12_2over3over2>::value_type; + using polynomial = math::polynomial; + using cache_type = math::detail::prepared_composition_coefficients; + std::vector powers(9, polynomial(4, value::one())); + powers[0] = {value::one()}; + const cache_type cache(powers, 1 << 20, 8); + BOOST_REQUIRE(cache.available()); + BOOST_CHECK_EQUAL(cache.simd_lanes(), 0); + typename cache_type::workspace scratch; + std::array output; + cache.evaluate_blocks(output, polynomial(9, value::one()), 0, 1, scratch); + BOOST_CHECK(output[0] == polynomial({value(9), value(8), value(8), value(8)})); +} + BOOST_AUTO_TEST_SUITE_END()