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 6a92fe87f..f2579d80a 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
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <algorithm>
 #include <array>
 #include <cstddef>
 #include <iterator>
@@ -27,6 +28,131 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast {
         static constexpr std::size_t product_coordinate_count = 54;
         using prepared_product_operand = std::array<base_backend_type, product_coordinate_count>;
 
+        // BN254 has v=w^2 and u=w^6-9, hence w^12-18*w^6+82=0. Evaluate
+        // degree-11 operands at 0,+1,-1,...,+11,-11: their degree-22 product is
+        // determined by 23 pointwise products, also after summing any number of operands.
+        struct rank23 {
+            static constexpr bool supported = Params::xi == std::array {9, 1} && Params::u_squared == -1 &&
+                                              Params::base_value_limb_count == 4 &&
+                                              Params::base_field_type::modulus_bits <= 254;
+            static constexpr std::size_t product_coordinate_count = 23;
+            using prepared_product_operand = std::array<base_backend_type, product_coordinate_count>;
+
+            template<typename Fp12Value>
+            static prepared_product_operand prepare_product_operand(const Fp12Value &x) {
+                static_assert(supported);
+                std::array<base_value_type, 12> flat;
+                for (std::size_t w = 0; w < 2; ++w) {
+                    for (std::size_t v = 0; v < 3; ++v) {
+                        const auto &a = x.data[w].data[v].data;
+                        flat[2 * v + w] = a[0] - (a[1].doubled().doubled().doubled() + a[1]);
+                        flat[6 + 2 * v + w] = a[1];
+                    }
+                }
+                prepared_product_operand result;
+                result[0] = flat[0].data.backend().base_data();
+                // Reuse even/odd evaluations for each +/- pair. All arithmetic is canonical
+                // Montgomery arithmetic; no integer coefficient growth enters the dot kernel.
+                for (std::size_t point = 1; point <= 11; ++point) {
+                    const base_value_type t(point), t_squared(point * point);
+                    auto even = flat[10], odd = flat[11];
+                    for (std::size_t j = 5; j > 0; --j) {
+                        if (point != 1) {
+                            even *= t_squared;
+                            odd *= t_squared;
+                        }
+                        even += flat[2 * (j - 1)];
+                        odd += flat[2 * (j - 1) + 1];
+                    }
+                    if (point != 1) {
+                        odd *= t;
+                    }
+                    result[2 * point - 1] = (even + odd).data.backend().base_data();
+                    result[2 * point] = (even - odd).data.backend().base_data();
+                }
+                return result;
+            }
+
+            // Input/output: [coordinate][digit][lane], canonical residues scaled by S=2^260.
+            // Fold interpolation, reduction modulo w^12-18*w^6+82, and the tower conversion
+            // into a fixed 12-by-23 map. The existing bounded dot kernel applies each row;
+            // only the caller's twelve final coordinates need conversion back to R=2^256.
+            template<typename Operations>
+            static void recombine_products52(limb *result, const limb *products) {
+                static_assert(supported);
+                using radix = montgomery_radix52<Params>;
+                constexpr std::size_t count = product_coordinate_count;
+                static const auto weights = [] {
+                    std::array<limb, 12 * count * radix::digits> packed;
+                    std::array<base_value_type, count> points;
+                    points[0] = base_value_type::zero();
+                    for (std::size_t i = 1; i < count; ++i) {
+                        points[i] = base_value_type((i + 1) / 2);
+                        if (i % 2 == 0) {
+                            points[i] = -points[i];
+                        }
+                    }
+                    // Fixed-size Lagrange basis construction keeps algebra independent of the
+                    // higher-level math polynomial module. Initialization happens once per kernel.
+                    for (std::size_t i = 0; i < count; ++i) {
+                        std::array<base_value_type, count> basis {};
+                        basis[0] = base_value_type::one();
+                        auto denominator = base_value_type::one();
+                        std::size_t degree = 0;
+                        for (std::size_t j = 0; j < count; ++j) {
+                            if (i == j) {
+                                continue;
+                            }
+                            basis[degree + 1] = basis[degree];
+                            for (std::size_t k = degree; k > 0; --k) {
+                                basis[k] = basis[k - 1] - points[j] * basis[k];
+                            }
+                            basis[0] *= -points[j];
+                            denominator *= points[i] - points[j];
+                            ++degree;
+                        }
+                        const auto inverse = denominator.inversed();
+                        for (auto &coefficient : basis) {
+                            coefficient *= inverse;
+                        }
+                        // Descending reduction also handles degrees 18..22, whose first
+                        // reduction still contributes to degrees at least twelve.
+                        for (std::size_t k = count - 1; k >= 12; --k) {
+                            basis[k - 6] += base_value_type(18) * basis[k];
+                            basis[k - 12] -= base_value_type(82) * basis[k];
+                        }
+                        for (std::size_t w = 0; w < 2; ++w) {
+                            for (std::size_t v = 0; v < 3; ++v) {
+                                const auto &imag = basis[6 + 2 * v + w];
+                                const std::array<base_value_type, 2> tower = {
+                                    basis[2 * v + w] + base_value_type(9) * imag, imag};
+                                for (std::size_t u = 0; u < 2; ++u) {
+                                    const auto digits = radix::from_montgomery64(tower[u].data.backend().base_data());
+                                    const std::size_t offset = ((6 * w + 2 * v + u) * count + i) * radix::digits;
+                                    std::copy(digits.begin(), digits.end(), packed.begin() + offset);
+                                }
+                            }
+                        }
+                    }
+                    return packed;
+                }();
+                static constexpr auto indices = [] {
+                    std::array<std::size_t, count> result {};
+                    for (std::size_t i = 0; i < count; ++i) {
+                        result[i] = i;
+                    }
+                    return result;
+                }();
+                // Each interpolation dot has only 23 canonical terms, within the existing
+                // 256-term IFMA bound. Its Montgomery reduction preserves the S scaling.
+                for (std::size_t c = 0; c < 12; ++c) {
+                    montgomery_inner_product52<Params, Operations>(result + c * radix::digits * Operations::lanes,
+                                                                   weights.data() + c * count * radix::digits, products,
+                                                                   indices.data(), count);
+                }
+            }
+        };
+
         // Packed input-side Fp2 value: two canonical Fp coefficients.
         struct fp2_base {
             limb_array data;
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 82ab40e01..6f18bdf8d 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
@@ -567,6 +567,31 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast {
         }
     };
 
+    // Lane-wise conditional modulus subtraction for the radix-2^52 dot-product reduction.
+    template<Fp12FastParams Params, typename Operations>
+    struct montgomery_radix52_simd {
+        using radix = montgomery_radix52<Params>;
+        using value_type = std::array<typename Operations::value_type, radix::digits>;
+
+        // x may be noncanonical, but must have normalized digits and x<S.
+        static void subtract_modulus(value_type &x) {
+            const auto zero = Operations::broadcast(0);
+            const auto one = Operations::broadcast(1);
+            const auto mask = Operations::broadcast(radix::mask);
+            value_type difference;
+            auto borrow = zero;
+            for (std::size_t i = 0; i < radix::digits; ++i) {
+                const auto subtrahend = Operations::add(Operations::broadcast(radix::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 < radix::digits; ++i) {
+                x[i] = Operations::select(keep, x[i], difference[i]);
+            }
+        }
+    };
+
     // 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.
@@ -578,26 +603,12 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast {
         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<vector, 5> 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<vector, 5> &x) {
-            std::array<vector, 5> 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<vector, 10> t;
@@ -639,7 +650,7 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast {
             }
             reduced[4] = t[9];
             for (unsigned i = 0; i < 4; ++i) {
-                subtract(reduced);
+                montgomery_radix52_simd<Params, Operations>::subtract_modulus(reduced);
             }
             auto carry = zero;
             for (std::size_t i = 0; i < 5; ++i) {
@@ -647,7 +658,7 @@ namespace nil::crypto3::algebra::fields::detail::fp12_fast {
                 total[i] = Operations::bit_and(sum, mask);
                 carry = Operations::shift52(sum);
             }
-            subtract(total);
+            montgomery_radix52_simd<Params, Operations>::subtract_modulus(total);
             first += terms;
         }
         for (std::size_t i = 0; i < 5; ++i) {
diff --git a/libs/algebra/test/fp12.cpp b/libs/algebra/test/fp12.cpp
index 342ef634b..34eb795a3 100644
--- a/libs/algebra/test/fp12.cpp
+++ b/libs/algebra/test/fp12.cpp
@@ -292,6 +292,200 @@ namespace {
         }
     }
 
+    template<typename Operations>
+    void check_rank23_recombination() {
+        namespace fast = fields::detail::fp12_fast;
+        using field = alt_bn128_254_fp12;
+        using value = field::value_type;
+        using params = field::extension_policy::fast_params;
+        using kernel = typename fast::fp12_fast<params>::rank23;
+        using base = typename fast::fp12_fast<params>::base_value_type;
+        using radix = fast::montgomery_radix52<params>;
+        constexpr std::size_t lanes = Operations::lanes, count = kernel::product_coordinate_count;
+        static_assert(count == 23);
+        // Independent interpolation oracle: evaluate each Lagrange polynomial directly at
+        // the tower element w using generic Fp12 multiplication. No flat-basis reduction or
+        // production interpolation weights are used to construct the expected result.
+        std::array<base, count> points;
+        points[0] = base::zero();
+        for (std::size_t i = 1; i < count; ++i) {
+            points[i] = base((i + 1) / 2);
+            if (i % 2 == 0) {
+                points[i] = -points[i];
+            }
+        }
+        std::array<value, count> weights;
+        for (std::size_t i = 0; i < count; ++i) {
+            weights[i] = value::one();
+            auto denominator = base::one();
+            for (std::size_t j = 0; j < count; ++j) {
+                if (i != j) {
+                    auto factor = value::zero();
+                    factor.coordinate(6) = base::one();    // w
+                    factor.coordinate(0) = -points[j];
+                    weights[i] = generic_fp12_mul<field>(weights[i], factor);
+                    denominator *= points[i] - points[j];
+                }
+            }
+            auto inverse = value::zero();
+            inverse.coordinate(0) = denominator.inversed();
+            weights[i] = generic_fp12_mul<field>(weights[i], inverse);
+        }
+        auto maximum = radix::modulus();
+        --maximum[0];
+        boost::random::mt19937 rng(0x2312);
+        // Every interpolation basis vector, maximum raw residues, all 144 tower basis
+        // products, boundary pairs, and random products; each lane carries independent data.
+        for (std::size_t pattern = 0; pattern < 26 + 144 + 36 + random_samples; ++pattern) {
+            std::array<fast::limb, count * 5 * lanes> products {};
+            std::array<value, lanes> expected;
+            expected.fill(value::zero());
+            for (std::size_t lane = 0; lane < lanes; ++lane) {
+                std::array<base, count> scalar;
+                if (pattern < 26) {
+                    for (std::size_t i = 0; i < count; ++i) {
+                        if (pattern == 24 || (pattern == 25 && (i + lane) % 2)) {
+                            radix::to_montgomery64(scalar[i].data.backend().base_data(), maximum);
+                        } else {
+                            scalar[i] = i == pattern ? base(lane + 1) : base::zero();
+                        }
+                        auto factor = value::zero();
+                        factor.coordinate(0) = scalar[i];
+                        expected[lane] += generic_fp12_mul<field>(weights[i], factor);
+                    }
+                } else {
+                    value left, right;
+                    if (pattern < 26 + 144) {
+                        left = value::zero();
+                        right = value::zero();
+                        left.coordinate((pattern - 26) / 12) = base(lane + 1);
+                        right.coordinate((pattern - 26) % 12) = base(lane + 2);
+                    } else if (pattern < 26 + 144 + 36) {
+                        left = boundary_fp12<field>((pattern - 26 - 144) / 6);
+                        right = boundary_fp12<field>((pattern - 26 - 144) % 6);
+                    } else {
+                        left = random_fp12<field>(rng);
+                        right = random_fp12<field>(rng);
+                    }
+                    const auto a = kernel::prepare_product_operand(left);
+                    const auto b = kernel::prepare_product_operand(right);
+                    for (std::size_t i = 0; i < count; ++i) {
+                        const auto &modulus = base::modulus_params.get_mod_obj().get_mod();
+                        BOOST_CHECK(eval_lt(a[i], modulus));
+                        BOOST_CHECK(eval_lt(b[i], modulus));
+                        base x, y;
+                        x.data.backend().base_data() = a[i];
+                        y.data.backend().base_data() = b[i];
+                        scalar[i] = x * y;
+                    }
+                    expected[lane] = generic_fp12_mul<field>(left, right);
+                }
+                for (std::size_t i = 0; i < count; ++i) {
+                    const auto digits = radix::from_montgomery64(scalar[i].data.backend().base_data());
+                    for (std::size_t digit = 0; digit < 5; ++digit) {
+                        products[(i * 5 + digit) * lanes + lane] = digits[digit];
+                    }
+                }
+            }
+            std::array<fast::limb, 12 * 5 * lanes + 2> combined;
+            combined.fill(0xAABBCCDDEEFF0011);
+            kernel::template recombine_products52<Operations>(combined.data() + 1, products.data());
+            BOOST_CHECK_EQUAL(combined.front(), 0xAABBCCDDEEFF0011);
+            BOOST_CHECK_EQUAL(combined.back(), 0xAABBCCDDEEFF0011);
+            for (std::size_t lane = 0; lane < lanes; ++lane) {
+                value actual;
+                for (std::size_t c = 0; c < 12; ++c) {
+                    typename radix::value_type packed;
+                    for (std::size_t digit = 0; digit < 5; ++digit) {
+                        packed[digit] = combined[1 + (c * 5 + digit) * lanes + lane];
+                        BOOST_CHECK_LE(packed[digit], radix::mask);
+                    }
+                    BOOST_CHECK(std::lexicographical_compare(packed.rbegin(), packed.rend(), radix::modulus().rbegin(),
+                                                             radix::modulus().rend()));
+                    radix::to_montgomery64(actual.coordinate(c).data.backend().base_data(), packed);
+                }
+                BOOST_CHECK_EQUAL(actual, expected[lane]);
+            }
+        }
+    }
+
+    template<typename Operations>
+    void check_rank23_fp12_dots() {
+        namespace fast = fields::detail::fp12_fast;
+        using value = alt_bn128_254_fp12::value_type;
+        using params = alt_bn128_254_fp12::extension_policy::fast_params;
+        using kernel = typename fast::fp12_fast<params>::rank23;
+        using radix = fast::montgomery_radix52<params>;
+        constexpr std::size_t lanes = Operations::lanes, terms = 514, coordinates = kernel::product_coordinate_count;
+        constexpr std::size_t panel_stride = 2 * terms * 5 * lanes;
+        boost::random::mt19937 rng(0x5212D07);
+        std::vector<value> left(terms);
+        std::vector<std::array<value, lanes>> right(terms);
+        std::vector<fast::limb> packed_left(coordinates * terms * 5), packed_right(coordinates * panel_stride, 0);
+        std::array<std::size_t, terms> indices;
+        for (std::size_t term = 0; term < terms; ++term) {
+            // The second half cancels the first across the 256-term reduction boundaries.
+            left[term] = term >= 257   ? -left[term - 257] :
+                         term % 3 == 0 ? boundary_fp12<alt_bn128_254_fp12>(term % 6) :
+                                         random_fp12<alt_bn128_254_fp12>(rng);
+            indices[term] = 2 * (terms - 1 - term);
+            const auto prepared_left = kernel::prepare_product_operand(left[term]);
+            for (std::size_t c = 0; c < coordinates; ++c) {
+                const auto digits = radix::from_montgomery64(prepared_left[c]);
+                std::copy(digits.begin(), digits.end(), packed_left.begin() + (c * terms + term) * 5);
+            }
+            for (std::size_t lane = 0; lane < lanes; ++lane) {
+                // Zero first/last lanes also model empty coefficients and a padded panel tail.
+                right[term][lane] = lane == 0 || lane + 1 == lanes ? value::zero() :
+                                    term >= 257                    ? right[term - 257][lane] :
+                                                                     random_fp12<alt_bn128_254_fp12>(rng);
+                const auto prepared_right = kernel::prepare_product_operand(right[term][lane]);
+                for (std::size_t c = 0; c < coordinates; ++c) {
+                    const auto digits = radix::from_montgomery64(prepared_right[c]);
+                    for (std::size_t digit = 0; digit < 5; ++digit) {
+                        packed_right[c * panel_stride + (indices[term] * 5 + digit) * lanes + lane] = digits[digit];
+                    }
+                }
+            }
+        }
+        std::array<value, lanes> expected;
+        expected.fill(value::zero());
+        std::size_t verified_terms = 0;
+        for (std::size_t count : {0, 1, 123, 255, 256, 257, 513, 514}) {
+            for (; verified_terms < count; ++verified_terms) {
+                for (std::size_t lane = 0; lane < lanes; ++lane) {
+                    expected[lane] +=
+                        generic_fp12_mul<alt_bn128_254_fp12>(left[verified_terms], right[verified_terms][lane]);
+                }
+            }
+            std::array<fast::limb, coordinates * 5 * lanes> products;
+            std::array<fast::limb, 12 * 5 * lanes> combined;
+            for (std::size_t c = 0; c < coordinates; ++c) {
+                fast::montgomery_inner_product52<params, Operations>(
+                    products.data() + c * 5 * lanes, packed_left.data() + c * terms * 5,
+                    packed_right.data() + c * panel_stride, indices.data(), count);
+            }
+            kernel::template recombine_products52<Operations>(combined.data(), products.data());
+            for (std::size_t lane = 0; lane < lanes; ++lane) {
+                value actual;
+                for (std::size_t c = 0; c < 12; ++c) {
+                    typename radix::value_type packed;
+                    for (std::size_t digit = 0; digit < 5; ++digit) {
+                        packed[digit] = combined[(c * 5 + digit) * lanes + lane];
+                        BOOST_CHECK_LE(packed[digit], radix::mask);
+                    }
+                    BOOST_CHECK(std::lexicographical_compare(packed.rbegin(), packed.rend(), radix::modulus().rbegin(),
+                                                             radix::modulus().rend()));
+                    radix::to_montgomery64(actual.coordinate(c).data.backend().base_data(), packed);
+                }
+                BOOST_CHECK_EQUAL(actual, expected[lane]);
+                if (count == terms) {
+                    BOOST_CHECK_EQUAL(actual, value::zero());
+                }
+            }
+        }
+    }
+
 }    // namespace
 
 BOOST_AUTO_TEST_SUITE(fp12_tests)
@@ -619,7 +813,21 @@ BOOST_AUTO_TEST_CASE(radix52_four_and_eight_lane_arithmetic_matches_reference) {
     check_radix52_dots<scalar_ifma52_operations<8>>();
 }
 
+BOOST_AUTO_TEST_CASE(radix52_rank23_matches_generic_fp12_in_four_and_eight_lanes) {
+    check_rank23_recombination<scalar_ifma52_operations<4>>();
+    check_rank23_recombination<scalar_ifma52_operations<8>>();
+    check_rank23_fp12_dots<scalar_ifma52_operations<4>>();
+    check_rank23_fp12_dots<scalar_ifma52_operations<8>>();
+}
+
 #if defined(CRYPTO3_HAS_AVX512IFMA)
+BOOST_AUTO_TEST_CASE(radix52_ifma_rank23_matches_generic_fp12) {
+    check_rank23_recombination<fields::detail::fp12_fast::ifma52_operations<4>>();
+    check_rank23_recombination<fields::detail::fp12_fast::ifma52_operations<8>>();
+    check_rank23_fp12_dots<fields::detail::fp12_fast::ifma52_operations<4>>();
+    check_rank23_fp12_dots<fields::detail::fp12_fast::ifma52_operations<8>>();
+}
+
 BOOST_AUTO_TEST_CASE(radix52_ifma_matches_reference_in_every_lane) {
     check_radix52_dots<fields::detail::fp12_fast::ifma52_operations<4>>();
     check_radix52_dots<fields::detail::fp12_fast::ifma52_operations<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 a2ae2080e..957a5898c 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
@@ -54,6 +54,7 @@ namespace nil::crypto3::math {
             using params = typename Value::field_type::extension_policy::fast_params;
             using kernel = algebra::fields::detail::fp12_fast::fp12_fast<params>;
             using prepared_type = typename kernel::prepared_product_operand;
+            using simd_kernel = typename kernel::rank23;
             using limb = algebra::fields::detail::fp12_fast::limb;
 
         public:
@@ -85,7 +86,7 @@ namespace nil::crypto3::math {
                     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) {
+                if constexpr (simd_kernel::supported) {
                     // 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 &&
@@ -156,7 +157,7 @@ namespace nil::crypto3::math {
                                  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 constexpr (simd_kernel::supported) {
                     if (simd_lanes_ == 4) {
                         evaluate_simd_blocks<4>(output, outer, first_block, block_count, scratch);
                         return;
@@ -218,7 +219,7 @@ namespace nil::crypto3::math {
             bool prepare_simd(const std::vector<Polynomial> &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<params>;
-                constexpr std::size_t coordinates = kernel::product_coordinate_count;
+                constexpr std::size_t coordinates = simd_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());
@@ -234,7 +235,7 @@ namespace nil::crypto3::math {
                         if (coefficient >= powers[exponent].size()) {
                             continue;
                         }
-                        const auto prepared = kernel::prepare_product_operand(powers[exponent][coefficient]);
+                        const auto prepared = simd_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;
@@ -255,7 +256,7 @@ namespace nil::crypto3::math {
                                       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<params>;
-                constexpr std::size_t coordinates = kernel::product_coordinate_count;
+                constexpr std::size_t coordinates = simd_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);
@@ -274,7 +275,7 @@ namespace nil::crypto3::math {
                     }
                     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]);
+                        const auto prepared = simd_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(),
@@ -287,9 +288,8 @@ namespace nil::crypto3::math {
                 }
                 for (std::size_t first = 0; any_active && first < width_;) {
                     const std::size_t active_lanes = std::min(Lanes, width_ - first);
-                    std::array<std::array<std::array<typename kernel::base_value_type, coordinates>, Lanes>,
-                               block_tile_size>
-                        products;
+                    // Keep the dot products in canonical radix-52 form until tower recombination.
+                    std::array<std::array<limb, coordinates * 5 * 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;
@@ -299,25 +299,28 @@ namespace nil::crypto3::math {
                             if (indices.empty()) {
                                 continue;
                             }
-                            std::array<limb, 5 * Lanes> reduced;
                             fast::montgomery_inner_product52<params, fast::ifma52_operations<Lanes>>(
-                                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);
-                            }
+                                products[block].data() + coordinate * 5 * Lanes,
+                                scratch.simd_scalars[block].data() + coordinate * indices.size() * 5, right,
+                                indices.data(), indices.size());
                         }
                     }
                     for (std::size_t block = 0; block < block_count; ++block) {
                         if (!scratch.exponents[block].empty()) {
+                            std::array<limb, 12 * 5 * Lanes> combined;
+                            simd_kernel::template recombine_products52<fast::ifma52_operations<Lanes>>(
+                                combined.data(), products[block].data());
                             for (std::size_t lane = 0; lane < active_lanes; ++lane) {
-                                output[block][first + lane] +=
-                                    kernel::template recombine_products<Value>(products[block][lane]);
+                                Value value;
+                                for (std::size_t coordinate = 0; coordinate < 12; ++coordinate) {
+                                    typename radix::value_type packed;
+                                    for (std::size_t digit = 0; digit < 5; ++digit) {
+                                        packed[digit] = combined[(coordinate * 5 + digit) * Lanes + lane];
+                                    }
+                                    radix::to_montgomery64(value.coordinate(coordinate).data.backend().base_data(),
+                                                           packed);
+                                }
+                                output[block][first + lane] += value;
                             }
                         }
                     }
@@ -344,7 +347,8 @@ 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
+     * Supported Fp12 coefficients may additionally cache prepared base-field forms (23 for BN254 SIMD,
+     * 54 for the scalar tower kernel). 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.
diff --git a/libs/math/test/polynomial_composition.cpp b/libs/math/test/polynomial_composition.cpp
index 468106878..b45c71c19 100644
--- a/libs/math/test/polynomial_composition.cpp
+++ b/libs/math/test/polynomial_composition.cpp
@@ -664,7 +664,7 @@ BOOST_AUTO_TEST_CASE(prepared_simd_panels_match_scalar_at_lane_and_term_boundari
             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);
+                    ((width + lanes - 1) / lanes) * lanes * 23 * 17 * 5 * sizeof(std::uint64_t);
                 BOOST_CHECK_EQUAL(candidate.bytes(), expected);
             }
 #else