From 8bcdeaac4f375d658e009eaaac90403829f5e450 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Fri, 28 Aug 2026 18:15:10 -0400 Subject: [PATCH 1/5] geometric interpolation --- libs/math/docs/geometric_lagrange.md | 74 +++++++++- .../domains/geometric_sequence_domain.hpp | 127 +++++++++++++++++- .../polynomial/operations/basis_change.hpp | 14 +- libs/math/test/geometric_sequence_domain.cpp | 88 +++++++++++- 4 files changed, 289 insertions(+), 14 deletions(-) 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; From 619ec14e19ea5f2962de8b69bacc64d590be4337 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Mon, 31 Aug 2026 17:23:10 -0400 Subject: [PATCH 2/5] x norm reconstruction --- libs/math/docs/polynomial_recovery.md | 95 ++++++--- .../polynomial_x_norm_reconstruction.hpp | 196 ++++++++++++++++++ libs/math/test/CMakeLists.txt | 1 + .../test/polynomial_x_norm_reconstruction.cpp | 189 +++++++++++++++++ 4 files changed, 452 insertions(+), 29 deletions(-) create mode 100644 libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp create mode 100644 libs/math/test/polynomial_x_norm_reconstruction.cpp diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md index 4df8901ac..d584bec0e 100644 --- a/libs/math/docs/polynomial_recovery.md +++ b/libs/math/docs/polynomial_recovery.md @@ -2,17 +2,18 @@ @tableofcontents -Crypto3.Math provides two complementary recovery primitives: +Crypto3.Math provides three complementary recovery facilities: -* square testing and square roots in a finite polynomial quotient field; and -* bounded rational reconstruction from a residue modulo a polynomial. +* square testing and square roots in a finite polynomial quotient field; +* bounded rational reconstruction from a residue modulo a polynomial; and +* one-call reconstruction of the fixed `X`-norm representation of an irreducible polynomial. Together they can recover representations of irreducible factors by the polynomial norm form P(X)^2 - X * Q(X)^2. -The library exposes the generic arithmetic operations. Applications remain responsible for selecting factors, degree -bounds, random sources, and any policy for combining or rejecting recovered representations. +The one-call operation owns the degree bounds and exact normalization for this norm equation. The caller supplies the +polynomial arithmetic context and the coefficient-field generator. ### Header map @@ -22,6 +23,7 @@ bounds, random sources, and any policy for combining or rejecting recovered repr | Coefficient-field square-root helpers | `` | | Quotient-field square testing and roots | `` | | Bounded rational reconstruction | `` | +| One-call `X`-norm reconstruction | `` | ## Field orders and coefficient square roots @@ -64,9 +66,9 @@ characterize squares if the quotient has zero divisors. `square_root_mod` uses Tonelli-Shanks in `K[X]/(B)`. Repeated calls should share two immutable precomputations: -1. `polynomial_divisor_context` caches the divisor and the inverse needed for modular reduction. -2. `polynomial_square_root_context` caches the decomposition - `order(K)^d - 1 = odd_order * 2^two_adicity` and a suitable quadratic nonresidue raised to `odd_order`. +* `polynomial_divisor_context` caches the divisor and the inverse needed for modular reduction. +* `polynomial_square_root_context` caches the decomposition + `order(K)^d - 1 = odd_order * 2^two_adicity` and a suitable quadratic nonresidue raised to `odd_order`. ```cpp #include @@ -146,25 +148,60 @@ reconstructs with bounds `degree(P) <= 0` and `degree(Q) <= 2` as P = 1, Q = 2 + 2X + X^2. -## Recovering a polynomial norm representation +## One-call `X`-norm reconstruction -The square-root and reconstruction APIs fit together as follows. For one irreducible factor `B` of degree `d`: +`recover_polynomial_x_norm_representation` composes quotient-field square roots, bounded rational reconstruction, and +coefficient-field normalization to recover polynomials `P` and `Q` satisfying -1. Represent the indeterminate by `x = {0, 1}` and test whether it is a square modulo `B`. -2. If it is square, compute `R` such that `R^2 = X mod B`. -3. Rationally reconstruct `R` with + P^2 - X * Q^2 = g. - maximum_numerator_degree = floor(d / 2), - maximum_denominator_degree = floor((d - 1) / 2). +The name explicitly identifies `X` as the quadratic element: the represented element is `P + Q * sqrt(X)`, whose norm +is `P^2 - X * Q^2`. The input `g` must be a canonical nonconstant irreducible polynomial. Irreducibility is a caller +precondition and is not tested. - This produces `P = R * Q mod B`. -4. Squaring the congruence gives +```cpp +#include + +// arithmetic_context and coefficient_generator are caller-owned. +auto representation = math::recover_polynomial_x_norm_representation( + g, arithmetic_context, coefficient_generator); + +if (representation) { + polynomial_type exact_norm = math::evaluate_polynomial_x_norm( + *representation, arithmetic_context); + // exact_norm == g +} +``` + +The coefficient generator returns coefficient-field values. The recovery operation adapts those values into canonical +degree-below-`degree(g)` representatives of `K[X]/(g)`. The generator remains caller-owned and must eventually supply +coefficients forming a quotient-field nonsquare. - P^2 - X * Q^2 = 0 mod B. +On success, the optional contains `polynomial_x_norm_representation`. Its `p` and `q` members hold the +two recovered coefficient polynomials. - The degree bounds make the left side have degree at most `d`, so it is a scalar multiple of `B`. +The function reduces `{0, 1}` modulo `g`, including when `g` is linear, and recovers `R` with `R^2 = X mod g`. It then +reconstructs `P = R * Q mod g` with + + degree(P) <= floor(degree(g) / 2), + degree(Q) <= floor((degree(g) - 1) / 2). + +These bounds imply + + P^2 - X * Q^2 = lambda * g + +for a coefficient-field scalar `lambda`. A nonzero square `lambda` is removed by scaling both outputs by +`sqrt(lambda^-1)`. The normalized norm is evaluated again and compared exactly with `g` before success is returned. +Both polynomial squares use the caller's arithmetic context; multiplication by `X` is a coefficient shift. + +| Outcome | Contract | +|---|---| +| Representation returned | The degree bounds hold and `evaluate_polynomial_x_norm(result, context) == g`. | +| No value returned | `X` is nonsquare modulo `g`, bounded reconstruction fails, or `lambda` is zero or nonsquare. | +| `std::invalid_argument` | The input is empty, noncanonical, zero, or constant, or a composed API contract is violated. | +| `std::logic_error` | An operation reported success but a required modular, scalar-multiple, or final exact identity is inconsistent. | -This is the local representation of `B` by the norm from adjoining a square root of `X`. Representations compose +This is the local representation of `g` by the norm from adjoining a square root of `X`. Representations compose multiplicatively: (P1^2 - X Q1^2) * (P2^2 - X Q2^2) @@ -172,8 +209,8 @@ multiplicatively: Consequently, a caller can factor a target polynomial, recover eligible irreducible factors independently, account for multiplicities and the scalar leading coefficient, and combine the local representations. Crypto3.Math deliberately -keeps that application-level policy separate from the generic factorization, quotient-field square-root, and bounded -rational-reconstruction primitives. +keeps that application-level policy separate from the generic factorization, quotient-field square-root, bounded +rational-reconstruction, and one-call `X`-norm reconstruction facilities. ### Worked example over F7 @@ -211,16 +248,16 @@ Then the scalar is corrected exactly: P'^2 - X Q'^2 = H. -This example performs the same local steps required for a higher-degree factor: factor, test whether `X` is square in -the quotient, compute its root, reconstruct bounded `P` and `Q`, and normalize the remaining coefficient-field -scalar. Combining several eligible factors then uses the multiplicative identity above. +The one-call API performs this square test, quotient-field square root, bounded reconstruction, scalar normalization, +and final exact verification using the supplied polynomial context and coefficient generator. Combining several +eligible factors then uses the multiplicative identity above. ## Reuse and performance -The [polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) describes these contexts in detail. Reuse -one polynomial arithmetic context throughout a recovery, one divisor context for every operation modulo the same `B`, -and one square-root context for repeated roots modulo that divisor. This avoids rebuilding divisor inverses, -multiplicative-group decompositions, nonresidue powers, backend plans, and scratch storage. +The [polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) describes these contexts in detail. The +one-call API reuses the caller's polynomial arithmetic context throughout recovery. It constructs one divisor context +and one square-root context and reuses them for all operations within that call. Direct users of the lower-level APIs +can retain those contexts across repeated operations modulo the same divisor. Let `Q` be the coefficient-field order, `d` the irreducible divisor degree, and `s` the two-adicity of `Q^d - 1`. The current implementations have the rough bounds shown below. diff --git a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp new file mode 100644 index 000000000..9143d2bc7 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp @@ -0,0 +1,196 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#ifndef CRYPTO3_MATH_POLYNOMIAL_X_NORM_RECONSTRUCTION_HPP +#define CRYPTO3_MATH_POLYNOMIAL_X_NORM_RECONSTRUCTION_HPP + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace nil::crypto3::math { + + /** Coefficients of P + Q * sqrt(X), represented by the polynomials P and Q. */ + template + struct polynomial_x_norm_representation { + Polynomial p; + Polynomial q; + }; + + /** + * Evaluate the polynomial norm of P + Q * sqrt(X): + * + * (P + Q * sqrt(X)) * (P - Q * sqrt(X)) = P^2 - X * Q^2. + * + * Both squares use the caller-owned arithmetic context. Multiplication by X is a coefficient shift and therefore + * does not require a polynomial product. + * + * @throws std::invalid_argument if P or Q is empty or noncanonical. + */ + template + typename Backend::polynomial_type evaluate_polynomial_x_norm( + const polynomial_x_norm_representation &representation, + polynomial_arithmetic::polynomial_context &arithmetic_context) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename polynomial_type::value_type; + + const auto is_canonical = [](const polynomial_type &polynomial) { + return !polynomial.empty() && + (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); + }; + if (!is_canonical(representation.p) || !is_canonical(representation.q)) { + throw std::invalid_argument("polynomial X-norm evaluation requires canonical nonempty inputs"); + } + + polynomial_type p_squared; + polynomial_type q_squared; + polynomial_type x_q_squared; + polynomial_type norm; + arithmetic_context.square(p_squared, representation.p); + arithmetic_context.square(q_squared, representation.q); + shift_left(x_q_squared, q_squared, 1); + subtraction(norm, p_squared, x_q_squared); + return norm; + } + + /** + * Recover P and Q satisfying + * + * P^2 - X * Q^2 = g, + * + * with degree(P) at most floor(degree(g) / 2) and degree(Q) at most + * floor((degree(g) - 1) / 2). The coefficient generator supplies field elements used to find a quadratic + * nonresidue in K[X]/(g); it remains owned by the caller. + * + * Irreducibility of g is a caller precondition and is not tested. The coefficient generator must eventually produce + * coefficients forming a nonsquare canonical representative of degree below degree(g). + * + * @return a normalized representation whose evaluated norm is exactly g; no value if X is nonsquare modulo g, + * bounded rational reconstruction fails, or the resulting nonzero scalar multiple of g cannot be + * normalized. + * @throws std::invalid_argument if g is empty, noncanonical, zero, or constant, or another documented precondition + * of a composed polynomial operation is violated. + * @throws std::logic_error if a composed operation reports success but its resulting identities are inconsistent. + */ + template + requires algebra::FieldValue && + std::constructible_from && + requires(typename Backend::polynomial_type &polynomial, Generator &generator, + const typename Backend::polynomial_type::value_type &value) { + polynomial[0] = generator(); + { value.is_square() } -> std::convertible_to; + } + std::optional> + recover_polynomial_x_norm_representation(const typename Backend::polynomial_type &g, + polynomial_arithmetic::polynomial_context &arithmetic_context, + Generator &coefficient_generator) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename polynomial_type::value_type; + using representation_type = polynomial_x_norm_representation; + + if (g.empty() || (g.size() > 1 && g[g.size() - 1] == value_type {})) { + throw std::invalid_argument("polynomial X-norm recovery requires a canonical nonempty polynomial"); + } + if (g.size() == 1) { + throw std::invalid_argument("polynomial X-norm recovery requires a nonconstant polynomial"); + } + + const std::size_t degree = g.size() - 1; + const std::size_t inverse_precision = std::max(1, degree - 1); + polynomial_divisor_context divisor_context(g, inverse_precision, arithmetic_context); + + // Work with the canonical representative of X in K[X]/(g). Reduction is necessary when g is linear. + const polynomial_type x = {value_type::zero(), value_type::one()}; + polynomial_type x_mod_g; + remainder(x_mod_g, x, divisor_context, arithmetic_context); + if (!is_square_mod(x_mod_g, divisor_context, arithmetic_context)) { + return std::nullopt; + } + + // Adapt caller-generated field coefficients into canonical representatives of K[X]/(g). No random source or + // polynomial backend is constructed here. + auto quotient_representative_generator = [&]() { + polynomial_type representative(degree); + for (std::size_t index = 0; index < degree; ++index) { + representative[index] = coefficient_generator(); + } + condense(representative); + return representative; + }; + polynomial_square_root_context square_root_context(divisor_context, arithmetic_context, + quotient_representative_generator); + + // Recover R with R^2 = X mod g. A failure here contradicts the successful square test above. + polynomial_type root; + if (!square_root_mod(root, x_mod_g, square_root_context, arithmetic_context)) { + throw std::logic_error("polynomial X-norm recovery failed after X was reported square modulo g"); + } + + // Reconstruct P = R * Q mod g with the unique degree bounds required by the norm equation. + polynomial_type p; + polynomial_type q; + if (!rational_reconstruct(p, q, root, g, degree / 2, (degree - 1) / 2, arithmetic_context)) { + return std::nullopt; + } + + // The two modular identities imply P^2 - X * Q^2 = lambda * g. The degree bounds ensure that lambda is a + // scalar. A nonzero remainder or nonconstant quotient would contradict those successful operations. + representation_type representation {std::move(p), std::move(q)}; + const polynomial_type norm = evaluate_polynomial_x_norm(representation, arithmetic_context); + polynomial_type scalar_quotient; + polynomial_type scalar_remainder; + divrem(scalar_quotient, scalar_remainder, norm, divisor_context, arithmetic_context); + if (!is_zero(scalar_remainder) || scalar_quotient.size() != 1) { + throw std::logic_error("polynomial X-norm reconstruction produced an inconsistent scalar multiple"); + } + + const value_type lambda = scalar_quotient[0]; + if (lambda.is_zero() || !lambda.is_square()) { + return std::nullopt; + } + + // Multiplying P and Q by sqrt(lambda^-1) changes their norm from lambda * g to exactly g. + const value_type normalization = algebra::fields::sqrt_known_square(lambda.inversed()); + scalar_multiplication(representation.p, representation.p, normalization); + scalar_multiplication(representation.q, representation.q, normalization); + + if (evaluate_polynomial_x_norm(representation, arithmetic_context) != g) { + throw std::logic_error("normalized polynomial X-norm representation failed exact verification"); + } + return std::move(representation); + } + +} // namespace nil::crypto3::math + +#endif // CRYPTO3_MATH_POLYNOMIAL_X_NORM_RECONSTRUCTION_HPP diff --git a/libs/math/test/CMakeLists.txt b/libs/math/test/CMakeLists.txt index 21c5b09f3..94c805953 100644 --- a/libs/math/test/CMakeLists.txt +++ b/libs/math/test/CMakeLists.txt @@ -55,6 +55,7 @@ set(TESTS_NAMES "polynomial_exponentiation" "polynomial_rational_reconstruction" "polynomial_square_root" + "polynomial_x_norm_reconstruction" "polynomial_composition" "polynomial_frobenius" "square_free_factorization" diff --git a/libs/math/test/polynomial_x_norm_reconstruction.cpp b/libs/math/test/polynomial_x_norm_reconstruction.cpp new file mode 100644 index 000000000..b999f0eda --- /dev/null +++ b/libs/math/test/polynomial_x_norm_reconstruction.cpp @@ -0,0 +1,189 @@ +//---------------------------------------------------------------------------// +// Copyright (c) 2026 +// +// MIT License +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +//---------------------------------------------------------------------------// + +#define BOOST_TEST_MODULE polynomial_x_norm_reconstruction_test + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace { + namespace fields = nil::crypto3::algebra::fields; + namespace math = nil::crypto3::math; + namespace polynomial_arithmetic = math::polynomial_arithmetic; + + using field_type = fields::babybear; + using value_type = field_type::value_type; + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + + using fq_field_type = fields::alt_bn128_base_field<254>; + using fq_value_type = fq_field_type::value_type; + using fq12_field_type = fields::fp12_2over3over2>; + using fq12_value_type = fq12_field_type::value_type; + using fq12_schoolbook_backend = polynomial_arithmetic::schoolbook_backend; + using fq12_mixed_radix_backend = polynomial_arithmetic::mixed_radix_backend; + + value_type first_quadratic_non_residue() { + value_type candidate(2); + while (candidate.is_square()) { + candidate = candidate + value_type::one(); + } + return candidate; + } + + fq12_value_type fq12_scalar(std::size_t value) { + fq12_value_type result = fq12_value_type::zero(); + result.coordinate(0) = fq_value_type(value); + return result; + } + + template + void check_bn254_fq12_linear_recovery(polynomial_arithmetic::polynomial_context &arithmetic_context, + std::size_t seed) { + using extension_polynomial_type = typename Backend::polynomial_type; + + const extension_polynomial_type g = {fq12_scalar(9), -fq12_scalar(4)}; + boost::random::mt19937 rng(seed); + auto coefficient_generator = [&] { return nil::crypto3::algebra::random_element(rng); }; + + const auto result = + math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + BOOST_REQUIRE(result.has_value()); + BOOST_CHECK_EQUAL(result->p.size(), 1); + BOOST_CHECK_EQUAL(result->q.size(), 1); + BOOST_CHECK(result->p[0] * result->p[0] == fq12_scalar(9)); + BOOST_CHECK(result->q[0] * result->q[0] == fq12_scalar(4)); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, arithmetic_context) == g); + } +} // namespace + +BOOST_AUTO_TEST_SUITE(polynomial_x_norm_reconstruction_test_suite) + +BOOST_AUTO_TEST_CASE(recovers_and_normalizes_an_irreducible_quadratic_with_the_stated_degree_bounds) { + polynomial_arithmetic::polynomial_context arithmetic_context; + const math::polynomial_x_norm_representation original = { + polynomial_type {value_type::one(), value_type::one()}, polynomial_type {value_type(9)}}; + const polynomial_type g = math::evaluate_polynomial_x_norm(original, arithmetic_context); + + // The discriminant of g is 81 * 77. Since 81 is square and 77 is nonsquare, g is irreducible. + BOOST_REQUIRE(!value_type(77).is_square()); + BOOST_REQUIRE(g == polynomial_type({value_type::one(), value_type::zero() - value_type(79), value_type::one()})); + + boost::random::mt19937 rng(0x584E4F52); + auto coefficient_generator = [&] { return nil::crypto3::algebra::random_element(rng); }; + const auto result = + math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + + BOOST_REQUIRE(result.has_value()); + BOOST_CHECK_LE(result->p.size() - 1, 1); + BOOST_CHECK_LE(result->q.size() - 1, 0); + // Rational reconstruction makes Q monic, initially producing lambda = 1/81. Exact normalization restores Q^2 = 81. + BOOST_CHECK(result->q[0] != value_type::one()); + BOOST_CHECK(result->q[0] * result->q[0] == value_type(81)); + BOOST_CHECK(math::evaluate_polynomial_x_norm(*result, arithmetic_context) == g); +} + +BOOST_AUTO_TEST_CASE(returns_no_value_when_x_is_nonsquare_modulo_the_irreducible_polynomial) { + const value_type non_residue = first_quadratic_non_residue(); + BOOST_REQUIRE((value_type::zero() - value_type::one()).is_square()); + const polynomial_type g = {value_type::zero() - non_residue, value_type::zero(), value_type::one()}; + polynomial_arithmetic::polynomial_context arithmetic_context; + std::size_t generator_calls = 0; + auto coefficient_generator = [&] { + ++generator_calls; + return value_type::zero(); + }; + + const auto result = + math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + BOOST_CHECK(!result.has_value()); + BOOST_CHECK_EQUAL(generator_calls, 0); +} + +BOOST_AUTO_TEST_CASE(returns_no_value_when_the_scalar_multiple_cannot_be_normalized) { + const value_type non_residue = first_quadratic_non_residue(); + const value_type minus_one = value_type::zero() - value_type::one(); + BOOST_REQUIRE(minus_one.is_square()); + + // X is one modulo non_residue * (X - 1), but the reconstructed norm is multiplied by + // lambda = -non_residue^-1. This lambda is nonsquare and cannot be removed by scaling P and Q. + const polynomial_type g = {value_type::zero() - non_residue, non_residue}; + polynomial_arithmetic::polynomial_context arithmetic_context; + auto coefficient_generator = [&] { return non_residue; }; + + const auto result = + math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + BOOST_CHECK(!result.has_value()); +} + +BOOST_AUTO_TEST_CASE(recovers_a_linear_bn254_fq12_norm_with_the_schoolbook_backend) { + polynomial_arithmetic::polynomial_context arithmetic_context; + check_bn254_fq12_linear_recovery(arithmetic_context, 0xF0125001); +} + +BOOST_AUTO_TEST_CASE(recovers_a_linear_bn254_fq12_norm_with_the_mixed_radix_backend) { + polynomial_arithmetic::polynomial_context arithmetic_context { + fq12_mixed_radix_backend(3)}; + check_bn254_fq12_linear_recovery(arithmetic_context, 0xF0125002); +} + +BOOST_AUTO_TEST_CASE(rejects_malformed_zero_and_constant_inputs) { + polynomial_arithmetic::polynomial_context arithmetic_context; + auto coefficient_generator = [] { return value_type::one(); }; + + polynomial_type empty; + empty.get_storage().clear(); + BOOST_CHECK_THROW( + math::recover_polynomial_x_norm_representation(empty, arithmetic_context, coefficient_generator), + std::invalid_argument); + + polynomial_type noncanonical(2); + noncanonical[0] = value_type::one(); + noncanonical[1] = value_type::zero(); + BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation(noncanonical, arithmetic_context, + coefficient_generator), + std::invalid_argument); + + BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation( + polynomial_type {value_type::zero()}, arithmetic_context, coefficient_generator), + std::invalid_argument); + BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation( + polynomial_type {value_type::one()}, arithmetic_context, coefficient_generator), + std::invalid_argument); +} + +BOOST_AUTO_TEST_SUITE_END() From 6f6fabd2d1a236dd0961f221198240e8fd61eebd Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Mon, 31 Aug 2026 17:32:53 -0400 Subject: [PATCH 3/5] docs --- libs/math/docs/polynomial_recovery.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md index d584bec0e..82dd6c431 100644 --- a/libs/math/docs/polynomial_recovery.md +++ b/libs/math/docs/polynomial_recovery.md @@ -6,13 +6,13 @@ Crypto3.Math provides three complementary recovery facilities: * square testing and square roots in a finite polynomial quotient field; * bounded rational reconstruction from a residue modulo a polynomial; and -* one-call reconstruction of the fixed `X`-norm representation of an irreducible polynomial. +* recovery of the fixed `X`-norm representation of an irreducible polynomial. Together they can recover representations of irreducible factors by the polynomial norm form P(X)^2 - X * Q(X)^2. -The one-call operation owns the degree bounds and exact normalization for this norm equation. The caller supplies the +The recovery operation owns the degree bounds and exact normalization for this norm equation. The caller supplies the polynomial arithmetic context and the coefficient-field generator. ### Header map @@ -23,7 +23,7 @@ polynomial arithmetic context and the coefficient-field generator. | Coefficient-field square-root helpers | `` | | Quotient-field square testing and roots | `` | | Bounded rational reconstruction | `` | -| One-call `X`-norm reconstruction | `` | +| `X`-norm reconstruction | `` | ## Field orders and coefficient square roots @@ -148,7 +148,7 @@ reconstructs with bounds `degree(P) <= 0` and `degree(Q) <= 2` as P = 1, Q = 2 + 2X + X^2. -## One-call `X`-norm reconstruction +## `X`-norm reconstruction `recover_polynomial_x_norm_representation` composes quotient-field square roots, bounded rational reconstruction, and coefficient-field normalization to recover polynomials `P` and `Q` satisfying @@ -210,7 +210,7 @@ multiplicatively: Consequently, a caller can factor a target polynomial, recover eligible irreducible factors independently, account for multiplicities and the scalar leading coefficient, and combine the local representations. Crypto3.Math deliberately keeps that application-level policy separate from the generic factorization, quotient-field square-root, bounded -rational-reconstruction, and one-call `X`-norm reconstruction facilities. +rational-reconstruction, and `X`-norm reconstruction facilities. ### Worked example over F7 @@ -248,14 +248,14 @@ Then the scalar is corrected exactly: P'^2 - X Q'^2 = H. -The one-call API performs this square test, quotient-field square root, bounded reconstruction, scalar normalization, +The recovery API performs this square test, quotient-field square root, bounded reconstruction, scalar normalization, and final exact verification using the supplied polynomial context and coefficient generator. Combining several eligible factors then uses the multiplicative identity above. ## Reuse and performance The [polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) describes these contexts in detail. The -one-call API reuses the caller's polynomial arithmetic context throughout recovery. It constructs one divisor context +recovery API reuses the caller's polynomial arithmetic context throughout recovery. It constructs one divisor context and one square-root context and reuses them for all operations within that call. Direct users of the lower-level APIs can retain those contexts across repeated operations modulo the same divisor. From 12710cab705180288576c73bc19ea3e5988bc4b7 Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Mon, 31 Aug 2026 18:06:25 -0400 Subject: [PATCH 4/5] renaming --- libs/math/docs/polynomial_recovery.md | 6 ++-- .../polynomial_x_norm_reconstruction.hpp | 7 +++-- .../test/polynomial_x_norm_reconstruction.cpp | 30 +++++++++---------- 3 files changed, 22 insertions(+), 21 deletions(-) diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md index 82dd6c431..373f131de 100644 --- a/libs/math/docs/polynomial_recovery.md +++ b/libs/math/docs/polynomial_recovery.md @@ -150,8 +150,8 @@ reconstructs with bounds `degree(P) <= 0` and `degree(Q) <= 2` as ## `X`-norm reconstruction -`recover_polynomial_x_norm_representation` composes quotient-field square roots, bounded rational reconstruction, and -coefficient-field normalization to recover polynomials `P` and `Q` satisfying +`recover_irreducible_polynomial_x_norm_representation` composes quotient-field square roots, bounded rational +reconstruction, and coefficient-field normalization to recover polynomials `P` and `Q` satisfying P^2 - X * Q^2 = g. @@ -163,7 +163,7 @@ precondition and is not tested. #include // arithmetic_context and coefficient_generator are caller-owned. -auto representation = math::recover_polynomial_x_norm_representation( +auto representation = math::recover_irreducible_polynomial_x_norm_representation( g, arithmetic_context, coefficient_generator); if (representation) { diff --git a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp index 9143d2bc7..a7782c3a7 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp @@ -112,9 +112,10 @@ namespace nil::crypto3::math { { value.is_square() } -> std::convertible_to; } std::optional> - recover_polynomial_x_norm_representation(const typename Backend::polynomial_type &g, - polynomial_arithmetic::polynomial_context &arithmetic_context, - Generator &coefficient_generator) { + recover_irreducible_polynomial_x_norm_representation( + const typename Backend::polynomial_type &g, + polynomial_arithmetic::polynomial_context &arithmetic_context, + Generator &coefficient_generator) { using polynomial_type = typename Backend::polynomial_type; using value_type = typename polynomial_type::value_type; using representation_type = polynomial_x_norm_representation; diff --git a/libs/math/test/polynomial_x_norm_reconstruction.cpp b/libs/math/test/polynomial_x_norm_reconstruction.cpp index b999f0eda..676f16229 100644 --- a/libs/math/test/polynomial_x_norm_reconstruction.cpp +++ b/libs/math/test/polynomial_x_norm_reconstruction.cpp @@ -80,8 +80,8 @@ namespace { boost::random::mt19937 rng(seed); auto coefficient_generator = [&] { return nil::crypto3::algebra::random_element(rng); }; - const auto result = - math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, arithmetic_context, + coefficient_generator); BOOST_REQUIRE(result.has_value()); BOOST_CHECK_EQUAL(result->p.size(), 1); BOOST_CHECK_EQUAL(result->q.size(), 1); @@ -105,8 +105,8 @@ BOOST_AUTO_TEST_CASE(recovers_and_normalizes_an_irreducible_quadratic_with_the_s boost::random::mt19937 rng(0x584E4F52); auto coefficient_generator = [&] { return nil::crypto3::algebra::random_element(rng); }; - const auto result = - math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, arithmetic_context, + coefficient_generator); BOOST_REQUIRE(result.has_value()); BOOST_CHECK_LE(result->p.size() - 1, 1); @@ -128,8 +128,8 @@ BOOST_AUTO_TEST_CASE(returns_no_value_when_x_is_nonsquare_modulo_the_irreducible return value_type::zero(); }; - const auto result = - math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, arithmetic_context, + coefficient_generator); BOOST_CHECK(!result.has_value()); BOOST_CHECK_EQUAL(generator_calls, 0); } @@ -145,8 +145,8 @@ BOOST_AUTO_TEST_CASE(returns_no_value_when_the_scalar_multiple_cannot_be_normali polynomial_arithmetic::polynomial_context arithmetic_context; auto coefficient_generator = [&] { return non_residue; }; - const auto result = - math::recover_polynomial_x_norm_representation(g, arithmetic_context, coefficient_generator); + const auto result = math::recover_irreducible_polynomial_x_norm_representation(g, arithmetic_context, + coefficient_generator); BOOST_CHECK(!result.has_value()); } @@ -167,21 +167,21 @@ BOOST_AUTO_TEST_CASE(rejects_malformed_zero_and_constant_inputs) { polynomial_type empty; empty.get_storage().clear(); - BOOST_CHECK_THROW( - math::recover_polynomial_x_norm_representation(empty, arithmetic_context, coefficient_generator), - std::invalid_argument); + BOOST_CHECK_THROW(math::recover_irreducible_polynomial_x_norm_representation( + empty, arithmetic_context, coefficient_generator), + std::invalid_argument); polynomial_type noncanonical(2); noncanonical[0] = value_type::one(); noncanonical[1] = value_type::zero(); - BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation(noncanonical, arithmetic_context, - coefficient_generator), + BOOST_CHECK_THROW(math::recover_irreducible_polynomial_x_norm_representation( + noncanonical, arithmetic_context, coefficient_generator), std::invalid_argument); - BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation( + BOOST_CHECK_THROW(math::recover_irreducible_polynomial_x_norm_representation( polynomial_type {value_type::zero()}, arithmetic_context, coefficient_generator), std::invalid_argument); - BOOST_CHECK_THROW(math::recover_polynomial_x_norm_representation( + BOOST_CHECK_THROW(math::recover_irreducible_polynomial_x_norm_representation( polynomial_type {value_type::one()}, arithmetic_context, coefficient_generator), std::invalid_argument); } From e76aa50252e7eac0a47c50cd4d2b59825a37a13e Mon Sep 17 00:00:00 2001 From: Riccardo Abbate Date: Mon, 31 Aug 2026 18:17:19 -0400 Subject: [PATCH 5/5] poly_x_norm_representation --- libs/math/docs/polynomial_recovery.md | 6 +- .../polynomial_x_norm_reconstruction.hpp | 66 +++++++++++++++++-- .../test/polynomial_x_norm_reconstruction.cpp | 38 +++++++++++ 3 files changed, 101 insertions(+), 9 deletions(-) diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md index 373f131de..3cf7d4b0b 100644 --- a/libs/math/docs/polynomial_recovery.md +++ b/libs/math/docs/polynomial_recovery.md @@ -201,12 +201,14 @@ Both polynomial squares use the caller's arithmetic context; multiplication by ` | `std::invalid_argument` | The input is empty, noncanonical, zero, or constant, or a composed API contract is violated. | | `std::logic_error` | An operation reported success but a required modular, scalar-multiple, or final exact identity is inconsistent. | -This is the local representation of `g` by the norm from adjoining a square root of `X`. Representations compose -multiplicatively: +This is the local representation of `g` by the norm from adjoining a square root of `X`. +`multiply_polynomial_x_norm_representations(left, right, context)` composes two representations using (P1^2 - X Q1^2) * (P2^2 - X Q2^2) = (P1 P2 + X Q1 Q2)^2 - X * (P1 Q2 + P2 Q1)^2. +All four polynomial products use the supplied arithmetic context, while multiplication by `X` is a coefficient shift. + Consequently, a caller can factor a target polynomial, recover eligible irreducible factors independently, account for multiplicities and the scalar leading coefficient, and combine the local representations. Crypto3.Math deliberately keeps that application-level policy separate from the generic factorization, quotient-field square-root, bounded diff --git a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp index a7782c3a7..77ebf2472 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_x_norm_reconstruction.hpp @@ -48,6 +48,21 @@ namespace nil::crypto3::math { Polynomial q; }; + namespace detail { + + template + bool is_canonical_polynomial_x_norm_representation( + const polynomial_x_norm_representation &representation) { + using value_type = typename Polynomial::value_type; + const auto is_canonical = [](const Polynomial &polynomial) { + return !polynomial.empty() && + (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); + }; + return is_canonical(representation.p) && is_canonical(representation.q); + } + + } // namespace detail + /** * Evaluate the polynomial norm of P + Q * sqrt(X): * @@ -63,13 +78,7 @@ namespace nil::crypto3::math { const polynomial_x_norm_representation &representation, polynomial_arithmetic::polynomial_context &arithmetic_context) { using polynomial_type = typename Backend::polynomial_type; - using value_type = typename polynomial_type::value_type; - - const auto is_canonical = [](const polynomial_type &polynomial) { - return !polynomial.empty() && - (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); - }; - if (!is_canonical(representation.p) || !is_canonical(representation.q)) { + if (!detail::is_canonical_polynomial_x_norm_representation(representation)) { throw std::invalid_argument("polynomial X-norm evaluation requires canonical nonempty inputs"); } @@ -84,6 +93,49 @@ namespace nil::crypto3::math { return norm; } + /** + * Multiply two polynomial X-norm representations using + * + * P = P1 * P2 + X * Q1 * Q2, + * Q = P1 * Q2 + Q1 * P2. + * + * The returned representation has norm equal to the product of the input norms. Every polynomial product uses the + * caller-owned arithmetic context; multiplication by X is a coefficient shift. + * + * @throws std::invalid_argument if an input polynomial is empty or noncanonical. + */ + template + polynomial_x_norm_representation multiply_polynomial_x_norm_representations( + const polynomial_x_norm_representation &left, + const polynomial_x_norm_representation &right, + polynomial_arithmetic::polynomial_context &arithmetic_context) { + using polynomial_type = typename Backend::polynomial_type; + using representation_type = polynomial_x_norm_representation; + + if (!detail::is_canonical_polynomial_x_norm_representation(left) || + !detail::is_canonical_polynomial_x_norm_representation(right)) { + throw std::invalid_argument("polynomial X-norm multiplication requires canonical nonempty inputs"); + } + + polynomial_type p_product; + polynomial_type q_product; + polynomial_type shifted_q_product; + polynomial_type result_p; + arithmetic_context.multiply(p_product, left.p, right.p); + arithmetic_context.multiply(q_product, left.q, right.q); + shift_left(shifted_q_product, q_product, 1); + addition(result_p, p_product, shifted_q_product); + + polynomial_type left_p_right_q; + polynomial_type left_q_right_p; + polynomial_type result_q; + arithmetic_context.multiply(left_p_right_q, left.p, right.q); + arithmetic_context.multiply(left_q_right_p, left.q, right.p); + addition(result_q, left_p_right_q, left_q_right_p); + + return representation_type {std::move(result_p), std::move(result_q)}; + } + /** * Recover P and Q satisfying * diff --git a/libs/math/test/polynomial_x_norm_reconstruction.cpp b/libs/math/test/polynomial_x_norm_reconstruction.cpp index 676f16229..89d6674e2 100644 --- a/libs/math/test/polynomial_x_norm_reconstruction.cpp +++ b/libs/math/test/polynomial_x_norm_reconstruction.cpp @@ -93,6 +93,44 @@ namespace { BOOST_AUTO_TEST_SUITE(polynomial_x_norm_reconstruction_test_suite) +BOOST_AUTO_TEST_CASE(multiplies_x_norm_representations_and_preserves_the_exact_norm_product) { + polynomial_arithmetic::polynomial_context arithmetic_context; + const math::polynomial_x_norm_representation left = { + polynomial_type {value_type::one(), value_type(2)}, polynomial_type {value_type(3)}}; + const math::polynomial_x_norm_representation right = { + polynomial_type {value_type(4), value_type(5)}, polynomial_type {value_type(6), value_type(7)}}; + + const auto product = + math::multiply_polynomial_x_norm_representations(left, right, arithmetic_context); + BOOST_CHECK(product.p == polynomial_type({value_type(4), value_type(31), value_type(31)})); + BOOST_CHECK(product.q == polynomial_type({value_type(18), value_type(34), value_type(14)})); + + const polynomial_type left_norm = math::evaluate_polynomial_x_norm(left, arithmetic_context); + const polynomial_type right_norm = math::evaluate_polynomial_x_norm(right, arithmetic_context); + polynomial_type expected_norm_product; + arithmetic_context.multiply(expected_norm_product, left_norm, right_norm); + BOOST_CHECK(math::evaluate_polynomial_x_norm(product, arithmetic_context) == expected_norm_product); + + const math::polynomial_x_norm_representation zero = {polynomial_type {value_type::zero()}, + polynomial_type {value_type::zero()}}; + const auto zero_product = + math::multiply_polynomial_x_norm_representations(zero, right, arithmetic_context); + BOOST_CHECK(zero_product.p == polynomial_type({value_type::zero()})); + BOOST_CHECK(zero_product.q == polynomial_type({value_type::zero()})); +} + +BOOST_AUTO_TEST_CASE(x_norm_representation_multiplication_rejects_malformed_inputs) { + polynomial_arithmetic::polynomial_context arithmetic_context; + const math::polynomial_x_norm_representation valid = {polynomial_type {value_type::one()}, + polynomial_type {value_type::zero()}}; + math::polynomial_x_norm_representation malformed = valid; + malformed.q.get_storage().clear(); + + BOOST_CHECK_THROW( + math::multiply_polynomial_x_norm_representations(valid, malformed, arithmetic_context), + std::invalid_argument); +} + BOOST_AUTO_TEST_CASE(recovers_and_normalizes_an_irreducible_quadratic_with_the_stated_degree_bounds) { polynomial_arithmetic::polynomial_context arithmetic_context; const math::polynomial_x_norm_representation original = {