Skip to content

Repository navigation

coefficient combination optimization #180

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
odraccir1 merged 3 commits into from
Oct 21, 2026
Merged

coefficient combination optimization #180

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
95cfd7e
Improve coefficient combination
odraccir1 Oct 21, 2026
152360e
Merge branch 'master' into ra/coefficient_combination_optimization
odraccir1 Oct 21, 2026
1a4dd93
fixing P3
odraccir1 Oct 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 126 additions & 0 deletions libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/fp12_fast.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#pragma once

#include <algorithm>
#include <array>
#include <cstddef>
#include <iterator>
Expand Down Expand Up @@ -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;
Expand Down
43 changes: 27 additions & 16 deletions libs/algebra/include/nil/crypto3/algebra/fields/detail/element/fp12_fast/limb_ops.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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;
Expand Down Expand Up @@ -639,15 +650,15 @@ 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) {
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);
montgomery_radix52_simd<Params, Operations>::subtract_modulus(total);
first += terms;
}
for (std::size_t i = 0; i < 5; ++i) {
Expand Down
Loading
Loading