diff --git a/libs/math/docs/polynomial_arithmetic.md b/libs/math/docs/polynomial_arithmetic.md index 10cd0d454..5409edb80 100644 --- a/libs/math/docs/polynomial_arithmetic.md +++ b/libs/math/docs/polynomial_arithmetic.md @@ -130,6 +130,23 @@ math::polynomial_divisor_context divisor_context( The context is immutable after construction. Operations receiving it assume it represents the intended divisor; when another precomputation also depends on `B`, callers must keep those contexts paired with the same polynomial. +### Fixed-operand FFT preparation + +For mixed-radix division, construction also prepares the FFT of the reversed-inverse prefix and one representation +of the divisor. These spectra own their data and belong to the divisor snapshot; changing the source polynomial +does not change them. The cache is bounded to at most two spectra and never grows during division. Constant divisors +and configurations selecting long division need no spectra. + +The backend exposes optional `prepare_low_product` / `try_multiply_low_prepared` operations through the arithmetic +context. A preparation records the fixed prefix, precision and transform length. Reuse requires the same precision +and a compatible transform; otherwise the divisor context uses ordinary `multiply_low`. Changing quotient lengths +does not build more cache entries. Backends such as schoolbook need only implement the three required operations. + +The second spectrum is either the original low-product divisor prefix or, when a strictly shorter supported +transform exists, the complete divisor folded for cyclic remainder reconstruction. Both are never retained together. +Preparation pays one forward FFT per nonzero fixed operand up front, plus spectrum storage and allocation. Repeated +compatible operations then avoid those forward FFTs. This trades setup time and memory for lower repeated cost. + ## Polynomial division `divrem` computes canonical `quotient` and `remainder` satisfying @@ -149,6 +166,35 @@ reversal turns division into the truncated product The divisor context must have inverse precision of at least `k` on this fast path. Long division does not use the precomputed inverse and therefore does not require that precision. +### Short cyclic remainder reconstruction + +Once the exact quotient is known, the identity `A = Q*B + R` gives, for a supported FFT order `M >= degree(B)`: + + R = (A mod (X^M - 1)) - (Q*B mod (X^M - 1)). + +The high coefficients wrap identically in both terms and cancel. Since `degree(R) < degree(B) <= M`, the difference +is the exact remainder. The implementation folds the dividend and the **complete** divisor, including its leading +coefficient. In particular, that coefficient wraps into cell zero when `M == degree(B)`. + +The mixed-radix backend uses this path for `degree(A) <= 2*degree(B) - 2`, the range from multiplying reduced residues. +It selects the smallest supported order at least `degree(B)` and uses it only when shorter than the ordinary +low-product transform. Outside that range, or when the active backend cannot reuse the prepared order, division +uses its original reconstruction. Quotient computation and inverse-precision validation are unchanged. + +`prepare_cyclic_remainder` and `try_cyclic_remainder` are optional backend/context operations. The latter requires +the exact quotient for the same divisor used in preparation; it returns `false` without changing output when the +degree range or transform is incompatible. Its output may alias the dividend or quotient. The separate preparation +type prevents confusing a cyclic convolution with a low product. + +Generic `multiply_low(A, B, k)` continues to mean `A*B mod X^k`: an ordinary k-point cyclic convolution cannot +replace it. Both reconstruction methods use the existing FFT kernels and preserve canonical results. + +For dense `squaremod` calls using inverse-based division and compatible nonzero preparations, squaring, quotient +calculation and remainder reconstruction each need one forward and one inverse FFT. Cyclic reconstruction shortens +only the remainder transforms; squaring and quotient calculation retain their existing transform lengths. Setup +performs the fixed-operand forward transforms once. These operation counts do not directly predict execution time: +pointwise products, folding, allocation and other division work also contribute. + Two convenience operations share the same dispatch: * `remainder` returns only `dividend mod B`. diff --git a/libs/math/include/nil/crypto3/math/polynomial/arithmetic/polynomial_division.hpp b/libs/math/include/nil/crypto3/math/polynomial/arithmetic/polynomial_division.hpp index 326eee707..6c676bc9a 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/arithmetic/polynomial_division.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/arithmetic/polynomial_division.hpp @@ -36,6 +36,54 @@ namespace nil::crypto3::math { namespace detail { + template + concept SupportsPreparedLowProduct = + requires(polynomial_arithmetic::polynomial_context &context, + typename Backend::polynomial_type &output, const typename Backend::polynomial_type &input, + std::size_t count) { + { context.prepare_low_product(input, count, count).has_value() } -> std::same_as; + { + context.try_multiply_low_prepared(output, input, *context.prepare_low_product(input, count, count), + count) + } -> std::same_as; + }; + + // Keep the optional spectra out of backends that only implement the required multiplication interface. + template + struct divisor_prepared_products { }; + + template + struct divisor_prepared_products { + using preparation_type = + decltype(std::declval &>().prepare_low_product( + std::declval(), std::size_t {}, std::size_t {})); + preparation_type inverse; + preparation_type divisor; + }; + + template + concept SupportsPreparedCyclicRemainder = + requires(polynomial_arithmetic::polynomial_context &context, + typename Backend::polynomial_type &output, const typename Backend::polynomial_type &input, + std::size_t count) { + { context.prepare_cyclic_remainder(input, count).has_value() } -> std::same_as; + { + context.try_cyclic_remainder(output, input, input, *context.prepare_cyclic_remainder(input, count)) + } -> std::same_as; + }; + + template + struct divisor_prepared_cyclic_remainder { }; + + template + struct divisor_prepared_cyclic_remainder { + using preparation_type = + decltype(std::declval &>().prepare_cyclic_remainder( + std::declval(), + std::size_t {})); + preparation_type value; + }; + inline bool use_basecase_division(const polynomial_arithmetic::polynomial_context_options &options, std::size_t divisor_coefficient_count, std::size_t quotient_coefficient_count) { @@ -52,6 +100,14 @@ namespace nil::crypto3::math { * rev(B)^-1 modulo X^inverse_precision. A later division may reuse this inverse when its quotient has at most * inverse_precision coefficients. * + * When supported by the backend, also prepare the two fixed low-product operands for a quotient with + * inverse_precision coefficients. These preparations belong to this divisor snapshot and never grow or change + * during division. Each product falls back to ordinary multiplication when its precision or transform length + * is incompatible with its preparation. + * No preparations are needed when that quotient size selects long division, or when B is constant. + * The modulus preparation uses a shorter cyclic FFT when supported and beneficial; otherwise it uses the + * original low-product FFT. At most one of these modulus spectra is retained, alongside the inverse spectrum. + * * @pre inverse_precision is positive. */ template @@ -82,6 +138,20 @@ namespace nil::crypto3::math { polynomial_type reversed_divisor(divisor_); reverse(reversed_divisor, reversed_divisor.size()); inverse_series(reversed_divisor_inverse_, reversed_divisor, inverse_precision_, arithmetic_context); + if constexpr (detail::SupportsPreparedLowProduct) { + if (degree() != 0 && + !detail::use_basecase_division(arithmetic_context.options(), divisor_.size(), inverse_precision_)) { + prepared_.inverse = arithmetic_context.prepare_low_product(reversed_divisor_inverse_, + inverse_precision_, inverse_precision_); + if constexpr (detail::SupportsPreparedCyclicRemainder) { + cyclic_.value = arithmetic_context.prepare_cyclic_remainder(divisor_, inverse_precision_); + if (cyclic_.value) { + return; + } + } + prepared_.divisor = arithmetic_context.prepare_low_product(divisor_, inverse_precision_, degree()); + } + } } const polynomial_type &divisor() const { @@ -100,10 +170,52 @@ namespace nil::crypto3::math { return inverse_precision_; } + // Multiply by the stored reversed inverse, retaining exactly coefficient_count low coefficients. The + // preparation must match this precision and the selected FFT length; otherwise use the ordinary product. + void multiply_low_reversed_inverse(polynomial_type &output, const polynomial_type &input, + std::size_t coefficient_count, + polynomial_arithmetic::polynomial_context &context) const { + if constexpr (detail::SupportsPreparedLowProduct) { + if (prepared_.inverse && + context.try_multiply_low_prepared(output, input, *prepared_.inverse, coefficient_count)) { + return; + } + } + context.multiply_low(output, input, reversed_divisor_inverse_, coefficient_count); + } + + // Only degree(B) low coefficients of Q*B are needed for remainder reconstruction. The leading coefficient + // of B cannot affect this prefix. Keep the same low-product operation whether or not preparation succeeds. + void multiply_low_divisor(polynomial_type &output, const polynomial_type &input, + polynomial_arithmetic::polynomial_context &context) const { + if constexpr (detail::SupportsPreparedLowProduct) { + if (prepared_.divisor && + context.try_multiply_low_prepared(output, input, *prepared_.divisor, degree())) { + return; + } + } + context.multiply_low(output, input, divisor_, degree()); + } + + // With exact quotient Q, cyclic reconstruction cancels the wrapped high coefficients of V and Q*B. + // The backend checks its bounded degree range and current transform lengths before changing output. + bool try_reconstruct_remainder(polynomial_type &output, const polynomial_type ÷nd, + const polynomial_type "ient, + polynomial_arithmetic::polynomial_context &context) const { + if constexpr (detail::SupportsPreparedCyclicRemainder) { + if (cyclic_.value) { + return context.try_cyclic_remainder(output, dividend, quotient, *cyclic_.value); + } + } + return false; + } + private: polynomial_type divisor_; polynomial_type reversed_divisor_inverse_; std::size_t inverse_precision_; + [[no_unique_address]] detail::divisor_prepared_products prepared_; + [[no_unique_address]] detail::divisor_prepared_cyclic_remainder cyclic_; }; /** @@ -116,6 +228,8 @@ namespace nil::crypto3::math { * * After recovering Q, only the first d coefficients of A - Q * B are needed because the remainder has degree less * than d. If n < d, the quotient is zero and the dividend is returned unchanged as the remainder. + * For n <= 2*d - 2, a backend may reconstruct R using a shorter cyclic convolution instead: folding A and + * Q*B modulo X^M - 1 with M >= d cancels their wrapped high coefficients and leaves exactly R. * * For small divisors or quotients, the arithmetic context selects quadratic long division instead. Its two * inclusive coefficient-count cutoffs are independently configurable, and setting either cutoff to zero disables @@ -182,8 +296,8 @@ namespace nil::crypto3::math { } polynomial_type reversed_quotient; - arithmetic_context.multiply_low(reversed_quotient, reversed_dividend, - divisor_context.reversed_divisor_inverse(), quotient_size); + divisor_context.multiply_low_reversed_inverse(reversed_quotient, reversed_dividend, quotient_size, + arithmetic_context); reversed_quotient.resize(quotient_size, value_type {}); reverse(reversed_quotient, quotient_size); condense(reversed_quotient); @@ -193,9 +307,9 @@ namespace nil::crypto3::math { if (divisor_degree == 0) { remainder_result.resize(1); remainder_result[0] = value_type {}; - } else { - arithmetic_context.multiply_low(remainder_result, quotient_result, divisor_context.divisor(), - divisor_degree); + } else if (!divisor_context.try_reconstruct_remainder(remainder_result, dividend, quotient_result, + arithmetic_context)) { + divisor_context.multiply_low_divisor(remainder_result, quotient_result, arithmetic_context); // multiply_low returns canonical output; restore the complete prefix before coefficient-wise subtraction. remainder_result.resize(divisor_degree, value_type {}); for (std::size_t i = 0; i < divisor_degree; ++i) { diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/mixed_radix_backend.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/mixed_radix_backend.hpp index b70b5041d..9631e68e1 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/backends/mixed_radix_backend.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/mixed_radix_backend.hpp @@ -27,6 +27,8 @@ #include #include +#include +#include #include #include #include @@ -44,8 +46,8 @@ namespace nil::crypto3::math::polynomial_arithmetic { * contains each product. It reuses its workspace across operations and is therefore intended * for sequential use. Concurrent operations require separate backend instances. * - * The configured transform must contain the complete product of the coefficients - * used by an operation. multiply_low discards input coefficients that cannot + * For multiply, square and multiply_low, the configured transform must contain the complete product of the + * coefficients used by the operation. multiply_low discards input coefficients that cannot * affect the requested low coefficients before computing that product. */ template @@ -54,6 +56,68 @@ namespace nil::crypto3::math::polynomial_arithmetic { using value_type = ValueType; using polynomial_type = math::polynomial; + /** + * Owned, immutable preparation of a fixed operand for multiply_low. The spectrum contains the FFT of + * the first min(fixed.size(), coefficient_count) coefficients, with enough zero padding for a linear + * product. It does not refer to the source polynomial, backend, FFT plan, or scratch storage. + */ + class prepared_low_product_type { + public: + std::size_t coefficient_count() const { + return coefficient_count_; + } + + // Zero means the requested fixed prefix is zero, so no transform was needed. + std::size_t transform_size() const { + return spectrum_.size(); + } + + std::size_t storage_bytes() const { + return spectrum_.capacity() * sizeof(value_type); + } + + private: + friend class mixed_radix_backend; + + prepared_low_product_type(std::size_t coefficient_count, std::size_t prefix_size, bool zero, + std::vector spectrum = {}) : + coefficient_count_(coefficient_count), prefix_size_(prefix_size), zero_(zero), + spectrum_(std::move(spectrum)) { + } + + std::size_t coefficient_count_; + std::size_t prefix_size_; + bool zero_; + std::vector spectrum_; + }; + + // Owns the FFT of the complete divisor modulo X^M - 1, including its leading coefficient. + // This preparation is separate from multiply_low, whose modulus is a power of X. + class prepared_cyclic_remainder_type { + public: + std::size_t divisor_degree() const { + return divisor_degree_; + } + + std::size_t transform_size() const { + return spectrum_.size(); + } + + std::size_t storage_bytes() const { + return spectrum_.capacity() * sizeof(value_type); + } + + private: + friend class mixed_radix_backend; + + prepared_cyclic_remainder_type(std::size_t divisor_degree, std::vector spectrum) : + divisor_degree_(divisor_degree), spectrum_(std::move(spectrum)) { + } + + std::size_t divisor_degree_; + std::vector spectrum_; + }; + explicit mixed_radix_backend(std::size_t transform_size) { // Every divisor of a valid transform order also has a root of unity. Caching those plans lets algorithms // whose operands grow in stages use a proportionally sized transform at each stage instead of repeatedly @@ -113,7 +177,169 @@ namespace nil::crypto3::math::polynomial_arithmetic { multiply_prefixes(output, left, left_size, right, right_size, result_size); } + /** + * Prepare a fixed right operand for a low product with a variable operand of variable_coefficient_count + * coefficients. This uses the same prefixes and smallest fitting transform as multiply_low. + * Returns nullopt if a positive precision has no variable coefficients or no plan can hold the complete + * prefix product. Zero precision and a zero fixed prefix need no FFT and produce a zero preparation. + * + * @pre fixed is a nonempty canonical coefficient polynomial. + */ + std::optional prepare_low_product(const polynomial_type &fixed, + std::size_t variable_coefficient_count, + std::size_t coefficient_count) { + const std::size_t prefix_size = std::min(fixed.size(), coefficient_count); + if (coefficient_count == 0) { + return prepared_low_product_type(coefficient_count, prefix_size, true); + } + const std::size_t variable_size = std::min(variable_coefficient_count, coefficient_count); + if (variable_size == 0) { + return std::nullopt; + } + if (math::is_zero(fixed.begin(), fixed.begin() + prefix_size)) { + return prepared_low_product_type(coefficient_count, prefix_size, true); + } + if (variable_size > std::numeric_limits::max() - (prefix_size - 1)) { + return std::nullopt; + } + const auto *plan = find_plan(variable_size + prefix_size - 1); + if (!plan) { + return std::nullopt; + } + + std::vector spectrum(fixed.begin(), fixed.begin() + prefix_size); + plan->fft(spectrum, workspace_); + compact_spectrum(spectrum); + return prepared_low_product_type(coefficient_count, prefix_size, false, std::move(spectrum)); + } + + /** + * Compute left * fixed mod X^coefficient_count using an existing preparation. Returns false, leaving + * output unchanged, if the precision differs or the current backend would select a different FFT length. + * Different variable-operand lengths and backend instances are compatible when that selected length + * agrees. Zero products need no matching FFT plan. Successful output is canonical and may alias left; + * the preparation is never modified. + * + * @pre left is a nonempty coefficient polynomial. Trailing zero coefficients are permitted, as in the + * fixed-length reversal of a dividend prefix during division. + */ + bool try_multiply_low_prepared(polynomial_type &output, const polynomial_type &left, + const prepared_low_product_type &prepared, std::size_t coefficient_count) { + if (coefficient_count != prepared.coefficient_count_) { + return false; + } + if (prepared.zero_) { + set_zero(output); + return true; + } + // A moved-from nonzero preparation no longer owns a spectrum. + if (prepared.spectrum_.empty()) { + return false; + } + const std::size_t left_size = std::min(left.size(), coefficient_count); + if (math::is_zero(left.begin(), left.begin() + left_size)) { + set_zero(output); + return true; + } + if (left_size > std::numeric_limits::max() - (prepared.prefix_size_ - 1)) { + return false; + } + const std::size_t product_size = left_size + prepared.prefix_size_ - 1; + const auto *plan = find_plan(product_size); + if (!plan || plan->size() != prepared.transform_size()) { + return false; + } + + polynomial_type transformed(left.begin(), left.begin() + left_size); + plan->fft(transformed.get_storage(), workspace_); + finish_product(output, transformed, prepared.spectrum_, std::min(coefficient_count, product_size), *plan); + return true; + } + + /** + * Prepare remainder reconstruction for products of reduced residues: degree(V) <= 2*n - 2, where + * n = degree(divisor). Select the smallest supported M >= n, and return nullopt unless it is shorter + * than the existing low-product transform for the expected quotient. Quotient lengths above n - 1 + * are capped during preparation; division outside this range uses its existing reconstruction. + * The complete divisor is folded modulo X^M - 1, so its leading coefficient wraps when M == n. + * + * @pre divisor is a nonempty canonical coefficient polynomial. + */ + std::optional prepare_cyclic_remainder(const polynomial_type &divisor, + std::size_t quotient_coefficient_count) { + const std::size_t n = divisor.size() - 1; + if (n < 2 || math::is_zero(divisor.begin(), divisor.begin() + n)) { + // A constant/linear divisor has no nontrivial reduced product to accelerate. For a monomial + // divisor the existing low-product remainder is zero and needs no transforms at all. + return std::nullopt; + } + const auto *plan = cyclic_remainder_plan(n, std::min(quotient_coefficient_count, n - 1)); + if (!plan) { + return std::nullopt; + } + std::vector spectrum; + fold_coefficients(spectrum, divisor, plan->size()); + plan->fft(spectrum, workspace_); + compact_spectrum(spectrum); + return prepared_cyclic_remainder_type(n, std::move(spectrum)); + } + + /** + * Reconstruct R from V = Q*H + R using a prepared H. For n = degree(H) and its prepared length M >= n, + * R = fold_M(V) - (Q*H mod (X^M - 1)): the wrapped high coefficients cancel because degree(R) < n. + * Returns false without changing output outside degree(V) <= 2*n - 2, if a compatible shorter plan is + * unavailable, or if the preparation has been moved from. Successful output is canonical and may alias + * either input. The preparation remains unchanged. + * + * @pre dividend and quotient are nonempty canonical polynomials, and quotient is the exact quotient + * of dividend by the divisor used to construct prepared. + */ + bool try_cyclic_remainder(polynomial_type &output, const polynomial_type ÷nd, + const polynomial_type "ient, const prepared_cyclic_remainder_type &prepared) { + const std::size_t n = prepared.divisor_degree_; + if (prepared.spectrum_.empty() || dividend.size() <= n || dividend.size() - n > n - 1 || + quotient.size() != dividend.size() - n) { + return false; + } + const auto *plan = cyclic_remainder_plan(n, quotient.size()); + if (!plan || plan->size() != prepared.transform_size()) { + return false; + } + + polynomial_type transformed, product, result; + fold_coefficients(transformed.get_storage(), quotient, plan->size()); + plan->fft(transformed.get_storage(), workspace_); + finish_product(product, transformed, prepared.spectrum_, plan->size(), *plan); + + // dividend is V and product holds cyclic Q*H. Folding V is essential: its wrapped coefficients + // cancel those of Q*H. Keep all M cells of the difference, including the tail from n through M - 1. + fold_coefficients(result.get_storage(), dividend, plan->size()); + for (std::size_t i = 0; i < product.size(); ++i) { + result[i] -= product[i]; + } + // Condense rather than truncate to n: an incorrect uncancelled tail must remain observable in tests. + math::condense(result); + output = std::move(result); + return true; + } + private: + static void compact_spectrum(std::vector &spectrum) { + // FFT exchanges buffers with workspace_. A preceding larger transform can therefore leave this + // smaller spectrum with a large allocation. Detach that excess capacity once during preparation. + if (spectrum.capacity() > spectrum.size()) { + std::vector(spectrum.begin(), spectrum.end()).swap(spectrum); + } + } + + static void fold_coefficients(std::vector &output, const polynomial_type &input, + std::size_t length) { + output.assign(length, value_type::zero()); + for (std::size_t i = 0; i < input.size(); ++i) { + output[i % length] += input[i]; + } + } + static std::vector divisor_transform_sizes(std::size_t transform_size) { if (transform_size == 0) { throw std::invalid_argument("mixed_radix_backend: expected transform size > 0"); @@ -148,17 +374,34 @@ namespace nil::crypto3::math::polynomial_arithmetic { output.assign(1, value_type::zero()); } - const mixed_radix_fft_plan &plan_for(std::size_t result_size) const { + const mixed_radix_fft_plan *find_plan(std::size_t result_size) const { // Plans are sorted by order, so the first plan large enough for the product introduces the least padding. const auto plan = std::lower_bound(plans_.begin(), plans_.end(), result_size, [](const mixed_radix_fft_plan &candidate, std::size_t size) { return candidate.size() < size; }); - if (plan == plans_.end()) { + return plan == plans_.end() ? nullptr : &*plan; + } + + const mixed_radix_fft_plan &plan_for(std::size_t result_size) const { + const auto *plan = find_plan(result_size); + if (!plan) { throw std::invalid_argument("mixed_radix_backend: transform is too small for the product"); } return *plan; } + const mixed_radix_fft_plan *cyclic_remainder_plan(std::size_t n, + std::size_t quotient_size) const { + if (n < 2 || quotient_size == 0 || quotient_size > n - 1 || + quotient_size > std::numeric_limits::max() - (n - 1)) { + return nullptr; + } + const auto *short_plan = find_plan(n); + // Within this range Q has quotient_size coefficients and the existing low product uses n from H. + const auto *linear_plan = find_plan(quotient_size + n - 1); + return short_plan && linear_plan && short_plan->size() < linear_plan->size() ? short_plan : nullptr; + } + void multiply_prefixes(polynomial_type &output, const polynomial_type &left, std::size_t left_size, const polynomial_type &right, std::size_t right_size, std::size_t result_size) { const mixed_radix_fft_plan &plan = plan_for(left_size + right_size - 1); @@ -168,6 +411,12 @@ namespace nil::crypto3::math::polynomial_arithmetic { plan.fft(transformed_left.get_storage(), workspace_); plan.fft(transformed_right.get_storage(), workspace_); + finish_product(output, transformed_left, transformed_right.get_storage(), result_size, plan); + } + + void finish_product(polynomial_type &output, polynomial_type &transformed_left, + const std::vector &transformed_right, std::size_t result_size, + const mixed_radix_fft_plan &plan) { for (std::size_t i = 0; i < plan.size(); ++i) { transformed_left[i] = transformed_left[i] * transformed_right[i]; } diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp index 3551f2845..0d9fc2b9b 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/polynomial_backend.hpp @@ -119,6 +119,47 @@ namespace nil::crypto3::math::polynomial_arithmetic { backend_.multiply_low(output, left, right, coefficient_count); } + // Optional capabilities: backends without prepared low products still satisfy PolynomialBackend. + // Forward the backend's owned preparation without adding another context or exposing its scratch storage. + auto prepare_low_product(const polynomial_type &fixed, std::size_t variable_coefficient_count, + std::size_t coefficient_count) + requires requires(backend_type &backend, const polynomial_type &operand, std::size_t count) { + backend.prepare_low_product(operand, count, count); + } + { + return backend_.prepare_low_product(fixed, variable_coefficient_count, coefficient_count); + } + + template + bool try_multiply_low_prepared(polynomial_type &output, const polynomial_type &left, + const PreparedOperand &prepared, std::size_t coefficient_count) + requires requires(backend_type &backend, polynomial_type &result, const polynomial_type &operand, + const PreparedOperand &fixed, std::size_t count) { + { backend.try_multiply_low_prepared(result, operand, fixed, count) } -> std::same_as; + } + { + return backend_.try_multiply_low_prepared(output, left, prepared, coefficient_count); + } + + auto prepare_cyclic_remainder(const polynomial_type &divisor, std::size_t quotient_coefficient_count) + requires requires(backend_type &backend, const polynomial_type &operand, std::size_t count) { + backend.prepare_cyclic_remainder(operand, count); + } + { + return backend_.prepare_cyclic_remainder(divisor, quotient_coefficient_count); + } + + template + bool try_cyclic_remainder(polynomial_type &output, const polynomial_type ÷nd, + const polynomial_type "ient, const PreparedRemainder &prepared) + requires requires(backend_type &backend, polynomial_type &result, const polynomial_type &operand, + const PreparedRemainder &fixed) { + { backend.try_cyclic_remainder(result, operand, operand, fixed) } -> std::same_as; + } + { + return backend_.try_cyclic_remainder(output, dividend, quotient, prepared); + } + const options_type &options() const { return options_; } diff --git a/libs/math/test/polynomial_backend.cpp b/libs/math/test/polynomial_backend.cpp index b1de72a77..9bd9515cb 100644 --- a/libs/math/test/polynomial_backend.cpp +++ b/libs/math/test/polynomial_backend.cpp @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -64,6 +65,47 @@ namespace { static_assert(polynomial_arithmetic::PolynomialBackend); static_assert(polynomial_arithmetic::PolynomialBackend); + template + concept HasPreparedLowProduct = requires(Context &context, const typename Context::polynomial_type &fixed) { + context.prepare_low_product(fixed, std::size_t(2), std::size_t(3)); + }; + + template + concept UsesPreparedLowProduct = + requires(Context &context, typename Context::polynomial_type &output, + const typename Context::polynomial_type &input, const Prepared &prepared) { + { context.try_multiply_low_prepared(output, input, prepared, std::size_t(3)) } -> std::same_as; + }; + + using fq_context = polynomial_arithmetic::polynomial_context; + using schoolbook_context = + polynomial_arithmetic::polynomial_context>; + static_assert(HasPreparedLowProduct); + static_assert(!HasPreparedLowProduct); + static_assert(UsesPreparedLowProduct); + static_assert(!UsesPreparedLowProduct); + static_assert(!UsesPreparedLowProduct); + + template + concept HasPreparedCyclicRemainder = requires(Context &context, const typename Context::polynomial_type &divisor) { + context.prepare_cyclic_remainder(divisor, std::size_t(5)); + }; + + template + concept UsesPreparedCyclicRemainder = + requires(Context &context, typename Context::polynomial_type &output, + const typename Context::polynomial_type &input, const Prepared &prepared) { + { context.try_cyclic_remainder(output, input, input, prepared) } -> std::same_as; + }; + + static_assert(HasPreparedCyclicRemainder); + static_assert(!HasPreparedCyclicRemainder); + static_assert(UsesPreparedCyclicRemainder); + static_assert( + !UsesPreparedCyclicRemainder); + static_assert(!UsesPreparedCyclicRemainder); + static_assert(!UsesPreparedLowProduct); + template math::polynomial expected_low_product(const math::polynomial &product, std::size_t coefficient_count) { @@ -142,6 +184,69 @@ namespace { return value; } + template + void check_prepared_low_products(Backend backend, const typename Backend::polynomial_type &left, + const typename Backend::polynomial_type &fixed) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename Backend::value_type; + polynomial_arithmetic::polynomial_context context(std::move(backend)); + polynomial_arithmetic::schoolbook_backend reference; + for (std::size_t count = 0; count <= left.size() + fixed.size() + 1; ++count) { + const auto prepared = context.prepare_low_product(fixed, left.size(), count); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->coefficient_count(), count); + polynomial_type variable = left; + for (int repeat = 0; repeat < 3; ++repeat) { + polynomial_type full_product, ordinary, result; + reference.multiply(full_product, variable, fixed); + const auto expected = expected_low_product(full_product, count); + context.multiply_low(ordinary, variable, fixed, count); + BOOST_CHECK(ordinary == expected); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, variable, *prepared, count)); + BOOST_CHECK(result == expected); + result = variable; + BOOST_REQUIRE(context.try_multiply_low_prepared(result, result, *prepared, count)); + BOOST_CHECK(result == expected); + // Different values with the same prefix length reuse the same immutable spectrum. + variable[0] += value_type::one(); + } + } + } + + template + void check_cyclic_reconstruction(Backend backend, const typename Backend::polynomial_type &divisor, + typename Backend::polynomial_type quotient, + const typename Backend::polynomial_type &remainder, std::size_t cyclic_length) { + using polynomial_type = typename Backend::polynomial_type; + using value_type = typename Backend::value_type; + polynomial_arithmetic::polynomial_context context(std::move(backend)); + polynomial_arithmetic::schoolbook_backend reference; + const auto prepared = context.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->divisor_degree(), divisor.size() - 1); + BOOST_CHECK_EQUAL(prepared->transform_size(), cyclic_length); + for (std::size_t repeat = 0; repeat < 3; ++repeat) { + polynomial_type dividend, expected_quotient, expected_remainder, result; + reference.multiply(dividend, quotient, divisor); + math::addition(dividend, dividend, remainder); + math::division(expected_quotient, expected_remainder, dividend, divisor); + BOOST_CHECK(expected_quotient == quotient); + BOOST_CHECK(expected_remainder == remainder); + BOOST_REQUIRE(context.try_cyclic_remainder(result, dividend, quotient, *prepared)); + BOOST_CHECK(result == expected_remainder); + // Reconstruction only condenses the M-cell difference; it never truncates to degree(H). Therefore + // this bound checks that every coefficient from degree(H) through M - 1 actually cancelled. + BOOST_CHECK_LT(result.size(), divisor.size()); + result = dividend; + BOOST_REQUIRE(context.try_cyclic_remainder(result, result, quotient, *prepared)); + BOOST_CHECK(result == expected_remainder); + result = quotient; + BOOST_REQUIRE(context.try_cyclic_remainder(result, dividend, result, *prepared)); + BOOST_CHECK(result == expected_remainder); + quotient[0] += value_type::one(); + } + } + } // namespace BOOST_AUTO_TEST_SUITE(polynomial_backend_test_suite) @@ -184,6 +289,218 @@ BOOST_AUTO_TEST_CASE(fq12_backends_conform) { } } +BOOST_AUTO_TEST_CASE(prepared_low_products_match_schoolbook_and_ordinary_products) { + check_prepared_low_products(fq_mixed_radix_backend(18), fq_mixed_radix_backend::polynomial_type {1, 2, 0, 3}, + fq_mixed_radix_backend::polynomial_type {0, 4, 5}); + check_prepared_low_products( + fq12_mixed_radix_backend(18), + fq12_mixed_radix_backend::polynomial_type {fq12_value(1), fq12_value_type::zero(), fq12_value(13)}, + fq12_mixed_radix_backend::polynomial_type {fq12_value_type::zero(), fq12_value(25), fq12_value(37)}); +} + +BOOST_AUTO_TEST_CASE(prepared_low_products_handle_zero_and_truncated_prefixes) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + fq_context context {fq_mixed_radix_backend(3)}; + const polynomial_type input {1, 2, 3}; + const polynomial_type zero {fq_value_type::zero()}; + const polynomial_type zero_prefix {0, 0, 7}; + polynomial_type result {fq_value_type(9)}; + + // The full product needs five coefficients, but the selected prefixes need only a size-three FFT. + auto prepared = context.prepare_low_product(input, input.size(), 2); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 3); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *prepared, 2)); + BOOST_CHECK(result == polynomial_type({1, 4})); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, zero, *prepared, 2)); + BOOST_CHECK(result == zero); + + for (const auto &fixed : {zero, zero_prefix}) { + const auto zero_prepared = context.prepare_low_product(fixed, input.size(), 2); + BOOST_REQUIRE(zero_prepared.has_value()); + BOOST_CHECK_EQUAL(zero_prepared->transform_size(), 0); + BOOST_CHECK_EQUAL(zero_prepared->storage_bytes(), 0); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *zero_prepared, 2)); + BOOST_CHECK(result == zero); + } + prepared = context.prepare_low_product(input, input.size(), 0); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 0); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *prepared, 0)); + BOOST_CHECK(result == zero); + + // A nonzero prefix may end in zero coefficients; final output must still be canonical. + const auto trailing_zero = context.prepare_low_product(polynomial_type {2, 0, 7}, 2, 2); + BOOST_REQUIRE(trailing_zero.has_value()); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, polynomial_type {3, 0, 1}, *trailing_zero, 2)); + BOOST_CHECK(result == polynomial_type({fq_value_type(6)})); +} + +BOOST_AUTO_TEST_CASE(prepared_low_products_validate_precision_and_transform_shape) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + fq_context context {fq_mixed_radix_backend(18)}; + const polynomial_type fixed {2, 0, 1}; + const auto prepared = context.prepare_low_product(fixed, 2, 10); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 6); + + // Several variable-prefix lengths select the same transform and may share this spectrum. + polynomial_arithmetic::schoolbook_backend reference; + for (const auto &left : {polynomial_type {1, 3}, polynomial_type {1, 3, 4}, polynomial_type {1, 3, 4, 5}}) { + polynomial_type result, expected; + reference.multiply(expected, left, fixed); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, left, *prepared, 10)); + BOOST_CHECK(result == expected); + } + + // A different precision cannot reuse a spectrum containing a different fixed prefix. A different + // selected length also cannot reuse it, even if the configured maximum can hold the product. + polynomial_type result {fq_value_type(11)}; + BOOST_CHECK(!context.try_multiply_low_prepared(result, polynomial_type {1, 3}, *prepared, 9)); + BOOST_CHECK(result == polynomial_type({fq_value_type(11)})); + for (const auto &input : {polynomial_type {fq_value_type::one()}, polynomial_type {1, 2, 3, 4, 5}}) { + result = input; + BOOST_CHECK(!context.try_multiply_low_prepared(result, result, *prepared, 10)); + BOOST_CHECK(result == input); + } + fq_context different_plan {fq_mixed_radix_backend(9)}; + BOOST_CHECK(!different_plan.try_multiply_low_prepared(result, polynomial_type {1, 3}, *prepared, 10)); + fq_context too_small {fq_mixed_radix_backend(3)}; + BOOST_CHECK(!too_small.try_multiply_low_prepared(result, polynomial_type {1, 3}, *prepared, 10)); + BOOST_CHECK(!too_small.prepare_low_product(fixed, 3, 3).has_value()); + BOOST_CHECK(!context.prepare_low_product(fixed, 0, 3).has_value()); + BOOST_CHECK(!context + .prepare_low_product(fixed, std::numeric_limits::max(), + std::numeric_limits::max()) + .has_value()); +} + +BOOST_AUTO_TEST_CASE(prepared_low_products_own_the_fixed_operand_and_survive_context_changes) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + const polynomial_type original {2, 0, 1}; + const polynomial_type input {1, 3}; + const auto prepared = [&] { + fq_context temporary {fq_mixed_radix_backend(6)}; + polynomial_type source = original; + auto result = temporary.prepare_low_product(source, input.size(), 10); + source[0] = fq_value_type(99); + return result; + }(); + BOOST_REQUIRE(prepared.has_value()); + auto copy = *prepared; + const auto moved = std::move(copy); + fq_context context {fq_mixed_radix_backend(18)}; + for (const auto *fixed : {&*prepared, &moved}) { + polynomial_type result, expected; + context.multiply_low(expected, input, original, 10); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *fixed, 10)); + BOOST_CHECK(result == expected); + } + const auto changed = context.prepare_low_product(polynomial_type {99, 0, 1}, input.size(), 10); + BOOST_REQUIRE(changed.has_value()); + polynomial_type result; + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *changed, 10)); + BOOST_CHECK(result == polynomial_type({99, 297, 1, 3})); + BOOST_REQUIRE(context.try_multiply_low_prepared(result, input, *prepared, 10)); + BOOST_CHECK(result == polynomial_type({2, 6, 1, 3})); +} + +BOOST_AUTO_TEST_CASE(prepared_low_products_do_not_retain_larger_workspace_allocations) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + fq_mixed_radix_backend backend(134); + polynomial_type large(60, fq_value_type::one()), result; + backend.square(result, large); // Populate the size-134 workspace first. + const polynomial_type small(16, fq_value_type::one()); + const auto prepared = backend.prepare_low_product(small, small.size(), small.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 67); + // The long-lived spectrum owns a fresh size-67 allocation instead of the size-134 scratch buffer. + BOOST_CHECK_LT(prepared->storage_bytes(), 134 * sizeof(fq_value_type)); + BOOST_REQUIRE(backend.try_multiply_low_prepared(result, small, *prepared, small.size())); + BOOST_REQUIRE_EQUAL(result.size(), small.size()); + for (std::size_t i = 0; i < small.size(); ++i) { + BOOST_CHECK(result[i] == fq_value_type(i + 1)); + } +} + +BOOST_AUTO_TEST_CASE(cyclic_reconstruction_folds_the_dividend_and_complete_divisor) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + // M == degree(H): H's nonunit leading coefficient must wrap into cell zero. + check_cyclic_reconstruction(fq_mixed_radix_backend(18), polynomial_type {2, 3, 5, 7, 11, 13, 17}, + polynomial_type {19, 23, 29, 31, 37}, polynomial_type {41, 43}, 6); + // M > degree(H): the extra coefficient in the folded difference must cancel to zero. + check_cyclic_reconstruction(fq_mixed_radix_backend(18), polynomial_type {2, 3, 5, 7, 11, 1}, + polynomial_type {13, 17, 19, 23}, polynomial_type {fq_value_type::zero()}, 6); + // Folding H can itself give zero, as for H = X^6 - 1 and M = 6. + check_cyclic_reconstruction(fq_mixed_radix_backend(18), polynomial_type {-fq_value_type::one(), 0, 0, 0, 0, 0, 1}, + polynomial_type {0, 0, 0, 0, 1}, polynomial_type {0, 0, 1}, 6); + + using extension_polynomial = fq12_mixed_radix_backend::polynomial_type; + check_cyclic_reconstruction(fq12_mixed_radix_backend(18), + extension_polynomial {fq12_value(1), fq12_value(13), fq12_value(25), fq12_value(37), + fq12_value(49), fq12_value(61)}, + extension_polynomial {fq12_value(73), fq12_value(85), fq12_value(97), fq12_value(109)}, + extension_polynomial {fq12_value(121), fq12_value(133)}, 6); +} + +BOOST_AUTO_TEST_CASE(cyclic_preparation_and_reuse_require_a_supported_shorter_transform_and_bounded_degree) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + fq_context context {fq_mixed_radix_backend(18)}; + const polynomial_type divisor {2, 3, 5, 7, 11, 13, 17}; + const auto prepared = context.prepare_cyclic_remainder(divisor, 5); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(divisor, 0).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(divisor, 1).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(polynomial_type {fq_value_type(7)}, 5).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(polynomial_type {7, 11}, 5).has_value()); + BOOST_CHECK(!context.prepare_cyclic_remainder(polynomial_type {0, 0, 0, 0, 0, 0, 1}, 5).has_value()); + fq_context same_length {fq_mixed_radix_backend(67)}; + BOOST_CHECK(!same_length.prepare_cyclic_remainder(divisor, 5).has_value()); + fq_context too_small {fq_mixed_radix_backend(3)}; + BOOST_CHECK(!too_small.prepare_cyclic_remainder(divisor, 5).has_value()); + + polynomial_arithmetic::schoolbook_backend reference; + for (std::size_t quotient_size : {1, 5, 6}) { + polynomial_type quotient(quotient_size, fq_value_type::one()), dividend, output {fq_value_type(99)}; + reference.multiply(dividend, quotient, divisor); + if (quotient_size == 5) { + // The current backend lacks the prepared order-six plan; the preparation remains valid elsewhere. + BOOST_CHECK(!same_length.try_cyclic_remainder(output, dividend, quotient, *prepared)); + BOOST_CHECK(output == polynomial_type({fq_value_type(99)})); + BOOST_CHECK(!too_small.try_cyclic_remainder(output, dividend, quotient, *prepared)); + BOOST_REQUIRE(context.try_cyclic_remainder(output, dividend, quotient, *prepared)); + BOOST_CHECK(math::is_zero(output)); + } else { + // A constant quotient has no shorter remainder transform; size six exceeds degree(V) <= 2*n - 2. + BOOST_CHECK(!context.try_cyclic_remainder(output, dividend, quotient, *prepared)); + BOOST_CHECK(output == polynomial_type({fq_value_type(99)})); + output = dividend; + BOOST_CHECK(!context.try_cyclic_remainder(output, output, quotient, *prepared)); + BOOST_CHECK(output == dividend); + } + } + polynomial_type small {0, 0, 1}, zero {fq_value_type::zero()}, output {fq_value_type(99)}; + BOOST_CHECK(!context.try_cyclic_remainder(output, small, zero, *prepared)); + BOOST_CHECK(output == polynomial_type({fq_value_type(99)})); +} + +BOOST_AUTO_TEST_CASE(short_cyclic_preparations_release_larger_workspace_allocations) { + using polynomial_type = fq_mixed_radix_backend::polynomial_type; + fq_mixed_radix_backend backend(134); + polynomial_type large(60, fq_value_type::one()), result; + backend.square(result, large); + const polynomial_type divisor(65, fq_value_type::one()), quotient(63, fq_value_type::one()); + const auto prepared = backend.prepare_cyclic_remainder(divisor, quotient.size()); + BOOST_REQUIRE(prepared.has_value()); + BOOST_CHECK_EQUAL(prepared->transform_size(), 67); + BOOST_CHECK_LT(prepared->storage_bytes(), 134 * sizeof(fq_value_type)); + polynomial_type dividend; + polynomial_arithmetic::schoolbook_backend {}.multiply(dividend, quotient, divisor); + dividend[0] += fq_value_type(19); + BOOST_REQUIRE(backend.try_cyclic_remainder(result, dividend, quotient, *prepared)); + BOOST_CHECK(result == polynomial_type({fq_value_type(19)})); +} + BOOST_AUTO_TEST_CASE(scalar_multiplication_supports_base_field_scalars) { using polynomial_type = math::polynomial; diff --git a/libs/math/test/polynomial_division.cpp b/libs/math/test/polynomial_division.cpp index 709d3c1a4..9154ff129 100644 --- a/libs/math/test/polynomial_division.cpp +++ b/libs/math/test/polynomial_division.cpp @@ -98,6 +98,124 @@ namespace { return value; } + struct prepared_product_counts { + std::size_t preparations = 0; + std::size_t prepared_products = 0; + std::size_t ordinary_products = 0; + std::size_t squares = 0; + std::size_t last_low_transform_size = 0; + std::size_t cyclic_preparations = 0; + std::size_t cyclic_remainders = 0; + std::size_t cyclic_transform_size = 0; + }; + + // Count calls at the backend boundary, as in the exponentiation tests; all arithmetic uses the real backend. + template + struct counting_mixed_radix_backend : polynomial_arithmetic::mixed_radix_backend { + using base_type = polynomial_arithmetic::mixed_radix_backend; + using polynomial_type = typename base_type::polynomial_type; + using prepared_type = typename base_type::prepared_low_product_type; + using cyclic_type = typename base_type::prepared_cyclic_remainder_type; + + counting_mixed_radix_backend(std::size_t order, prepared_product_counts &counts) : + base_type(order), counts(counts) { + } + + auto prepare_low_product(const polynomial_type &fixed, std::size_t variable_count, std::size_t count) { + ++counts.preparations; + auto prepared = base_type::prepare_low_product(fixed, variable_count, count); + counts.last_low_transform_size = prepared ? prepared->transform_size() : 0; + return prepared; + } + + bool try_multiply_low_prepared(polynomial_type &output, const polynomial_type &input, + const prepared_type &prepared, std::size_t count) { + const bool used = base_type::try_multiply_low_prepared(output, input, prepared, count); + counts.prepared_products += used; + return used; + } + + void multiply_low(polynomial_type &output, const polynomial_type &left, const polynomial_type &right, + std::size_t count) { + ++counts.ordinary_products; + base_type::multiply_low(output, left, right, count); + } + + void square(polynomial_type &output, const polynomial_type &input) { + ++counts.squares; + base_type::square(output, input); + } + + auto prepare_cyclic_remainder(const polynomial_type &divisor, std::size_t quotient_count) { + ++counts.cyclic_preparations; + auto prepared = base_type::prepare_cyclic_remainder(divisor, quotient_count); + counts.cyclic_transform_size = prepared ? prepared->transform_size() : 0; + return prepared; + } + + bool try_cyclic_remainder(polynomial_type &output, const polynomial_type ÷nd, + const polynomial_type "ient, const cyclic_type &prepared) { + const bool used = base_type::try_cyclic_remainder(output, dividend, quotient, prepared); + counts.cyclic_remainders += used; + return used; + } + + prepared_product_counts &counts; + }; + + template + void check_reused_division(const typename Backend::polynomial_type ÷nd, + const math::polynomial_divisor_context &divisor, + polynomial_arithmetic::polynomial_context &context) { + using polynomial_type = typename Backend::polynomial_type; + polynomial_type expected_quotient, expected_remainder, quotient, remainder, reconstructed; + math::division(expected_quotient, expected_remainder, dividend, divisor.divisor()); + math::divrem(quotient, remainder, dividend, divisor, context); + BOOST_CHECK(quotient == expected_quotient); + BOOST_CHECK(remainder == expected_remainder); + BOOST_CHECK(remainder.size() < divisor.divisor().size()); + polynomial_arithmetic::schoolbook_backend {}.multiply(reconstructed, quotient, + divisor.divisor()); + math::addition(reconstructed, reconstructed, remainder); + BOOST_CHECK(reconstructed == dividend); + } + + template + void check_cached_squaremod(const math::polynomial &divisor, math::polynomial input, + std::size_t maximum_order = 58, std::size_t cyclic_order = 0) { + using backend_type = counting_mixed_radix_backend; + using polynomial_type = typename backend_type::polynomial_type; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context context {backend_type(maximum_order, counts)}; + math::polynomial_divisor_context fixed(divisor, divisor.size() - 2, context); + BOOST_CHECK_EQUAL(counts.preparations, cyclic_order ? 1 : 2); + BOOST_CHECK_EQUAL(counts.cyclic_transform_size, cyclic_order); + if (cyclic_order != 0) { + BOOST_CHECK_EQUAL(counts.last_low_transform_size, maximum_order); + } + + for (std::size_t repeat = 0; repeat < 3; ++repeat) { + polynomial_type square, quotient, expected, result; + polynomial_arithmetic::schoolbook_backend {}.square(square, input); + math::division(quotient, expected, square, divisor); + counts = {}; + math::squaremod(result, input, fixed, context); + BOOST_CHECK(result == expected); + BOOST_CHECK_EQUAL(counts.squares, 1); + BOOST_CHECK_EQUAL(counts.prepared_products, cyclic_order ? 1 : 2); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, cyclic_order ? 1 : 0); + BOOST_CHECK_EQUAL(counts.ordinary_products, 0); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + // One square and two operations using fixed spectra use six transforms. With cyclic reconstruction, + // only the final two transforms use the shorter cyclic_order; squaring and quotient work stay full size. + polynomial_type alias = input; + math::squaremod(alias, alias, fixed, context); + BOOST_CHECK(alias == expected); + check_reused_division(square, fixed, context); + input[0] += ValueType::one(); + } + } + } // namespace BOOST_AUTO_TEST_SUITE(polynomial_division_test_suite) @@ -451,4 +569,221 @@ BOOST_AUTO_TEST_CASE(fast_divrem_supports_extension_coefficients_and_base_field_ BOOST_CHECK(mixed_radix_remainder == expected_remainder); } +BOOST_AUTO_TEST_CASE(squaremod_reuses_fixed_products_for_prime_and_extension_coefficients) { + using polynomial_type = math::polynomial; + polynomial_type factor(11, fq_value_type::one()), input(11, fq_value_type::one()), divisor; + for (std::size_t i = 0; i < 10; ++i) { + factor[i] = fq_value_type(i + 2); + input[i] = fq_value_type(2 * i + 3); + } + // A reducible degree-11 modulus, first monic and then nonmonic. + polynomial_arithmetic::schoolbook_backend {}.multiply(divisor, factor, polynomial_type {1, 1}); + check_cached_squaremod(divisor, input); + for (auto &coefficient : divisor) { + coefficient *= fq_value_type(7); + } + check_cached_squaremod(divisor, input); + + math::polynomial extension_divisor(12, fq12_value_type::one()); + math::polynomial extension_input(11, fq12_value_type::one()); + for (std::size_t i = 0; i < extension_input.size(); ++i) { + extension_divisor[i] = fq12_value(12 * i + 1); + extension_input[i] = fq12_value(12 * i + 145); + } + check_cached_squaremod(extension_divisor, extension_input); + extension_divisor.back() = fq12_value(289); + check_cached_squaremod(extension_divisor, extension_input); +} + +BOOST_AUTO_TEST_CASE(prepared_division_falls_back_for_changed_quotient_lengths_without_growing_the_cache) { + using backend_type = counting_mixed_radix_backend<>; + using polynomial_type = typename backend_type::polynomial_type; + const polynomial_type divisor {2, 3, 5, 7, 11, 13, 1}; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context context(backend_type(18, counts), options); + const math::polynomial_divisor_context fixed(divisor, 5, context); + BOOST_CHECK_EQUAL(counts.preparations, 1); + BOOST_CHECK_EQUAL(counts.cyclic_preparations, 1); + + for (const auto &remainder : {polynomial_type {fq_value_type::zero()}, polynomial_type {19, 23}}) { + // Quotient size five reuses the inverse and cyclic modulus. Sizes four and two reuse only the cyclic + // modulus. At size one the ordinary remainder transform is equally short, so both products fall back. + for (std::size_t size : {5, 4, 2, 1, 5}) { + polynomial_type quotient(size, fq_value_type::zero()), dividend; + quotient.back() = fq_value_type(17); + polynomial_arithmetic::schoolbook_backend {}.multiply(dividend, quotient, divisor); + math::addition(dividend, dividend, remainder); + counts = {}; + check_reused_division(dividend, fixed, context); + BOOST_CHECK_EQUAL(counts.prepared_products, size == 5 ? 1 : 0); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, size > 1 ? 1 : 0); + BOOST_CHECK_EQUAL(counts.ordinary_products, size == 5 ? 0 : (size == 1 ? 2 : 1)); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + + polynomial_type expected_quotient, expected_remainder, aliased = dividend; + math::divrem(aliased, expected_remainder, aliased, fixed, context); + BOOST_CHECK(aliased == quotient); + BOOST_CHECK(expected_remainder == remainder); + aliased = dividend; + math::divrem(expected_quotient, aliased, aliased, fixed, context); + BOOST_CHECK(expected_quotient == quotient); + BOOST_CHECK(aliased == remainder); + } + } + for (const auto &small : {polynomial_type {fq_value_type::zero()}, polynomial_type {0, 0, 1}}) { + counts = {}; + check_reused_division(small, fixed, context); + BOOST_CHECK_EQUAL(counts.prepared_products + counts.ordinary_products + counts.preparations + + counts.cyclic_preparations + counts.cyclic_remainders, + 0); + } +} + +BOOST_AUTO_TEST_CASE(prepared_divisor_snapshots_survive_copy_move_and_arithmetic_context_changes) { + using backend_type = counting_mixed_radix_backend<>; + using polynomial_type = typename backend_type::polynomial_type; + const polynomial_type original {2, 3, 5, 7, 11, 13, 1}; + const polynomial_type dividend {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + const auto fixed = [&] { + polynomial_arithmetic::polynomial_context temporary(backend_type(18, counts), options); + polynomial_type source = original; + math::polynomial_divisor_context result(source, 5, temporary); + source[0] = fq_value_type(99); + return result; + }(); + auto copy = fixed; + const auto moved = std::move(copy); + for (std::size_t order : {522, 67}) { + polynomial_arithmetic::polynomial_context context(backend_type(order, counts), options); + for (const auto *snapshot : {&fixed, &moved}) { + counts = {}; + check_reused_division(dividend, *snapshot, context); + BOOST_CHECK(snapshot->divisor() == original); + BOOST_CHECK_EQUAL(counts.prepared_products, order == 522 ? 1 : 0); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, order == 522 ? 1 : 0); + BOOST_CHECK_EQUAL(counts.ordinary_products, order == 522 ? 0 : 2); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + } + polynomial_type different = original; + different[0] = fq_value_type(99); + const math::polynomial_divisor_context changed(different, 5, context); + check_reused_division(dividend, changed, context); + check_reused_division(dividend, fixed, context); + } +} + +BOOST_AUTO_TEST_CASE(prepared_division_handles_sparse_reversals_and_zero_fixed_prefixes) { + using backend_type = counting_mixed_radix_backend<>; + using polynomial_type = typename backend_type::polynomial_type; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context context(backend_type(18, counts), options); + // Reversing the top five coefficients of X^10 + X^4 produces [1, 0, 0, 0, 0]. Preserve this padded + // prefix length when selecting the quotient FFT. The inverse is constant for both moduli. + const polynomial_type dividend {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1}; + for (const auto &divisor : {polynomial_type {1, 0, 0, 0, 0, 0, 1}, polynomial_type {0, 0, 0, 0, 0, 0, 1}}) { + const math::polynomial_divisor_context fixed(divisor, 5, context); + counts = {}; + check_reused_division(dividend, fixed, context); + BOOST_CHECK_EQUAL(counts.prepared_products + counts.cyclic_remainders, 2); + BOOST_CHECK_EQUAL(counts.ordinary_products, 0); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + } +} + +BOOST_AUTO_TEST_CASE(unavailable_or_unused_preparations_preserve_existing_division_behavior) { + using backend_type = counting_mixed_radix_backend<>; + using polynomial_type = typename backend_type::polynomial_type; + const polynomial_type divisor {2, 3, 5}; + const polynomial_type dividend {7, 11, 13, 17, 19}; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context context(backend_type(9, counts), options); + // Newton inversion to precision seven fits order nine, but preparing its full length-seven inverse product + // would need order thirteen or larger. The optional preparation must not make construction fail. + const math::polynomial_divisor_context partial(divisor, 7, context); + BOOST_CHECK_EQUAL(counts.preparations, 2); + counts = {}; + check_reused_division(dividend, partial, context); + BOOST_CHECK_EQUAL(counts.ordinary_products, 1); + BOOST_CHECK_EQUAL(counts.prepared_products, 1); + BOOST_CHECK_EQUAL(counts.preparations, 0); + + polynomial_arithmetic::polynomial_context basecase(backend_type(9, counts)); + counts = {}; + const math::polynomial_divisor_context unprepared(divisor, 3, basecase); + BOOST_CHECK_EQUAL(counts.preparations, 0); + counts = {}; + check_reused_division(dividend, unprepared, basecase); + BOOST_CHECK_EQUAL(counts.ordinary_products + counts.prepared_products, 0); + // A context constructed under different cutoffs still works when Newton division is subsequently requested. + check_reused_division(dividend, unprepared, context); + BOOST_CHECK_EQUAL(counts.ordinary_products, 2); + BOOST_CHECK_EQUAL(counts.preparations, 0); + counts = {}; + const math::polynomial_divisor_context constant(polynomial_type {fq_value_type(7)}, 3, context); + BOOST_CHECK_EQUAL(counts.preparations, 0); +} + +BOOST_AUTO_TEST_CASE(squaremod_uses_full_quotient_transforms_and_short_cyclic_reconstruction) { + // degree(H) = 64 selects M = 67. A square of a reduced dense residue and each quotient product select 134. + math::polynomial divisor(65, fq_value_type::one()), input(64, fq_value_type::one()); + for (std::size_t i = 0; i < input.size(); ++i) { + divisor[i] = fq_value_type(i + 2); + input[i] = fq_value_type(2 * i + 3); + } + check_cached_squaremod(divisor, input, 134, 67); + + using polynomial_type = math::polynomial; + polynomial_type factor(64, fq12_value_type::one()), extension_input(64, fq12_value_type::one()), extension_divisor; + for (std::size_t i = 0; i < factor.size() - 1; ++i) { + factor[i] = fq12_value(12 * i + 1); + extension_input[i] = fq12_value(12 * i + 769); + } + // Include a reducible Fp12 modulus, both monic and nonmonic, with coefficients spanning the extension field. + polynomial_arithmetic::schoolbook_backend {}.multiply( + extension_divisor, factor, polynomial_type {fq12_value(1537), fq12_value_type::one()}); + check_cached_squaremod(extension_divisor, extension_input, 134, 67); + for (auto &coefficient : extension_divisor) { + coefficient *= fq12_value(1549); + } + check_cached_squaremod(extension_divisor, extension_input, 134, 67); +} + +BOOST_AUTO_TEST_CASE(cyclic_reconstruction_falls_back_outside_the_reduced_product_degree_range) { + using backend_type = counting_mixed_radix_backend<>; + using polynomial_type = typename backend_type::polynomial_type; + const polynomial_type divisor {2, 3, 5, 7, 11, 13, 1}; + prepared_product_counts counts; + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context context(backend_type(18, counts), options); + // A larger inverse precision still prepares the cyclic modulus for the bounded range. General division + // with degree(V) > 2*degree(H) - 2 must retain the ordinary reconstruction and its precision checks. + const math::polynomial_divisor_context fixed(divisor, 10, context); + BOOST_CHECK_EQUAL(counts.cyclic_transform_size, 6); + for (std::size_t size : {5, 6, 8, 5}) { + polynomial_type quotient(size, fq_value_type::one()), dividend; + polynomial_arithmetic::schoolbook_backend {}.multiply(dividend, quotient, divisor); + dividend[0] += fq_value_type(17); + counts = {}; + check_reused_division(dividend, fixed, context); + BOOST_CHECK_EQUAL(counts.cyclic_remainders, size == 5 ? 1 : 0); + BOOST_CHECK_EQUAL(counts.ordinary_products, size == 5 ? 1 : 2); + BOOST_CHECK_EQUAL(counts.preparations + counts.cyclic_preparations, 0); + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/libs/math/test/polynomial_frobenius.cpp b/libs/math/test/polynomial_frobenius.cpp index 0f754f5b0..87d1877aa 100644 --- a/libs/math/test/polynomial_frobenius.cpp +++ b/libs/math/test/polynomial_frobenius.cpp @@ -35,6 +35,7 @@ #include #include +#include #include namespace { @@ -171,6 +172,44 @@ BOOST_AUTO_TEST_CASE(frobenius_modulo_a_nonzero_constant_is_zero) { BOOST_CHECK(result == polynomial_type({fq_value_type::zero()})); } +BOOST_AUTO_TEST_CASE(mixed_radix_frobenius_reuses_cyclic_reduction_and_cached_composition) { + using mixed_backend = polynomial_arithmetic::mixed_radix_backend; + using reference_backend = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename mixed_backend::polynomial_type; + polynomial_type divisor(7, fq12_value_type::one()), input(6, fq12_value_type::one()); + for (std::size_t i = 0; i < divisor.size(); ++i) { + divisor[i] = fq12_value(12 * i + 1); + } + for (std::size_t i = 0; i < input.size(); ++i) { + input[i] = fq12_value(12 * i + 85); + } + // Force inverse-based division at this bounded size. For degree-six H and a dense reduced product, + // the quotient uses order nine, while reconstruction uses order six instead of the low product's 18. + polynomial_arithmetic::polynomial_context_options options; + options.basecase_divisor_coefficient_cutoff = 0; + options.basecase_quotient_coefficient_cutoff = 0; + polynomial_arithmetic::polynomial_context mixed(mixed_backend(18), options); + polynomial_arithmetic::polynomial_context reference; + const math::polynomial_divisor_context reference_divisor(divisor, 5, reference); + const math::polynomial_frobenius_context frobenius(divisor, mixed); + const auto order = fields::field_order(); + const polynomial_type x {fq12_value_type::zero(), fq12_value_type::one()}; + polynomial_type expected, result; + math::powmod(expected, x, order, reference_divisor, reference); + BOOST_CHECK(frobenius.x_to_field_order() == expected); + math::powmod(expected, input, order, reference_divisor, reference); + math::frobenius_map(result, input, frobenius, mixed); + BOOST_CHECK(result == expected); + auto alias = input; + math::frobenius_map(alias, alias, frobenius, mixed); + BOOST_CHECK(alias == expected); + // Reuse the same composition cache with another input. Frobenius fixes the added constant one. + input[0] += fq12_value_type::one(); + expected[0] += fq12_value_type::one(); + math::frobenius_map(result, input, frobenius, mixed); + BOOST_CHECK(result == expected); +} + BOOST_AUTO_TEST_CASE(frobenius_context_rejects_invalid_precomputation_inputs) { using backend_type = polynomial_arithmetic::schoolbook_backend; using polynomial_type = typename backend_type::polynomial_type;