diff --git a/libs/math/docs/geometric_lagrange.md b/libs/math/docs/geometric_lagrange.md index 338dc9212..5cf375c1f 100644 --- a/libs/math/docs/geometric_lagrange.md +++ b/libs/math/docs/geometric_lagrange.md @@ -1,4 +1,4 @@ -# Exact geometric-domain Lagrange weights {#math_geometric_lagrange} +# Exact geometric domains and interpolation {#math_geometric_lagrange} @tableofcontents @@ -7,7 +7,8 @@ x_i = r^i, 0 <= i < m, where `r` is the field's configured geometric generator. The domain provides a linear-time field-element API for -evaluating every Lagrange basis polynomial at an arbitrary point. +evaluating every Lagrange basis polynomial at an arbitrary point and a context-based API for exact geometric +interpolation. The relevant headers are: @@ -16,6 +17,9 @@ The relevant headers are: | Montgomery batch inversion | `` | | Evaluation-domain interface | `` | | Exact geometric domain | `` | +| Polynomial context | `` | +| Schoolbook backend | `` | +| Mixed-radix backend | `` | ## Domain precomputation @@ -30,6 +34,13 @@ For the Lagrange path, the constructor precomputes: Z(X) = product(X - x_i, i = 0 .. m - 1). +For exact geometric interpolation, it also precomputes: + +* the triangular powers `T_i = r^(i * (i - 1) / 2)` and their inverses; +* the denominator products `D_0 = 1` and `D_i = product(1 - r^j, j = 1 .. i)`; +* the inverses `1 / D_i`; and +* the fixed interpolation kernel `T_i / D_i`. + For general points, constructing every `w_i` from pairwise differences would take quadratic work. Geometric points instead satisfy the recurrence @@ -53,6 +64,59 @@ If `r` has exact order `m`, the final denominator vanishes and the implementatio Construction therefore takes `O(m)` field operations, one field inversion, and `O(m)` stored field elements. +## Exact geometric interpolation + +The concrete geometric domain provides a non-virtual, compile-time backend-selected API: + +```cpp +template +typename Backend::polynomial_type interpolate( + const std::vector& evaluations, + polynomial_arithmetic::polynomial_context& context +) const; +``` + +The input must contain exactly one evaluation for each point `1, r, ..., r^(m-1)`. Any other count throws +`std::invalid_argument`. The domain points remain `FieldType::value_type`, while the evaluations and returned +coefficients may belong to a compatible extension field such as Fq12. + +Interpolation proceeds in seven stages: + +1. Validate that the input contains exactly `m` evaluations. +2. Build `scaled[i] = evaluations[i] * (-1)^i / D_i` and embed the fixed kernel `T_i / D_i` into the backend's + coefficient field. +3. Multiply `scaled` by the fixed kernel through the supplied context and retain coefficients `0` through `m-1`. +4. Recover each Newton coefficient by multiplying convolution coefficient `i` by `1 / T_i`, then form the dynamic + Newton input by multiplying it by `D_i`. +5. Form the fixed Newton-to-monomial kernel `(-1)^i * T_i / D_i` and embed it in reverse order. +6. Multiply the reversed fixed kernel by the dynamic Newton input through the same supplied context. +7. Set output coefficient `i` to product coefficient `m - 1 + i` multiplied by `1 / D_i`, then remove trailing + zeros. + +The second product is transposed multiplication: the fixed Newton-to-monomial kernel is reversed, not the dynamic +Newton input. + +A schoolbook context requires no transform configuration: + +```cpp +using backend_type = polynomial_arithmetic::schoolbook_backend; +polynomial_arithmetic::polynomial_context context; + +const auto coefficients = domain.interpolate(evaluations, context); +``` + +A mixed-radix context must use a valid transform order supporting at least `2 * m - 1` product coefficients: + +```cpp +using backend_type = polynomial_arithmetic::mixed_radix_backend; +polynomial_arithmetic::polynomial_context context {backend_type(transform_order)}; + +const auto coefficients = domain.interpolate(evaluations, context); +``` + +Both products use the caller's context. Interpolation does not construct a backend, invoke the legacy transform, or +perform field inversions. + ## Evaluating all weights The combined overload returns both the weights and `Z(t)` without repeating the product: @@ -89,9 +153,11 @@ a compatible extension field. | Operation | Rough current cost | |---|---| | Construct a size-`m` domain | `O(m)` field operations and one inversion | +| Interpolate `m` evaluations | Two backend products, `O(m)` additional field operations, and no inversions | | Evaluate all weights at one off-domain point | `O(m)` field operations and one inversion | | Evaluate at a domain point | `O(m)` work and no inversion | | Evaluate all weights at `k` independent points | `O(k * m)` field operations and at most `k` inversions | -Returning `m` weights already requires linear output work. All scratch storage used by an evaluation is local, while -the domain precomputation is immutable, so concurrent weight evaluations on one domain do not share mutable scratch. +Returning `m` weights already requires linear output work. All scratch storage used by an evaluation or interpolation +is local, while the domain precomputation is immutable. Concurrent interpolation calls may share one domain, but each +call requires a separate polynomial context because backends may reuse mutable plans, caches, or scratch storage. diff --git a/libs/math/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp b/libs/math/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp index a6bed5a53..ee0a78d30 100644 --- a/libs/math/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp +++ b/libs/math/include/nil/crypto3/math/domains/geometric_sequence_domain.hpp @@ -32,6 +32,7 @@ #include #include +#include #include #include @@ -53,8 +54,14 @@ namespace nil { * and one field inversion. Therefore, after constructing one reusable domain, evaluating the weights at k * points takes O(m * k), not quadratic work per point. * - * This complexity guarantee applies to constructor precomputation and the field-element Lagrange - * overloads. The legacy transforms and powers-based overload are documented separately below. + * The concrete interpolate overload reconstructs a degree-below-m polynomial from exactly one evaluation + * at each domain point. Its evaluations and coefficients may belong to a compatible extension field. It + * reuses cached geometric factors, performs no interpolation-time field inversions, and routes both + * polynomial products through the caller's reusable polynomial context. + * + * These guarantees apply to constructor precomputation, the field-element Lagrange overloads, and the + * context-based interpolation overload. The legacy transforms and powers-based overload are documented + * separately below. */ template class geometric_sequence_domain : public evaluation_domain { @@ -65,6 +72,10 @@ namespace nil { field_value_type generator; std::vector geometric_sequence; std::vector geometric_triangular_sequence; + std::vector inverse_geometric_triangular_sequence; + std::vector interpolation_denominator_products; + std::vector inverse_interpolation_denominator_products; + std::vector interpolation_kernel; std::vector barycentric_weights; polynomial vanishing_polynomial; @@ -72,13 +83,14 @@ namespace nil { * Here m is the exact domain size: there are m points x_0, ..., x_(m-1), and transforms consume * exactly m coefficients or evaluations. * - * Build the reusable field data once, in four linear stages: + * Build the reusable field data once, in five linear stages: * 1. Generate r^i and r^(i(i-1)/2), while validating that the domain points are distinct. * 2. Batch-invert a packed denominator vector using one field inversion. Its first entry is r, * and entry i > 0 is 1 - r^i, so the result supplies both r^-1 and every inverse needed by * the recurrences below. - * 3. Derive the barycentric weights by recurrence. - * 4. Derive the coefficients of Z(X) = product_(i=0)^(m-1) (X - r^i) by recurrence. + * 3. Derive the geometric interpolation factors by recurrence. + * 4. Derive the barycentric weights by recurrence. + * 5. Derive the coefficients of Z(X) = product_(i=0)^(m-1) (X - r^i) by recurrence. * * The completed object is immutable and constructor-only scratch remains local, so concurrent * Lagrange evaluations share only read-only state. @@ -87,6 +99,10 @@ namespace nil { generator(fields::arithmetic_params::geometric_generator), geometric_sequence(m, field_value_type::zero()), geometric_triangular_sequence(m, field_value_type::zero()), + inverse_geometric_triangular_sequence(m, field_value_type::zero()), + interpolation_denominator_products(m, field_value_type::zero()), + inverse_interpolation_denominator_products(m, field_value_type::zero()), + interpolation_kernel(m, field_value_type::zero()), barycentric_weights(m, field_value_type::zero()), vanishing_polynomial(m + 1, field_value_type::zero()) { if (generator.is_zero()) { @@ -127,6 +143,29 @@ namespace nil { inverse_geometric_sequence[i] = inverse_geometric_sequence[i - 1] * inverse_denominators[0]; } + /* + * Cache every base-field factor used by geometric interpolation. With + * + * D_i = product_(j=1)^i (1 - r^j), + * + * inverse_denominators supplies the recurrence for 1 / D_i without any further field + * inversions. The interpolation kernel is reused for both products; the second product embeds + * it in reverse order with alternating signs. + */ + inverse_geometric_triangular_sequence[0] = field_value_type::one(); + interpolation_denominator_products[0] = field_value_type::one(); + inverse_interpolation_denominator_products[0] = field_value_type::one(); + interpolation_kernel[0] = field_value_type::one(); + for (std::size_t i = 1; i < m; ++i) { + inverse_geometric_triangular_sequence[i] = + inverse_geometric_triangular_sequence[i - 1] * inverse_geometric_sequence[i - 1]; + interpolation_denominator_products[i] = + interpolation_denominator_products[i - 1] * denominators[i]; + inverse_interpolation_denominator_products[i] = + inverse_interpolation_denominator_products[i - 1] * inverse_denominators[i]; + interpolation_kernel[i] = + geometric_triangular_sequence[i] * inverse_interpolation_denominator_products[i]; + } /* * Let Z(X) = product_(j=0)^(m-1) (X - x_j) be the domain's vanishing polynomial and let * Z'(X) be its formal derivative. At a domain point, @@ -199,6 +238,84 @@ namespace nil { evaluation_domain(validate_size(m)), precomputation_(m) { } + /** + * Interpolate one value at each geometric domain point using the caller's multiplication context. + * Domain-dependent factors remain in the base field while polynomial coefficients may belong to an + * extension field. Both convolutions use the same compile-time-selected backend instance. + * + * The implementation proceeds in seven stages: + * 1. Validate that there is exactly one evaluation per domain point. + * 2. Scale the evaluations and embed the fixed evaluation-to-Newton kernel. + * 3. Convolve them to recover the scaled Newton coefficients. + * 4. Remove the geometric triangular scaling and prepare the dynamic Newton input. + * 5. Embed the fixed Newton-to-monomial kernel in signed reverse order. + * 6. Convolve the reversed fixed kernel with the dynamic Newton input. + * 7. Extract and scale the transposed-product coefficients into canonical monomial form. + */ + template + typename Backend::polynomial_type + interpolate(const std::vector &evaluations, + polynomial_arithmetic::polynomial_context &context) const { + using polynomial_type = typename Backend::polynomial_type; + using coefficient_type = typename polynomial_type::value_type; + + // 1. Validate the exact evaluation count. + if (evaluations.size() != this->m) { + throw std::invalid_argument("geometric: expected one evaluation per domain point"); + } + + // 2. Scale the evaluations and embed the fixed evaluation-to-Newton kernel. + polynomial_type scaled_evaluations(this->m, coefficient_type::zero()); + polynomial_type embedded_interpolation_kernel(this->m, coefficient_type::zero()); + for (std::size_t i = 0; i < this->m; ++i) { + const field_value_type evaluation_factor = + i % 2 == 0 ? precomputation_.inverse_interpolation_denominator_products[i] : + -precomputation_.inverse_interpolation_denominator_products[i]; + scaled_evaluations[i] = evaluations[i] * evaluation_factor; + embedded_interpolation_kernel[i] = + coefficient_type::one() * precomputation_.interpolation_kernel[i]; + } + condense(scaled_evaluations); + + // 3. Recover the scaled Newton coefficients with the first convolution. + polynomial_type interpolation_convolution; + context.multiply(interpolation_convolution, scaled_evaluations, embedded_interpolation_kernel); + interpolation_convolution.resize(this->m, coefficient_type::zero()); + + // 4. Remove the triangular scaling and prepare the dynamic Newton input. + polynomial_type newton_input(this->m, coefficient_type::zero()); + for (std::size_t i = 0; i < this->m; ++i) { + const coefficient_type newton_coefficient = + interpolation_convolution[i] * precomputation_.inverse_geometric_triangular_sequence[i]; + newton_input[i] = newton_coefficient * precomputation_.interpolation_denominator_products[i]; + } + condense(newton_input); + + // 5. Embed the fixed Newton-to-monomial kernel in signed reverse order. + polynomial_type reversed_newton_basis_kernel(this->m, coefficient_type::zero()); + for (std::size_t i = 0; i < this->m; ++i) { + const std::size_t kernel_index = this->m - 1 - i; + const field_value_type kernel_value = kernel_index % 2 == 0 ? + precomputation_.interpolation_kernel[kernel_index] : + -precomputation_.interpolation_kernel[kernel_index]; + reversed_newton_basis_kernel[i] = coefficient_type::one() * kernel_value; + } + + // 6. Perform the transposed product, reversing the fixed kernel rather than the dynamic input. + polynomial_type monomial_convolution; + context.multiply(monomial_convolution, reversed_newton_basis_kernel, newton_input); + monomial_convolution.resize(2 * this->m - 1, coefficient_type::zero()); + + // 7. Extract, scale, and normalize the monomial coefficients. + polynomial_type result(this->m, coefficient_type::zero()); + for (std::size_t i = 0; i < this->m; ++i) { + result[i] = monomial_convolution[this->m - 1 + i] * + precomputation_.inverse_interpolation_denominator_products[i]; + } + condense(result); + return result; + } + /* * TODO: These legacy transforms are independent of the field-element Lagrange evaluation below. * They still perform individual inversions and use the generic polynomial-multiplication backend, diff --git a/libs/math/include/nil/crypto3/math/polynomial/operations/basis_change.hpp b/libs/math/include/nil/crypto3/math/polynomial/operations/basis_change.hpp index e21c9001d..48cce2552 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/operations/basis_change.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/operations/basis_change.hpp @@ -286,7 +286,19 @@ namespace nil { z[i] = -z[i]; } - w = transpose_multiplication(n - 1, w, u); + /* + * This transposed product reverses the fixed Newton-basis kernel u, not the dynamic coefficients w. + * Keep w as the algebraic multiplication operand so this legacy path continues to support + * extension-field and group-valued coefficients scaled by base-field elements. + */ + std::vector reversed_u(u); + reverse(reversed_u, n); + std::vector product; + multiplication(product, w, reversed_u); + product.resize(2 * n - 1, value_type::zero()); + for (std::size_t i = 0; i < n; ++i) { + w[i] = product[n - 1 + i]; + } for (std::size_t i = 0; i < n; i++) { a[i] = w[i] * z[i]; diff --git a/libs/math/test/geometric_sequence_domain.cpp b/libs/math/test/geometric_sequence_domain.cpp index 8bd080271..49c9100b0 100644 --- a/libs/math/test/geometric_sequence_domain.cpp +++ b/libs/math/test/geometric_sequence_domain.cpp @@ -31,21 +31,30 @@ #include #include +#include #include #include +#include #include #include #include +#include +#include +#include using namespace nil::crypto3; namespace { using bn254_fq = algebra::fields::alt_bn128<254>; - - template - ValueType evaluate(const Coefficients &coefficients, const ValueType &point) { - ValueType result = ValueType::zero(); + using bn254_fq12 = algebra::fields::fp12_2over3over2; + using fq_value_type = bn254_fq::value_type; + using fq12_value_type = bn254_fq12::value_type; + + template + typename Coefficients::value_type evaluate(const Coefficients &coefficients, const PointType &point) { + using coefficient_type = typename Coefficients::value_type; + coefficient_type result = coefficient_type::zero(); for (auto it = coefficients.rbegin(); it != coefficients.rend(); ++it) { result = result * point + *it; } @@ -68,6 +77,43 @@ namespace { return result; } + fq12_value_type fq12_value(std::size_t first_coordinate) { + fq12_value_type value = fq12_value_type::zero(); + for (std::size_t i = 0; i < bn254_fq12::arity; ++i) { + value.coordinate(i) = fq_value_type(first_coordinate + i); + } + return value; + } + + template + void check_backend_aware_interpolation(math::polynomial_arithmetic::polynomial_context &context) { + using polynomial_type = typename Backend::polynomial_type; + + constexpr std::size_t domain_size = 5; + const math::geometric_sequence_domain domain(domain_size); + const std::vector> coefficient_cases = { + {fq12_value(1), fq12_value(13), fq12_value(25), fq12_value(37), fq12_value(49)}, + {fq12_value(7), fq12_value(19), fq12_value(31), fq12_value_type::zero(), fq12_value_type::zero()}, + std::vector(domain_size, fq12_value_type::zero())}; + + for (const std::vector &coefficients : coefficient_cases) { + std::vector evaluations(domain_size, fq12_value_type::zero()); + for (std::size_t i = 0; i < domain_size; ++i) { + evaluations[i] = evaluate(coefficients, domain.get_domain_element(i)); + } + + polynomial_type expected(coefficients.begin(), coefficients.end()); + math::condense(expected); + const polynomial_type actual = domain.interpolate(evaluations, context); + BOOST_CHECK(actual == expected); + } + + const std::vector too_few(domain_size - 1, fq12_value_type::zero()); + const std::vector too_many(domain_size + 1, fq12_value_type::zero()); + BOOST_CHECK_THROW(domain.interpolate(too_few, context), std::invalid_argument); + BOOST_CHECK_THROW(domain.interpolate(too_many, context), std::invalid_argument); + } + } // namespace BOOST_AUTO_TEST_SUITE(geometric_sequence_domain_test_suite) @@ -248,6 +294,40 @@ BOOST_AUTO_TEST_CASE(lagrange_weights_interpolate_and_are_unit_vectors_on_the_do } } +BOOST_AUTO_TEST_CASE(backend_aware_interpolation_supports_fq12_coefficients) { + using schoolbook_backend = math::polynomial_arithmetic::schoolbook_backend; + using mixed_radix_backend = math::polynomial_arithmetic::mixed_radix_backend; + + BOOST_TEST_CONTEXT("schoolbook") { + math::polynomial_arithmetic::polynomial_context context; + check_backend_aware_interpolation(context); + } + BOOST_TEST_CONTEXT("mixed radix") { + // A size-five interpolation performs two length-five products, each requiring nine coefficients. + math::polynomial_arithmetic::polynomial_context context {mixed_radix_backend(9)}; + check_backend_aware_interpolation(context); + } +} + +BOOST_AUTO_TEST_CASE(legacy_inverse_fft_retains_the_correct_newton_basis_orientation) { + using field_type = algebra::fields::alt_bn128_scalar_field<254>; + using value_type = field_type::value_type; + + constexpr std::size_t domain_size = 5; + math::geometric_sequence_domain domain(domain_size); + math::evaluation_domain &abstract_domain = domain; + const std::vector coefficients = {value_type(3u), value_type(5u), value_type(7u), value_type(11u), + value_type(13u)}; + std::vector evaluations(domain_size, value_type::zero()); + for (std::size_t i = 0; i < domain_size; ++i) { + evaluations[i] = evaluate(coefficients, domain.get_domain_element(i)); + } + + abstract_domain.inverse_fft(evaluations); + + BOOST_CHECK(evaluations == coefficients); +} + BOOST_AUTO_TEST_CASE(add_poly_z_adds_the_scaled_vanishing_polynomial) { using value_type = bn254_fq::value_type;