Cache fixed FFT operands and shorten polynomial modular reduction

libs/math/docs/polynomial_arithmetic.md

@@ -130,6 +130,23 @@ math::polynomial_divisor_context<backend_type> 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`.

libs/math/include/nil/crypto3/math/polynomial/arithmetic/polynomial_division.hpp

@@ -36,6 +36,54 @@
 namespace nil::crypto3::math {
 
     namespace detail {
+        template<typename Backend>
+        concept SupportsPreparedLowProduct =
+            requires(polynomial_arithmetic::polynomial_context<Backend> &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<bool>;
+                {
+                    context.try_multiply_low_prepared(output, input, *context.prepare_low_product(input, count, count),
+                                                      count)
+                } -> std::same_as<bool>;
+            };
+
+        // Keep the optional spectra out of backends that only implement the required multiplication interface.
+        template<typename Backend>
+        struct divisor_prepared_products { };
+
+        template<SupportsPreparedLowProduct Backend>
+        struct divisor_prepared_products<Backend> {
+            using preparation_type =
+                decltype(std::declval<polynomial_arithmetic::polynomial_context<Backend> &>().prepare_low_product(
+                    std::declval<const typename Backend::polynomial_type &>(), std::size_t {}, std::size_t {}));
+            preparation_type inverse;
+            preparation_type divisor;
+        };
+
+        template<typename Backend>
+        concept SupportsPreparedCyclicRemainder =
+            requires(polynomial_arithmetic::polynomial_context<Backend> &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<bool>;
+                {
+                    context.try_cyclic_remainder(output, input, input, *context.prepare_cyclic_remainder(input, count))
+                } -> std::same_as<bool>;
+            };
+
+        template<typename Backend>
+        struct divisor_prepared_cyclic_remainder { };
+
+        template<SupportsPreparedCyclicRemainder Backend>
+        struct divisor_prepared_cyclic_remainder<Backend> {
+            using preparation_type =
+                decltype(std::declval<polynomial_arithmetic::polynomial_context<Backend> &>().prepare_cyclic_remainder(
+                    std::declval<const typename Backend::polynomial_type &>(),
+                    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<polynomial_arithmetic::PolynomialBackend Backend>
@@ -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<backend_type>) {
+                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<backend_type>) {
+                        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<backend_type> &context) const {
+            if constexpr (detail::SupportsPreparedLowProduct<backend_type>) {
+                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<backend_type> &context) const {
+            if constexpr (detail::SupportsPreparedLowProduct<backend_type>) {
+                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 &dividend,
+                                       const polynomial_type &quotient,
+                                       polynomial_arithmetic::polynomial_context<backend_type> &context) const {
+            if constexpr (detail::SupportsPreparedCyclicRemainder<backend_type>) {
+                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<backend_type> prepared_;
+        [[no_unique_address]] detail::divisor_prepared_cyclic_remainder<backend_type> 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) {