diff --git a/libs/math/docs/polynomial_arithmetic.md b/libs/math/docs/polynomial_arithmetic.md index 5409edb80..9f0832b57 100644 --- a/libs/math/docs/polynomial_arithmetic.md +++ b/libs/math/docs/polynomial_arithmetic.md @@ -10,13 +10,13 @@ The main dependency chain is: | Layer | Main types and operations | |---|---| -| Multiplication | `PolynomialBackend`, `schoolbook_backend`, `mixed_radix_backend` | +| Multiplication | `PolynomialBackend`, `schoolbook_backend`, `mixed_radix_backend`, `bn254_small_prime_fft_backend` | | Arithmetic lifetime | `polynomial_context` and `polynomial_context_options` | | Division | `inverse_series`, `polynomial_divisor_context`, `divrem`, `remainder`, `exact_division` | | GCD | `gcd`, with optional internal half-GCD reduction | -| Quotient ring | `mulmod`, `squaremod`, `powmod` | +| Quotient ring | `mulmod`, `squaremod`, `powmod`, `powmod_x` | | Modular composition | `compose_mod_reference`, `polynomial_composition_precomputation`, `compose_mod` | -| Frobenius | `polynomial_frobenius_context`, `frobenius_map` | +| Frobenius | `polynomial_frobenius_context`, `frobenius_map`, `frobenius_product`, `frobenius_weighted_sum` | ### Header map @@ -25,13 +25,19 @@ The main dependency chain is: | Representation concepts | `` | | Backend concept and arithmetic context | `` | | Schoolbook and mixed-radix backends | ``, `` | +| BN254 Fp12 small-prime backend | `` | | Power-series inversion | `` | | Divisor context and division | `` | | GCD | `` | | Quotient-ring multiplication | `` | | Quotient-ring exponentiation | `` | | Modular composition | `` | -| Frobenius maps | `` | +| Frobenius maps, products, and weighted sums | `` | + +Internal headers live beside their owning subsystem. Relative to `nil/crypto3/math/`, the small-prime backend's +parameters, arithmetic, transforms, conversion, CRT export, and pointwise kernels are grouped under +`polynomial/backends/detail/small_prime_fft/`. The prepared composition helper lives in +`polynomial/quotient_ring/detail/`. General helpers shared across math components remain in `detail/`. ## Multiplication backends @@ -47,22 +53,79 @@ backend.multiply_low(output, left, right, coefficient_count); Inputs and outputs use ascending coefficient order and canonical representation. `multiply_low` computes the product modulo `X^coefficient_count`; a zero coefficient count produces the canonical zero polynomial. -Crypto3.Math currently provides two implementations: +Crypto3.Math currently provides three implementations: * `schoolbook_backend` uses direct quadratic coefficient multiplication. It is the reference backend and is usually preferable for small operands or products with one very short operand. * `mixed_radix_backend` uses roots from `RootFieldType` to transform `ValueType` coefficients. Its constructor takes an explicit maximum transform order, caches plans for that order's divisors, and reuses transform scratch storage. The configured order must cover every requested product. +* `bn254_small_prime_fft_backend` multiplies BN254 Fp12 coefficient polynomials through auxiliary small-prime + transforms and exact signed CRT reconstruction. It takes a power-of-two maximum transform length and caches all + smaller power-of-two plans. Scalar and linear operands use direct coefficient arithmetic. This separation permits, for example, roots in a base field to transform values in an extension field. Backend selection is explicit; higher-level algorithms do not silently replace the caller's backend. +### BN254 Fp12 small-prime products + +This backend has a fixed coefficient type, available as `bn254_small_prime_fft_backend::value_type`. It converts +the twelve base-field coordinates of each coefficient into residues modulo auxiliary 50-bit primes. Each prime has +twelve transform channels; pointwise operations apply the BN254 tower identities in that prime's quotient ring. +Inverse transforms and signed CRT reconstruction return exact BN254 Fp12 coefficients. + +For base-field modulus `p` and at most `t` overlapping polynomial terms, the conservative coordinate bound is + + C = 102 * t * (p - 1)^2. + +The product of the selected primes must exceed `2*C`. Construction selects the shortest sufficient catalogue prefix +for the configured maximum length; all smaller plans share that prefix, coefficient converter, and CRT data. The +prime count is therefore fixed for a backend instance rather than selected again for every call. The 32,768-point +configuration uses eleven primes. CRT export uses a certified direct reduction into Montgomery representation when +possible and an exact reconstruction fallback otherwise. + +The FMA kernels use integer-valued binary64 residues with proved bounds. They require the supported floating-point +compiler settings and round-to-nearest mode; other supported settings or rounding modes select exact integer +arithmetic. The calling thread's rounding mode must remain unchanged throughout an operation. AVX2/FMA and AVX-512 +transform kernels depend on the compiler target, and IFMA acceleration also depends on the project's feature +configuration. A binary must run on a CPU supporting its compiled instruction set. This backend executes serially +and provides no constant-time guarantee. + +The constructor is + +```cpp +bn254_small_prime_fft_backend(max_transform_size, + max_spectrum_bytes, + max_matrix_spectrum_bytes); +``` + +Only `max_transform_size` is required. Ordinary FFT products need a transform containing the full linear product of +the retained input prefixes. `max_spectrum_bytes` defaults to the maximum `size_t` and bounds **one** spectrum at +the configured maximum length. A spectrum contains `prime_count * 12 * transform_length` doubles. Another working +spectrum, twiddle tables, output coefficients, and independently prepared operands require additional memory. +`max_matrix_spectrum_bytes` defaults to 256 MiB and bounds additional shared input spectra for optional matrix +products; zero disables that capability. Neither setting is a process-wide memory limit. + +For example, an explicitly configured context can be constructed as follows: + +```cpp +#include + +namespace pa = nil::crypto3::math::polynomial_arithmetic; +using backend_type = pa::bn254_small_prime_fft_backend; + +pa::polynomial_context arithmetic_context(backend_type(32768)); +``` + +Construction checks transform support, CRT capacity, and the spectrum limit before allocating plans. Product calls +reuse spectrum capacity, and dedicated squaring converts and transforms only one operand. Results are assembled +before replacing the output, preserving aliasing support and leaving output unchanged on failure. + ## Polynomial contexts `polynomial_arithmetic::polynomial_context` owns one backend and the algorithm-selection parameters used by division, GCD, composition, factorization, and recovery. Keeping the context alive lets a stateful backend reuse plans, -configuration, and scratch storage: +configuration, and scratch storage. Continuing with the small-prime backend above: ```cpp namespace pa = nil::crypto3::math::polynomial_arithmetic; @@ -71,7 +134,7 @@ pa::polynomial_context_options options; options.gcd_half_gcd_cutoff = 0; pa::polynomial_context arithmetic_context( - backend_type{}, options); + backend_type(32768), options); ``` The context forwards `multiply`, `square`, and `multiply_low` to its backend. It is intended for sequential reuse. @@ -83,9 +146,12 @@ Separate contexts are required for concurrent calls when the backend owns mutabl |---|---:|---| | `basecase_divisor_coefficient_cutoff` | 10 | Use quadratic long division when the divisor has at most this many coefficients. | | `basecase_quotient_coefficient_cutoff` | 2 | Use quadratic long division when the quotient has at most this many coefficients. | -| `half_gcd_basecase_cutoff` | 30 | Below this size, recursive half-GCD constructs its transformation iteratively. | +| `half_gcd_basecase_cutoff` | 30 | At or below this coefficient count, recursive half-GCD constructs its transformation iteratively; zero disables this base case. | | `gcd_half_gcd_cutoff` | 0 | Use half-GCD when the smaller operand reaches this size; zero disables half-GCD. | -| `modular_composition_cached_power_limit` | maximum `size_t` | Limit the number of Brent-Kung baby powers retained in memory. | +| `modular_composition_cached_power_limit` | maximum `size_t` | Cap the number of Brent-Kung baby powers; constructing a composition cache requires a positive limit. | +| `distinct_degree_screening_limit` | 0 | Screen initial irreducible-factor degrees before the full Kaltofen-Shoup baby table; zero disables screening. | +| `modular_composition_prepared_byte_limit` | 0 | Bound the additional prepared Fp12 coefficient payload per composition cache; zero disables preparation. | +| `modular_composition_simd_lanes` | 8 | Select four or eight IFMA lanes when supported; zero selects scalar prepared dots. The byte budget still controls whether preparation is enabled. | The two division cutoffs are independent: satisfying either enabled cutoff selects long division. Setting one of them to zero disables that criterion. Half-GCD is disabled by default because its crossover depends strongly on the @@ -132,19 +198,21 @@ another precomputation also depends on `B`, callers must keep those contexts pai ### Fixed-operand FFT preparation -For mixed-radix division, construction also prepares the FFT of the reversed-inverse prefix and one representation +With either FFT backend, construction can also prepare 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. +and a compatible transform; the small-prime backend also checks the catalogue prefix. 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 +Preparation pays one forward FFT per channel of each nonzero fixed operand up front, plus spectrum storage and +allocation. For the small-prime backend, the channels are the selected primes times twelve coordinates. Repeated compatible operations then avoid those forward FFTs. This trades setup time and memory for lower repeated cost. ## Polynomial division @@ -176,8 +244,8 @@ The high coefficients wrap identically in both terms and cancel. Since `degree(R 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 +Both FFT backends use this path for `degree(A) <= 2*degree(B) - 2`, the range from multiplying reduced residues. +Each 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. @@ -190,10 +258,11 @@ Generic `multiply_low(A, B, k)` continues to mean `A*B mod X^k`: an ordinary k-p 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 +calculation and remainder reconstruction each need one forward and one inverse FFT per channel. 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. +pointwise products, folding, allocation and other division work also contribute. For the small-prime backend, +coefficient conversion and CRT export also contribute to each product. Two convenience operations share the same dispatch: @@ -217,7 +286,14 @@ pair through Half-GCD recursively derives a product of these two-by-two polynomial matrices from the high coefficient halves, then applies it to the complete inputs. A reduction lowers the second polynomial to roughly half the original first size. The public `gcd` operation selects this internal path only when `gcd_half_gcd_cutoff` is nonzero and the smaller operand -reaches the cutoff. Below `half_gcd_basecase_cutoff`, it constructs the same transformation iteratively. +reaches the cutoff. At or below `half_gcd_basecase_cutoff`, it constructs the same transformation iteratively. + +Short-quotient Euclidean steps use direct division without constructing a reciprocal context. Recursive half-GCD +reuses the already transformed high halves when reconstructing the full pair. Backends may also expose +`try_multiply_matrix` for a row-major 2-by-2 polynomial matrix times a 2-by-`Columns` matrix. The small-prime +backend supports one or two columns, sharing input spectra and reconstructing each sum of products once when its +shape, CRT bound, and workspace budget permit. Outputs must be distinct and pointers nonnull; outputs may alias +inputs. A declined operation leaves all outputs unchanged, and half-GCD uses ordinary products and additions. Half-GCD reduces the number of sequential Euclidean divisions but introduces polynomial-matrix products and temporary storage. It is therefore not unconditionally faster; its cutoff should be tuned for the active backend and coefficient @@ -231,13 +307,16 @@ quotient to be a field. * `mulmod` computes `left * right mod B`. * `squaremod` uses the backend's dedicated square operation and reduces the result. * `powmod` performs binary exponentiation, reducing the base and every intermediate product. +* `powmod_x` computes `X^exponent mod B` with the same binary schedule, replacing multiplication by the base with + a coefficient shift and reduction of its single possible leading term. Outputs are canonical and may alias their inputs. If `B` has degree `d`, multiplying representatives already reduced modulo `B` requires at most `d - 1` inverse coefficients. Reducing an initially unreduced operand may require more. Modulo a nonzero constant, every result is the zero polynomial because the quotient is the zero ring. -`powmod` accepts built-in integer exponents and compatible multiprecision integer types. Exponents must be -nonnegative; exponent zero returns the quotient-ring identity. +Both power functions accept built-in integer exponents and compatible multiprecision integer types. Exponents must +be nonnegative; exponent zero returns the quotient-ring identity. For a positive `L`-bit exponent, `powmod_x` performs +`L - 1` modular squares and at most `L - 1` linear multiply-by-X passes. It handles nonmonic and linear divisors. ## Brent-Kung modular composition @@ -262,11 +341,32 @@ coefficient count, balances precomputation against giant-step multiplications. Construct and reuse the precomputation when composing several outer polynomials with the same inner polynomial and divisor. The one-off overload constructs it internally. `modular_composition_cached_power_limit` caps `k`, reducing -memory from cached powers at the cost of more giant-step multiplications. +memory from cached powers at the cost of more giant-step multiplications. Zero is rejected by the Brent-Kung APIs; +it does not select Horner composition automatically. `compose_mod_reference` remains available explicitly. + +### Prepared coefficient combinations + +Supported Fp12 coefficient types can additionally prepare the linear forms used by the block dot products. Set +`modular_composition_prepared_byte_limit` before constructing the composition cache to enable this preparation. +Its default is zero. `prepared_coefficient_bytes()` reports the retained payload. The budget applies to each cache +separately and excludes the ordinary baby powers, giant step, temporary block outputs, and per-call scratch. -The current coefficientwise formation of block values costs `O(L * degree(B))` field operations and can dominate at -large comparable degrees. Brent-Kung reduces modular multiplications; it does not remove the quadratic work needed to -form these linear combinations. +The scalar prepared kernel uses 54 base-field product coordinates from the Fp12 tower. BN254 builds with +`CRYPTO3_HAS_AVX512IFMA` can use a 23-coordinate evaluation/interpolation representation with four or eight lanes. +`modular_composition_simd_lanes` accepts 0, 4, or 8; zero keeps the scalar prepared kernel. Unsupported SIMD builds +retain the scalar prepared path when it fits the budget. Unsupported coefficient types or a budget too small for +any preparation retain the ordinary coefficient loop. + +Prepared evaluation batches up to four adjacent blocks, prepares each selected outer coefficient once per block, +and accumulates base-field dot products before recombination. Blocks with fewer than eight nonzero nonconstant +terms use the ordinary loop. SIMD preparation additionally requires a cached coefficient width of at least four. +All paths retain descending Horner order and exact field results. Prepared data is owned by the precomputation and +read-only during composition; copies own independent data, and scratch belongs to the call. This optimization is +independent of the polynomial multiplication backend and can be used with schoolbook, mixed-radix, or small-prime products. + +The formation of block values still costs `O(L * degree(B))` field operations and can dominate at large comparable +degrees. Brent-Kung reduces modular multiplications; it does not remove the quadratic work needed to form these linear +combinations. Prepared dots and SIMD improve its constants without changing that asymptotic bound. ## Iterated Frobenius maps @@ -275,9 +375,9 @@ Let the finite coefficient field `K` contain `Q` elements. Every coefficient sat A(X)^Q mod B = A(X^Q mod B) mod B. -`polynomial_frobenius_context` computes and stores `X^Q mod B`, owns the required divisor context, and builds a -Brent-Kung precomputation for composition with that value. The context can then apply many Frobenius maps without -repeating the field-order exponentiation or rebuilding baby powers. +`polynomial_frobenius_context` computes `X^Q mod B` using `powmod_x`, owns the required divisor context, and +builds a Brent-Kung precomputation for composition with that value. The context can then apply many Frobenius maps +without repeating the field-order exponentiation or rebuilding baby powers. ```cpp math::polynomial_frobenius_context frobenius_context( @@ -293,6 +393,37 @@ One map raises a quotient-ring element to the `Q`-th power. The iterated overloa need not be irreducible. This is a polynomial quotient-ring Frobenius operation and is distinct from any specialized coordinate-level Frobenius implementation supplied by an extension-field element type. +### Frobenius products and weighted sums + +Two operations compute aggregates of conjugates by binary doubling. Let `F` be a power of the coefficient-field +Frobenius, so `F` fixes coefficients, and supply the canonical reduced image `F(X) = X^(Q^s) mod B`: + +```cpp +math::frobenius_product(output, input, term_count, image, + divisor_context, arithmetic_context); +math::frobenius_weighted_sum(output, input, weight, term_count, image, + divisor_context, arithmetic_context); +``` + +The product is `input * F(input) * ... * F^(term_count - 1)(input)`. For the weighted sum, define + + P_0 = 1, + P_i = weight * F(weight) * ... * F^(i - 1)(weight). + +The result is `sum(P_i * F^i(input), i = 0 .. term_count - 1)`. An empty product is the quotient-ring identity; +an empty sum is zero. A singleton is the reduced input. Outputs are canonical and may alias any polynomial input. +Inputs must be canonical, and the supplied image must be reduced. Its algebraic identity is a caller precondition; +the functions do not repeat the exponentiation to verify it. Irreducibility of `B` is not required. + +These aggregates use `O(log(term_count))` compositions and modular products, with composition-cache construction +included at each needed doubling. Only the initial cache and one doubling cache are live at once. The iterated +`frobenius_map` overload above still applies one map per iteration; it does not use this aggregate algorithm. + +Square-root recovery also uses an internal characteristic-power product for extension coefficients. If `Q = p^e`, +raising to `p` requires a coefficient Frobenius twist as well as composition with `X^p`. The public aggregate APIs +fix the coefficient field; passing `X^p` to them in place of `X^Q` is generally incorrect when `e > 1`. See +[polynomial recovery](@ref math_polynomial_recovery) for the characteristic-sized square-root seed. + ## Context ownership and reuse The arithmetic context normally has the longest lifetime. Divisor, composition, Frobenius, and square-root contexts @@ -309,6 +440,7 @@ coefficient-field order. These are rough bounds for the current implementation: |---|---| | Schoolbook multiplication or squaring | `O(n^2)` | | Mixed-radix multiplication | `O(N * sum(radices))`, where `N` is the selected transform size | +| Small-prime multiplication | `O(r * N * log N)` small-residue transform work, plus conversion, tower pointwise work, and CRT/export; `r` is the selected prime count and there are twelve coordinate channels per prime | | Newton series inversion | `O(M(n))` | | Divisor-context construction | `O(M(n))` | | Long division | `O(n^2)` | @@ -317,11 +449,15 @@ coefficient-field order. These are rough bounds for the current implementation: | Half-GCD | `O(M(n) * log n)` | | `mulmod` or `squaremod` | `O(M(n))` with fast division; `O(n^2)` with schoolbook arithmetic | | `powmod` | `O(M(n) * log e)` with fast division | +| `powmod_x` | `O(M(n) * log e)` with fast division, with linear shifts/reductions for set bits | | Reference modular composition | `O(n * M(n))` | | Current Brent-Kung composition | `O(sqrt(n) * M(n) + n^2)` | | Frobenius-context construction | `O(M(n) * log Q + sqrt(n) * M(n))` | | One cached Frobenius map | `O(sqrt(n) * M(n) + n^2)` | +| Frobenius product or weighted sum of `t` terms | `O(log t)` compositions, composition-cache preparations, and modular products for `t > 1` | -The `n^2` term in Brent-Kung and Frobenius maps is the current coefficient-by-coefficient block-formation phase. -Mixed-radix `multiply_low` truncates its inputs but still transforms their complete prefix product, so it has the same -rough asymptotic cost as multiplying those prefixes in full. +The `n^2` term in Brent-Kung and Frobenius maps is the block-formation phase, including when prepared dots are enabled. +Both FFT backends' `multiply_low` operations truncate their inputs but still transform their complete prefix product, +so they have the same rough asymptotic cost as multiplying those prefixes in full. They do not currently use truncated +transforms. The transform-only cost for the small-prime backend is not a complete product-cost estimate: conversion +and exact export can be substantial and must be included in measurements. diff --git a/libs/math/docs/polynomial_factorization.md b/libs/math/docs/polynomial_factorization.md index bd5e3f029..ed4163b46 100644 --- a/libs/math/docs/polynomial_factorization.md +++ b/libs/math/docs/polynomial_factorization.md @@ -91,9 +91,26 @@ Cantor-Zassenhaus repeats randomized split trials until a nontrivial factor is f an arbitrary retry limit. A generator that does not adequately explore the coefficient field can therefore prevent progress. Callers that require a time or work limit should enforce it through their generator or at a higher level. -For large products, callers can instead construct a `mixed_radix_backend` with an adequate -transform order and place it in the polynomial context. This permits base-field roots to transform extension-field -coefficients. Backend choice changes the arithmetic implementation, not the factorization API or result. +For large products, callers can construct a `mixed_radix_backend` with an adequate transform +order, or a `bn254_small_prime_fft_backend` for BN254 Fp12 coefficients, and place it in the polynomial context. +The latter uses auxiliary-prime transforms and exact CRT export, so its transform lengths do not depend on roots +available in the coefficient field. Both backends reuse plans and can prepare fixed divisor operands. Backend choice +changes the arithmetic implementation without changing the factorization API or its reconstruction contract. See +the [backend and memory contracts](@ref math_polynomial_arithmetic) before selecting the maximum transform size. + +Several independent options affect the work in this pipeline: + +| Option | Affected work | Default | +|---|---|---| +| `gcd_half_gcd_cutoff` | Recursive half-GCD during square-free decomposition, screening, and splitting | 0: disabled | +| `distinct_degree_screening_limit` | Early degree checks before the full baby table | 0: disabled | +| `modular_composition_cached_power_limit` | Baby powers retained by each Brent-Kung composition cache | Maximum `size_t`: permits the usual square-root block size | +| `modular_composition_prepared_byte_limit` | Additional prepared Fp12 coefficient combinations | 0: disabled | +| `modular_composition_simd_lanes` | Four- or eight-lane IFMA for eligible prepared combinations | 8; still requires a positive byte budget and build support | + +Selecting the small-prime backend does not enable half-GCD or prepared coefficient combinations automatically. +Choose their cutoffs and budgets from representative complete operations, including preparation and memory costs. +Composition budgets apply to each cache, and simultaneous workers need separate arithmetic contexts and generators. ## Staged factorization @@ -114,11 +131,19 @@ The factor that triggers the stop is included in `result.factors`, and `result.c `complete == true`. Square-free and distinct-degree factorization also expose staged callbacks, allowing a caller to stop at the earliest useful stage rather than computing a complete irreducible factorization. -Complete factorization splits each distinct-degree group and invokes its irreducible-factor callbacks as soon as the -group is discovered. A callback stop also stops distinct-degree discovery, so later groups and square-free parts are -not processed. Factors retain their square-free multiplicities and the same group/splitting order. The shared -generator is consumed in the same order, including when callbacks use it, and callback or generator exceptions -propagate to the caller. +Complete factorization finishes square-free decomposition before invoking any irreducible-factor callback. It then +splits each distinct-degree group and invokes those callbacks as soon as the group is discovered. A callback stop +skips the remaining distinct-degree discovery and equal-degree splitting, including work on later square-free parts. +Factors retain their square-free multiplicities and group/splitting order. The caller-owned generator +is shared with callbacks, so their sampling affects its state for subsequent splitting. Callback and generator +exceptions propagate to the caller. Reproducibility assumes the same implementation, configuration, and callback +behavior; changing a recovery algorithm can change sampling even with the same initial seed. + +The public factorization APIs collect the factors they report. Internally, complete factorization releases each +distinct-degree group after splitting and any required exact division instead of retaining all processed groups. +Its internal callback path also allows norm recovery to consume irreducible factors without collecting a second +history. Square-free parts, active Frobenius precomputation, and pending equal-degree work remain available for the +algorithm; this does not make its total memory use constant. By default, each nonlinear square-free part prepares the initial Frobenius power and full baby-step table before discovering groups. To screen small degrees first, set a numeric limit on the existing arithmetic options: @@ -145,7 +170,7 @@ both rejected and suitable inputs before choosing a limit. Passing a bounded scr recovery will succeed. Splitting and callback work run while the active Frobenius precomputation remains alive, which can increase peak -memory. Splitting copies the emitted group because distinct-degree refinement or screening still needs its stored +memory. Splitting copies the emitted group because distinct-degree refinement or screening still needs its original polynomial after the callback returns. No default screening bound is selected from the behavior of one candidate. ## Preconditions and current restrictions @@ -189,5 +214,8 @@ backend. | Complete factorization | The sum of the square-free, Kaltofen-Shoup, and equal-degree stages | The Cantor-Zassenhaus bounds are expected bounds for uniform random samples; a caller-supplied generator has no fixed -retry limit. Each Frobenius map uses the current Brent-Kung implementation and therefore includes its quadratic -coefficient-combination phase. Staged callbacks stop paying these costs once the caller has obtained enough factors. +retry limit. Initial `X^Q mod B` construction uses `powmod_x`; subsequent Frobenius maps use Brent-Kung composition. +Its coefficient-combination phase remains quadratic even when prepared dots and SIMD reduce the constant cost. +Staged callbacks stop paying these costs once the caller has obtained enough factors. If a callback performs norm +recovery, that additional work is outside these factorization-only bounds; the +[recovery guide](@ref math_polynomial_recovery) describes odd- and even-degree root paths and factor combination. diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md index cf5550bcb..6272d51f9 100644 --- a/libs/math/docs/polynomial_recovery.md +++ b/libs/math/docs/polynomial_recovery.md @@ -23,6 +23,7 @@ polynomial arithmetic context and the coefficient-field generator. | Field orders and multiplicative-group decomposition | `` | | Coefficient-field square-root helpers | `` | | Quotient-field square testing and roots | `` | +| Frobenius products and weighted sums | `` | | Bounded rational reconstruction | `` | | `X`-norm reconstruction | `` | @@ -45,6 +46,9 @@ criterion. `sqrt()` requires a square input and asserts that precondition; calle test first. The implementation uses the field's multiplicative-group decomposition and Tonelli-Shanks rather than converting the value to a generic integer representation. +`algebra::fields::sqrt_known_square(value)` is the helper used after a square test or a mathematical invariant has +established that the input is square. Its precondition is the same; it is not a square-classification API. + ## Quotient-field square testing Let `K` be a finite field and let `B` be irreducible of degree `d`. The quotient `K[X]/(B)` is a field with @@ -65,6 +69,8 @@ characterize squares if the quotient has zero divisors. ## Quotient-field square roots +### General residues: Tonelli-Shanks + `square_root_mod` uses Tonelli-Shanks in `K[X]/(B)`. Repeated calls should share two immutable precomputations: * `polynomial_divisor_context` caches the divisor and the inverse needed for modular reduction. @@ -106,6 +112,68 @@ The square-root operation supports odd characteristic. It returns `false` and st square; zero is returned as its own root. Output may alias input. The coefficient field may itself be an extension field, including BN254 Fq12. +### Specialized roots of X + +Norm recovery needs a square root of `x = X mod B`. Write `d = degree(B)` and `Q = order(K)`. Two specialized APIs +avoid constructing a quotient-field nonresidue and avoid exponentiations whose bit length grows with `d`: + +| API | Divisor degree | Root construction | +|---|---|---| +| `square_root_mod_x_odd_degree(output, divisor_context, arithmetic_context)` | Positive odd degree | A coefficient-field norm root and a product of Frobenius conjugates; no generator. | +| `square_root_mod_x_even_degree(output, divisor_context, arithmetic_context, generator)` | Positive even degree | A weighted Frobenius trace, followed by a coefficient-field root; samples only after a zero trace. | + +Both require an irreducible divisor over an odd-characteristic field with coefficient square testing and +`sqrt_known_square` support. They support nonmonic divisors. The norm + + N = (-1)^d * B(0) / leading_coefficient(B) + +decides whether `x` is square. A nonsquare returns `false` and sets output to zero; success returns a canonical +reduced root. Irreducibility is a caller precondition and is not rechecked. Degree `d - 1` divisor-inverse coefficients +suffice for `d > 1`. Compositions require a positive `modular_composition_cached_power_limit`. Invalid degree parity, +characteristic two, an invalid generated representative, or insufficient arithmetic capacity is rejected. + +For odd `d = 1`, `x` is already a coefficient-field element, so the answer is its scalar root, including zero when +`B` is a multiple of `X`. This case needs neither composition nor sampling. For larger odd degrees, let + + b = x^((Q + 1) / 2), + D = b^Q * b^(Q^3) * ... * b^(Q^(d - 2)). + +Then `D^2 = N/x`, so `x * D / sqrt(N)` is a root. The product is evaluated by binary doubling using +`O(log d)` compositions and modular products after seed construction. No random coefficients are consumed. + +For even degrees, let `F(u) = u^Q`, `lambda = x^((Q - 1) / 2)`, and + + P_0 = 1, + P_i = lambda * F(lambda) * ... * F^(i - 1)(lambda), + W(c) = sum(P_i * F^i(c), i = 0 .. d - 1). + +For a root `r` of `x`, `r * W(c) = Trace(r*c)` belongs to the coefficient field. When `W(c)` is nonzero, +`C = x * W(c)^2` is a nonzero coefficient-field square, and `x * W(c) / sqrt(C)` is a root of `x`. +Each weighted sum uses `O(log d)` compositions and modular products. + +The even-degree API tries `c = 1` first without calling the generator. If the weighted trace is zero, it asks the +generator for another canonical polynomial of degree below `d`. It retains `lambda` and `X^Q` across retries; +each weighted-sum call constructs its own composition caches. Uniform quotient-field sampling gives zero trace with +probability `1/Q`. There is no fixed retry limit: the generator must eventually supply a nonzero weighted trace, and +generator exceptions propagate. This low-level generator returns whole polynomial representatives; the norm-recovery +APIs below instead accept a coefficient generator and adapt it internally. + +### Shared seed for extension coefficients + +Both specialized methods need `lambda = x^((Q - 1) / 2)` and `X^Q mod B`. When the coefficient field has order +`Q = p^e`, `e > 1`, and its value type supplies `Frobenius_map`, the seed starts with + + a = x^((p - 1) / 2), + x^p = x * a^2, + lambda = a * a^p * ... * a^(p^(e - 1)), + x^Q = x * lambda^2. + +The quotient-ring binary exponentiation therefore uses the characteristic-sized exponent `(p - 1)/2`. The conjugate +product uses `O(log e)` compositions and products, with the required coefficient Frobenius twists. Prime fields and +coefficient types without that map use `powmod_x` with `(Q - 1)/2` directly. Multiplication by `x` uses a linear shift/reduction. +The two methods share this seed implementation; ordinary `polynomial_frobenius_context` construction still computes +its initial `X^Q` using `powmod_x` with exponent `Q`. + ## Bounded rational reconstruction Given a reduced residue `R` modulo a nonconstant polynomial `B`, `rational_reconstruct` searches for polynomials `P` @@ -178,14 +246,27 @@ if (representation) { } ``` -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. +The root path is selected before constructing any Tonelli-Shanks context: + +| Input and configuration | Selected root path | Coefficient-generator use for this root | +|---|---|---| +| Linear `g` | Scalar root through `square_root_mod_x_odd_degree` | None | +| Odd degree greater than one, positive composition cache limit | Frobenius product | None | +| Even degree, positive composition cache limit | Weighted Frobenius trace | Only after the initial `c = 1` trace is zero | +| Degree greater than one, zero composition cache limit | Tonelli-Shanks fallback | Quotient-field nonresidue search | + +The coefficient generator remains caller-owned. When sampling is needed, recovery draws `degree(g)` field values +per representative and removes trailing zero coefficients. Trace retries must eventually yield a nonzero weighted +trace; the Tonelli-Shanks fallback must eventually find a nonsquare. The selected root and normalized `(P, Q)` need +not match those produced by another root algorithm. All successful paths satisfy the same exact norm identity. + +The zero-cache fallback is available in irreducible recovery, where the factor is already supplied. When general +recovery needs complete factorization, its Frobenius compositions still require a positive cache limit. On success, the optional contains `polynomial_x_norm_representation`. Its `p` and `q` members hold the two recovered coefficient polynomials. -The function reduces `{0, 1}` modulo `g`, including when `g` is linear, and recovers `R` with `R^2 = X mod g`. It then +The selected path recovers `R` with `R^2 = X mod g`, interpreting `X` as its residue also for linear `g`. It then reconstructs `P = R * Q mod g` with degree(P) <= floor(degree(g) / 2), @@ -255,7 +336,14 @@ Complete factorization writes the remaining input as where `c` is the leading coefficient and each `g` is monic and irreducible. An even factor power `g^(2r)` has the immediate representation `(g^r, 0)`. For an odd power `g^(2r+1)`, irreducible recovery first obtains `(P_g, Q_g)` for the unpaired copy of `g`; scaling both components by `g^r` then represents the complete factor power. Factor recovery -stops as soon as an odd-multiplicity factor cannot be represented. +stops as soon as an odd-multiplicity factor cannot be represented. It consumes each irreducible factor without +retaining a collected factor list; only the resulting representations are saved for the balanced combination below. + +Multiplicity parity and irreducible-degree parity serve different purposes. Even **multiplicity** gives the direct +square representation above, regardless of the factor's degree. For an unpaired copy, the factor's **degree** chooses +the odd-degree product or even-degree trace root algorithm. Recovery runs inside the irreducible-factor callback; +it is not deferred until every factor has passed its eligibility checks. See +[staged factorization](@ref math_polynomial_factorization) for screening and callback ordering. The factor-power representations are combined with `multiply_polynomial_x_norm_representations` in a balanced product tree. Finally, both components are scaled by `sqrt(c)` to incorporate the factorization's leading coefficient. A @@ -267,15 +355,18 @@ When `-1` is square, separate normalization loses no representations: every moni scalar obstructions in individual factors and the leading scalar can cancel in `H`; the algorithm does not exploit such cancellation. -The coefficient generator is shared by complete factorization and odd-factor recovery. It must meet the factorization -generator contract and must allow each square-root context to find a quotient-field nonsquare. All polynomial products -and squares use the caller's compile-time-selected arithmetic context; no polynomial backend or random engine is -constructed internally. +The coefficient generator is shared by complete factorization and recovery of odd-multiplicity factors. It must meet +the factorization generator contract and support any trace retries or fallback nonresidue searches described above. +Odd-degree root recovery consumes no samples, but equal-degree factorization can still do so. Changing the root path +can therefore change the generator state seen by later factorization work and the resulting factor order. A seeded +run is reproducible with the same implementation and configuration; identical sampling across different algorithms +is not part of the contract. All polynomial products and squares use the caller's compile-time-selected arithmetic +context; no polynomial backend or random engine is constructed internally. | Outcome | Contract | |---|---| | Representation returned | `evaluate_polynomial_x_norm(result, context) == H` exactly. | -| No value returned | A special-case or necessary square test fails, an odd factor is unrepresentable, bounded recovery fails, or scalar normalization is impossible. | +| No value returned | A special-case or necessary square test fails, an odd-multiplicity factor is unrepresentable, bounded recovery fails, or scalar normalization is impossible. | | `std::invalid_argument` | `H` is empty or noncanonical, or a composed API contract is violated. | | `std::logic_error` | Completed factorization, recovery, combination, normalization, or final verification is internally inconsistent. | @@ -330,24 +421,36 @@ 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 -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. +recovery API reuses the caller's polynomial arithmetic context throughout recovery and constructs one divisor context +per irreducible recovery. Only the Tonelli-Shanks fallback constructs a `polynomial_square_root_context`. Direct users +of that lower-level API can retain both contexts across repeated roots modulo the same divisor. + +The specialized paths construct composition caches for the shared seed and the degree-dependent product or trace. +Each binary chain retains only its initial cache and one doubling cache at a time, and caches are released between +phases. Optional prepared coefficient payloads obey `modular_composition_prepared_byte_limit` **per cache**. Ordinary +cached powers, backend scratch, and concurrent operations require additional memory; this option is not a total +memory limit. Measure seed construction, degree-dependent work, and complete recovery together when assessing speed. -Let `Q` be the coefficient-field order, `d` the irreducible divisor degree, and `s` the two-adicity of `Q^d - 1`. +Let `Q = p^e` 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. -The table counts quotient-ring multiplications and Euclidean steps. The cost of each underlying polynomial -multiplication, reduction, and division is given in the +The table counts quotient-ring multiplications, compositions, and Euclidean steps. Composition includes preparing +caches where needed; it is not a constant-time field operation. The cost of the underlying arithmetic is given in the [polynomial arithmetic cost overview](@ref math_polynomial_arithmetic) and depends on the selected backend. | Operation | Rough current bound | |---|---| | General `is_square_mod` | `O(d * log Q)` quotient-ring multiplications | | `is_square_mod(X)` norm shortcut | One coefficient-field square test | -| Square-root context construction | Expected `O(d * log Q)` quotient-ring multiplications | -| `square_root_mod` | `O(d * log Q + s^2)` quotient-ring multiplications | +| Tonelli-Shanks context construction | Expected `O(d * log Q)` quotient-ring multiplications with uniform sampling | +| General `square_root_mod` | `O(d * log Q + s^2)` quotient-ring multiplications | +| Shared seed with coefficient Frobenius support | `O(log p)` quotient-ring products plus `O(log e)` compositions/products and coefficient twists | +| Shared seed without coefficient Frobenius support | `O(log Q)` quotient-ring products | +| Specialized odd-degree root, after its seed | `O(log d)` compositions/products, coefficient-field root/normalization, and linear finalization; linear divisors need only the scalar root | +| Specialized even-degree root, after its seed | `O(log d)` compositions/products per trace attempt, coefficient-field root/normalization, and linear finalization | | Rational reconstruction | Up to `O(d)` sequential Euclidean steps; `O(d^2)` with schoolbook short-by-long products | -The square-root context has no deterministic bound with a caller-supplied generator because it searches until it -finds a nonresidue. Rational reconstruction does not currently use half-GCD acceleration. +Tonelli-Shanks nonresidue search and even-degree trace retries have no deterministic attempt bound for an arbitrary +caller-supplied generator. Odd-degree specialized roots are deterministic and require no quotient-field sampling. +Rational reconstruction does not currently use half-GCD acceleration. These operation bounds describe the algorithms; +they do not imply a fixed speedup independent of degree, backend, cache configuration, or CPU. diff --git a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp b/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp deleted file mode 100644 index 27256ca4e..000000000 --- a/libs/math/include/nil/crypto3/math/detail/bn254_small_prime_fft.hpp +++ /dev/null @@ -1,2161 +0,0 @@ -//---------------------------------------------------------------------------// -// 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_DETAIL_BN254_SMALL_PRIME_FFT_HPP -#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include - -#ifndef CRYPTO3_SMALL_PRIME_FFT_RADIX -#define CRYPTO3_SMALL_PRIME_FFT_RADIX 4 -#endif -static_assert(CRYPTO3_SMALL_PRIME_FFT_RADIX == 4 || CRYPTO3_SMALL_PRIME_FFT_RADIX == 8 || - CRYPTO3_SMALL_PRIME_FFT_RADIX == 16, - "small-prime FFT stage radix must be 4, 8 or 16"); - -#if defined(CRYPTO3_HAS_AVX512IFMA) && defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_AVAILABLE) && \ - !defined(CRYPTO3_SMALL_PRIME_CRT_DISABLE_IFMA) -#define CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED -#endif - -namespace nil::crypto3::math::detail { - - // Independent coefficient lanes, using the existing IFMA operation policy. A CRT digit is - // below 2^50 and each five-digit export weight is below p. With at most twelve primes, - // every unnormalized column is below 24*2^52 and the complete sum is below 2^54*p. - template - struct small_prime_fft_crt_accumulator { - using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; - using vector = typename Operations::value_type; - std::array columns; - - small_prime_fft_crt_accumulator() { - columns.fill(Operations::broadcast(0)); - } - - void add(const typename radix::value_type &weight, vector digits) { - for (std::size_t limb = 0; limb < weight.size(); ++limb) { - const auto w = Operations::broadcast(weight[limb]); - columns[limb] = Operations::multiply_add_low(columns[limb], digits, w); - columns[limb + 1] = Operations::multiply_add_high(columns[limb + 1], digits, w); - } - } - - // The certificate reciprocal has up to 61 bits, so it cannot be an IFMA operand - // without splitting. Compute floor(d*r/2^50) exactly from four 52-bit split products. - static vector certificate_term(vector digits, std::uint64_t reciprocal) { - const auto zero = Operations::broadcast(0); - const auto low = Operations::broadcast(reciprocal & radix::mask); - const auto high = Operations::broadcast(reciprocal >> 52); - const auto a = Operations::multiply_add_low(zero, digits, low); - const auto b = Operations::multiply_add_high(zero, digits, low); - const auto c = Operations::multiply_add_low(zero, digits, high); - const auto d = Operations::multiply_add_high(zero, digits, high); - return Operations::add( - Operations::add(Operations::template shift_right<50>(a), Operations::template shift_left<2>(b)), - Operations::add(Operations::template shift_left<2>(c), Operations::template shift_left<54>(d))); - } - - struct certificate { - vector quotients; - vector accepted; - }; - - static certificate certify(vector lower, std::uint64_t scaled_bound, std::size_t prime_count) { - constexpr std::uint64_t scale = std::uint64_t(1) << 60; - const auto zero = Operations::broadcast(0); - const auto one = Operations::broadcast(1); - const auto quotient = - Operations::template shift_right<60>(Operations::add(lower, Operations::broadcast(scale / 2))); - const auto center = Operations::template shift_left<60>(quotient); - const auto bound = Operations::broadcast(scaled_bound); - const auto upper = Operations::add(lower, Operations::broadcast(2 * prime_count)); - const auto lower_ok = - Operations::select(Operations::less(Operations::add(lower, bound), center), zero, one); - const auto upper_ok = - Operations::select(Operations::less(Operations::add(center, bound), upper), zero, one); - return {quotient, Operations::bit_and(lower_ok, upper_ok)}; - } - - // Add the selected centering correction and divide by S=2^260. The caller prepares - // weights and corrections with R*S (R=2^256), so the result is already R-Montgomery. - std::array reduce(const std::array &correction) const { - const auto zero = Operations::broadcast(0); - const auto mask = Operations::broadcast(radix::mask); - const auto inverse = Operations::broadcast(radix::p_dash()); - std::array t; - t.fill(zero); - std::copy(columns.begin(), columns.end(), t.begin()); - for (std::size_t i = 0; i < correction.size(); ++i) - t[i] = Operations::add(t[i], correction[i]); - for (std::size_t i = 0; i < 9; ++i) { - t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); - t[i] = Operations::bit_and(t[i], mask); - } - // Each cancellation adds at most five split products to a column. All columns - // stay below 2^57; the initial sum plus correction is < 2^54*p < S*p. - for (std::size_t i = 0; i < 5; ++i) { - const auto m = Operations::multiply_add_low(zero, t[i], inverse); - for (std::size_t j = 0; j < 5; ++j) { - const auto p = Operations::broadcast(radix::modulus()[j]); - t[i + j] = Operations::multiply_add_low(t[i + j], m, p); - t[i + j + 1] = Operations::multiply_add_high(t[i + j + 1], m, p); - } - t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); - } - std::array reduced; - for (std::size_t i = 5; i < 9; ++i) { - reduced[i - 5] = Operations::bit_and(t[i], mask); - t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); - } - reduced[4] = t[9]; - // REDC(T)<2*p::subtract_modulus(reduced); - return reduced; - } - }; - - /** - * Convert canonical BN254 Fp12 coordinates to a prefix of the small-prime catalogue. - * Coordinate order is Crypto3's [1,u,v,uv,v^2,uv^2,w,uw,vw,uvw,v^2w,uv^2w]. - * Each output cell is an exact integer in [0,q), congruent to the canonical base-field - * coordinate modulo q. In particular, a base-field -1 is lifted to p-1 before reduction. - * - * Precomputed powers are reusable across operands and sizes. Conversion extracts each - * coordinate once for all primes, using fixed-size scratch and caller-owned output storage. - * Kernel selection occurs once per call; the calling thread's rounding mode must remain - * unchanged until it returns. This conversion has no constant-time guarantee. - */ - class bn254_small_prime_fft_converter { - static constexpr std::size_t chunk_bits = 52; - static constexpr std::size_t chunk_count = 5; - static constexpr std::size_t block_size = 32; - using wide_type = boost::multiprecision::uint128_modular_t; - - struct prime_context { - small_prime_fft_arithmetic arithmetic; - std::array powers; - - explicit prime_context(std::uint64_t modulus) : arithmetic(modulus) { - std::uint64_t power = 1; - for (auto &weight : powers) { - power = ((wide_type(power) << chunk_bits) % modulus).convert_to(); - // Center the weights so that |weight| < q/2. These conversions are exact - // under every rounding mode, and setup products fit in 102 bits. - const auto centered = static_cast(power) - - (power > modulus / 2 ? static_cast(modulus) : 0); - weight = static_cast(centered); - } - } - }; - - public: - using field_type = algebra::fields::fp12_2over3over2>; - using value_type = field_type::value_type; - - explicit bn254_small_prime_fft_converter(std::size_t prime_count) { - if (prime_count == 0 || prime_count > small_prime_fft_primes.size()) { - throw std::invalid_argument("bn254_small_prime_fft_converter: invalid catalogue prefix"); - } - primes_.reserve(prime_count); - for (std::size_t i = 0; i < prime_count; ++i) { - primes_.emplace_back(small_prime_fft_primes[i].modulus); - } - } - - std::size_t prime_count() const { - return primes_.size(); - } - - /** - * Write output[(prime*12 + coordinate)*stride + coefficient], padding each channel - * with zeros from coefficients.size() to stride. Pass a subspan to convert a prefix. - * stride need not be a power of two; the transform planner supplies that restriction. - * - * Require coefficients.size() <= stride and output.size() == prime_count()*12*stride. - * Check dimensions without overflow before writing any output; throw invalid_argument - * on mismatch. Empty input writes zeros, including when reusing a previous output buffer. - */ - void convert(std::span output, std::span coefficients, std::size_t stride) const { - const std::size_t channels = field_type::arity * primes_.size(); - if (channels == 0 || coefficients.size() > stride || stride > output.size() / channels || - output.size() != channels * stride) { - throw std::invalid_argument("bn254_small_prime_fft_converter: invalid buffer dimensions"); - } - if (output.empty()) { - return; - } - const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { - return prime.arithmetic.fma_supported(); - }); - if (use_fma) { - convert_impl(output, coefficients, stride); - } else { - convert_impl(output, coefficients, stride); - } - } - - private: - template - void convert_impl(std::span output, std::span coefficients, - std::size_t stride) const { - for (std::size_t channel = 0; channel < field_type::arity * primes_.size(); ++channel) { - auto padding = output.subspan(channel * stride + coefficients.size(), stride - coefficients.size()); - std::fill(padding.begin(), padding.end(), 0.0); - } - - std::array, block_size> chunks; - for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { - for (std::size_t offset = 0; offset < coefficients.size();) { - const std::size_t count = std::min(block_size, coefficients.size() - offset); - for (std::size_t i = 0; i < count; ++i) { - std::array words {}; - using base_backend_type = typename Params::base_field_type::integral_type::backend_type; - if constexpr (base_backend_type::limb_bits == 64) { - namespace fp12_fast = algebra::fields::detail::fp12_fast; - // The stored value is a*R mod p, with R = 2^256. Zero-extend it - // into the low half of montgomery; one reduction recovers a in [0,p). - typename Params::limb_array montgomery {}; - const auto stored = fp12_fast::load_limbs( - coefficients[offset + i].coordinate(coordinate).data.backend().base_data()); - std::copy(stored.begin(), stored.end(), montgomery.begin()); - std::array canonical; - fp12_fast::montgomery_reduce(canonical.data(), montgomery); - // pack only extracts five 52-bit chunks; it does not rescale them. - // Packing the canonical value preserves the lift of field -1 to p-1. - words = fp12_fast::montgomery_radix52::pack(canonical.data()); - } else { - const auto canonical = coefficients[offset + i].coordinate(coordinate).to_integral(); - boost::multiprecision::export_bits(canonical, words.begin(), chunk_bits, false); - } - for (std::size_t j = 0; j < chunk_count; ++j) { - chunks[i][j] = static_cast(words[j]); - } - } - for (std::size_t prime = 0; prime < primes_.size(); ++prime) { - const auto &context = primes_[prime]; - const auto &arithmetic = context.arithmetic; - const double modulus = static_cast(arithmetic.modulus()); - double *destination = - output.data() + (prime * field_type::arity + coordinate) * stride + offset; - for (std::size_t i = 0; i < count; ++i) { - double sum; - if constexpr (UseFma) { - sum = arithmetic.reduce_fma(chunks[i][0]); - } else { - sum = arithmetic.reduce_integer(chunks[i][0]); - } - // Each chunk is < 2^52 and each weight has magnitude < q/2, so the - // products are < 4*q^2 (q > 2^49). The initial remainder is < q; - // adding four terms of magnitude < 1.5*q keeps every partial sum - // inside (-7*q,7*q), hence below 2^53 and exactly representable. - for (std::size_t j = 1; j < chunk_count; ++j) { - if constexpr (UseFma) { - sum += arithmetic.multiply_fma(chunks[i][j], context.powers[j - 1]); - } else { - sum += arithmetic.multiply_integer(chunks[i][j], context.powers[j - 1]); - } - } - if constexpr (UseFma) { - sum = arithmetic.reduce_fma(sum); - } else { - sum = arithmetic.reduce_integer(sum); - } - destination[i] = sum < 0 ? sum + modulus : sum; - } - } - offset += count; - } - } - } - - std::vector primes_; - }; - - /** - * Reconstruct signed integer coordinates from a prefix of the small-prime catalogue, then - * reduce them into BN254 Fp12. The supplied bound C must cover the original integer expression, - * with M > 2*C for the prime product M. Centering must precede reduction into the target field: - * the CRT representative M-1 denotes -1, not the field element (M-1) mod p. - * - * Cofactors, inverses and Montgomery export weights are prepared once. An integer interval - * certificate permits direct export when it proves both centering and the coefficient bound; - * uncertain cases use exact wide reconstruction. Both paths use fixed-size scratch and - * caller-owned output, with no arbitrary-precision setup work in the coefficient loop. - * This operation has no constant-time guarantee. - */ - class bn254_small_prime_fft_reconstructor { - using base_field_type = algebra::fields::alt_bn128_base_field<254>; - static constexpr std::size_t block_size = 32; - static constexpr std::size_t correction_count = std::bit_width(small_prime_fft_primes.size() - 1); - // Each CRT term is < M, so their sum is < prime_count*M < 12*2^600 < 2^604. - static constexpr unsigned crt_bits = 50 * small_prime_fft_primes.size() + correction_count; - using integer_type = boost::multiprecision::number>; - static constexpr unsigned export_split = 384; - using export_low_type = - boost::multiprecision::number>; - using export_type = boost::multiprecision::uint512_modular_t; - static_assert(50 * small_prime_fft_primes.size() <= 600, "CRT export folding requires M < 2^600"); - using fast_params = bn254_small_prime_fft_converter::field_type::extension_policy::fast_params; - using direct_weight_type = std::array; -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - using radix52 = algebra::fields::detail::fp12_fast::montgomery_radix52; -#endif - static constexpr bool direct_export_supported = base_field_type::integral_type::backend_type::limb_bits == 64; - static constexpr unsigned certificate_bits = 60; - static constexpr std::uint64_t certificate_scale = std::uint64_t(1) << certificate_bits; - static_assert(std::all_of(small_prime_fft_primes.begin(), small_prime_fft_primes.end(), [](const auto &prime) { - return prime.modulus > (std::uint64_t(1) << 49) && prime.modulus < (std::uint64_t(1) << 50); - })); - - struct prime_context { - small_prime_fft_arithmetic arithmetic; - integer_type cofactor; - double inverse; - std::uint64_t certificate_reciprocal; - direct_weight_type export_weight; -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - radix52::value_type export_weight52; -#endif - - prime_context(std::uint64_t modulus, const integer_type &cofactor, const direct_weight_type &weight) : - arithmetic(modulus), cofactor(cofactor), export_weight(weight) { - using wide_type = boost::multiprecision::uint128_modular_t; - const auto residue = (cofactor % modulus).convert_to(); - inverse = static_cast(boost::multiprecision::inverse_mod(wide_type(residue), wide_type(modulus)) - .convert_to()); - certificate_reciprocal = - ((wide_type(1) << (certificate_bits + 50)) / modulus).convert_to(); - } - }; - - public: - using field_type = bn254_small_prime_fft_converter::field_type; - using value_type = field_type::value_type; - - bn254_small_prime_fft_reconstructor(const small_prime_fft_crt_params ¶ms, - const boost::multiprecision::cpp_int &coefficient_bound) { - if (params.prime_count == 0 || params.prime_count > small_prime_fft_primes.size() || - coefficient_bound < 0) { - throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid CRT parameters"); - } - boost::multiprecision::cpp_int product = 1; - for (std::size_t i = 0; i < params.prime_count; ++i) { - product *= small_prime_fft_primes[i].modulus; - } - if (params.modulus != product || 2 * coefficient_bound >= product) { - throw std::invalid_argument( - "bn254_small_prime_fft_reconstructor: inconsistent modulus or insufficient capacity"); - } - - // The checks above establish positivity and capacity before conversion to fixed precision. - using setup_backend = integer_type::backend_type::cpp_int_type; - multiples_[0].backend().from_cpp_int(setup_backend(product.backend())); - bound_.backend().from_cpp_int(setup_backend(coefficient_bound.backend())); - half_modulus_ = multiples_[0] >> 1; - const boost::multiprecision::cpp_int sum_bound = params.prime_count * product - 1; - constexpr unsigned limb_bits = integer_type::backend_type::limb_bits; - accumulator_limbs_ = (boost::multiprecision::msb(sum_bound) + 1 + limb_bits - 1) / limb_bits; - export_radix_ = (export_type(1) << export_split) % export_type(base_field_type::modulus); - for (std::size_t bit = 1; bit < correction_count; ++bit) { - multiples_[bit] = multiples_[bit - 1] << 1; - } - primes_.reserve(params.prime_count); - const boost::multiprecision::cpp_int characteristic(base_field_type::modulus.backend().to_cpp_int()); - const boost::multiprecision::cpp_int r_squared = - (boost::multiprecision::cpp_int(1) << 512) % characteristic; - scaled_bound_ = ((coefficient_bound << certificate_bits) / product).convert_to(); - for (std::size_t i = 0; i < params.prime_count; ++i) { - const auto modulus = small_prime_fft_primes[i].modulus; - const boost::multiprecision::cpp_int weight = (product / modulus) * r_squared % characteristic; - direct_weight_type limbs {}; - boost::multiprecision::export_bits(weight, limbs.begin(), 64, false); - primes_.emplace_back(modulus, multiples_[0] / modulus, limbs); -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - // S/R=16. Export after radix-52 REDC must retain the field's R=2^256 scaling. - const boost::multiprecision::cpp_int vector_weight = (weight * 16) % characteristic; - limbs.fill(0); - boost::multiprecision::export_bits(vector_weight, limbs.begin(), 64, false); - primes_.back().export_weight52 = radix52::pack(limbs.data()); -#endif - } - // The centering quotient is in [0,prime_count]. Store canonical -k*M*R^2 mod p. - for (std::size_t k = 0; k <= params.prime_count; ++k) { - boost::multiprecision::cpp_int correction = (k * product * r_squared) % characteristic; - if (correction != 0) { - correction = characteristic - correction; - } - boost::multiprecision::export_bits(correction, export_corrections_[k].begin(), 64, false); -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - const boost::multiprecision::cpp_int vector_correction = (correction * 16) % characteristic; - direct_weight_type limbs {}; - boost::multiprecision::export_bits(vector_correction, limbs.begin(), 64, false); - const auto packed = radix52::pack(limbs.data()); - for (std::size_t digit = 0; digit < packed.size(); ++digit) - export_corrections52_[digit][k] = packed[digit]; -#endif - } - } - - std::size_t prime_count() const { - return primes_.size(); - } - - /** - * Combine an inverse-transform scale with (M/q)^-1, returning a centered residue. - * Require an integer-valued scale of magnitude < 2^53; invalid prime indices throw. - */ - double weighted_normalization(std::size_t prime, double scale) const { - const auto &context = primes_.at(prime); - const double q = static_cast(context.arithmetic.modulus()); - double result = - context.arithmetic.multiply_integer(context.arithmetic.reduce_integer(scale), context.inverse); - if (result > q / 2) { - result -= q; - } else if (result < -q / 2) { - result += q; - } - return result; - } - - /** - * Read residues[(prime*12 + coordinate)*stride + coefficient] for each requested output. - * Require output.size() <= stride and residues.size() == prime_count()*12*stride. - * Check dimensions without overflow before writing any output; throw invalid_argument on - * mismatch. An empty output does no work; entries after output.size() in each channel are - * not read. Input residues must already include any inverse-transform normalization. - * - * @pre Each consumed residue is an integer-valued double of magnitude < 2^53. The thread's - * rounding mode must remain unchanged during the call; kernel selection occurs once. - * - * Throw overflow_error if a centered coordinate exceeds the declared bound; earlier output - * coordinates may already have been written. This check cannot detect an incorrect original - * bound when the integer expression has wrapped by M, so the caller must justify that bound. - */ - void reconstruct(std::span output, std::span residues, std::size_t stride) const { - reconstruct_dispatch(output, residues, stride); - } - - /** - * Reconstruct residues already multiplied by (M/q)^-1 modulo each q. All dimensions, - * residue bounds and failure semantics are those of reconstruct(). The backend can fuse - * this weight with 1/N in the pointwise pass because the inverse transform is linear. - */ - void reconstruct_weighted(std::span output, std::span residues, - std::size_t stride) const { - reconstruct_dispatch(output, residues, stride); - } - - private: - template - void reconstruct_dispatch(std::span output, std::span residues, - std::size_t stride) const { - const std::size_t channels = field_type::arity * primes_.size(); - if (channels == 0 || output.size() > stride || stride > residues.size() / channels || - residues.size() != channels * stride) { - throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid buffer dimensions"); - } - if (output.empty()) { - return; - } - const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { - return prime.arithmetic.fma_supported(); - }); - if (use_fma) { -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - reconstruct_ifma(output, residues, stride); - return; -#endif - reconstruct_impl(output, residues, stride); - } else { - reconstruct_impl(output, residues, stride); - } - } - - void add_cofactor_product(integer_type &sum, const integer_type &cofactor, std::uint64_t digit) const { - using limb_type = boost::multiprecision::limb_type; - using wide_type = boost::multiprecision::double_limb_type; - constexpr unsigned limb_bits = integer_type::backend_type::limb_bits; - if constexpr (limb_bits >= 50) { - auto *destination = sum.backend().limbs(); - const auto *source = cofactor.backend().limbs(); - wide_type carry = 0; - // A digit is < 2^50. Each product plus the accumulator limb and carry fits - // in a double limb. The complete sum is < prime_count*M, so its final carry - // fits in the active width prepared by the constructor; higher limbs stay zero. - for (std::size_t limb = 0; limb < accumulator_limbs_; ++limb) { - const wide_type value = wide_type(source[limb]) * digit + destination[limb] + carry; - destination[limb] = static_cast(value); - carry = value >> limb_bits; - } - assert(carry == 0); - assert((destination[sum.backend().size() - 1] & ~integer_type::backend_type::upper_limb_mask) == 0); - } else { - // A 50-bit digit spans more than one limb on a 32-bit limb backend. - sum += cofactor * digit; - } - } - - base_field_type::value_type export_magnitude(const integer_type &magnitude) const { - // Centering and the coefficient-bound check have already established magnitude < M < 2^600. - // Write magnitude = low + high*2^384 and replace 2^384 by its cached residue modulo p. - // The fold is < 2^384 + 2^216*2^254 < 2^471, within the existing BN254 Barrett reducer's - // range. This avoids its general-division fallback while retaining the field's normal - // Montgomery conversion. The fixed-width low conversion intentionally discards high bits. - const export_type high(magnitude >> export_split); - export_type folded {export_low_type(magnitude)}; - if (high != 0) { - folded += high * export_radix_; - } - return base_field_type::value_type(folded); - } - - base_field_type::value_type export_exact(integer_type magnitude) const { - // Each canonical CRT term is < M. Reduce their sum to [0,M) without division. - for (std::size_t bit = std::bit_width(primes_.size() - 1); bit-- > 0;) { - if (magnitude >= multiples_[bit]) { - magnitude -= multiples_[bit]; - } - } - const bool negative = magnitude > half_modulus_; - if (negative) { - magnitude = multiples_[0] - magnitude; - } - if (magnitude > bound_) { - throw std::overflow_error("bn254_small_prime_fft_reconstructor: coefficient exceeds CRT bound"); - } - const auto value = export_magnitude(magnitude); - return negative ? -value : value; - } - - static void add_export_product(fast_params::limb_array &sum, const direct_weight_type &weight, - std::uint64_t digit) { - using wide_type = algebra::fields::detail::fp12_fast::wide_limb; - wide_type carry = 0; - for (std::size_t limb = 0; limb < weight.size(); ++limb) { - const wide_type value = wide_type(weight[limb]) * digit + sum[limb] + carry; - sum[limb] = static_cast(value); - carry = value >> 64; - } - // At most twelve digits < 2^50 times weights < p, plus one correction < p. - // The total is < 2^54*p < 2^308, so the fifth limb cannot overflow. - sum[weight.size()] += static_cast(carry); - } - - bool certify_export(std::uint64_t lower, std::size_t "ient) const { - // For d < q < 2^50 and a=floor(2^110/q), floor(d*a/2^50) <= 2^60*d/q - // < floor(d*a/2^50)+2. Summing encloses 2^60*sum(d/q) in [lower,upper]. - // With at most twelve primes, all quantities below fit in uint64_t. - const std::uint64_t upper = lower + 2 * primes_.size(); - quotient = (lower + certificate_scale / 2) >> certificate_bits; - const std::uint64_t center = std::uint64_t(quotient) << certificate_bits; - // Containment in k +/- C/M proves the bound AND the unique centering quotient, - // since C/M < 1/2. Using floor(2^60*C/M) makes this test conservative. - return lower + scaled_bound_ >= center && upper <= center + scaled_bound_; - } - - template - base_field_type::value_type export_direct(typename Params::limb_array &sum, std::size_t quotient) const { - add_export_product(sum, export_corrections_[quotient], 1); - // sum = (sum(d_i*M_i)-k*M)*R^2 mod p, and sum < 2^54*p < R*p for R=2^256. - // One existing REDC therefore yields a canonical Montgomery coordinate directly, - // without materializing the wide CRT integer, folding, or a second field conversion. - direct_weight_type reduced; - algebra::fields::detail::fp12_fast::montgomery_reduce(reduced.data(), sum); - base_field_type::value_type result; - std::copy(reduced.begin(), reduced.end(), result.data.backend().base_data().limbs()); - return result; - } - - template - void reconstruct_impl(std::span output, std::span residues, - std::size_t stride) const { - for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { - for (std::size_t offset = 0; offset < output.size();) { - const std::size_t count = std::min(block_size, output.size() - offset); - using accumulator_type = std::conditional_t; - std::array sums {}; - std::array lower {}; - std::array, small_prime_fft_primes.size()> saved_digits; - for (std::size_t prime = 0; prime < primes_.size(); ++prime) { - const auto &context = primes_[prime]; - const auto &arithmetic = context.arithmetic; - const double modulus = static_cast(arithmetic.modulus()); - const double *source = - residues.data() + (prime * field_type::arity + coordinate) * stride + offset; - const auto accumulate = [&](std::size_t i, std::uint64_t digit) { - if constexpr (Direct) { - using wide_type = algebra::fields::detail::fp12_fast::wide_limb; - saved_digits[prime][i] = digit; - lower[i] += static_cast( - (wide_type(digit) * context.certificate_reciprocal) >> 50); - add_export_product(sums[i], context.export_weight, digit); - } else { - add_cofactor_product(sums[i], context.cofactor, digit); - } - }; - std::size_t i = 0; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { - if (count >= 4) { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic); - const __m256d inverse = _mm256_set1_pd(context.inverse); - const __m256d q = _mm256_set1_pd(modulus); - for (; count - i >= 4; i += 4) { - __m256d digits = packed.reduce_fma(_mm256_loadu_pd(source + i)); - if constexpr (!Weighted) { - digits = packed.multiply_fma(digits, inverse); - } - const __m256d negative = _mm256_cmp_pd(digits, _mm256_setzero_pd(), _CMP_LT_OQ); - digits = _mm256_add_pd(digits, _mm256_and_pd(negative, q)); - std::array canonical; - _mm256_storeu_pd(canonical.data(), digits); - for (std::size_t lane = 0; lane < canonical.size(); ++lane) { - accumulate(i + lane, static_cast(canonical[lane])); - } - } - } - } -#endif - for (; i < count; ++i) { - double digit; - if constexpr (UseFma) { - digit = arithmetic.reduce_fma(source[i]); - if constexpr (!Weighted) { - digit = arithmetic.multiply_fma(digit, context.inverse); - } - } else { - digit = arithmetic.reduce_integer(source[i]); - if constexpr (!Weighted) { - digit = arithmetic.multiply_integer(digit, context.inverse); - } - } - // Both operands have magnitude < q, so the product remainder is in (-q,q). - // Canonical digits give 0 <= digit*(M/q) < M and bound the complete sum. - const auto canonical = static_cast(digit < 0 ? digit + modulus : digit); - accumulate(i, canonical); - } - } - for (std::size_t i = 0; i < count; ++i) { - auto &destination = output[offset + i].coordinate(coordinate); - if constexpr (Direct) { - std::size_t quotient; - if (certify_export(lower[i], quotient)) { - destination = export_direct(sums[i], quotient); - } else { - integer_type exact_sum = 0; - for (std::size_t prime = 0; prime < primes_.size(); ++prime) { - add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][i]); - } - destination = export_exact(exact_sum); - } - } else { - destination = export_exact(sums[i]); - } - } - offset += count; - } - } - } - -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - template - void reconstruct_ifma(std::span output, std::span residues, - std::size_t stride) const { - using operations = algebra::fields::detail::fp12_fast::ifma52_operations<8>; - for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { - for (std::size_t offset = 0; offset < output.size();) { - const std::size_t count = std::min(operations::lanes, output.size() - offset); - const auto active = static_cast<__mmask8>((1u << count) - 1); - using accumulator_type = small_prime_fft_crt_accumulator; - accumulator_type accumulator; - auto lower = operations::broadcast(0); - std::array, small_prime_fft_primes.size()> - saved_digits; - for (std::size_t prime = 0; prime < primes_.size(); ++prime) { - const auto &context = primes_[prime]; - const small_prime_fft_arithmetic::wide_fma_kernel packed(context.arithmetic); - const double *source = - residues.data() + (prime * field_type::arity + coordinate) * stride + offset; - auto digits = packed.reduce_fma(_mm512_maskz_loadu_pd(active, source)); - if constexpr (!Weighted) { - digits = packed.multiply_fma(digits, _mm512_set1_pd(context.inverse)); - } - const auto negative = _mm512_cmp_pd_mask(digits, _mm512_setzero_pd(), _CMP_LT_OQ); - digits = _mm512_mask_add_pd(digits, negative, digits, - _mm512_set1_pd(static_cast(context.arithmetic.modulus()))); - const operations::value_type canonical {_mm512_cvttpd_epu64(digits)}; - operations::store(saved_digits[prime].data(), canonical); - accumulator.add(context.export_weight52, canonical); - lower = operations::add( - lower, accumulator_type::certificate_term(canonical, context.certificate_reciprocal)); - } - const auto certificate = accumulator_type::certify(lower, scaled_bound_, primes_.size()); - const auto accepted = operations::less(operations::broadcast(0), certificate.accepted) & active; - std::array, 5> reduced; - if (accepted != 0) { - std::array corrections; - for (std::size_t digit = 0; digit < corrections.size(); ++digit) { - // Every quotient is in [0,prime_count], including uncertified lanes. - corrections[digit] = {_mm512_i64gather_epi64(certificate.quotients.data, - export_corrections52_[digit].data(), 8)}; - } - const auto values = accumulator.reduce(corrections); - for (std::size_t digit = 0; digit < values.size(); ++digit) - operations::store(reduced[digit].data(), values[digit]); - } - for (std::size_t lane = 0; lane < count; ++lane) { - auto &destination = output[offset + lane].coordinate(coordinate); - if ((accepted & (1u << lane)) != 0) { - // Numeric repacking only: the R*S weights already account for S/R. - base_field_type::value_type value; - auto *words = value.data.backend().base_data().limbs(); - words[0] = reduced[0][lane] | (reduced[1][lane] << 52); - words[1] = (reduced[1][lane] >> 12) | (reduced[2][lane] << 40); - words[2] = (reduced[2][lane] >> 24) | (reduced[3][lane] << 28); - words[3] = (reduced[3][lane] >> 36) | (reduced[4][lane] << 16); - destination = value; - } else { - integer_type exact_sum = 0; - for (std::size_t prime = 0; prime < primes_.size(); ++prime) { - add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][lane]); - } - destination = export_exact(exact_sum); - } - } - offset += count; - } - } - } -#endif - - std::vector primes_; - std::array multiples_; - integer_type half_modulus_; - integer_type bound_; - std::size_t accumulator_limbs_; - export_type export_radix_; - std::uint64_t scaled_bound_; - std::array export_corrections_ {}; -#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) - std::array, 5> export_corrections52_ {}; -#endif - }; - - /** - * In-place, power-of-two transforms for one prime from the catalogue. The forward transform - * maps natural-order coefficients to evaluations at omega^k in bit-reversed index order. - * Both directions use this order, so convolution needs no permutation passes. - * - * A reusable table of max(1,N/2) centered twiddles serves both directions. Radix-4 recursion - * fuses pairs of layers, with radix-2 leaves for odd depths and special cases for unit twiddles. - * AVX2/FMA builds fuse 8-, 16- and 32-point leaves; AVX-512 adds paired 32-/64-point - * leaves and eight-lane stages. CRYPTO3_SMALL_PRIME_FFT_RADIX selects 4 (default), 8 - * or 16 for large SIMD stages; the scalar/integer fallback keeps the radix-4 schedule. - * Transform calls allocate no buffers; recursion uses O(log N) stack space. - * - * @pre Inputs are integer-valued doubles in (-3*q,3*q). The calling thread's rounding mode - * must remain unchanged during a call. Each call selects FMA or integer arithmetic once. - * The forward output is in (-3*q,3*q); the normalized inverse output is in [0,q). - * This operation has no constant-time guarantee. - */ - class small_prime_fft_plan { - public: - small_prime_fft_plan(std::size_t prime_index, std::size_t size) : - arithmetic_(checked_prime(prime_index, size).modulus), size_(size), - twiddles_(std::max(std::size_t(1), size / 2), 1.0) { - using wide_type = boost::multiprecision::uint128_modular_t; - const auto &prime = small_prime_fft_primes[prime_index]; - wide_type root = prime.root; - const unsigned depth = std::countr_zero(size_); - for (unsigned level = prime.two_adicity; level > depth; --level) { - root = root * root % prime.modulus; - } - omega_ = root.convert_to(); - size_inverse_ = centered( - static_cast(boost::multiprecision::inverse_mod(wide_type(size_), wide_type(prime.modulus)) - .convert_to())); - - // At a layer with 2^d blocks, block j uses omega_(2^(d+1))^bitreverse(j,d). - // Each new table segment is the preceding prefix times a primitive root of twice - // the previous order: [1], [i], [w8,i*w8], [w16,i*w16,w8*w16,i*w8*w16], ... . - for (std::size_t half = size_ / 4; half != 0; half >>= 1) { - twiddles_[half] = centered(static_cast(root.convert_to())); - root = root * root % prime.modulus; - } - for (std::size_t half = 1; half < size_ / 2; half <<= 1) { - const double step = twiddles_[half]; - for (std::size_t i = 1; i < half; ++i) { - twiddles_[half + i] = centered(arithmetic_.multiply_integer(twiddles_[i], step)); - } - } - } - - std::size_t size() const { - return size_; - } - - std::uint64_t modulus() const { - return arithmetic_.modulus(); - } - - std::uint64_t omega() const { - return omega_; - } - - /** Exact centered representative of 1/N modulo q, also usable to fuse scaling into pointwise work. */ - double size_inverse() const { - return size_inverse_; - } - - /** Require exactly N entries; invalid dimensions throw before any writes. N=1 is the identity. */ - void fft(std::span values) const { - check_buffer(values); - if (arithmetic_.fma_supported()) { - forward_recursive(values.data(), size_, 0); - } else { - forward_recursive(values.data(), size_, 0); - } - } - - /** - * Require exactly N entries in the forward spectrum order. By default, scale by 1/N and - * return canonical residues. With normalize=false, return N times the inverse in (-2*q,2*q), - * allowing the caller to fuse the scaling elsewhere. Invalid dimensions throw before writes. - */ - void inverse_fft(std::span values, bool normalize = true) const { - check_buffer(values); - if (arithmetic_.fma_supported()) { - inverse_impl(values, normalize); - } else { - inverse_impl(values, normalize); - } - } - - /** - * Unscaled inverse for pointwise products whose integer representatives are in (-q,q). - * The caller has already folded 1/N into its pointwise scale. This stricter input bound - * permits omitting the first leaf reductions; the output still lies in (-2*q,2*q). - * Dimension errors throw before writes, as in inverse_fft(). - */ - void inverse_fft_pointwise(std::span values) const { - check_buffer<1>(values); - if (arithmetic_.fma_supported()) { - inverse_impl(values, false); - } else { - inverse_impl(values, false); - } - } - - private: - static const small_prime_fft_prime &checked_prime(std::size_t prime_index, std::size_t size) { - if (prime_index >= small_prime_fft_primes.size() || !std::has_single_bit(size) || - std::countr_zero(size) > small_prime_fft_primes[prime_index].two_adicity) { - throw std::invalid_argument("small_prime_fft_plan: invalid prime index or transform size"); - } - return small_prime_fft_primes[prime_index]; - } - - template - void check_buffer(std::span values) const { - if (values.size() != size_ || twiddles_.empty()) { - throw std::invalid_argument("small_prime_fft_plan: invalid buffer dimensions or moved-from plan"); - } -#ifndef NDEBUG - const double limit = Bound * static_cast(modulus()); - for (double value : values) { - assert(std::isfinite(value) && std::trunc(value) == value && std::abs(value) < limit); - } -#endif - } - - double centered(double value) const { - const double q = static_cast(modulus()); - return value > q / 2 ? value - q : (value < -q / 2 ? value + q : value); - } - - BOOST_FORCEINLINE double inverse_twiddle(std::size_t index) const { - if (index == 0) { - return 1.0; - } - // Within [b,2*b), reflected exponents add to half the root order, hence w^-1=-w_reflected. - const std::size_t base = std::bit_floor(index); - return -twiddles_[base + (2 * base - 1 - index)]; - } - - // Keep arithmetic selection outside the butterfly loops, with compile-time dispatch inside them. - template - BOOST_FORCEINLINE static Value reduce(const Arithmetic &arithmetic, Value value) { - if constexpr (UseFma) { - return arithmetic.reduce_fma(value); - } else { - return arithmetic.reduce_integer(value); - } - } - - template - BOOST_FORCEINLINE static Value multiply(const Arithmetic &arithmetic, Value left, Value right) { - if constexpr (UseFma) { - return arithmetic.multiply_fma(left, right); - } else { - return arithmetic.multiply_integer(left, right); - } - } - - template - BOOST_FORCEINLINE static void multiply_pair(const Arithmetic &arithmetic, Value &a, Value &b, Value wa, - Value wb) { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma && requires { arithmetic.multiply_fma_pair(a, b, wa, wb); }) { - arithmetic.multiply_fma_pair(a, b, wa, wb); - } else -#endif - { - a = multiply(arithmetic, a, wa); - b = multiply(arithmetic, b, wb); - } - } - - template - BOOST_FORCEINLINE static void forward_butterfly(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, - Value &v3, Value w2, Value w, Value iw) { -#if defined(__clang__) -#pragma clang fp reassociate(off) -#pragma clang fp contract(off) -#endif - const Value a0 = reduce(arithmetic, v0); - const Value a1 = v1; - Value a2 = v2; - Value a3 = v3; - if constexpr (UnitTwiddle) { - a2 = reduce(arithmetic, a2); - a3 = reduce(arithmetic, a3); - } else { - multiply_pair(arithmetic, a2, a3, w2, w2); - } - const Value even_sum = a0 + a2; - const Value even_difference = a0 - a2; - Value odd_sum = a1 + a3; - Value odd_difference = a1 - a3; - if constexpr (UnitTwiddle) { - odd_sum = reduce(arithmetic, odd_sum); - odd_difference = multiply(arithmetic, odd_difference, iw); - } else { - multiply_pair(arithmetic, odd_sum, odd_difference, w, iw); - } - // Input magnitudes < 3*q and weights < q/2 imply |a2|,|a3| < q and - // |a1 +/- a3| < 4*q. Every product is < 2*q^2, so the odd terms are < q. - // Thus all outputs are < 3*q, with every integer sum exact below 2^53. - v0 = even_sum + odd_sum; - v1 = even_sum - odd_sum; - v2 = even_difference + odd_difference; - v3 = even_difference - odd_difference; - } - -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - BOOST_FORCEINLINE static void transpose4(__m256d &v0, __m256d &v1, __m256d &v2, __m256d &v3) { - const __m256d a = _mm256_unpacklo_pd(v0, v1); - const __m256d b = _mm256_unpackhi_pd(v0, v1); - const __m256d c = _mm256_unpacklo_pd(v2, v3); - const __m256d d = _mm256_unpackhi_pd(v2, v3); - v0 = _mm256_permute2f128_pd(a, c, 0x20); - v1 = _mm256_permute2f128_pd(b, d, 0x20); - v2 = _mm256_permute2f128_pd(a, c, 0x31); - v3 = _mm256_permute2f128_pd(b, d, 0x31); - } - - // Aligned table indices, with unaligned loads: the owning vector needs no special allocator. - template - BOOST_FORCEINLINE __m256d twiddle4(std::size_t index) const { - if constexpr (!Inverse) { - return _mm256_loadu_pd(twiddles_.data() + index); - } else { - if (index == 0) { - return _mm256_setr_pd(1.0, inverse_twiddle(1), inverse_twiddle(2), inverse_twiddle(3)); - } - // Four consecutive entries stay in one [b,2*b) segment. Reflect, reverse, negate. - const std::size_t base = std::bit_floor(index); - const __m256d reflected = _mm256_loadu_pd(twiddles_.data() + 3 * base - index - 4); - return _mm256_sub_pd(_mm256_setzero_pd(), _mm256_permute4x64_pd(reflected, 0x1b)); - } - } - - template - BOOST_FORCEINLINE void leaf_twiddles(std::size_t node, __m256d &w2, __m256d &w, __m256d &iw) const { - w2 = twiddle4(4 * node); - const __m256d low = twiddle4(8 * node); - const __m256d high = twiddle4(8 * node + 4); - w = _mm256_permute4x64_pd(_mm256_unpacklo_pd(low, high), 0xd8); - iw = _mm256_permute4x64_pd(_mm256_unpackhi_pd(low, high), 0xd8); - } - - template - BOOST_FORCEINLINE void forward16(const small_prime_fft_arithmetic::packed_fma_kernel &packed, __m256d &v0, - __m256d &v1, __m256d &v2, __m256d &v3, std::size_t node) const { - forward_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(twiddles_[node]), - _mm256_set1_pd(twiddles_[2 * node]), - _mm256_set1_pd(twiddles_[2 * node + 1])); - // The four child transforms now run across SIMD lanes, without scalar leaves or reloads. - transpose4(v0, v1, v2, v3); - __m256d w2, w, iw; - leaf_twiddles(node, w2, w, iw); - forward_butterfly(packed, v0, v1, v2, v3, w2, w, iw); - transpose4(v0, v1, v2, v3); - } - - template - BOOST_FORCEINLINE void inverse16(const small_prime_fft_arithmetic::packed_fma_kernel &packed, __m256d &v0, - __m256d &v1, __m256d &v2, __m256d &v3, std::size_t node) const { - transpose4(v0, v1, v2, v3); - __m256d w2, w, iw; - leaf_twiddles(node, w2, w, iw); - inverse_butterfly(packed, v0, v1, v2, v3, w2, w, iw, !BoundedInput); - transpose4(v0, v1, v2, v3); - inverse_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(inverse_twiddle(node)), - _mm256_set1_pd(inverse_twiddle(2 * node)), - _mm256_set1_pd(inverse_twiddle(2 * node + 1)), false); - } - - template - void base16(double *values, std::size_t node) const { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - __m256d v0 = _mm256_loadu_pd(values); - __m256d v1 = _mm256_loadu_pd(values + 4); - __m256d v2 = _mm256_loadu_pd(values + 8); - __m256d v3 = _mm256_loadu_pd(values + 12); - if constexpr (Inverse) { - inverse16(packed, v0, v1, v2, v3, node); - } else { - forward16(packed, v0, v1, v2, v3, node); - } - _mm256_storeu_pd(values, v0); - _mm256_storeu_pd(values + 4, v1); - _mm256_storeu_pd(values + 8, v2); - _mm256_storeu_pd(values + 12, v3); - } - - template - BOOST_FORCEINLINE static void radix2_packed(const Arithmetic &packed, Value &left, Value &right, Value weight) { - if constexpr (Inverse) { - // Both subtransforms already return magnitudes < 2*q. No input reductions are needed. - const Value sum = left + right; - const Value difference = left - right; - left = packed.reduce_fma(sum); - right = UnitTwiddle ? packed.reduce_fma(difference) : packed.multiply_fma(difference, weight); - } else { - const Value a = packed.reduce_fma(left); - const Value b = UnitTwiddle ? packed.reduce_fma(right) : packed.multiply_fma(right, weight); - left = a + b; - right = a - b; - } - } - - template - BOOST_FORCEINLINE static void stage4(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, Value &v3, - const double *weights) { - const Value w2 = Arithmetic::broadcast(weights[0]); - const Value w = Arithmetic::broadcast(weights[1]); - const Value iw = Arithmetic::broadcast(weights[2]); - if constexpr (Inverse) { - inverse_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw, false); - } else { - forward_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw); - } - } - - template - BOOST_FORCEINLINE void stage_weights(double *weights, std::size_t node) const { - weights[0] = Inverse ? inverse_twiddle(node) : twiddles_[node]; - weights[1] = Inverse ? inverse_twiddle(2 * node) : twiddles_[2 * node]; - weights[2] = Inverse ? inverse_twiddle(2 * node + 1) : twiddles_[2 * node + 1]; - } - - template - BOOST_FORCEINLINE void stage4(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, Value &v3, - std::size_t node) const { - double weights[3]; - stage_weights(weights, node); - stage4(arithmetic, v0, v1, v2, v3, weights); - } - - template - BOOST_FORCEINLINE static void fused_outer(const Arithmetic &arithmetic, Value *v, const double *weights) { - if constexpr (Radix == 8) { - const Value weight = Arithmetic::broadcast(weights[0]); - radix2_packed(arithmetic, v[0], v[4], weight); - radix2_packed(arithmetic, v[1], v[5], weight); - radix2_packed(arithmetic, v[2], v[6], weight); - radix2_packed(arithmetic, v[3], v[7], weight); - } else { - stage4(arithmetic, v[0], v[4], v[8], v[12], weights); - stage4(arithmetic, v[1], v[5], v[9], v[13], weights); - stage4(arithmetic, v[2], v[6], v[10], v[14], weights); - stage4(arithmetic, v[3], v[7], v[11], v[15], weights); - } - } - - // Fuse three/four binary layers without intermediate array passes. Every constituent - // butterfly retains the forward <3*q and inverse <2*q bounds; no wider lazy range is used. - // UnitTwiddle is valid only at node zero, hence only the first child inherits it. - template - BOOST_FORCEINLINE static void fused_stage(const Arithmetic &arithmetic, Value *v, const double *weights) { - static_assert(Radix == 8 || Radix == 16); - if constexpr (!Inverse) { - fused_outer(arithmetic, v, weights); - } - constexpr std::size_t outer_weights = Radix == 8 ? 1 : 3; - stage4(arithmetic, v[0], v[1], v[2], v[3], weights + outer_weights); - stage4(arithmetic, v[4], v[5], v[6], v[7], weights + outer_weights + 3); - if constexpr (Radix == 16) { - stage4(arithmetic, v[8], v[9], v[10], v[11], weights + outer_weights + 6); - stage4(arithmetic, v[12], v[13], v[14], v[15], weights + outer_weights + 9); - } - if constexpr (Inverse) { - fused_outer(arithmetic, v, weights); - } - } - - template - void large_radix(double *values, std::size_t stride, std::size_t node) const { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - using arithmetic_type = small_prime_fft_arithmetic::wide_fma_kernel; -#else - using arithmetic_type = small_prime_fft_arithmetic::packed_fma_kernel; -#endif - const arithmetic_type arithmetic(arithmetic_); - // Dispatch requires stride >=16, a power of two divisible by the SIMD width. - assert(stride >= 16 && stride % arithmetic_type::lanes == 0); - // The node is fixed throughout this coefficient loop. Resolve reflections and - // table lookups once, retaining scalar weights to avoid forcing Radix-1 live vectors. - double weights[Radix - 1]; - if constexpr (Radix == 8) { - weights[0] = Inverse ? inverse_twiddle(node) : twiddles_[node]; - } else { - stage_weights(weights, node); - } - constexpr std::size_t outer_weights = Radix == 8 ? 1 : 3; -#pragma GCC unroll 4 - for (std::size_t child = 0; child < Radix / 4; ++child) { - stage_weights(weights + outer_weights + 3 * child, Radix / 4 * node + child); - } - for (std::size_t i = 0; i < stride; i += arithmetic_type::lanes) { - arithmetic_type::value_type v[Radix]; -#pragma GCC unroll 16 - for (std::size_t part = 0; part < Radix; ++part) { - v[part] = arithmetic_type::load(values + part * stride + i); - } - fused_stage(arithmetic, v, weights); -#pragma GCC unroll 16 - for (std::size_t part = 0; part < Radix; ++part) { - arithmetic_type::store(values + part * stride + i, v[part]); - } - } - } - - template - void base32(double *values, std::size_t node) const { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - __m256d v0 = _mm256_loadu_pd(values); - __m256d v1 = _mm256_loadu_pd(values + 4); - __m256d v2 = _mm256_loadu_pd(values + 8); - __m256d v3 = _mm256_loadu_pd(values + 12); - __m256d v4 = _mm256_loadu_pd(values + 16); - __m256d v5 = _mm256_loadu_pd(values + 20); - __m256d v6 = _mm256_loadu_pd(values + 24); - __m256d v7 = _mm256_loadu_pd(values + 28); - const __m256d weight = _mm256_set1_pd(Inverse ? inverse_twiddle(node) : twiddles_[node]); - // One radix-2 layer and two register-resident 16-point transforms cover all five layers. - if constexpr (Inverse) { - inverse16(packed, v0, v1, v2, v3, 2 * node); - inverse16(packed, v4, v5, v6, v7, 2 * node + 1); - } - radix2_packed(packed, v0, v4, weight); - radix2_packed(packed, v1, v5, weight); - radix2_packed(packed, v2, v6, weight); - radix2_packed(packed, v3, v7, weight); - if constexpr (!Inverse) { - forward16(packed, v0, v1, v2, v3, 2 * node); - forward16(packed, v4, v5, v6, v7, 2 * node + 1); - } - _mm256_storeu_pd(values, v0); - _mm256_storeu_pd(values + 4, v1); - _mm256_storeu_pd(values + 8, v2); - _mm256_storeu_pd(values + 12, v3); - _mm256_storeu_pd(values + 16, v4); - _mm256_storeu_pd(values + 20, v5); - _mm256_storeu_pd(values + 24, v6); - _mm256_storeu_pd(values + 28, v7); - } - - template - void base8(double *values, std::size_t node) const { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - const __m256d low = _mm256_loadu_pd(values); - const __m256d high = _mm256_loadu_pd(values + 4); - __m256d v0, v1, v2, v3; - if constexpr (Inverse) { - __m256d a = _mm256_permute4x64_pd(_mm256_unpacklo_pd(low, high), 0xd8); - __m256d b = _mm256_permute4x64_pd(_mm256_unpackhi_pd(low, high), 0xd8); - if constexpr (!BoundedInput) { - a = packed.reduce_fma(a); - b = packed.reduce_fma(b); - } - const __m256d sum = packed.reduce_fma(a + b); - const __m256d difference = packed.multiply_fma(a - b, twiddle4(4 * node)); - const __m256d even = _mm256_unpacklo_pd(sum, difference); - const __m256d odd = _mm256_unpackhi_pd(sum, difference); - v0 = _mm256_permute2f128_pd(even, even, 0x00); - v1 = _mm256_permute2f128_pd(odd, odd, 0x00); - v2 = _mm256_permute2f128_pd(even, even, 0x11); - v3 = _mm256_permute2f128_pd(odd, odd, 0x11); - inverse_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(inverse_twiddle(node)), - _mm256_set1_pd(inverse_twiddle(2 * node)), - _mm256_set1_pd(inverse_twiddle(2 * node + 1)), false); - _mm256_storeu_pd(values, _mm256_permute2f128_pd(v0, v1, 0x20)); - _mm256_storeu_pd(values + 4, _mm256_permute2f128_pd(v2, v3, 0x20)); - } else { - // Two independent columns occupy the lower lanes; the duplicated upper lanes are discarded. - v0 = _mm256_permute2f128_pd(low, low, 0x00); - v1 = _mm256_permute2f128_pd(low, low, 0x11); - v2 = _mm256_permute2f128_pd(high, high, 0x00); - v3 = _mm256_permute2f128_pd(high, high, 0x11); - forward_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(twiddles_[node]), - _mm256_set1_pd(twiddles_[2 * node]), - _mm256_set1_pd(twiddles_[2 * node + 1])); - const __m256d a = packed.reduce_fma( - _mm256_permute2f128_pd(_mm256_unpacklo_pd(v0, v1), _mm256_unpacklo_pd(v2, v3), 0x20)); - const __m256d b = packed.multiply_fma( - _mm256_permute2f128_pd(_mm256_unpackhi_pd(v0, v1), _mm256_unpackhi_pd(v2, v3), 0x20), - twiddle4(4 * node)); - const __m256d even = _mm256_unpacklo_pd(a + b, a - b); - const __m256d odd = _mm256_unpackhi_pd(a + b, a - b); - _mm256_storeu_pd(values, _mm256_permute2f128_pd(even, odd, 0x20)); - _mm256_storeu_pd(values + 4, _mm256_permute2f128_pd(even, odd, 0x31)); - } - } - - template - void basecase(double *values, std::size_t length, std::size_t node) const { - if (length == 16) { - base16(values, node); - } else if (length == 32) { - base32(values, node); - } else { - base8(values, node); - } - } - -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - BOOST_FORCEINLINE static __m512d join4(__m256d low, __m256d high) { - return _mm512_insertf64x4(_mm512_castpd256_pd512(low), high, 1); - } - - // Two independent 4x4 transposes, one in each 256-bit half. The stored spectrum - // remains bit-reversed; these permutations only change the lanes during a leaf. - BOOST_FORCEINLINE static void transpose4(__m512d &v0, __m512d &v1, __m512d &v2, __m512d &v3) { - const __m512d a = _mm512_unpacklo_pd(v0, v1); - const __m512d b = _mm512_unpackhi_pd(v0, v1); - const __m512d c = _mm512_unpacklo_pd(v2, v3); - const __m512d d = _mm512_unpackhi_pd(v2, v3); - const __m512i low = _mm512_setr_epi64(0, 1, 8, 9, 4, 5, 12, 13); - const __m512i high = _mm512_setr_epi64(2, 3, 10, 11, 6, 7, 14, 15); - v0 = _mm512_permutex2var_pd(a, low, c); - v1 = _mm512_permutex2var_pd(b, low, d); - v2 = _mm512_permutex2var_pd(a, high, c); - v3 = _mm512_permutex2var_pd(b, high, d); - } - - template - BOOST_FORCEINLINE __m512d paired_twiddle(std::size_t low, std::size_t high) const { - return join4(_mm256_set1_pd(Inverse ? inverse_twiddle(low) : twiddles_[low]), - _mm256_set1_pd(Inverse ? inverse_twiddle(high) : twiddles_[high])); - } - - // Keep the complete leaf in one function: outlining this step forces live vectors - // through memory across the call. Recursive transform boundaries remain uninlined. - template - BOOST_FORCEINLINE void paired_outer(const small_prime_fft_arithmetic::wide_fma_kernel &arithmetic, __m512d &v0, - __m512d &v1, __m512d &v2, __m512d &v3, std::size_t node) const { - const __m512d w2 = paired_twiddle(node, node + 1); - const __m512d w = paired_twiddle(2 * node, 2 * node + 2); - const __m512d iw = paired_twiddle(2 * node + 1, 2 * node + 3); - if constexpr (Inverse) { - inverse_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw, false); - } else { - forward_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw); - } - } - - template - BOOST_FORCEINLINE void base16_pair(const small_prime_fft_arithmetic::wide_fma_kernel &arithmetic, __m512d &v0, - __m512d &v1, __m512d &v2, __m512d &v3, std::size_t node) const { - // Group matching rows of two adjacent 16-point transforms in the two vector halves. - const __m512d a = _mm512_shuffle_f64x2(v0, v2, 0x44); - const __m512d b = _mm512_shuffle_f64x2(v0, v2, 0xee); - const __m512d c = _mm512_shuffle_f64x2(v1, v3, 0x44); - const __m512d d = _mm512_shuffle_f64x2(v1, v3, 0xee); - v0 = a; - v1 = b; - v2 = c; - v3 = d; - if constexpr (!Inverse) { - paired_outer(arithmetic, v0, v1, v2, v3, node); - } - transpose4(v0, v1, v2, v3); - __m256d lo_w2, lo_w, lo_iw, hi_w2, hi_w, hi_iw; - leaf_twiddles(node, lo_w2, lo_w, lo_iw); - leaf_twiddles(node + 1, hi_w2, hi_w, hi_iw); - if constexpr (Inverse) { - inverse_butterfly(arithmetic, v0, v1, v2, v3, join4(lo_w2, hi_w2), join4(lo_w, hi_w), - join4(lo_iw, hi_iw), !BoundedInput); - } else { - forward_butterfly(arithmetic, v0, v1, v2, v3, join4(lo_w2, hi_w2), join4(lo_w, hi_w), - join4(lo_iw, hi_iw)); - } - transpose4(v0, v1, v2, v3); - if constexpr (Inverse) { - paired_outer(arithmetic, v0, v1, v2, v3, node); - } - // Restore contiguous bit-reversed order in the same four input/output vectors. - const __m512d out0 = _mm512_shuffle_f64x2(v0, v1, 0x44); - const __m512d out1 = _mm512_shuffle_f64x2(v2, v3, 0x44); - const __m512d out2 = _mm512_shuffle_f64x2(v0, v1, 0xee); - const __m512d out3 = _mm512_shuffle_f64x2(v2, v3, 0xee); - v0 = out0; - v1 = out1; - v2 = out2; - v3 = out3; - } - - template - void base32_wide(double *values, std::size_t node) const { - const small_prime_fft_arithmetic::wide_fma_kernel arithmetic(arithmetic_); - __m512d v0 = _mm512_loadu_pd(values); - __m512d v1 = _mm512_loadu_pd(values + 8); - __m512d v2 = _mm512_loadu_pd(values + 16); - __m512d v3 = _mm512_loadu_pd(values + 24); - const __m512d weight = _mm512_set1_pd(Inverse ? inverse_twiddle(node) : twiddles_[node]); - if constexpr (Inverse) { - base16_pair(arithmetic, v0, v1, v2, v3, 2 * node); - } - radix2_packed(arithmetic, v0, v2, weight); - radix2_packed(arithmetic, v1, v3, weight); - if constexpr (!Inverse) { - base16_pair(arithmetic, v0, v1, v2, v3, 2 * node); - } - _mm512_storeu_pd(values, v0); - _mm512_storeu_pd(values + 8, v1); - _mm512_storeu_pd(values + 16, v2); - _mm512_storeu_pd(values + 24, v3); - } - - template - void base64_wide(double *values, std::size_t node) const { - const small_prime_fft_arithmetic::wide_fma_kernel arithmetic(arithmetic_); - __m512d v[8]; -#pragma GCC unroll 8 - for (std::size_t i = 0; i < 8; ++i) { - v[i] = _mm512_loadu_pd(values + 8 * i); - } - if constexpr (!Inverse) { - stage4(arithmetic, v[0], v[2], v[4], v[6], node); - stage4(arithmetic, v[1], v[3], v[5], v[7], node); - } - base16_pair(arithmetic, v[0], v[1], v[2], v[3], 4 * node); - base16_pair(arithmetic, v[4], v[5], v[6], v[7], 4 * node + 2); - if constexpr (Inverse) { - stage4(arithmetic, v[0], v[2], v[4], v[6], node); - stage4(arithmetic, v[1], v[3], v[5], v[7], node); - } -#pragma GCC unroll 8 - for (std::size_t i = 0; i < 8; ++i) { - _mm512_storeu_pd(values + 8 * i, v[i]); - } - } - - template - void wide_basecase(double *values, std::size_t length, std::size_t node) const { - if (length == 32) { - if (node == 0) { - base32_wide(values, node); - } else { - base32_wide(values, node); - } - } else { - if (node == 0) { - base64_wide(values, node); - } else { - base64_wide(values, node); - } - } - } -#endif -#endif - - template - void forward_radix4(double *values, std::size_t quarter, std::size_t node, std::size_t first, - std::size_t last) const { - const double w2 = twiddles_[node]; - const double w = twiddles_[2 * node]; - const double iw = twiddles_[2 * node + 1]; - std::size_t i = first; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - if constexpr (UseFma) { - if (quarter >= 8) { - const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); - const __m512d packed_w2 = _mm512_set1_pd(w2); - const __m512d packed_w = _mm512_set1_pd(w); - const __m512d packed_iw = _mm512_set1_pd(iw); -#pragma GCC unroll 2 - for (; last - i >= 8; i += 8) { - __m512d v0 = _mm512_loadu_pd(values + i); - __m512d v1 = _mm512_loadu_pd(values + i + quarter); - __m512d v2 = _mm512_loadu_pd(values + i + 2 * quarter); - __m512d v3 = _mm512_loadu_pd(values + i + 3 * quarter); - forward_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw); - _mm512_storeu_pd(values + i, v0); - _mm512_storeu_pd(values + i + quarter, v1); - _mm512_storeu_pd(values + i + 2 * quarter, v2); - _mm512_storeu_pd(values + i + 3 * quarter, v3); - } - } - } -#endif -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { - if (quarter >= 4) { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - const __m256d packed_w2 = _mm256_set1_pd(w2); - const __m256d packed_w = _mm256_set1_pd(w); - const __m256d packed_iw = _mm256_set1_pd(iw); -#pragma GCC unroll 2 - for (; last - i >= 4; i += 4) { - __m256d v0 = _mm256_loadu_pd(values + i); - __m256d v1 = _mm256_loadu_pd(values + i + quarter); - __m256d v2 = _mm256_loadu_pd(values + i + 2 * quarter); - __m256d v3 = _mm256_loadu_pd(values + i + 3 * quarter); - forward_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw); - _mm256_storeu_pd(values + i, v0); - _mm256_storeu_pd(values + i + quarter, v1); - _mm256_storeu_pd(values + i + 2 * quarter, v2); - _mm256_storeu_pd(values + i + 3 * quarter, v3); - } - } - } -#endif - for (; i < last; ++i) { - forward_butterfly(arithmetic_, values[i], values[i + quarter], - values[i + 2 * quarter], values[i + 3 * quarter], w2, w, iw); - } - } - -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - template - void four_layer_block(double *values, std::size_t length, std::size_t node) const { - const std::size_t quarter = length / 4; - const std::size_t sixteenth = length / 16; - // Sixteen strips of 64 doubles occupy 8 KiB. Complete both radix-4 passes on - // each strip group while it is hot, before recursing into its sixteen children. - // This only changes visitation order; every butterfly retains the same inputs and bounds. - constexpr std::size_t strip_size = 64; - const auto pass = [&](double *base, std::size_t stride, std::size_t j, std::size_t first, - std::size_t last) { - if constexpr (Inverse) { - if (j == 0) { - inverse_radix4(base, stride, j, first, last); - } else { - inverse_radix4(base, stride, j, first, last); - } - } else { - if (j == 0) { - forward_radix4(base, stride, j, first, last); - } else { - forward_radix4(base, stride, j, first, last); - } - } - }; - for (std::size_t first = 0; first < sixteenth; first += strip_size) { - const std::size_t last = std::min(first + strip_size, sixteenth); - if constexpr (Inverse) { - for (std::size_t child = 0; child < 4; ++child) { - pass(values + child * quarter, sixteenth, 4 * node + child, first, last); - } - } - for (std::size_t part = 0; part < 4; ++part) { - pass(values, quarter, node, first + part * sixteenth, last + part * sixteenth); - } - if constexpr (!Inverse) { - for (std::size_t child = 0; child < 4; ++child) { - pass(values + child * quarter, sixteenth, 4 * node + child, first, last); - } - } - } - } -#endif - - template - void forward_recursive(double *values, std::size_t length, std::size_t node) const { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - if (length == 32 || length == 64) { - wide_basecase(values, length, node); - return; - } -#endif - if (length >= 8 && length <= 32) { - if (node == 0) { - basecase(values, length, node); - } else { - basecase(values, length, node); - } - return; - } -#if CRYPTO3_SMALL_PRIME_FFT_RADIX != 4 - constexpr std::size_t radix = CRYPTO3_SMALL_PRIME_FFT_RADIX; - if (length >= 16 * radix) { - const auto stride = length / radix; - if (node == 0) { - large_radix(values, stride, node); - } else { - large_radix(values, stride, node); - } - for (std::size_t child = 0; child < radix; ++child) { - forward_recursive(values + child * stride, stride, radix * node + child); - } - return; - } -#endif - if (length >= 4096) { - four_layer_block(values, length, node); - for (std::size_t child = 0; child < 16; ++child) { - forward_recursive(values + child * (length / 16), length / 16, 16 * node + child); - } - return; - } - } -#endif - if (length == 1) { - return; - } - if (length == 2) { - const double left = reduce(arithmetic_, values[0]); - const double right = node == 0 ? reduce(arithmetic_, values[1]) : - multiply(arithmetic_, values[1], twiddles_[node]); - values[0] = left + right; - values[1] = left - right; - return; - } - const std::size_t quarter = length / 4; - if (node == 0) { - forward_radix4(values, quarter, node, 0, quarter); - } else { - forward_radix4(values, quarter, node, 0, quarter); - } - if (quarter > 1) { - for (std::size_t child = 0; child < 4; ++child) { - forward_recursive(values + child * quarter, quarter, 4 * node + child); - } - } - } - - template - BOOST_FORCEINLINE static void inverse_butterfly(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, - Value &v3, Value w2, Value w, Value iw, bool leaf) { -#if defined(__clang__) -#pragma clang fp reassociate(off) -#pragma clang fp contract(off) -#endif - // Leaf inputs may be forward residues of magnitude < 3*q. Reduce them before - // adding four entries; non-leaf inputs already satisfy the inverse bound < 2*q. - const Value b0 = leaf ? reduce(arithmetic, v0) : v0; - const Value b1 = leaf ? reduce(arithmetic, v1) : v1; - const Value b2 = leaf ? reduce(arithmetic, v2) : v2; - const Value b3 = leaf ? reduce(arithmetic, v3) : v3; - const Value even_sum = b0 + b1; - const Value odd_sum = b2 + b3; - Value left = b0 - b1; - Value right = b2 - b3; - if constexpr (UnitTwiddle) { - left = reduce(arithmetic, left); - right = multiply(arithmetic, right, iw); - } else { - multiply_pair(arithmetic, left, right, w, iw); - } - // Four inputs < 2*q give sums < 8*q < 2^53. Times a centered weight this is - // < 4*q^2, so the a2 remainder is < 1.5*q. Left and right are each < q; - // their sum is < 2*q and their weighted difference is < q. The invariant closes. - v0 = reduce(arithmetic, even_sum + odd_sum); - v1 = left + right; - v2 = even_sum - odd_sum; - v3 = left - right; - if constexpr (UnitTwiddle) { - v2 = reduce(arithmetic, v2); - v3 = reduce(arithmetic, v3); - } else { - multiply_pair(arithmetic, v2, v3, w2, w2); - } - } - - template - void inverse_radix4(double *values, std::size_t quarter, std::size_t node, std::size_t first, - std::size_t last) const { - const double w2 = inverse_twiddle(node); - const double w = inverse_twiddle(2 * node); - const double iw = inverse_twiddle(2 * node + 1); - std::size_t i = first; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - if constexpr (UseFma) { - if (quarter >= 8) { - const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); - const __m512d packed_w2 = _mm512_set1_pd(w2); - const __m512d packed_w = _mm512_set1_pd(w); - const __m512d packed_iw = _mm512_set1_pd(iw); -#pragma GCC unroll 2 - for (; last - i >= 8; i += 8) { - __m512d v0 = _mm512_loadu_pd(values + i); - __m512d v1 = _mm512_loadu_pd(values + i + quarter); - __m512d v2 = _mm512_loadu_pd(values + i + 2 * quarter); - __m512d v3 = _mm512_loadu_pd(values + i + 3 * quarter); - inverse_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw, - false); - _mm512_storeu_pd(values + i, v0); - _mm512_storeu_pd(values + i + quarter, v1); - _mm512_storeu_pd(values + i + 2 * quarter, v2); - _mm512_storeu_pd(values + i + 3 * quarter, v3); - } - } - } -#endif -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { - if (quarter >= 4) { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - const __m256d packed_w2 = _mm256_set1_pd(w2); - const __m256d packed_w = _mm256_set1_pd(w); - const __m256d packed_iw = _mm256_set1_pd(iw); -#pragma GCC unroll 2 - for (; last - i >= 4; i += 4) { - __m256d v0 = _mm256_loadu_pd(values + i); - __m256d v1 = _mm256_loadu_pd(values + i + quarter); - __m256d v2 = _mm256_loadu_pd(values + i + 2 * quarter); - __m256d v3 = _mm256_loadu_pd(values + i + 3 * quarter); - inverse_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw, - false); - _mm256_storeu_pd(values + i, v0); - _mm256_storeu_pd(values + i + quarter, v1); - _mm256_storeu_pd(values + i + 2 * quarter, v2); - _mm256_storeu_pd(values + i + 3 * quarter, v3); - } - } - } -#endif - for (; i < last; ++i) { - inverse_butterfly(arithmetic_, values[i], values[i + quarter], - values[i + 2 * quarter], values[i + 3 * quarter], w2, w, iw, - quarter == 1 && !BoundedInput); - } - } - - template - void inverse_recursive(double *values, std::size_t length, std::size_t node) const { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - if (length == 32 || length == 64) { - wide_basecase(values, length, node); - return; - } -#endif - if (length >= 8 && length <= 32) { - if (node == 0) { - basecase(values, length, node); - } else { - basecase(values, length, node); - } - return; - } -#if CRYPTO3_SMALL_PRIME_FFT_RADIX != 4 - constexpr std::size_t radix = CRYPTO3_SMALL_PRIME_FFT_RADIX; - if (length >= 16 * radix) { - const auto stride = length / radix; - for (std::size_t child = 0; child < radix; ++child) { - inverse_recursive(values + child * stride, stride, radix * node + child); - } - if (node == 0) { - large_radix(values, stride, node); - } else { - large_radix(values, stride, node); - } - return; - } -#endif - if (length >= 4096) { - for (std::size_t child = 0; child < 16; ++child) { - inverse_recursive(values + child * (length / 16), length / 16, - 16 * node + child); - } - four_layer_block(values, length, node); - return; - } - } -#endif - if (length == 1) { - if constexpr (!BoundedInput) { - values[0] = reduce(arithmetic_, values[0]); - } - return; - } - if (length == 2) { - const double left = BoundedInput ? values[0] : reduce(arithmetic_, values[0]); - const double right = BoundedInput ? values[1] : reduce(arithmetic_, values[1]); - values[0] = reduce(arithmetic_, left + right); - values[1] = node == 0 ? reduce(arithmetic_, left - right) : - multiply(arithmetic_, left - right, inverse_twiddle(node)); - return; - } - const std::size_t quarter = length / 4; - if (quarter > 1) { - for (std::size_t child = 0; child < 4; ++child) { - inverse_recursive(values + child * quarter, quarter, 4 * node + child); - } - } - if (node == 0) { - inverse_radix4(values, quarter, node, 0, quarter); - } else { - inverse_radix4(values, quarter, node, 0, quarter); - } - } - - template - void inverse_impl(std::span values, bool normalize) const { - inverse_recursive(values.data(), size_, 0); - if (normalize) { - const double q = static_cast(modulus()); - std::size_t i = 0; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - if constexpr (UseFma) { - const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); - const __m512d scale = _mm512_set1_pd(size_inverse_); - const __m512d packed_q = _mm512_set1_pd(q); - for (; values.size() - i >= 8; i += 8) { - __m512d value = packed.multiply_fma(_mm512_loadu_pd(values.data() + i), scale); - const auto negative = _mm512_cmp_pd_mask(value, _mm512_setzero_pd(), _CMP_LT_OQ); - value = _mm512_mask_add_pd(value, negative, value, packed_q); - _mm512_storeu_pd(values.data() + i, value); - } - } -#endif -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { - if (values.size() >= 4) { - const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); - const __m256d scale = _mm256_set1_pd(size_inverse_); - const __m256d packed_q = _mm256_set1_pd(q); - for (; values.size() - i >= 4; i += 4) { - __m256d value = packed.multiply_fma(_mm256_loadu_pd(values.data() + i), scale); - const __m256d negative = _mm256_cmp_pd(value, _mm256_setzero_pd(), _CMP_LT_OQ); - value = _mm256_add_pd(value, _mm256_and_pd(negative, packed_q)); - _mm256_storeu_pd(values.data() + i, value); - } - } - } -#endif - for (; i < values.size(); ++i) { - double &value = values[i]; - value = multiply(arithmetic_, value, size_inverse_); - value = value < 0 ? value + q : value; - } - } - } - - small_prime_fft_arithmetic arithmetic_; - std::size_t size_; - std::uint64_t omega_; - double size_inverse_; - std::vector twiddles_; - }; - - /** - * Pointwise products in the BN254 tower u^2=-1, v^3=9+u, w^2=v, over one small modulus q. - * This is a quotient ring; irreducibility of the tower polynomials modulo q is not required. - * Coordinate channels have the same order as the converter, with one prime per call. - * - * Reuse the field library's Karatsuba tower identities and quadratic-extension squaring - * identity with scalar, four-lane or eight-lane residue kernels. The existing field elements - * bind their base arithmetic to the target modulus and cannot represent these small-prime intermediates. - * - * @pre Input cells are integer-valued doubles of magnitude < 2^53. scale is an integer with - * |scale| < q/2, normally the FFT plan's centered 1/N. Buffers are disjoint or identical. - * The thread's rounding mode stays unchanged during the call. Output cells are in (-q,q). - * No buffers are allocated and no constant-time guarantee is provided. - */ -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) -// The packed kernels use vector coordinates in the existing tower arrays. -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wignored-attributes" -#endif - class bn254_small_prime_fft_pointwise { - - template - struct kernel { - using quadratic_type = std::array; - using cubic_type = std::array; - using tower_type = std::array; - - const Arithmetic &arithmetic; - - Value reduce(Value value) const { - if constexpr (UseFma) { - return arithmetic.reduce_fma(value); - } else { - return arithmetic.reduce_integer(value); - } - } - - Value multiply(Value a, Value b) const { - if constexpr (UseFma) { - return arithmetic.multiply_fma(a, b); - } else { - return arithmetic.multiply_integer(a, b); - } - } - - BOOST_FORCEINLINE void multiply_pair(Value &a, Value &b, Value wa, Value wb) const { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma && requires { arithmetic.multiply_fma_pair(a, b, wa, wb); }) { - arithmetic.multiply_fma_pair(a, b, wa, wb); - } else -#endif - { - a = multiply(a, wa); - b = multiply(b, wb); - } - } - - template - Value add(Value a, Value b) const { - if constexpr (Reduce) { - return reduce(a + b); - } else { - return a + b; - } - } - - // Three reduced coordinates give an exact integer of magnitude < 3*q. - // Callers can defer this reduction to the next bounded sum or final scaling. - template - Value subtract(Value a, Value b, Value c) const { - if constexpr (Reduce) { - return reduce(a - b - c); - } else { - return a - b - c; - } - } - - template - std::array add(const std::array &a, const std::array &b) const { - std::array result; - for (std::size_t i = 0; i < N; ++i) { - result[i] = add(a[i], b[i]); - } - return result; - } - - template - std::array subtract(const std::array &a, const std::array &b, - const std::array &c) const { - std::array result; - for (std::size_t i = 0; i < N; ++i) { - result[i] = subtract(a[i], b[i], c[i]); - } - return result; - } - - // Quadratic inputs and returned coordinates are in (-q,q). Pair sums are < 2*q; - // the Karatsuba product is < 4*q^2 and its remainder < 1.5*q. Correcting by the - // two diagonal products gives a sum < 3.5*q, hence exact in binary64 before reduction. - BOOST_FORCEINLINE quadratic_type multiply(const quadratic_type &a, const quadratic_type &b) const { -#if defined(__clang__) -#pragma clang fp reassociate(off) -#pragma clang fp contract(off) -#endif - // Keep the two independent products in registers and interleave their FMA chains. - Value diagonal0 = a[0]; - Value diagonal1 = a[1]; - multiply_pair(diagonal0, diagonal1, b[0], b[1]); - const Value cross = multiply(a[0] + a[1], b[0] + b[1]); - return {reduce(diagonal0 - diagonal1), reduce(cross - diagonal0 - diagonal1)}; - } - - quadratic_type multiply_xi(const quadratic_type &a) const { -#if defined(__clang__) -#pragma clang fp reassociate(off) -#pragma clang fp contract(off) -#endif - // (9+u)*(a0+a1*u) = (9*a0-a1) + (a0+9*a1)*u. Reduce 8*a first: - // |8*a| < 8*q < 2^53, whereas forming 9*a directly could lose an integer bit. - return {reduce(reduce(8 * a[0]) + a[0] - a[1]), reduce(reduce(8 * a[1]) + a[1] + a[0])}; - } - - quadratic_type add_multiply_xi(const quadratic_type &a, const quadratic_type &b) const { -#if defined(__clang__) -#pragma clang fp reassociate(off) -#pragma clang fp contract(off) -#endif - // a + (9+u)*b, with |a| < 3*q and |b| < q. After reducing 8*b, - // each four-term sum is < 6*q < 2^53, so only its final reduction is needed. - return {reduce(reduce(8 * b[0]) + b[0] - b[1] + a[0]), reduce(reduce(8 * b[1]) + b[1] + b[0] + a[1])}; - } - - cubic_type multiply_v(const cubic_type &a) const { - return {multiply_xi(a[2]), a[0], a[1]}; - } - - cubic_type multiply(const cubic_type &a, const cubic_type &b) const { - const auto diagonal0 = multiply(a[0], b[0]); - const auto diagonal1 = multiply(a[1], b[1]); - const auto diagonal2 = multiply(a[2], b[2]); - // cross12 is reduced before multiplication by 8 inside multiply_xi. - // The other cross terms stay < 3*q until these final sums; the cubic result is < q. - // Finish each output pair before forming the next cross term to limit live temporaries. - cubic_type result; - { - const auto cross12 = subtract(multiply(add(a[1], a[2]), add(b[1], b[2])), diagonal1, diagonal2); - result[0] = add_multiply_xi(diagonal0, cross12); - } - { - const auto cross01 = - subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0, diagonal1); - result[1] = add_multiply_xi(cross01, diagonal2); - } - { - const auto cross02 = - subtract(multiply(add(a[0], a[2]), add(b[0], b[2])), diagonal0, diagonal2); - result[2] = add(cross02, diagonal1); - } - return result; - } - - // Both tower operations return coordinates of magnitude < 3*q. The final - // |scale| < q/2 gives a product < 1.5*q^2, whose FMA remainder is already in (-q,q). - // Keeping these sums unreduced is also safe when two tower products are added below. - tower_type multiply(const tower_type &a, const tower_type &b) const { - const auto diagonal0 = multiply(a[0], b[0]); - const auto diagonal1 = multiply(a[1], b[1]); - const auto cross = subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0, diagonal1); - return {add(diagonal0, multiply_v(diagonal1)), cross}; - } - - tower_type square(const tower_type &a) const { - // The quadratic-extension square needs two cubic products instead of three: - // c0 = (a0+a1)*(a0+v*a1) - a0*a1 - v*a0*a1; c1 = 2*a0*a1. - const auto cross = multiply(a[0], a[1]); - const auto combined = multiply(add(a[0], a[1]), add(a[0], multiply_v(a[1]))); - return {subtract(combined, cross, multiply_v(cross)), add(cross, cross)}; - } - - template - BOOST_FORCEINLINE void load(tower_type &a, tower_type &b, std::span left, - std::span right, std::size_t stride, std::size_t index) const { - constexpr std::size_t width = sizeof(Value) / sizeof(double); - for (std::size_t w = 0; w < 2; ++w) { - for (std::size_t v = 0; v < 3; ++v) { - for (std::size_t u = 0; u < 2; ++u) { - const std::size_t offset = (6 * w + 2 * v + u) * stride + index; - if constexpr (width == 1) { - a[w][v][u] = reduce(left[offset]); - if constexpr (!Square) { - b[w][v][u] = reduce(right[offset]); - } -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - } else { - a[w][v][u] = reduce(Arithmetic::load(left.data() + offset)); - if constexpr (!Square) { - b[w][v][u] = reduce(Arithmetic::load(right.data() + offset)); - } -#endif - } - } - } - } - } - }; - - public: - explicit bn254_small_prime_fft_pointwise(std::uint64_t modulus) : arithmetic_(modulus) { - } - - /** Overwrite left with scale*left*right. Both buffers must contain exactly 12*stride cells. */ - void multiply(std::span left, std::span right, std::size_t stride, double scale) const { - check_buffers(left, right, stride, scale); - if (arithmetic_.fma_supported()) { - apply(left, left, right, scale); - } else { - apply(left, left, right, scale); - } - } - - /** Overwrite values with scale*values^2, loading each tower element once. */ - void square(std::span values, std::size_t stride, double scale) const { - check_buffers(values, values, stride, scale); - if (arithmetic_.fma_supported()) { - apply(values, values, values, scale); - } else { - apply(values, values, values, scale); - } - } - - /** Compute scale*(a*b+c*d). Every buffer has 12*stride cells; output may equal any input. */ - void multiply_sum(std::span output, std::span a, std::span b, - std::span c, std::span d, std::size_t stride, - double scale) const { - check_buffers(output, a, stride, scale); - check_buffers(output, b, stride, scale); - check_buffers(output, c, stride, scale); - check_buffers(output, d, stride, scale); - if (arithmetic_.fma_supported()) { - apply(output, a, b, scale, c, d); - } else { - apply(output, a, b, scale, c, d); - } - } - - private: - void check_buffers(std::span left, std::span right, std::size_t stride, - double scale) const { - if (left.size() % 12 != 0 || left.size() / 12 != stride || left.size() != right.size()) { - throw std::invalid_argument("bn254_small_prime_fft_pointwise: invalid buffer dimensions"); - } - assert(std::isfinite(scale) && std::trunc(scale) == scale && - std::abs(scale) < static_cast(arithmetic_.modulus()) / 2); - } - - template - void apply(std::span output, std::span left, std::span right, double scale, - std::span extra_left = {}, std::span extra_right = {}) const { - // Derive the indexing bound from the buffer extent checked by the public entry points. - const std::size_t stride = output.size() / 12; - std::size_t offset = 0; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - if constexpr (UseFma) { -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) - using wide_kernel = small_prime_fft_arithmetic::wide_fma_kernel; - offset = stride - stride % wide_kernel::lanes; - if (offset != 0) { - const wide_kernel wide(arithmetic_); - apply_range(wide, output, left, right, extra_left, extra_right, 0, offset, - wide_kernel::broadcast(scale)); - } -#endif - using packed_kernel = small_prime_fft_arithmetic::packed_fma_kernel; - const std::size_t packed_end = stride - stride % packed_kernel::lanes; - if (offset < packed_end) { - const packed_kernel packed(arithmetic_); - apply_range(packed, output, left, right, extra_left, extra_right, offset, - packed_end, packed_kernel::broadcast(scale)); - offset = packed_end; - } - } -#endif - apply_range(arithmetic_, output, left, right, extra_left, extra_right, offset, stride, - scale); - } - - template - static void apply_range(const Arithmetic &arithmetic, std::span output, std::span left, - std::span right, std::span extra_left, - std::span extra_right, std::size_t begin, std::size_t end, Value scale) { - constexpr std::size_t width = sizeof(Value) / sizeof(double); - const std::size_t stride = output.size() / 12; - const kernel operations {arithmetic}; - using tower_type = typename kernel::tower_type; - for (std::size_t i = begin; i < end; i += width) { - tower_type a, b; - operations.template load(a, b, left, right, stride, i); - tower_type result; - if constexpr (Square) { - result = operations.square(a); - } else { - result = operations.multiply(a, b); - if constexpr (Sum) { - // Reuse the input slots after the first product. Load both remaining - // inputs before storing anything so output may still alias any input. - operations.template load(a, b, extra_left, extra_right, stride, i); - // Each product coordinate has magnitude < 3*q; their sum is exact below 6*q. - // Reduce that sum before the usual inverse-transform normalization. - result = operations.add(result, operations.multiply(a, b)); - } - } - for (std::size_t w = 0; w < 2; ++w) { - for (std::size_t v = 0; v < 3; ++v) { - // Fold inverse-transform normalization into this existing write pass, - // pairing the independent scales of the two quadratic coordinates. - const std::size_t offset = (6 * w + 2 * v) * stride + i; - Value first = result[w][v][0]; - Value second = result[w][v][1]; - operations.multiply_pair(first, second, scale, scale); - if constexpr (width == 1) { - output[offset] = first; - output[offset + stride] = second; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) - } else { - Arithmetic::store(output.data() + offset, first); - Arithmetic::store(output.data() + offset + stride, second); -#endif - } - } - } - } - } - - small_prime_fft_arithmetic arithmetic_; - }; -#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) -#pragma GCC diagnostic pop -#endif - -} // namespace nil::crypto3::math::detail - -#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/arithmetic/half_gcd.hpp b/libs/math/include/nil/crypto3/math/polynomial/arithmetic/half_gcd.hpp index b057b1465..e4ee92a50 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/arithmetic/half_gcd.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/arithmetic/half_gcd.hpp @@ -364,13 +364,7 @@ namespace nil::crypto3::math::detail { template void validate_half_gcd_inputs(const typename Backend::polynomial_type &first, const typename Backend::polynomial_type &second) { - using value_type = typename Backend::polynomial_type::value_type; - - const auto is_canonical = [](const typename Backend::polynomial_type &polynomial) { - return polynomial.size() != 0 && - (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); - }; - if (!is_canonical(first) || !is_canonical(second)) { + if (!is_canonical_polynomial(first) || !is_canonical_polynomial(second)) { throw std::invalid_argument("half_gcd_reduce: inputs must be canonical polynomials"); } if (is_zero(first)) { diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp index 6876af248..c6399e24b 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/bn254_small_prime_fft_backend.hpp @@ -36,7 +36,7 @@ #include #include -#include +#include #include #include #include diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft.hpp new file mode 100644 index 000000000..e0a8ecfa9 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft.hpp @@ -0,0 +1,35 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_BN254_SMALL_PRIME_FFT_HPP +#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP + +// Internal components of the BN254 small-prime polynomial backend. Keep this aggregate +// include for callers that use the complete conversion, transform and reconstruction pipeline. +#include +#include +#include +#include + +#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_conversion.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_conversion.hpp new file mode 100644 index 000000000..42423b78e --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_conversion.hpp @@ -0,0 +1,204 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_BN254_SMALL_PRIME_FFT_CONVERSION_HPP +#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_CONVERSION_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace nil::crypto3::math::detail { + + /** + * Convert canonical BN254 Fp12 coordinates to a prefix of the small-prime catalogue. + * Coordinate order is Crypto3's [1,u,v,uv,v^2,uv^2,w,uw,vw,uvw,v^2w,uv^2w]. + * Each output cell is an exact integer in [0,q), congruent to the canonical base-field + * coordinate modulo q. In particular, a base-field -1 is lifted to p-1 before reduction. + * + * Precomputed powers are reusable across operands and sizes. Conversion extracts each + * coordinate once for all primes, using fixed-size scratch and caller-owned output storage. + * Kernel selection occurs once per call; the calling thread's rounding mode must remain + * unchanged until it returns. This conversion has no constant-time guarantee. + */ + class bn254_small_prime_fft_converter { + static constexpr std::size_t chunk_bits = 52; + static constexpr std::size_t chunk_count = 5; + static constexpr std::size_t block_size = 32; + using wide_type = boost::multiprecision::uint128_modular_t; + + struct prime_context { + small_prime_fft_arithmetic arithmetic; + std::array powers; + + explicit prime_context(std::uint64_t modulus) : arithmetic(modulus) { + std::uint64_t power = 1; + for (auto &weight : powers) { + power = ((wide_type(power) << chunk_bits) % modulus).convert_to(); + // Center the weights so that |weight| < q/2. These conversions are exact + // under every rounding mode, and setup products fit in 102 bits. + const auto centered = static_cast(power) - + (power > modulus / 2 ? static_cast(modulus) : 0); + weight = static_cast(centered); + } + } + }; + + public: + using field_type = algebra::fields::fp12_2over3over2>; + using value_type = field_type::value_type; + + explicit bn254_small_prime_fft_converter(std::size_t prime_count) { + if (prime_count == 0 || prime_count > small_prime_fft_primes.size()) { + throw std::invalid_argument("bn254_small_prime_fft_converter: invalid catalogue prefix"); + } + primes_.reserve(prime_count); + for (std::size_t i = 0; i < prime_count; ++i) { + primes_.emplace_back(small_prime_fft_primes[i].modulus); + } + } + + std::size_t prime_count() const { + return primes_.size(); + } + + /** + * Write output[(prime*12 + coordinate)*stride + coefficient], padding each channel + * with zeros from coefficients.size() to stride. Pass a subspan to convert a prefix. + * stride need not be a power of two; the transform planner supplies that restriction. + * + * Require coefficients.size() <= stride and output.size() == prime_count()*12*stride. + * Check dimensions without overflow before writing any output; throw invalid_argument + * on mismatch. Empty input writes zeros, including when reusing a previous output buffer. + */ + void convert(std::span output, std::span coefficients, std::size_t stride) const { + const std::size_t channels = field_type::arity * primes_.size(); + if (channels == 0 || coefficients.size() > stride || stride > output.size() / channels || + output.size() != channels * stride) { + throw std::invalid_argument("bn254_small_prime_fft_converter: invalid buffer dimensions"); + } + if (output.empty()) { + return; + } + const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { + return prime.arithmetic.fma_supported(); + }); + if (use_fma) { + convert_impl(output, coefficients, stride); + } else { + convert_impl(output, coefficients, stride); + } + } + + private: + template + void convert_impl(std::span output, std::span coefficients, + std::size_t stride) const { + for (std::size_t channel = 0; channel < field_type::arity * primes_.size(); ++channel) { + auto padding = output.subspan(channel * stride + coefficients.size(), stride - coefficients.size()); + std::fill(padding.begin(), padding.end(), 0.0); + } + + std::array, block_size> chunks; + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + for (std::size_t offset = 0; offset < coefficients.size();) { + const std::size_t count = std::min(block_size, coefficients.size() - offset); + for (std::size_t i = 0; i < count; ++i) { + std::array words {}; + using base_backend_type = typename Params::base_field_type::integral_type::backend_type; + if constexpr (base_backend_type::limb_bits == 64) { + namespace fp12_fast = algebra::fields::detail::fp12_fast; + // The stored value is a*R mod p, with R = 2^256. Zero-extend it + // into the low half of montgomery; one reduction recovers a in [0,p). + typename Params::limb_array montgomery {}; + const auto stored = fp12_fast::load_limbs( + coefficients[offset + i].coordinate(coordinate).data.backend().base_data()); + std::copy(stored.begin(), stored.end(), montgomery.begin()); + std::array canonical; + fp12_fast::montgomery_reduce(canonical.data(), montgomery); + // pack only extracts five 52-bit chunks; it does not rescale them. + // Packing the canonical value preserves the lift of field -1 to p-1. + words = fp12_fast::montgomery_radix52::pack(canonical.data()); + } else { + const auto canonical = coefficients[offset + i].coordinate(coordinate).to_integral(); + boost::multiprecision::export_bits(canonical, words.begin(), chunk_bits, false); + } + for (std::size_t j = 0; j < chunk_count; ++j) { + chunks[i][j] = static_cast(words[j]); + } + } + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + const auto &context = primes_[prime]; + const auto &arithmetic = context.arithmetic; + const double modulus = static_cast(arithmetic.modulus()); + double *destination = + output.data() + (prime * field_type::arity + coordinate) * stride + offset; + for (std::size_t i = 0; i < count; ++i) { + double sum; + if constexpr (UseFma) { + sum = arithmetic.reduce_fma(chunks[i][0]); + } else { + sum = arithmetic.reduce_integer(chunks[i][0]); + } + // Each chunk is < 2^52 and each weight has magnitude < q/2, so the + // products are < 4*q^2 (q > 2^49). The initial remainder is < q; + // adding four terms of magnitude < 1.5*q keeps every partial sum + // inside (-7*q,7*q), hence below 2^53 and exactly representable. + for (std::size_t j = 1; j < chunk_count; ++j) { + if constexpr (UseFma) { + sum += arithmetic.multiply_fma(chunks[i][j], context.powers[j - 1]); + } else { + sum += arithmetic.multiply_integer(chunks[i][j], context.powers[j - 1]); + } + } + if constexpr (UseFma) { + sum = arithmetic.reduce_fma(sum); + } else { + sum = arithmetic.reduce_integer(sum); + } + destination[i] = sum < 0 ? sum + modulus : sum; + } + } + offset += count; + } + } + } + + std::vector primes_; + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_CONVERSION_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_crt.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_crt.hpp new file mode 100644 index 000000000..d92d43539 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_crt.hpp @@ -0,0 +1,623 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_BN254_SMALL_PRIME_FFT_CRT_HPP +#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_CRT_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#if defined(CRYPTO3_HAS_AVX512IFMA) && defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_AVAILABLE) && \ + !defined(CRYPTO3_SMALL_PRIME_CRT_DISABLE_IFMA) +#define CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED +#endif + +namespace nil::crypto3::math::detail { + + // Independent coefficient lanes, using the existing IFMA operation policy. A CRT digit is + // below 2^50 and each five-digit export weight is below p. With at most twelve primes, + // every unnormalized column is below 24*2^52 and the complete sum is below 2^54*p. + template + struct small_prime_fft_crt_accumulator { + using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; + using vector = typename Operations::value_type; + std::array columns; + + small_prime_fft_crt_accumulator() { + columns.fill(Operations::broadcast(0)); + } + + void add(const typename radix::value_type &weight, vector digits) { + for (std::size_t limb = 0; limb < weight.size(); ++limb) { + const auto w = Operations::broadcast(weight[limb]); + columns[limb] = Operations::multiply_add_low(columns[limb], digits, w); + columns[limb + 1] = Operations::multiply_add_high(columns[limb + 1], digits, w); + } + } + + // The certificate reciprocal has up to 61 bits, so it cannot be an IFMA operand + // without splitting. Compute floor(d*r/2^50) exactly from four 52-bit split products. + static vector certificate_term(vector digits, std::uint64_t reciprocal) { + const auto zero = Operations::broadcast(0); + const auto low = Operations::broadcast(reciprocal & radix::mask); + const auto high = Operations::broadcast(reciprocal >> 52); + const auto a = Operations::multiply_add_low(zero, digits, low); + const auto b = Operations::multiply_add_high(zero, digits, low); + const auto c = Operations::multiply_add_low(zero, digits, high); + const auto d = Operations::multiply_add_high(zero, digits, high); + return Operations::add( + Operations::add(Operations::template shift_right<50>(a), Operations::template shift_left<2>(b)), + Operations::add(Operations::template shift_left<2>(c), Operations::template shift_left<54>(d))); + } + + struct certificate { + vector quotients; + vector accepted; + }; + + static certificate certify(vector lower, std::uint64_t scaled_bound, std::size_t prime_count) { + constexpr std::uint64_t scale = std::uint64_t(1) << 60; + const auto zero = Operations::broadcast(0); + const auto one = Operations::broadcast(1); + const auto quotient = + Operations::template shift_right<60>(Operations::add(lower, Operations::broadcast(scale / 2))); + const auto center = Operations::template shift_left<60>(quotient); + const auto bound = Operations::broadcast(scaled_bound); + const auto upper = Operations::add(lower, Operations::broadcast(2 * prime_count)); + const auto lower_ok = + Operations::select(Operations::less(Operations::add(lower, bound), center), zero, one); + const auto upper_ok = + Operations::select(Operations::less(Operations::add(center, bound), upper), zero, one); + return {quotient, Operations::bit_and(lower_ok, upper_ok)}; + } + + // Add the selected centering correction and divide by S=2^260. The caller prepares + // weights and corrections with R*S (R=2^256), so the result is already R-Montgomery. + std::array reduce(const std::array &correction) const { + const auto zero = Operations::broadcast(0); + const auto mask = Operations::broadcast(radix::mask); + const auto inverse = Operations::broadcast(radix::p_dash()); + std::array t; + t.fill(zero); + std::copy(columns.begin(), columns.end(), t.begin()); + for (std::size_t i = 0; i < correction.size(); ++i) + t[i] = Operations::add(t[i], correction[i]); + for (std::size_t i = 0; i < 9; ++i) { + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + t[i] = Operations::bit_and(t[i], mask); + } + // Each cancellation adds at most five split products to a column. All columns + // stay below 2^57; the initial sum plus correction is < 2^54*p < S*p. + for (std::size_t i = 0; i < 5; ++i) { + const auto m = Operations::multiply_add_low(zero, t[i], inverse); + for (std::size_t j = 0; j < 5; ++j) { + const auto p = Operations::broadcast(radix::modulus()[j]); + t[i + j] = Operations::multiply_add_low(t[i + j], m, p); + t[i + j + 1] = Operations::multiply_add_high(t[i + j + 1], m, p); + } + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + std::array reduced; + for (std::size_t i = 5; i < 9; ++i) { + reduced[i - 5] = Operations::bit_and(t[i], mask); + t[i + 1] = Operations::add(t[i + 1], Operations::shift52(t[i])); + } + reduced[4] = t[9]; + // REDC(T)<2*p::subtract_modulus(reduced); + return reduced; + } + }; + + /** + * Reconstruct signed integer coordinates from a prefix of the small-prime catalogue, then + * reduce them into BN254 Fp12. The supplied bound C must cover the original integer expression, + * with M > 2*C for the prime product M. Centering must precede reduction into the target field: + * the CRT representative M-1 denotes -1, not the field element (M-1) mod p. + * + * Cofactors, inverses and Montgomery export weights are prepared once. An integer interval + * certificate permits direct export when it proves both centering and the coefficient bound; + * uncertain cases use exact wide reconstruction. Both paths use fixed-size scratch and + * caller-owned output, with no arbitrary-precision setup work in the coefficient loop. + * This operation has no constant-time guarantee. + */ + class bn254_small_prime_fft_reconstructor { + using base_field_type = algebra::fields::alt_bn128_base_field<254>; + static constexpr std::size_t block_size = 32; + static constexpr std::size_t correction_count = std::bit_width(small_prime_fft_primes.size() - 1); + // Each CRT term is < M, so their sum is < prime_count*M < 12*2^600 < 2^604. + static constexpr unsigned crt_bits = 50 * small_prime_fft_primes.size() + correction_count; + using integer_type = boost::multiprecision::number>; + static constexpr unsigned export_split = 384; + using export_low_type = + boost::multiprecision::number>; + using export_type = boost::multiprecision::uint512_modular_t; + static_assert(50 * small_prime_fft_primes.size() <= 600, "CRT export folding requires M < 2^600"); + using fast_params = bn254_small_prime_fft_converter::field_type::extension_policy::fast_params; + using direct_weight_type = std::array; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + using radix52 = algebra::fields::detail::fp12_fast::montgomery_radix52; +#endif + static constexpr bool direct_export_supported = base_field_type::integral_type::backend_type::limb_bits == 64; + static constexpr unsigned certificate_bits = 60; + static constexpr std::uint64_t certificate_scale = std::uint64_t(1) << certificate_bits; + static_assert(std::all_of(small_prime_fft_primes.begin(), small_prime_fft_primes.end(), [](const auto &prime) { + return prime.modulus > (std::uint64_t(1) << 49) && prime.modulus < (std::uint64_t(1) << 50); + })); + + struct prime_context { + small_prime_fft_arithmetic arithmetic; + integer_type cofactor; + double inverse; + std::uint64_t certificate_reciprocal; + direct_weight_type export_weight; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + radix52::value_type export_weight52; +#endif + + prime_context(std::uint64_t modulus, const integer_type &cofactor, const direct_weight_type &weight) : + arithmetic(modulus), cofactor(cofactor), export_weight(weight) { + using wide_type = boost::multiprecision::uint128_modular_t; + const auto residue = (cofactor % modulus).convert_to(); + inverse = static_cast(boost::multiprecision::inverse_mod(wide_type(residue), wide_type(modulus)) + .convert_to()); + certificate_reciprocal = + ((wide_type(1) << (certificate_bits + 50)) / modulus).convert_to(); + } + }; + + public: + using field_type = bn254_small_prime_fft_converter::field_type; + using value_type = field_type::value_type; + + bn254_small_prime_fft_reconstructor(const small_prime_fft_crt_params ¶ms, + const boost::multiprecision::cpp_int &coefficient_bound) { + if (params.prime_count == 0 || params.prime_count > small_prime_fft_primes.size() || + coefficient_bound < 0) { + throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid CRT parameters"); + } + boost::multiprecision::cpp_int product = 1; + for (std::size_t i = 0; i < params.prime_count; ++i) { + product *= small_prime_fft_primes[i].modulus; + } + if (params.modulus != product || 2 * coefficient_bound >= product) { + throw std::invalid_argument( + "bn254_small_prime_fft_reconstructor: inconsistent modulus or insufficient capacity"); + } + + // The checks above establish positivity and capacity before conversion to fixed precision. + using setup_backend = integer_type::backend_type::cpp_int_type; + multiples_[0].backend().from_cpp_int(setup_backend(product.backend())); + bound_.backend().from_cpp_int(setup_backend(coefficient_bound.backend())); + half_modulus_ = multiples_[0] >> 1; + const boost::multiprecision::cpp_int sum_bound = params.prime_count * product - 1; + constexpr unsigned limb_bits = integer_type::backend_type::limb_bits; + accumulator_limbs_ = (boost::multiprecision::msb(sum_bound) + 1 + limb_bits - 1) / limb_bits; + export_radix_ = (export_type(1) << export_split) % export_type(base_field_type::modulus); + for (std::size_t bit = 1; bit < correction_count; ++bit) { + multiples_[bit] = multiples_[bit - 1] << 1; + } + primes_.reserve(params.prime_count); + const boost::multiprecision::cpp_int characteristic(base_field_type::modulus.backend().to_cpp_int()); + const boost::multiprecision::cpp_int r_squared = + (boost::multiprecision::cpp_int(1) << 512) % characteristic; + scaled_bound_ = ((coefficient_bound << certificate_bits) / product).convert_to(); + for (std::size_t i = 0; i < params.prime_count; ++i) { + const auto modulus = small_prime_fft_primes[i].modulus; + const boost::multiprecision::cpp_int weight = (product / modulus) * r_squared % characteristic; + direct_weight_type limbs {}; + boost::multiprecision::export_bits(weight, limbs.begin(), 64, false); + primes_.emplace_back(modulus, multiples_[0] / modulus, limbs); +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + // S/R=16. Export after radix-52 REDC must retain the field's R=2^256 scaling. + const boost::multiprecision::cpp_int vector_weight = (weight * 16) % characteristic; + limbs.fill(0); + boost::multiprecision::export_bits(vector_weight, limbs.begin(), 64, false); + primes_.back().export_weight52 = radix52::pack(limbs.data()); +#endif + } + // The centering quotient is in [0,prime_count]. Store canonical -k*M*R^2 mod p. + for (std::size_t k = 0; k <= params.prime_count; ++k) { + boost::multiprecision::cpp_int correction = (k * product * r_squared) % characteristic; + if (correction != 0) { + correction = characteristic - correction; + } + boost::multiprecision::export_bits(correction, export_corrections_[k].begin(), 64, false); +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + const boost::multiprecision::cpp_int vector_correction = (correction * 16) % characteristic; + direct_weight_type limbs {}; + boost::multiprecision::export_bits(vector_correction, limbs.begin(), 64, false); + const auto packed = radix52::pack(limbs.data()); + for (std::size_t digit = 0; digit < packed.size(); ++digit) + export_corrections52_[digit][k] = packed[digit]; +#endif + } + } + + std::size_t prime_count() const { + return primes_.size(); + } + + /** + * Combine an inverse-transform scale with (M/q)^-1, returning a centered residue. + * Require an integer-valued scale of magnitude < 2^53; invalid prime indices throw. + */ + double weighted_normalization(std::size_t prime, double scale) const { + const auto &context = primes_.at(prime); + const double q = static_cast(context.arithmetic.modulus()); + double result = + context.arithmetic.multiply_integer(context.arithmetic.reduce_integer(scale), context.inverse); + if (result > q / 2) { + result -= q; + } else if (result < -q / 2) { + result += q; + } + return result; + } + + /** + * Read residues[(prime*12 + coordinate)*stride + coefficient] for each requested output. + * Require output.size() <= stride and residues.size() == prime_count()*12*stride. + * Check dimensions without overflow before writing any output; throw invalid_argument on + * mismatch. An empty output does no work; entries after output.size() in each channel are + * not read. Input residues must already include any inverse-transform normalization. + * + * @pre Each consumed residue is an integer-valued double of magnitude < 2^53. The thread's + * rounding mode must remain unchanged during the call; kernel selection occurs once. + * + * Throw overflow_error if a centered coordinate exceeds the declared bound; earlier output + * coordinates may already have been written. This check cannot detect an incorrect original + * bound when the integer expression has wrapped by M, so the caller must justify that bound. + */ + void reconstruct(std::span output, std::span residues, std::size_t stride) const { + reconstruct_dispatch(output, residues, stride); + } + + /** + * Reconstruct residues already multiplied by (M/q)^-1 modulo each q. All dimensions, + * residue bounds and failure semantics are those of reconstruct(). The backend can fuse + * this weight with 1/N in the pointwise pass because the inverse transform is linear. + */ + void reconstruct_weighted(std::span output, std::span residues, + std::size_t stride) const { + reconstruct_dispatch(output, residues, stride); + } + + private: + template + void reconstruct_dispatch(std::span output, std::span residues, + std::size_t stride) const { + const std::size_t channels = field_type::arity * primes_.size(); + if (channels == 0 || output.size() > stride || stride > residues.size() / channels || + residues.size() != channels * stride) { + throw std::invalid_argument("bn254_small_prime_fft_reconstructor: invalid buffer dimensions"); + } + if (output.empty()) { + return; + } + const bool use_fma = std::all_of(primes_.begin(), primes_.end(), [](const prime_context &prime) { + return prime.arithmetic.fma_supported(); + }); + if (use_fma) { +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + reconstruct_ifma(output, residues, stride); + return; +#endif + reconstruct_impl(output, residues, stride); + } else { + reconstruct_impl(output, residues, stride); + } + } + + void add_cofactor_product(integer_type &sum, const integer_type &cofactor, std::uint64_t digit) const { + using limb_type = boost::multiprecision::limb_type; + using wide_type = boost::multiprecision::double_limb_type; + constexpr unsigned limb_bits = integer_type::backend_type::limb_bits; + if constexpr (limb_bits >= 50) { + auto *destination = sum.backend().limbs(); + const auto *source = cofactor.backend().limbs(); + wide_type carry = 0; + // A digit is < 2^50. Each product plus the accumulator limb and carry fits + // in a double limb. The complete sum is < prime_count*M, so its final carry + // fits in the active width prepared by the constructor; higher limbs stay zero. + for (std::size_t limb = 0; limb < accumulator_limbs_; ++limb) { + const wide_type value = wide_type(source[limb]) * digit + destination[limb] + carry; + destination[limb] = static_cast(value); + carry = value >> limb_bits; + } + assert(carry == 0); + assert((destination[sum.backend().size() - 1] & ~integer_type::backend_type::upper_limb_mask) == 0); + } else { + // A 50-bit digit spans more than one limb on a 32-bit limb backend. + sum += cofactor * digit; + } + } + + base_field_type::value_type export_magnitude(const integer_type &magnitude) const { + // Centering and the coefficient-bound check have already established magnitude < M < 2^600. + // Write magnitude = low + high*2^384 and replace 2^384 by its cached residue modulo p. + // The fold is < 2^384 + 2^216*2^254 < 2^471, within the existing BN254 Barrett reducer's + // range. This avoids its general-division fallback while retaining the field's normal + // Montgomery conversion. The fixed-width low conversion intentionally discards high bits. + const export_type high(magnitude >> export_split); + export_type folded {export_low_type(magnitude)}; + if (high != 0) { + folded += high * export_radix_; + } + return base_field_type::value_type(folded); + } + + base_field_type::value_type export_exact(integer_type magnitude) const { + // Each canonical CRT term is < M. Reduce their sum to [0,M) without division. + for (std::size_t bit = std::bit_width(primes_.size() - 1); bit-- > 0;) { + if (magnitude >= multiples_[bit]) { + magnitude -= multiples_[bit]; + } + } + const bool negative = magnitude > half_modulus_; + if (negative) { + magnitude = multiples_[0] - magnitude; + } + if (magnitude > bound_) { + throw std::overflow_error("bn254_small_prime_fft_reconstructor: coefficient exceeds CRT bound"); + } + const auto value = export_magnitude(magnitude); + return negative ? -value : value; + } + + static void add_export_product(fast_params::limb_array &sum, const direct_weight_type &weight, + std::uint64_t digit) { + using wide_type = algebra::fields::detail::fp12_fast::wide_limb; + wide_type carry = 0; + for (std::size_t limb = 0; limb < weight.size(); ++limb) { + const wide_type value = wide_type(weight[limb]) * digit + sum[limb] + carry; + sum[limb] = static_cast(value); + carry = value >> 64; + } + // At most twelve digits < 2^50 times weights < p, plus one correction < p. + // The total is < 2^54*p < 2^308, so the fifth limb cannot overflow. + sum[weight.size()] += static_cast(carry); + } + + bool certify_export(std::uint64_t lower, std::size_t "ient) const { + // For d < q < 2^50 and a=floor(2^110/q), floor(d*a/2^50) <= 2^60*d/q + // < floor(d*a/2^50)+2. Summing encloses 2^60*sum(d/q) in [lower,upper]. + // With at most twelve primes, all quantities below fit in uint64_t. + const std::uint64_t upper = lower + 2 * primes_.size(); + quotient = (lower + certificate_scale / 2) >> certificate_bits; + const std::uint64_t center = std::uint64_t(quotient) << certificate_bits; + // Containment in k +/- C/M proves the bound AND the unique centering quotient, + // since C/M < 1/2. Using floor(2^60*C/M) makes this test conservative. + return lower + scaled_bound_ >= center && upper <= center + scaled_bound_; + } + + template + base_field_type::value_type export_direct(typename Params::limb_array &sum, std::size_t quotient) const { + add_export_product(sum, export_corrections_[quotient], 1); + // sum = (sum(d_i*M_i)-k*M)*R^2 mod p, and sum < 2^54*p < R*p for R=2^256. + // One existing REDC therefore yields a canonical Montgomery coordinate directly, + // without materializing the wide CRT integer, folding, or a second field conversion. + direct_weight_type reduced; + algebra::fields::detail::fp12_fast::montgomery_reduce(reduced.data(), sum); + base_field_type::value_type result; + std::copy(reduced.begin(), reduced.end(), result.data.backend().base_data().limbs()); + return result; + } + + template + void reconstruct_impl(std::span output, std::span residues, + std::size_t stride) const { + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + for (std::size_t offset = 0; offset < output.size();) { + const std::size_t count = std::min(block_size, output.size() - offset); + using accumulator_type = std::conditional_t; + std::array sums {}; + std::array lower {}; + std::array, small_prime_fft_primes.size()> saved_digits; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + const auto &context = primes_[prime]; + const auto &arithmetic = context.arithmetic; + const double modulus = static_cast(arithmetic.modulus()); + const double *source = + residues.data() + (prime * field_type::arity + coordinate) * stride + offset; + const auto accumulate = [&](std::size_t i, std::uint64_t digit) { + if constexpr (Direct) { + using wide_type = algebra::fields::detail::fp12_fast::wide_limb; + saved_digits[prime][i] = digit; + lower[i] += static_cast( + (wide_type(digit) * context.certificate_reciprocal) >> 50); + add_export_product(sums[i], context.export_weight, digit); + } else { + add_cofactor_product(sums[i], context.cofactor, digit); + } + }; + std::size_t i = 0; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { + if (count >= 4) { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic); + const __m256d inverse = _mm256_set1_pd(context.inverse); + const __m256d q = _mm256_set1_pd(modulus); + for (; count - i >= 4; i += 4) { + __m256d digits = packed.reduce_fma(_mm256_loadu_pd(source + i)); + if constexpr (!Weighted) { + digits = packed.multiply_fma(digits, inverse); + } + const __m256d negative = _mm256_cmp_pd(digits, _mm256_setzero_pd(), _CMP_LT_OQ); + digits = _mm256_add_pd(digits, _mm256_and_pd(negative, q)); + std::array canonical; + _mm256_storeu_pd(canonical.data(), digits); + for (std::size_t lane = 0; lane < canonical.size(); ++lane) { + accumulate(i + lane, static_cast(canonical[lane])); + } + } + } + } +#endif + for (; i < count; ++i) { + double digit; + if constexpr (UseFma) { + digit = arithmetic.reduce_fma(source[i]); + if constexpr (!Weighted) { + digit = arithmetic.multiply_fma(digit, context.inverse); + } + } else { + digit = arithmetic.reduce_integer(source[i]); + if constexpr (!Weighted) { + digit = arithmetic.multiply_integer(digit, context.inverse); + } + } + // Both operands have magnitude < q, so the product remainder is in (-q,q). + // Canonical digits give 0 <= digit*(M/q) < M and bound the complete sum. + const auto canonical = static_cast(digit < 0 ? digit + modulus : digit); + accumulate(i, canonical); + } + } + for (std::size_t i = 0; i < count; ++i) { + auto &destination = output[offset + i].coordinate(coordinate); + if constexpr (Direct) { + std::size_t quotient; + if (certify_export(lower[i], quotient)) { + destination = export_direct(sums[i], quotient); + } else { + integer_type exact_sum = 0; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][i]); + } + destination = export_exact(exact_sum); + } + } else { + destination = export_exact(sums[i]); + } + } + offset += count; + } + } + } + +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + template + void reconstruct_ifma(std::span output, std::span residues, + std::size_t stride) const { + using operations = algebra::fields::detail::fp12_fast::ifma52_operations<8>; + for (std::size_t coordinate = 0; coordinate < field_type::arity; ++coordinate) { + for (std::size_t offset = 0; offset < output.size();) { + const std::size_t count = std::min(operations::lanes, output.size() - offset); + const auto active = static_cast<__mmask8>((1u << count) - 1); + using accumulator_type = small_prime_fft_crt_accumulator; + accumulator_type accumulator; + auto lower = operations::broadcast(0); + std::array, small_prime_fft_primes.size()> + saved_digits; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + const auto &context = primes_[prime]; + const small_prime_fft_arithmetic::wide_fma_kernel packed(context.arithmetic); + const double *source = + residues.data() + (prime * field_type::arity + coordinate) * stride + offset; + auto digits = packed.reduce_fma(_mm512_maskz_loadu_pd(active, source)); + if constexpr (!Weighted) { + digits = packed.multiply_fma(digits, _mm512_set1_pd(context.inverse)); + } + const auto negative = _mm512_cmp_pd_mask(digits, _mm512_setzero_pd(), _CMP_LT_OQ); + digits = _mm512_mask_add_pd(digits, negative, digits, + _mm512_set1_pd(static_cast(context.arithmetic.modulus()))); + const operations::value_type canonical {_mm512_cvttpd_epu64(digits)}; + operations::store(saved_digits[prime].data(), canonical); + accumulator.add(context.export_weight52, canonical); + lower = operations::add( + lower, accumulator_type::certificate_term(canonical, context.certificate_reciprocal)); + } + const auto certificate = accumulator_type::certify(lower, scaled_bound_, primes_.size()); + const auto accepted = operations::less(operations::broadcast(0), certificate.accepted) & active; + std::array, 5> reduced; + if (accepted != 0) { + std::array corrections; + for (std::size_t digit = 0; digit < corrections.size(); ++digit) { + // Every quotient is in [0,prime_count], including uncertified lanes. + corrections[digit] = {_mm512_i64gather_epi64(certificate.quotients.data, + export_corrections52_[digit].data(), 8)}; + } + const auto values = accumulator.reduce(corrections); + for (std::size_t digit = 0; digit < values.size(); ++digit) + operations::store(reduced[digit].data(), values[digit]); + } + for (std::size_t lane = 0; lane < count; ++lane) { + auto &destination = output[offset + lane].coordinate(coordinate); + if ((accepted & (1u << lane)) != 0) { + // Numeric repacking only: the R*S weights already account for S/R. + base_field_type::value_type value; + auto *words = value.data.backend().base_data().limbs(); + words[0] = reduced[0][lane] | (reduced[1][lane] << 52); + words[1] = (reduced[1][lane] >> 12) | (reduced[2][lane] << 40); + words[2] = (reduced[2][lane] >> 24) | (reduced[3][lane] << 28); + words[3] = (reduced[3][lane] >> 36) | (reduced[4][lane] << 16); + destination = value; + } else { + integer_type exact_sum = 0; + for (std::size_t prime = 0; prime < primes_.size(); ++prime) { + add_cofactor_product(exact_sum, primes_[prime].cofactor, saved_digits[prime][lane]); + } + destination = export_exact(exact_sum); + } + } + offset += count; + } + } + } +#endif + + std::vector primes_; + std::array multiples_; + integer_type half_modulus_; + integer_type bound_; + std::size_t accumulator_limbs_; + export_type export_radix_; + std::uint64_t scaled_bound_; + std::array export_corrections_ {}; +#if defined(CRYPTO3_SMALL_PRIME_CRT_IFMA_SELECTED) + std::array, 5> export_corrections52_ {}; +#endif + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_CRT_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_pointwise.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_pointwise.hpp new file mode 100644 index 000000000..2a0d0be34 --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/bn254_small_prime_fft_pointwise.hpp @@ -0,0 +1,383 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_BN254_SMALL_PRIME_FFT_POINTWISE_HPP +#define CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_POINTWISE_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace nil::crypto3::math::detail { + + /** + * Pointwise products in the BN254 tower u^2=-1, v^3=9+u, w^2=v, over one small modulus q. + * This is a quotient ring; irreducibility of the tower polynomials modulo q is not required. + * Coordinate channels have the same order as the converter, with one prime per call. + * + * Reuse the field library's Karatsuba tower identities and quadratic-extension squaring + * identity with scalar, four-lane or eight-lane residue kernels. The existing field elements + * bind their base arithmetic to the target modulus and cannot represent these small-prime intermediates. + * + * @pre Input cells are integer-valued doubles of magnitude < 2^53. scale is an integer with + * |scale| < q/2, normally the FFT plan's centered 1/N. Buffers are disjoint or identical. + * The thread's rounding mode stays unchanged during the call. Output cells are in (-q,q). + * No buffers are allocated and no constant-time guarantee is provided. + */ +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) +// The packed kernels use vector coordinates in the existing tower arrays. +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wignored-attributes" +#endif + class bn254_small_prime_fft_pointwise { + + template + struct kernel { + using quadratic_type = std::array; + using cubic_type = std::array; + using tower_type = std::array; + + const Arithmetic &arithmetic; + + Value reduce(Value value) const { + if constexpr (UseFma) { + return arithmetic.reduce_fma(value); + } else { + return arithmetic.reduce_integer(value); + } + } + + Value multiply(Value a, Value b) const { + if constexpr (UseFma) { + return arithmetic.multiply_fma(a, b); + } else { + return arithmetic.multiply_integer(a, b); + } + } + + BOOST_FORCEINLINE void multiply_pair(Value &a, Value &b, Value wa, Value wb) const { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma && requires { arithmetic.multiply_fma_pair(a, b, wa, wb); }) { + arithmetic.multiply_fma_pair(a, b, wa, wb); + } else +#endif + { + a = multiply(a, wa); + b = multiply(b, wb); + } + } + + template + Value add(Value a, Value b) const { + if constexpr (Reduce) { + return reduce(a + b); + } else { + return a + b; + } + } + + // Three reduced coordinates give an exact integer of magnitude < 3*q. + // Callers can defer this reduction to the next bounded sum or final scaling. + template + Value subtract(Value a, Value b, Value c) const { + if constexpr (Reduce) { + return reduce(a - b - c); + } else { + return a - b - c; + } + } + + template + std::array add(const std::array &a, const std::array &b) const { + std::array result; + for (std::size_t i = 0; i < N; ++i) { + result[i] = add(a[i], b[i]); + } + return result; + } + + template + std::array subtract(const std::array &a, const std::array &b, + const std::array &c) const { + std::array result; + for (std::size_t i = 0; i < N; ++i) { + result[i] = subtract(a[i], b[i], c[i]); + } + return result; + } + + // Quadratic inputs and returned coordinates are in (-q,q). Pair sums are < 2*q; + // the Karatsuba product is < 4*q^2 and its remainder < 1.5*q. Correcting by the + // two diagonal products gives a sum < 3.5*q, hence exact in binary64 before reduction. + BOOST_FORCEINLINE quadratic_type multiply(const quadratic_type &a, const quadratic_type &b) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + // Keep the two independent products in registers and interleave their FMA chains. + Value diagonal0 = a[0]; + Value diagonal1 = a[1]; + multiply_pair(diagonal0, diagonal1, b[0], b[1]); + const Value cross = multiply(a[0] + a[1], b[0] + b[1]); + return {reduce(diagonal0 - diagonal1), reduce(cross - diagonal0 - diagonal1)}; + } + + quadratic_type multiply_xi(const quadratic_type &a) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + // (9+u)*(a0+a1*u) = (9*a0-a1) + (a0+9*a1)*u. Reduce 8*a first: + // |8*a| < 8*q < 2^53, whereas forming 9*a directly could lose an integer bit. + return {reduce(reduce(8 * a[0]) + a[0] - a[1]), reduce(reduce(8 * a[1]) + a[1] + a[0])}; + } + + quadratic_type add_multiply_xi(const quadratic_type &a, const quadratic_type &b) const { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + // a + (9+u)*b, with |a| < 3*q and |b| < q. After reducing 8*b, + // each four-term sum is < 6*q < 2^53, so only its final reduction is needed. + return {reduce(reduce(8 * b[0]) + b[0] - b[1] + a[0]), reduce(reduce(8 * b[1]) + b[1] + b[0] + a[1])}; + } + + cubic_type multiply_v(const cubic_type &a) const { + return {multiply_xi(a[2]), a[0], a[1]}; + } + + cubic_type multiply(const cubic_type &a, const cubic_type &b) const { + const auto diagonal0 = multiply(a[0], b[0]); + const auto diagonal1 = multiply(a[1], b[1]); + const auto diagonal2 = multiply(a[2], b[2]); + // cross12 is reduced before multiplication by 8 inside multiply_xi. + // The other cross terms stay < 3*q until these final sums; the cubic result is < q. + // Finish each output pair before forming the next cross term to limit live temporaries. + cubic_type result; + { + const auto cross12 = subtract(multiply(add(a[1], a[2]), add(b[1], b[2])), diagonal1, diagonal2); + result[0] = add_multiply_xi(diagonal0, cross12); + } + { + const auto cross01 = + subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0, diagonal1); + result[1] = add_multiply_xi(cross01, diagonal2); + } + { + const auto cross02 = + subtract(multiply(add(a[0], a[2]), add(b[0], b[2])), diagonal0, diagonal2); + result[2] = add(cross02, diagonal1); + } + return result; + } + + // Both tower operations return coordinates of magnitude < 3*q. The final + // |scale| < q/2 gives a product < 1.5*q^2, whose FMA remainder is already in (-q,q). + // Keeping these sums unreduced is also safe when two tower products are added below. + tower_type multiply(const tower_type &a, const tower_type &b) const { + const auto diagonal0 = multiply(a[0], b[0]); + const auto diagonal1 = multiply(a[1], b[1]); + const auto cross = subtract(multiply(add(a[0], a[1]), add(b[0], b[1])), diagonal0, diagonal1); + return {add(diagonal0, multiply_v(diagonal1)), cross}; + } + + tower_type square(const tower_type &a) const { + // The quadratic-extension square needs two cubic products instead of three: + // c0 = (a0+a1)*(a0+v*a1) - a0*a1 - v*a0*a1; c1 = 2*a0*a1. + const auto cross = multiply(a[0], a[1]); + const auto combined = multiply(add(a[0], a[1]), add(a[0], multiply_v(a[1]))); + return {subtract(combined, cross, multiply_v(cross)), add(cross, cross)}; + } + + template + BOOST_FORCEINLINE void load(tower_type &a, tower_type &b, std::span left, + std::span right, std::size_t stride, std::size_t index) const { + constexpr std::size_t width = sizeof(Value) / sizeof(double); + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + for (std::size_t u = 0; u < 2; ++u) { + const std::size_t offset = (6 * w + 2 * v + u) * stride + index; + if constexpr (width == 1) { + a[w][v][u] = reduce(left[offset]); + if constexpr (!Square) { + b[w][v][u] = reduce(right[offset]); + } +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + } else { + a[w][v][u] = reduce(Arithmetic::load(left.data() + offset)); + if constexpr (!Square) { + b[w][v][u] = reduce(Arithmetic::load(right.data() + offset)); + } +#endif + } + } + } + } + } + }; + + public: + explicit bn254_small_prime_fft_pointwise(std::uint64_t modulus) : arithmetic_(modulus) { + } + + /** Overwrite left with scale*left*right. Both buffers must contain exactly 12*stride cells. */ + void multiply(std::span left, std::span right, std::size_t stride, double scale) const { + check_buffers(left, right, stride, scale); + if (arithmetic_.fma_supported()) { + apply(left, left, right, scale); + } else { + apply(left, left, right, scale); + } + } + + /** Overwrite values with scale*values^2, loading each tower element once. */ + void square(std::span values, std::size_t stride, double scale) const { + check_buffers(values, values, stride, scale); + if (arithmetic_.fma_supported()) { + apply(values, values, values, scale); + } else { + apply(values, values, values, scale); + } + } + + /** Compute scale*(a*b+c*d). Every buffer has 12*stride cells; output may equal any input. */ + void multiply_sum(std::span output, std::span a, std::span b, + std::span c, std::span d, std::size_t stride, + double scale) const { + check_buffers(output, a, stride, scale); + check_buffers(output, b, stride, scale); + check_buffers(output, c, stride, scale); + check_buffers(output, d, stride, scale); + if (arithmetic_.fma_supported()) { + apply(output, a, b, scale, c, d); + } else { + apply(output, a, b, scale, c, d); + } + } + + private: + void check_buffers(std::span left, std::span right, std::size_t stride, + double scale) const { + if (left.size() % 12 != 0 || left.size() / 12 != stride || left.size() != right.size()) { + throw std::invalid_argument("bn254_small_prime_fft_pointwise: invalid buffer dimensions"); + } + assert(std::isfinite(scale) && std::trunc(scale) == scale && + std::abs(scale) < static_cast(arithmetic_.modulus()) / 2); + } + + template + void apply(std::span output, std::span left, std::span right, double scale, + std::span extra_left = {}, std::span extra_right = {}) const { + // Derive the indexing bound from the buffer extent checked by the public entry points. + const std::size_t stride = output.size() / 12; + std::size_t offset = 0; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + using wide_kernel = small_prime_fft_arithmetic::wide_fma_kernel; + offset = stride - stride % wide_kernel::lanes; + if (offset != 0) { + const wide_kernel wide(arithmetic_); + apply_range(wide, output, left, right, extra_left, extra_right, 0, offset, + wide_kernel::broadcast(scale)); + } +#endif + using packed_kernel = small_prime_fft_arithmetic::packed_fma_kernel; + const std::size_t packed_end = stride - stride % packed_kernel::lanes; + if (offset < packed_end) { + const packed_kernel packed(arithmetic_); + apply_range(packed, output, left, right, extra_left, extra_right, offset, + packed_end, packed_kernel::broadcast(scale)); + offset = packed_end; + } + } +#endif + apply_range(arithmetic_, output, left, right, extra_left, extra_right, offset, stride, + scale); + } + + template + static void apply_range(const Arithmetic &arithmetic, std::span output, std::span left, + std::span right, std::span extra_left, + std::span extra_right, std::size_t begin, std::size_t end, Value scale) { + constexpr std::size_t width = sizeof(Value) / sizeof(double); + const std::size_t stride = output.size() / 12; + const kernel operations {arithmetic}; + using tower_type = typename kernel::tower_type; + for (std::size_t i = begin; i < end; i += width) { + tower_type a, b; + operations.template load(a, b, left, right, stride, i); + tower_type result; + if constexpr (Square) { + result = operations.square(a); + } else { + result = operations.multiply(a, b); + if constexpr (Sum) { + // Reuse the input slots after the first product. Load both remaining + // inputs before storing anything so output may still alias any input. + operations.template load(a, b, extra_left, extra_right, stride, i); + // Each product coordinate has magnitude < 3*q; their sum is exact below 6*q. + // Reduce that sum before the usual inverse-transform normalization. + result = operations.add(result, operations.multiply(a, b)); + } + } + for (std::size_t w = 0; w < 2; ++w) { + for (std::size_t v = 0; v < 3; ++v) { + // Fold inverse-transform normalization into this existing write pass, + // pairing the independent scales of the two quadratic coordinates. + const std::size_t offset = (6 * w + 2 * v) * stride + i; + Value first = result[w][v][0]; + Value second = result[w][v][1]; + operations.multiply_pair(first, second, scale, scale); + if constexpr (width == 1) { + output[offset] = first; + output[offset + stride] = second; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + } else { + Arithmetic::store(output.data() + offset, first); + Arithmetic::store(output.data() + offset + stride, second); +#endif + } + } + } + } + } + + small_prime_fft_arithmetic arithmetic_; + }; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) +#pragma GCC diagnostic pop +#endif + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_BN254_SMALL_PRIME_FFT_POINTWISE_HPP diff --git a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_arithmetic.hpp similarity index 98% rename from libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp rename to libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_arithmetic.hpp index 72a530d14..f19635f24 100644 --- a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_arithmetic.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_arithmetic.hpp @@ -195,7 +195,7 @@ namespace nil::crypto3::math::detail { __m256d reduce_fma(__m256d value) const { assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double lanes[4]; + alignas(32) double lanes[4]; _mm256_storeu_pd(lanes, value); for (double lane : lanes) { assert(valid_integer(lane)); @@ -220,7 +220,7 @@ namespace nil::crypto3::math::detail { __m256d multiply_fma(__m256d left, __m256d right) const { assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double a[4], b[4]; + alignas(32) double a[4], b[4]; _mm256_storeu_pd(a, left); _mm256_storeu_pd(b, right); for (std::size_t i = 0; i < 4; ++i) { @@ -257,7 +257,7 @@ namespace nil::crypto3::math::detail { assert(&left0 != &left1); assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double a[8], b[8]; + alignas(32) double a[8], b[8]; _mm256_storeu_pd(a, left0); _mm256_storeu_pd(a + 4, left1); _mm256_storeu_pd(b, right0); @@ -329,7 +329,8 @@ namespace nil::crypto3::math::detail { __m512d reduce_fma(__m512d value) const { assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double lanes[8]; + // Keep debug scratch SIMD-aligned even when sanitizers change the stack layout. + alignas(64) double lanes[8]; _mm512_storeu_pd(lanes, value); for (double lane : lanes) { assert(valid_integer(lane)); @@ -354,7 +355,7 @@ namespace nil::crypto3::math::detail { __m512d multiply_fma(__m512d left, __m512d right) const { assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double a[8], b[8]; + alignas(64) double a[8], b[8]; _mm512_storeu_pd(a, left); _mm512_storeu_pd(b, right); for (std::size_t i = 0; i < 8; ++i) { @@ -391,7 +392,7 @@ namespace nil::crypto3::math::detail { assert(&left0 != &left1); assert(arithmetic_.fma_supported()); #ifndef NDEBUG - double a[16], b[16]; + alignas(64) double a[16], b[16]; _mm512_storeu_pd(a, left0); _mm512_storeu_pd(a + 8, left1); _mm512_storeu_pd(b, right0); diff --git a/libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_params.hpp similarity index 100% rename from libs/math/include/nil/crypto3/math/detail/small_prime_fft_params.hpp rename to libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_params.hpp diff --git a/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_plan.hpp b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_plan.hpp new file mode 100644 index 000000000..39fbba73b --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/backends/detail/small_prime_fft/small_prime_fft_plan.hpp @@ -0,0 +1,1090 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_SMALL_PRIME_FFT_PLAN_HPP +#define CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PLAN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#ifndef CRYPTO3_SMALL_PRIME_FFT_RADIX +#define CRYPTO3_SMALL_PRIME_FFT_RADIX 4 +#endif +static_assert(CRYPTO3_SMALL_PRIME_FFT_RADIX == 4 || CRYPTO3_SMALL_PRIME_FFT_RADIX == 8 || + CRYPTO3_SMALL_PRIME_FFT_RADIX == 16, + "small-prime FFT stage radix must be 4, 8 or 16"); + +namespace nil::crypto3::math::detail { + + /** + * In-place, power-of-two transforms for one prime from the catalogue. The forward transform + * maps natural-order coefficients to evaluations at omega^k in bit-reversed index order. + * Both directions use this order, so convolution needs no permutation passes. + * + * A reusable table of max(1,N/2) centered twiddles serves both directions. Radix-4 recursion + * fuses pairs of layers, with radix-2 leaves for odd depths and special cases for unit twiddles. + * AVX2/FMA builds fuse 8-, 16- and 32-point leaves; AVX-512 adds paired 32-/64-point + * leaves and eight-lane stages. CRYPTO3_SMALL_PRIME_FFT_RADIX selects 4 (default), 8 + * or 16 for large SIMD stages; the scalar/integer fallback keeps the radix-4 schedule. + * Transform calls allocate no buffers; recursion uses O(log N) stack space. + * + * @pre Inputs are integer-valued doubles in (-3*q,3*q). The calling thread's rounding mode + * must remain unchanged during a call. Each call selects FMA or integer arithmetic once. + * The forward output is in (-3*q,3*q); the normalized inverse output is in [0,q). + * This operation has no constant-time guarantee. + */ + class small_prime_fft_plan { + public: + small_prime_fft_plan(std::size_t prime_index, std::size_t size) : + arithmetic_(checked_prime(prime_index, size).modulus), size_(size), + twiddles_(std::max(std::size_t(1), size / 2), 1.0) { + using wide_type = boost::multiprecision::uint128_modular_t; + const auto &prime = small_prime_fft_primes[prime_index]; + wide_type root = prime.root; + const unsigned depth = std::countr_zero(size_); + for (unsigned level = prime.two_adicity; level > depth; --level) { + root = root * root % prime.modulus; + } + omega_ = root.convert_to(); + size_inverse_ = centered( + static_cast(boost::multiprecision::inverse_mod(wide_type(size_), wide_type(prime.modulus)) + .convert_to())); + + // At a layer with 2^d blocks, block j uses omega_(2^(d+1))^bitreverse(j,d). + // Each new table segment is the preceding prefix times a primitive root of twice + // the previous order: [1], [i], [w8,i*w8], [w16,i*w16,w8*w16,i*w8*w16], ... . + for (std::size_t half = size_ / 4; half != 0; half >>= 1) { + twiddles_[half] = centered(static_cast(root.convert_to())); + root = root * root % prime.modulus; + } + for (std::size_t half = 1; half < size_ / 2; half <<= 1) { + const double step = twiddles_[half]; + for (std::size_t i = 1; i < half; ++i) { + twiddles_[half + i] = centered(arithmetic_.multiply_integer(twiddles_[i], step)); + } + } + } + + std::size_t size() const { + return size_; + } + + std::uint64_t modulus() const { + return arithmetic_.modulus(); + } + + std::uint64_t omega() const { + return omega_; + } + + /** Exact centered representative of 1/N modulo q, also usable to fuse scaling into pointwise work. */ + double size_inverse() const { + return size_inverse_; + } + + /** Require exactly N entries; invalid dimensions throw before any writes. N=1 is the identity. */ + void fft(std::span values) const { + check_buffer(values); + if (arithmetic_.fma_supported()) { + forward_recursive(values.data(), size_, 0); + } else { + forward_recursive(values.data(), size_, 0); + } + } + + /** + * Require exactly N entries in the forward spectrum order. By default, scale by 1/N and + * return canonical residues. With normalize=false, return N times the inverse in (-2*q,2*q), + * allowing the caller to fuse the scaling elsewhere. Invalid dimensions throw before writes. + */ + void inverse_fft(std::span values, bool normalize = true) const { + check_buffer(values); + if (arithmetic_.fma_supported()) { + inverse_impl(values, normalize); + } else { + inverse_impl(values, normalize); + } + } + + /** + * Unscaled inverse for pointwise products whose integer representatives are in (-q,q). + * The caller has already folded 1/N into its pointwise scale. This stricter input bound + * permits omitting the first leaf reductions; the output still lies in (-2*q,2*q). + * Dimension errors throw before writes, as in inverse_fft(). + */ + void inverse_fft_pointwise(std::span values) const { + check_buffer<1>(values); + if (arithmetic_.fma_supported()) { + inverse_impl(values, false); + } else { + inverse_impl(values, false); + } + } + + private: + static const small_prime_fft_prime &checked_prime(std::size_t prime_index, std::size_t size) { + if (prime_index >= small_prime_fft_primes.size() || !std::has_single_bit(size) || + std::countr_zero(size) > small_prime_fft_primes[prime_index].two_adicity) { + throw std::invalid_argument("small_prime_fft_plan: invalid prime index or transform size"); + } + return small_prime_fft_primes[prime_index]; + } + + template + void check_buffer(std::span values) const { + if (values.size() != size_ || twiddles_.empty()) { + throw std::invalid_argument("small_prime_fft_plan: invalid buffer dimensions or moved-from plan"); + } +#ifndef NDEBUG + const double limit = Bound * static_cast(modulus()); + for (double value : values) { + assert(std::isfinite(value) && std::trunc(value) == value && std::abs(value) < limit); + } +#endif + } + + double centered(double value) const { + const double q = static_cast(modulus()); + return value > q / 2 ? value - q : (value < -q / 2 ? value + q : value); + } + + BOOST_FORCEINLINE double inverse_twiddle(std::size_t index) const { + if (index == 0) { + return 1.0; + } + // Within [b,2*b), reflected exponents add to half the root order, hence w^-1=-w_reflected. + const std::size_t base = std::bit_floor(index); + return -twiddles_[base + (2 * base - 1 - index)]; + } + + // Keep arithmetic selection outside the butterfly loops, with compile-time dispatch inside them. + template + BOOST_FORCEINLINE static Value reduce(const Arithmetic &arithmetic, Value value) { + if constexpr (UseFma) { + return arithmetic.reduce_fma(value); + } else { + return arithmetic.reduce_integer(value); + } + } + + template + BOOST_FORCEINLINE static Value multiply(const Arithmetic &arithmetic, Value left, Value right) { + if constexpr (UseFma) { + return arithmetic.multiply_fma(left, right); + } else { + return arithmetic.multiply_integer(left, right); + } + } + + template + BOOST_FORCEINLINE static void multiply_pair(const Arithmetic &arithmetic, Value &a, Value &b, Value wa, + Value wb) { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma && requires { arithmetic.multiply_fma_pair(a, b, wa, wb); }) { + arithmetic.multiply_fma_pair(a, b, wa, wb); + } else +#endif + { + a = multiply(arithmetic, a, wa); + b = multiply(arithmetic, b, wb); + } + } + + template + BOOST_FORCEINLINE static void forward_butterfly(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, + Value &v3, Value w2, Value w, Value iw) { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + const Value a0 = reduce(arithmetic, v0); + const Value a1 = v1; + Value a2 = v2; + Value a3 = v3; + if constexpr (UnitTwiddle) { + a2 = reduce(arithmetic, a2); + a3 = reduce(arithmetic, a3); + } else { + multiply_pair(arithmetic, a2, a3, w2, w2); + } + const Value even_sum = a0 + a2; + const Value even_difference = a0 - a2; + Value odd_sum = a1 + a3; + Value odd_difference = a1 - a3; + if constexpr (UnitTwiddle) { + odd_sum = reduce(arithmetic, odd_sum); + odd_difference = multiply(arithmetic, odd_difference, iw); + } else { + multiply_pair(arithmetic, odd_sum, odd_difference, w, iw); + } + // Input magnitudes < 3*q and weights < q/2 imply |a2|,|a3| < q and + // |a1 +/- a3| < 4*q. Every product is < 2*q^2, so the odd terms are < q. + // Thus all outputs are < 3*q, with every integer sum exact below 2^53. + v0 = even_sum + odd_sum; + v1 = even_sum - odd_sum; + v2 = even_difference + odd_difference; + v3 = even_difference - odd_difference; + } + +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + BOOST_FORCEINLINE static void transpose4(__m256d &v0, __m256d &v1, __m256d &v2, __m256d &v3) { + const __m256d a = _mm256_unpacklo_pd(v0, v1); + const __m256d b = _mm256_unpackhi_pd(v0, v1); + const __m256d c = _mm256_unpacklo_pd(v2, v3); + const __m256d d = _mm256_unpackhi_pd(v2, v3); + v0 = _mm256_permute2f128_pd(a, c, 0x20); + v1 = _mm256_permute2f128_pd(b, d, 0x20); + v2 = _mm256_permute2f128_pd(a, c, 0x31); + v3 = _mm256_permute2f128_pd(b, d, 0x31); + } + + // Aligned table indices, with unaligned loads: the owning vector needs no special allocator. + template + BOOST_FORCEINLINE __m256d twiddle4(std::size_t index) const { + if constexpr (!Inverse) { + return _mm256_loadu_pd(twiddles_.data() + index); + } else { + if (index == 0) { + return _mm256_setr_pd(1.0, inverse_twiddle(1), inverse_twiddle(2), inverse_twiddle(3)); + } + // Four consecutive entries stay in one [b,2*b) segment. Reflect, reverse, negate. + const std::size_t base = std::bit_floor(index); + const __m256d reflected = _mm256_loadu_pd(twiddles_.data() + 3 * base - index - 4); + return _mm256_sub_pd(_mm256_setzero_pd(), _mm256_permute4x64_pd(reflected, 0x1b)); + } + } + + template + BOOST_FORCEINLINE void leaf_twiddles(std::size_t node, __m256d &w2, __m256d &w, __m256d &iw) const { + w2 = twiddle4(4 * node); + const __m256d low = twiddle4(8 * node); + const __m256d high = twiddle4(8 * node + 4); + w = _mm256_permute4x64_pd(_mm256_unpacklo_pd(low, high), 0xd8); + iw = _mm256_permute4x64_pd(_mm256_unpackhi_pd(low, high), 0xd8); + } + + template + BOOST_FORCEINLINE void forward16(const small_prime_fft_arithmetic::packed_fma_kernel &packed, __m256d &v0, + __m256d &v1, __m256d &v2, __m256d &v3, std::size_t node) const { + forward_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(twiddles_[node]), + _mm256_set1_pd(twiddles_[2 * node]), + _mm256_set1_pd(twiddles_[2 * node + 1])); + // The four child transforms now run across SIMD lanes, without scalar leaves or reloads. + transpose4(v0, v1, v2, v3); + __m256d w2, w, iw; + leaf_twiddles(node, w2, w, iw); + forward_butterfly(packed, v0, v1, v2, v3, w2, w, iw); + transpose4(v0, v1, v2, v3); + } + + template + BOOST_FORCEINLINE void inverse16(const small_prime_fft_arithmetic::packed_fma_kernel &packed, __m256d &v0, + __m256d &v1, __m256d &v2, __m256d &v3, std::size_t node) const { + transpose4(v0, v1, v2, v3); + __m256d w2, w, iw; + leaf_twiddles(node, w2, w, iw); + inverse_butterfly(packed, v0, v1, v2, v3, w2, w, iw, !BoundedInput); + transpose4(v0, v1, v2, v3); + inverse_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(inverse_twiddle(node)), + _mm256_set1_pd(inverse_twiddle(2 * node)), + _mm256_set1_pd(inverse_twiddle(2 * node + 1)), false); + } + + template + void base16(double *values, std::size_t node) const { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + __m256d v0 = _mm256_loadu_pd(values); + __m256d v1 = _mm256_loadu_pd(values + 4); + __m256d v2 = _mm256_loadu_pd(values + 8); + __m256d v3 = _mm256_loadu_pd(values + 12); + if constexpr (Inverse) { + inverse16(packed, v0, v1, v2, v3, node); + } else { + forward16(packed, v0, v1, v2, v3, node); + } + _mm256_storeu_pd(values, v0); + _mm256_storeu_pd(values + 4, v1); + _mm256_storeu_pd(values + 8, v2); + _mm256_storeu_pd(values + 12, v3); + } + + template + BOOST_FORCEINLINE static void radix2_packed(const Arithmetic &packed, Value &left, Value &right, Value weight) { + if constexpr (Inverse) { + // Both subtransforms already return magnitudes < 2*q. No input reductions are needed. + const Value sum = left + right; + const Value difference = left - right; + left = packed.reduce_fma(sum); + right = UnitTwiddle ? packed.reduce_fma(difference) : packed.multiply_fma(difference, weight); + } else { + const Value a = packed.reduce_fma(left); + const Value b = UnitTwiddle ? packed.reduce_fma(right) : packed.multiply_fma(right, weight); + left = a + b; + right = a - b; + } + } + + template + BOOST_FORCEINLINE static void stage4(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, Value &v3, + const double *weights) { + const Value w2 = Arithmetic::broadcast(weights[0]); + const Value w = Arithmetic::broadcast(weights[1]); + const Value iw = Arithmetic::broadcast(weights[2]); + if constexpr (Inverse) { + inverse_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw, false); + } else { + forward_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw); + } + } + + template + BOOST_FORCEINLINE void stage_weights(double *weights, std::size_t node) const { + weights[0] = Inverse ? inverse_twiddle(node) : twiddles_[node]; + weights[1] = Inverse ? inverse_twiddle(2 * node) : twiddles_[2 * node]; + weights[2] = Inverse ? inverse_twiddle(2 * node + 1) : twiddles_[2 * node + 1]; + } + + template + BOOST_FORCEINLINE void stage4(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, Value &v3, + std::size_t node) const { + double weights[3]; + stage_weights(weights, node); + stage4(arithmetic, v0, v1, v2, v3, weights); + } + + template + BOOST_FORCEINLINE static void fused_outer(const Arithmetic &arithmetic, Value *v, const double *weights) { + if constexpr (Radix == 8) { + const Value weight = Arithmetic::broadcast(weights[0]); + radix2_packed(arithmetic, v[0], v[4], weight); + radix2_packed(arithmetic, v[1], v[5], weight); + radix2_packed(arithmetic, v[2], v[6], weight); + radix2_packed(arithmetic, v[3], v[7], weight); + } else { + stage4(arithmetic, v[0], v[4], v[8], v[12], weights); + stage4(arithmetic, v[1], v[5], v[9], v[13], weights); + stage4(arithmetic, v[2], v[6], v[10], v[14], weights); + stage4(arithmetic, v[3], v[7], v[11], v[15], weights); + } + } + + // Fuse three/four binary layers without intermediate array passes. Every constituent + // butterfly retains the forward <3*q and inverse <2*q bounds; no wider lazy range is used. + // UnitTwiddle is valid only at node zero, hence only the first child inherits it. + template + BOOST_FORCEINLINE static void fused_stage(const Arithmetic &arithmetic, Value *v, const double *weights) { + static_assert(Radix == 8 || Radix == 16); + if constexpr (!Inverse) { + fused_outer(arithmetic, v, weights); + } + constexpr std::size_t outer_weights = Radix == 8 ? 1 : 3; + stage4(arithmetic, v[0], v[1], v[2], v[3], weights + outer_weights); + stage4(arithmetic, v[4], v[5], v[6], v[7], weights + outer_weights + 3); + if constexpr (Radix == 16) { + stage4(arithmetic, v[8], v[9], v[10], v[11], weights + outer_weights + 6); + stage4(arithmetic, v[12], v[13], v[14], v[15], weights + outer_weights + 9); + } + if constexpr (Inverse) { + fused_outer(arithmetic, v, weights); + } + } + + template + void large_radix(double *values, std::size_t stride, std::size_t node) const { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + using arithmetic_type = small_prime_fft_arithmetic::wide_fma_kernel; +#else + using arithmetic_type = small_prime_fft_arithmetic::packed_fma_kernel; +#endif + const arithmetic_type arithmetic(arithmetic_); + // Dispatch requires stride >=16, a power of two divisible by the SIMD width. + assert(stride >= 16 && stride % arithmetic_type::lanes == 0); + // The node is fixed throughout this coefficient loop. Resolve reflections and + // table lookups once, retaining scalar weights to avoid forcing Radix-1 live vectors. + double weights[Radix - 1]; + if constexpr (Radix == 8) { + weights[0] = Inverse ? inverse_twiddle(node) : twiddles_[node]; + } else { + stage_weights(weights, node); + } + constexpr std::size_t outer_weights = Radix == 8 ? 1 : 3; +#pragma GCC unroll 4 + for (std::size_t child = 0; child < Radix / 4; ++child) { + stage_weights(weights + outer_weights + 3 * child, Radix / 4 * node + child); + } + for (std::size_t i = 0; i < stride; i += arithmetic_type::lanes) { + arithmetic_type::value_type v[Radix]; +#pragma GCC unroll 16 + for (std::size_t part = 0; part < Radix; ++part) { + v[part] = arithmetic_type::load(values + part * stride + i); + } + fused_stage(arithmetic, v, weights); +#pragma GCC unroll 16 + for (std::size_t part = 0; part < Radix; ++part) { + arithmetic_type::store(values + part * stride + i, v[part]); + } + } + } + + template + void base32(double *values, std::size_t node) const { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + __m256d v0 = _mm256_loadu_pd(values); + __m256d v1 = _mm256_loadu_pd(values + 4); + __m256d v2 = _mm256_loadu_pd(values + 8); + __m256d v3 = _mm256_loadu_pd(values + 12); + __m256d v4 = _mm256_loadu_pd(values + 16); + __m256d v5 = _mm256_loadu_pd(values + 20); + __m256d v6 = _mm256_loadu_pd(values + 24); + __m256d v7 = _mm256_loadu_pd(values + 28); + const __m256d weight = _mm256_set1_pd(Inverse ? inverse_twiddle(node) : twiddles_[node]); + // One radix-2 layer and two register-resident 16-point transforms cover all five layers. + if constexpr (Inverse) { + inverse16(packed, v0, v1, v2, v3, 2 * node); + inverse16(packed, v4, v5, v6, v7, 2 * node + 1); + } + radix2_packed(packed, v0, v4, weight); + radix2_packed(packed, v1, v5, weight); + radix2_packed(packed, v2, v6, weight); + radix2_packed(packed, v3, v7, weight); + if constexpr (!Inverse) { + forward16(packed, v0, v1, v2, v3, 2 * node); + forward16(packed, v4, v5, v6, v7, 2 * node + 1); + } + _mm256_storeu_pd(values, v0); + _mm256_storeu_pd(values + 4, v1); + _mm256_storeu_pd(values + 8, v2); + _mm256_storeu_pd(values + 12, v3); + _mm256_storeu_pd(values + 16, v4); + _mm256_storeu_pd(values + 20, v5); + _mm256_storeu_pd(values + 24, v6); + _mm256_storeu_pd(values + 28, v7); + } + + template + void base8(double *values, std::size_t node) const { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + const __m256d low = _mm256_loadu_pd(values); + const __m256d high = _mm256_loadu_pd(values + 4); + __m256d v0, v1, v2, v3; + if constexpr (Inverse) { + __m256d a = _mm256_permute4x64_pd(_mm256_unpacklo_pd(low, high), 0xd8); + __m256d b = _mm256_permute4x64_pd(_mm256_unpackhi_pd(low, high), 0xd8); + if constexpr (!BoundedInput) { + a = packed.reduce_fma(a); + b = packed.reduce_fma(b); + } + const __m256d sum = packed.reduce_fma(a + b); + const __m256d difference = packed.multiply_fma(a - b, twiddle4(4 * node)); + const __m256d even = _mm256_unpacklo_pd(sum, difference); + const __m256d odd = _mm256_unpackhi_pd(sum, difference); + v0 = _mm256_permute2f128_pd(even, even, 0x00); + v1 = _mm256_permute2f128_pd(odd, odd, 0x00); + v2 = _mm256_permute2f128_pd(even, even, 0x11); + v3 = _mm256_permute2f128_pd(odd, odd, 0x11); + inverse_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(inverse_twiddle(node)), + _mm256_set1_pd(inverse_twiddle(2 * node)), + _mm256_set1_pd(inverse_twiddle(2 * node + 1)), false); + _mm256_storeu_pd(values, _mm256_permute2f128_pd(v0, v1, 0x20)); + _mm256_storeu_pd(values + 4, _mm256_permute2f128_pd(v2, v3, 0x20)); + } else { + // Two independent columns occupy the lower lanes; the duplicated upper lanes are discarded. + v0 = _mm256_permute2f128_pd(low, low, 0x00); + v1 = _mm256_permute2f128_pd(low, low, 0x11); + v2 = _mm256_permute2f128_pd(high, high, 0x00); + v3 = _mm256_permute2f128_pd(high, high, 0x11); + forward_butterfly(packed, v0, v1, v2, v3, _mm256_set1_pd(twiddles_[node]), + _mm256_set1_pd(twiddles_[2 * node]), + _mm256_set1_pd(twiddles_[2 * node + 1])); + const __m256d a = packed.reduce_fma( + _mm256_permute2f128_pd(_mm256_unpacklo_pd(v0, v1), _mm256_unpacklo_pd(v2, v3), 0x20)); + const __m256d b = packed.multiply_fma( + _mm256_permute2f128_pd(_mm256_unpackhi_pd(v0, v1), _mm256_unpackhi_pd(v2, v3), 0x20), + twiddle4(4 * node)); + const __m256d even = _mm256_unpacklo_pd(a + b, a - b); + const __m256d odd = _mm256_unpackhi_pd(a + b, a - b); + _mm256_storeu_pd(values, _mm256_permute2f128_pd(even, odd, 0x20)); + _mm256_storeu_pd(values + 4, _mm256_permute2f128_pd(even, odd, 0x31)); + } + } + + template + void basecase(double *values, std::size_t length, std::size_t node) const { + if (length == 16) { + base16(values, node); + } else if (length == 32) { + base32(values, node); + } else { + base8(values, node); + } + } + +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + BOOST_FORCEINLINE static __m512d join4(__m256d low, __m256d high) { + return _mm512_insertf64x4(_mm512_castpd256_pd512(low), high, 1); + } + + // Two independent 4x4 transposes, one in each 256-bit half. The stored spectrum + // remains bit-reversed; these permutations only change the lanes during a leaf. + BOOST_FORCEINLINE static void transpose4(__m512d &v0, __m512d &v1, __m512d &v2, __m512d &v3) { + const __m512d a = _mm512_unpacklo_pd(v0, v1); + const __m512d b = _mm512_unpackhi_pd(v0, v1); + const __m512d c = _mm512_unpacklo_pd(v2, v3); + const __m512d d = _mm512_unpackhi_pd(v2, v3); + const __m512i low = _mm512_setr_epi64(0, 1, 8, 9, 4, 5, 12, 13); + const __m512i high = _mm512_setr_epi64(2, 3, 10, 11, 6, 7, 14, 15); + v0 = _mm512_permutex2var_pd(a, low, c); + v1 = _mm512_permutex2var_pd(b, low, d); + v2 = _mm512_permutex2var_pd(a, high, c); + v3 = _mm512_permutex2var_pd(b, high, d); + } + + template + BOOST_FORCEINLINE __m512d paired_twiddle(std::size_t low, std::size_t high) const { + return join4(_mm256_set1_pd(Inverse ? inverse_twiddle(low) : twiddles_[low]), + _mm256_set1_pd(Inverse ? inverse_twiddle(high) : twiddles_[high])); + } + + // Keep the complete leaf in one function: outlining this step forces live vectors + // through memory across the call. Recursive transform boundaries remain uninlined. + template + BOOST_FORCEINLINE void paired_outer(const small_prime_fft_arithmetic::wide_fma_kernel &arithmetic, __m512d &v0, + __m512d &v1, __m512d &v2, __m512d &v3, std::size_t node) const { + const __m512d w2 = paired_twiddle(node, node + 1); + const __m512d w = paired_twiddle(2 * node, 2 * node + 2); + const __m512d iw = paired_twiddle(2 * node + 1, 2 * node + 3); + if constexpr (Inverse) { + inverse_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw, false); + } else { + forward_butterfly(arithmetic, v0, v1, v2, v3, w2, w, iw); + } + } + + template + BOOST_FORCEINLINE void base16_pair(const small_prime_fft_arithmetic::wide_fma_kernel &arithmetic, __m512d &v0, + __m512d &v1, __m512d &v2, __m512d &v3, std::size_t node) const { + // Group matching rows of two adjacent 16-point transforms in the two vector halves. + const __m512d a = _mm512_shuffle_f64x2(v0, v2, 0x44); + const __m512d b = _mm512_shuffle_f64x2(v0, v2, 0xee); + const __m512d c = _mm512_shuffle_f64x2(v1, v3, 0x44); + const __m512d d = _mm512_shuffle_f64x2(v1, v3, 0xee); + v0 = a; + v1 = b; + v2 = c; + v3 = d; + if constexpr (!Inverse) { + paired_outer(arithmetic, v0, v1, v2, v3, node); + } + transpose4(v0, v1, v2, v3); + __m256d lo_w2, lo_w, lo_iw, hi_w2, hi_w, hi_iw; + leaf_twiddles(node, lo_w2, lo_w, lo_iw); + leaf_twiddles(node + 1, hi_w2, hi_w, hi_iw); + if constexpr (Inverse) { + inverse_butterfly(arithmetic, v0, v1, v2, v3, join4(lo_w2, hi_w2), join4(lo_w, hi_w), + join4(lo_iw, hi_iw), !BoundedInput); + } else { + forward_butterfly(arithmetic, v0, v1, v2, v3, join4(lo_w2, hi_w2), join4(lo_w, hi_w), + join4(lo_iw, hi_iw)); + } + transpose4(v0, v1, v2, v3); + if constexpr (Inverse) { + paired_outer(arithmetic, v0, v1, v2, v3, node); + } + // Restore contiguous bit-reversed order in the same four input/output vectors. + const __m512d out0 = _mm512_shuffle_f64x2(v0, v1, 0x44); + const __m512d out1 = _mm512_shuffle_f64x2(v2, v3, 0x44); + const __m512d out2 = _mm512_shuffle_f64x2(v0, v1, 0xee); + const __m512d out3 = _mm512_shuffle_f64x2(v2, v3, 0xee); + v0 = out0; + v1 = out1; + v2 = out2; + v3 = out3; + } + + template + void base32_wide(double *values, std::size_t node) const { + const small_prime_fft_arithmetic::wide_fma_kernel arithmetic(arithmetic_); + __m512d v0 = _mm512_loadu_pd(values); + __m512d v1 = _mm512_loadu_pd(values + 8); + __m512d v2 = _mm512_loadu_pd(values + 16); + __m512d v3 = _mm512_loadu_pd(values + 24); + const __m512d weight = _mm512_set1_pd(Inverse ? inverse_twiddle(node) : twiddles_[node]); + if constexpr (Inverse) { + base16_pair(arithmetic, v0, v1, v2, v3, 2 * node); + } + radix2_packed(arithmetic, v0, v2, weight); + radix2_packed(arithmetic, v1, v3, weight); + if constexpr (!Inverse) { + base16_pair(arithmetic, v0, v1, v2, v3, 2 * node); + } + _mm512_storeu_pd(values, v0); + _mm512_storeu_pd(values + 8, v1); + _mm512_storeu_pd(values + 16, v2); + _mm512_storeu_pd(values + 24, v3); + } + + template + void base64_wide(double *values, std::size_t node) const { + const small_prime_fft_arithmetic::wide_fma_kernel arithmetic(arithmetic_); + __m512d v[8]; +#pragma GCC unroll 8 + for (std::size_t i = 0; i < 8; ++i) { + v[i] = _mm512_loadu_pd(values + 8 * i); + } + if constexpr (!Inverse) { + stage4(arithmetic, v[0], v[2], v[4], v[6], node); + stage4(arithmetic, v[1], v[3], v[5], v[7], node); + } + base16_pair(arithmetic, v[0], v[1], v[2], v[3], 4 * node); + base16_pair(arithmetic, v[4], v[5], v[6], v[7], 4 * node + 2); + if constexpr (Inverse) { + stage4(arithmetic, v[0], v[2], v[4], v[6], node); + stage4(arithmetic, v[1], v[3], v[5], v[7], node); + } +#pragma GCC unroll 8 + for (std::size_t i = 0; i < 8; ++i) { + _mm512_storeu_pd(values + 8 * i, v[i]); + } + } + + template + void wide_basecase(double *values, std::size_t length, std::size_t node) const { + if (length == 32) { + if (node == 0) { + base32_wide(values, node); + } else { + base32_wide(values, node); + } + } else { + if (node == 0) { + base64_wide(values, node); + } else { + base64_wide(values, node); + } + } + } +#endif +#endif + + template + void forward_radix4(double *values, std::size_t quarter, std::size_t node, std::size_t first, + std::size_t last) const { + const double w2 = twiddles_[node]; + const double w = twiddles_[2 * node]; + const double iw = twiddles_[2 * node + 1]; + std::size_t i = first; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + if constexpr (UseFma) { + if (quarter >= 8) { + const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); + const __m512d packed_w2 = _mm512_set1_pd(w2); + const __m512d packed_w = _mm512_set1_pd(w); + const __m512d packed_iw = _mm512_set1_pd(iw); +#pragma GCC unroll 2 + for (; last - i >= 8; i += 8) { + __m512d v0 = _mm512_loadu_pd(values + i); + __m512d v1 = _mm512_loadu_pd(values + i + quarter); + __m512d v2 = _mm512_loadu_pd(values + i + 2 * quarter); + __m512d v3 = _mm512_loadu_pd(values + i + 3 * quarter); + forward_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw); + _mm512_storeu_pd(values + i, v0); + _mm512_storeu_pd(values + i + quarter, v1); + _mm512_storeu_pd(values + i + 2 * quarter, v2); + _mm512_storeu_pd(values + i + 3 * quarter, v3); + } + } + } +#endif +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { + if (quarter >= 4) { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + const __m256d packed_w2 = _mm256_set1_pd(w2); + const __m256d packed_w = _mm256_set1_pd(w); + const __m256d packed_iw = _mm256_set1_pd(iw); +#pragma GCC unroll 2 + for (; last - i >= 4; i += 4) { + __m256d v0 = _mm256_loadu_pd(values + i); + __m256d v1 = _mm256_loadu_pd(values + i + quarter); + __m256d v2 = _mm256_loadu_pd(values + i + 2 * quarter); + __m256d v3 = _mm256_loadu_pd(values + i + 3 * quarter); + forward_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw); + _mm256_storeu_pd(values + i, v0); + _mm256_storeu_pd(values + i + quarter, v1); + _mm256_storeu_pd(values + i + 2 * quarter, v2); + _mm256_storeu_pd(values + i + 3 * quarter, v3); + } + } + } +#endif + for (; i < last; ++i) { + forward_butterfly(arithmetic_, values[i], values[i + quarter], + values[i + 2 * quarter], values[i + 3 * quarter], w2, w, iw); + } + } + +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + template + void four_layer_block(double *values, std::size_t length, std::size_t node) const { + const std::size_t quarter = length / 4; + const std::size_t sixteenth = length / 16; + // Sixteen strips of 64 doubles occupy 8 KiB. Complete both radix-4 passes on + // each strip group while it is hot, before recursing into its sixteen children. + // This only changes visitation order; every butterfly retains the same inputs and bounds. + constexpr std::size_t strip_size = 64; + const auto pass = [&](double *base, std::size_t stride, std::size_t j, std::size_t first, + std::size_t last) { + if constexpr (Inverse) { + if (j == 0) { + inverse_radix4(base, stride, j, first, last); + } else { + inverse_radix4(base, stride, j, first, last); + } + } else { + if (j == 0) { + forward_radix4(base, stride, j, first, last); + } else { + forward_radix4(base, stride, j, first, last); + } + } + }; + for (std::size_t first = 0; first < sixteenth; first += strip_size) { + const std::size_t last = std::min(first + strip_size, sixteenth); + if constexpr (Inverse) { + for (std::size_t child = 0; child < 4; ++child) { + pass(values + child * quarter, sixteenth, 4 * node + child, first, last); + } + } + for (std::size_t part = 0; part < 4; ++part) { + pass(values, quarter, node, first + part * sixteenth, last + part * sixteenth); + } + if constexpr (!Inverse) { + for (std::size_t child = 0; child < 4; ++child) { + pass(values + child * quarter, sixteenth, 4 * node + child, first, last); + } + } + } + } +#endif + + template + void forward_recursive(double *values, std::size_t length, std::size_t node) const { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + if (length == 32 || length == 64) { + wide_basecase(values, length, node); + return; + } +#endif + if (length >= 8 && length <= 32) { + if (node == 0) { + basecase(values, length, node); + } else { + basecase(values, length, node); + } + return; + } +#if CRYPTO3_SMALL_PRIME_FFT_RADIX != 4 + constexpr std::size_t radix = CRYPTO3_SMALL_PRIME_FFT_RADIX; + if (length >= 16 * radix) { + const auto stride = length / radix; + if (node == 0) { + large_radix(values, stride, node); + } else { + large_radix(values, stride, node); + } + for (std::size_t child = 0; child < radix; ++child) { + forward_recursive(values + child * stride, stride, radix * node + child); + } + return; + } +#endif + if (length >= 4096) { + four_layer_block(values, length, node); + for (std::size_t child = 0; child < 16; ++child) { + forward_recursive(values + child * (length / 16), length / 16, 16 * node + child); + } + return; + } + } +#endif + if (length == 1) { + return; + } + if (length == 2) { + const double left = reduce(arithmetic_, values[0]); + const double right = node == 0 ? reduce(arithmetic_, values[1]) : + multiply(arithmetic_, values[1], twiddles_[node]); + values[0] = left + right; + values[1] = left - right; + return; + } + const std::size_t quarter = length / 4; + if (node == 0) { + forward_radix4(values, quarter, node, 0, quarter); + } else { + forward_radix4(values, quarter, node, 0, quarter); + } + if (quarter > 1) { + for (std::size_t child = 0; child < 4; ++child) { + forward_recursive(values + child * quarter, quarter, 4 * node + child); + } + } + } + + template + BOOST_FORCEINLINE static void inverse_butterfly(const Arithmetic &arithmetic, Value &v0, Value &v1, Value &v2, + Value &v3, Value w2, Value w, Value iw, bool leaf) { +#if defined(__clang__) +#pragma clang fp reassociate(off) +#pragma clang fp contract(off) +#endif + // Leaf inputs may be forward residues of magnitude < 3*q. Reduce them before + // adding four entries; non-leaf inputs already satisfy the inverse bound < 2*q. + const Value b0 = leaf ? reduce(arithmetic, v0) : v0; + const Value b1 = leaf ? reduce(arithmetic, v1) : v1; + const Value b2 = leaf ? reduce(arithmetic, v2) : v2; + const Value b3 = leaf ? reduce(arithmetic, v3) : v3; + const Value even_sum = b0 + b1; + const Value odd_sum = b2 + b3; + Value left = b0 - b1; + Value right = b2 - b3; + if constexpr (UnitTwiddle) { + left = reduce(arithmetic, left); + right = multiply(arithmetic, right, iw); + } else { + multiply_pair(arithmetic, left, right, w, iw); + } + // Four inputs < 2*q give sums < 8*q < 2^53. Times a centered weight this is + // < 4*q^2, so the a2 remainder is < 1.5*q. Left and right are each < q; + // their sum is < 2*q and their weighted difference is < q. The invariant closes. + v0 = reduce(arithmetic, even_sum + odd_sum); + v1 = left + right; + v2 = even_sum - odd_sum; + v3 = left - right; + if constexpr (UnitTwiddle) { + v2 = reduce(arithmetic, v2); + v3 = reduce(arithmetic, v3); + } else { + multiply_pair(arithmetic, v2, v3, w2, w2); + } + } + + template + void inverse_radix4(double *values, std::size_t quarter, std::size_t node, std::size_t first, + std::size_t last) const { + const double w2 = inverse_twiddle(node); + const double w = inverse_twiddle(2 * node); + const double iw = inverse_twiddle(2 * node + 1); + std::size_t i = first; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + if constexpr (UseFma) { + if (quarter >= 8) { + const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); + const __m512d packed_w2 = _mm512_set1_pd(w2); + const __m512d packed_w = _mm512_set1_pd(w); + const __m512d packed_iw = _mm512_set1_pd(iw); +#pragma GCC unroll 2 + for (; last - i >= 8; i += 8) { + __m512d v0 = _mm512_loadu_pd(values + i); + __m512d v1 = _mm512_loadu_pd(values + i + quarter); + __m512d v2 = _mm512_loadu_pd(values + i + 2 * quarter); + __m512d v3 = _mm512_loadu_pd(values + i + 3 * quarter); + inverse_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw, + false); + _mm512_storeu_pd(values + i, v0); + _mm512_storeu_pd(values + i + quarter, v1); + _mm512_storeu_pd(values + i + 2 * quarter, v2); + _mm512_storeu_pd(values + i + 3 * quarter, v3); + } + } + } +#endif +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { + if (quarter >= 4) { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + const __m256d packed_w2 = _mm256_set1_pd(w2); + const __m256d packed_w = _mm256_set1_pd(w); + const __m256d packed_iw = _mm256_set1_pd(iw); +#pragma GCC unroll 2 + for (; last - i >= 4; i += 4) { + __m256d v0 = _mm256_loadu_pd(values + i); + __m256d v1 = _mm256_loadu_pd(values + i + quarter); + __m256d v2 = _mm256_loadu_pd(values + i + 2 * quarter); + __m256d v3 = _mm256_loadu_pd(values + i + 3 * quarter); + inverse_butterfly(packed, v0, v1, v2, v3, packed_w2, packed_w, packed_iw, + false); + _mm256_storeu_pd(values + i, v0); + _mm256_storeu_pd(values + i + quarter, v1); + _mm256_storeu_pd(values + i + 2 * quarter, v2); + _mm256_storeu_pd(values + i + 3 * quarter, v3); + } + } + } +#endif + for (; i < last; ++i) { + inverse_butterfly(arithmetic_, values[i], values[i + quarter], + values[i + 2 * quarter], values[i + 3 * quarter], w2, w, iw, + quarter == 1 && !BoundedInput); + } + } + + template + void inverse_recursive(double *values, std::size_t length, std::size_t node) const { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + if (length == 32 || length == 64) { + wide_basecase(values, length, node); + return; + } +#endif + if (length >= 8 && length <= 32) { + if (node == 0) { + basecase(values, length, node); + } else { + basecase(values, length, node); + } + return; + } +#if CRYPTO3_SMALL_PRIME_FFT_RADIX != 4 + constexpr std::size_t radix = CRYPTO3_SMALL_PRIME_FFT_RADIX; + if (length >= 16 * radix) { + const auto stride = length / radix; + for (std::size_t child = 0; child < radix; ++child) { + inverse_recursive(values + child * stride, stride, radix * node + child); + } + if (node == 0) { + large_radix(values, stride, node); + } else { + large_radix(values, stride, node); + } + return; + } +#endif + if (length >= 4096) { + for (std::size_t child = 0; child < 16; ++child) { + inverse_recursive(values + child * (length / 16), length / 16, + 16 * node + child); + } + four_layer_block(values, length, node); + return; + } + } +#endif + if (length == 1) { + if constexpr (!BoundedInput) { + values[0] = reduce(arithmetic_, values[0]); + } + return; + } + if (length == 2) { + const double left = BoundedInput ? values[0] : reduce(arithmetic_, values[0]); + const double right = BoundedInput ? values[1] : reduce(arithmetic_, values[1]); + values[0] = reduce(arithmetic_, left + right); + values[1] = node == 0 ? reduce(arithmetic_, left - right) : + multiply(arithmetic_, left - right, inverse_twiddle(node)); + return; + } + const std::size_t quarter = length / 4; + if (quarter > 1) { + for (std::size_t child = 0; child < 4; ++child) { + inverse_recursive(values + child * quarter, quarter, 4 * node + child); + } + } + if (node == 0) { + inverse_radix4(values, quarter, node, 0, quarter); + } else { + inverse_radix4(values, quarter, node, 0, quarter); + } + } + + template + void inverse_impl(std::span values, bool normalize) const { + inverse_recursive(values.data(), size_, 0); + if (normalize) { + const double q = static_cast(modulus()); + std::size_t i = 0; +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX512_SELECTED) + if constexpr (UseFma) { + const small_prime_fft_arithmetic::wide_fma_kernel packed(arithmetic_); + const __m512d scale = _mm512_set1_pd(size_inverse_); + const __m512d packed_q = _mm512_set1_pd(q); + for (; values.size() - i >= 8; i += 8) { + __m512d value = packed.multiply_fma(_mm512_loadu_pd(values.data() + i), scale); + const auto negative = _mm512_cmp_pd_mask(value, _mm512_setzero_pd(), _CMP_LT_OQ); + value = _mm512_mask_add_pd(value, negative, value, packed_q); + _mm512_storeu_pd(values.data() + i, value); + } + } +#endif +#if defined(CRYPTO3_SMALL_PRIME_FFT_AVX2_SELECTED) + if constexpr (UseFma) { + if (values.size() >= 4) { + const small_prime_fft_arithmetic::packed_fma_kernel packed(arithmetic_); + const __m256d scale = _mm256_set1_pd(size_inverse_); + const __m256d packed_q = _mm256_set1_pd(q); + for (; values.size() - i >= 4; i += 4) { + __m256d value = packed.multiply_fma(_mm256_loadu_pd(values.data() + i), scale); + const __m256d negative = _mm256_cmp_pd(value, _mm256_setzero_pd(), _CMP_LT_OQ); + value = _mm256_add_pd(value, _mm256_and_pd(negative, packed_q)); + _mm256_storeu_pd(values.data() + i, value); + } + } + } +#endif + for (; i < values.size(); ++i) { + double &value = values[i]; + value = multiply(arithmetic_, value, size_inverse_); + value = value < 0 ? value + q : value; + } + } + } + + small_prime_fft_arithmetic arithmetic_; + std::size_t size_; + std::uint64_t omega_; + double size_inverse_; + std::vector twiddles_; + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_SMALL_PRIME_FFT_PLAN_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/factorization/complete_factorization.hpp b/libs/math/include/nil/crypto3/math/polynomial/factorization/complete_factorization.hpp index d122279e2..38920860d 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/factorization/complete_factorization.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/factorization/complete_factorization.hpp @@ -27,7 +27,6 @@ #include #include -#include #include #include @@ -35,6 +34,48 @@ namespace nil::crypto3::math { + namespace detail { + + /** + * Consume irreducible factors without retaining a result history. The callback receives each factor as an + * rvalue and may move it into owned storage; otherwise it is released when the callback returns. Its return + * value controls both equal-degree splitting and distinct-degree discovery. Square-free decomposition still + * finishes before any irreducible-factor callback, preserving arithmetic and shared-generator ordering. + */ + template + factorization_control + complete_factorization_impl(typename Backend::polynomial_type::value_type &leading_coefficient, + const typename Backend::polynomial_type &input, + polynomial_arithmetic::polynomial_context &arithmetic_context, + Generator &generator, FactorCallback &&factor_callback) { + using polynomial_type = typename Backend::polynomial_type; + + auto square_free_result = square_free_factorization(input, arithmetic_context); + leading_coefficient = square_free_result.leading_coefficient; + for (auto &square_free_factor : square_free_result.factors) { + const std::size_t multiplicity = square_free_factor.multiplicity; + const std::size_t block_size = kaltofen_shoup_block_size(square_free_factor.polynomial.size() - 1); + + const auto control = kaltofen_shoup_factor_monic_square_free( + nullptr, std::move(square_free_factor.polynomial), block_size, arithmetic_context, + [&](const distinct_degree_factor °ree_group) { + // Splitting consumes a copy; Kaltofen-Shoup still needs the original group for exact division + // if this callback continues. No previously processed groups are retained. + return factor_distinct_degree_group( + degree_group, arithmetic_context, generator, [&](polynomial_type &&factor) { + return factor_callback( + polynomial_factor {std::move(factor), multiplicity}); + }); + }); + if (control == factorization_control::stop_factorization) { + return control; + } + } + return factorization_control::continue_factorization; + } + + } // namespace detail + /** * Factor a polynomial into monic irreducible factors by composing the three finite-field factorization stages: * @@ -71,34 +112,14 @@ namespace nil::crypto3::math { using polynomial_type = typename Backend::polynomial_type; using result_type = polynomial_factorization_result; - auto square_free_result = square_free_factorization(input, arithmetic_context); - result_type result; - result.leading_coefficient = square_free_result.leading_coefficient; - for (auto &square_free_factor : square_free_result.factors) { - const std::size_t multiplicity = square_free_factor.multiplicity; - const std::size_t block_size = detail::kaltofen_shoup_block_size(square_free_factor.polynomial.size() - 1); - - std::vector> degree_groups; - const auto control = detail::kaltofen_shoup_factor_monic_square_free( - degree_groups, std::move(square_free_factor.polynomial), block_size, arithmetic_context, - [&](const distinct_degree_factor °ree_group) { - // Split each newly discovered degree group now and report its irreducible factors immediately. - // This lets the caller stop before Kaltofen-Shoup searches for later degree groups. - // The splitter takes its group by value, so passing degree_group intentionally copies it. - // Splitting consumes that copy; Kaltofen-Shoup still needs the original polynomial for exact - // division if this callback continues. - return detail::factor_distinct_degree_group( - degree_group, arithmetic_context, generator, [&](polynomial_type &&factor) { - result.factors.push_back({std::move(factor), multiplicity}); - return factor_callback(result.factors.back()); - }); - }); - if (control == factorization_control::stop_factorization) { - result.complete = false; - return result; - } - } + const auto control = + detail::complete_factorization_impl(result.leading_coefficient, input, arithmetic_context, + generator, [&](polynomial_factor &&factor) { + result.factors.push_back(std::move(factor)); + return factor_callback(result.factors.back()); + }); + result.complete = control == factorization_control::continue_factorization; return result; } diff --git a/libs/math/include/nil/crypto3/math/polynomial/factorization/kaltofen_shoup_distinct_degree_factorization.hpp b/libs/math/include/nil/crypto3/math/polynomial/factorization/kaltofen_shoup_distinct_degree_factorization.hpp index 309a221f5..ed9950dfa 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/factorization/kaltofen_shoup_distinct_degree_factorization.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/factorization/kaltofen_shoup_distinct_degree_factorization.hpp @@ -222,8 +222,9 @@ namespace nil::crypto3::math::detail { * After the lower degrees in this block have been removed, taking the GCD with coarse_block therefore extracts * exactly the factors of that degree. * - * The extracted groups are appended in increasing degree order. The callback is invoked immediately after each - * group is appended and may stop the split early. + * The extracted groups are reported in increasing degree order. If output is nonnull, each group is appended + * before invoking the callback. Otherwise the group is kept only until its callback and exact division finish. + * The callback borrows the group and may stop the split before that division. * * @return stop_factorization if the callback requests an early stop; continue_factorization otherwise. * @throws std::invalid_argument if first_factor_degree is zero, or if degree_count is zero or greater than the @@ -235,7 +236,7 @@ namespace nil::crypto3::math::detail { template requires DistinctDegreeFactorCallback factorization_control kaltofen_shoup_split_coarse_block( - std::vector> &output, + std::vector> *output, typename Backend::polynomial_type coarse_block, const typename Backend::polynomial_type &giant_step, std::size_t first_factor_degree, std::size_t degree_count, const kaltofen_shoup_frobenius_precomputation &precomputation, @@ -259,12 +260,13 @@ namespace nil::crypto3::math::detail { continue; } - output.push_back({std::move(factor), first_factor_degree + offset}); - if (factor_callback(output.back()) == factorization_control::stop_factorization) { + distinct_degree_factor group {std::move(factor), first_factor_degree + offset}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + if (factor_callback(reported) == factorization_control::stop_factorization) { return factorization_control::stop_factorization; } - factorization_exact_quotient(quotient, coarse_block, output.back().polynomial, arithmetic_context); + factorization_exact_quotient(quotient, coarse_block, reported.polynomial, arithmetic_context); coarse_block = std::move(quotient); } @@ -281,6 +283,9 @@ namespace nil::crypto3::math::detail { * a giant-step polynomial for every block while preserving increasing degree order and immediate callback stops. * Optional small-degree screening uses successive Frobenius differences before building the full baby table. * The powers and divisor context remain relative to the original input even after screened factors are removed. + * A nonnull output collects each group before invoking the callback. With a null output, processed groups are + * released after their callback and any remaining exact division. The callback borrows its group only for the + * duration of the call. * * Once the smallest unprocessed factor degree is greater than half the degree of the unclassified polynomial, * that polynomial is either constant or irreducible. Otherwise it would contain at least two factors of at least @@ -294,7 +299,7 @@ namespace nil::crypto3::math::detail { template requires DistinctDegreeFactorCallback factorization_control kaltofen_shoup_factor_monic_square_free( - std::vector> &output, + std::vector> *output, typename Backend::polynomial_type input, std::size_t block_size, polynomial_arithmetic::polynomial_context &arithmetic_context, FactorCallback &&factor_callback) { using polynomial_type = typename Backend::polynomial_type; @@ -305,8 +310,9 @@ namespace nil::crypto3::math::detail { } // A linear polynomial is already one degree-one factor, so no Frobenius precomputation is needed. if (input.size() == 2) { - output.push_back({std::move(input), 1}); - return factor_callback(output.back()); + distinct_degree_factor group {std::move(input), 1}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + return factor_callback(reported); } polynomial_frobenius_context frobenius_context(input, arithmetic_context); @@ -332,11 +338,12 @@ namespace nil::crypto3::math::detail { if (factor.size() > 1) { // All lower degrees have been removed, so this is an exact-degree group, not necessarily one // irreducible factor. Let the existing callback split it and apply multiplicity-aware recovery. - output.push_back({std::move(factor), degree}); - if (factor_callback(output.back()) == factorization_control::stop_factorization) + distinct_degree_factor group {std::move(factor), degree}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + if (factor_callback(reported) == factorization_control::stop_factorization) return factorization_control::stop_factorization; polynomial_type quotient; - factorization_exact_quotient(quotient, unclassified, output.back().polynomial, arithmetic_context); + factorization_exact_quotient(quotient, unclassified, reported.polynomial, arithmetic_context); unclassified = std::move(quotient); } } @@ -344,8 +351,9 @@ namespace nil::crypto3::math::detail { if (degree > (unclassified.size() - 1) / 2) { if (unclassified.size() > 1) { const auto remaining_degree = unclassified.size() - 1; - output.push_back({std::move(unclassified), remaining_degree}); - return factor_callback(output.back()); + distinct_degree_factor group {std::move(unclassified), remaining_degree}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + return factor_callback(reported); } return factorization_control::continue_factorization; } @@ -374,8 +382,9 @@ namespace nil::crypto3::math::detail { // Every irreducible factor in this block has degree at least first_factor_degree. If the block's total // degree is less than twice that bound, it cannot contain two factors and is itself irreducible. if (coarse_block_degree / 2 < first_factor_degree) { - output.push_back({std::move(coarse_block), coarse_block_degree}); - if (factor_callback(output.back()) == factorization_control::stop_factorization) { + distinct_degree_factor group {std::move(coarse_block), coarse_block_degree}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + if (factor_callback(reported) == factorization_control::stop_factorization) { return factorization_control::stop_factorization; } } else if (kaltofen_shoup_split_coarse_block(output, std::move(coarse_block), giant_step, @@ -396,8 +405,9 @@ namespace nil::crypto3::math::detail { if (unclassified.size() > 1) { const std::size_t irreducible_factor_degree = unclassified.size() - 1; - output.push_back({std::move(unclassified), irreducible_factor_degree}); - if (factor_callback(output.back()) == factorization_control::stop_factorization) { + distinct_degree_factor group {std::move(unclassified), irreducible_factor_degree}; + const auto &reported = output ? output->emplace_back(std::move(group)) : group; + if (factor_callback(reported) == factorization_control::stop_factorization) { return factorization_control::stop_factorization; } } @@ -438,7 +448,7 @@ namespace nil::crypto3::math { } const std::size_t block_size = detail::kaltofen_shoup_block_size(monic_input.size() - 1); - if (detail::kaltofen_shoup_factor_monic_square_free(result.factors, std::move(monic_input), block_size, + if (detail::kaltofen_shoup_factor_monic_square_free(&result.factors, std::move(monic_input), block_size, arithmetic_context, factor_callback) == factorization_control::stop_factorization) { result.complete = false; diff --git a/libs/math/include/nil/crypto3/math/polynomial/operations/basic_operations.hpp b/libs/math/include/nil/crypto3/math/polynomial/operations/basic_operations.hpp index 479ee1f69..5263fb5c3 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/operations/basic_operations.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/operations/basic_operations.hpp @@ -64,6 +64,16 @@ namespace nil { std::default_initializable && std::equality_comparable; + // Canonical coefficient storage is nonempty, with no trailing zero except for the zero polynomial [0]. + template + requires std::default_initializable && + std::equality_comparable + bool is_canonical_polynomial(const Polynomial &polynomial) { + using value_type = typename Polynomial::value_type; + return polynomial.size() != 0 && + (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); + } + } // namespace detail /** diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/detail/prepared_composition_coefficients.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/detail/prepared_composition_coefficients.hpp new file mode 100644 index 000000000..03ae91e6f --- /dev/null +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/detail/prepared_composition_coefficients.hpp @@ -0,0 +1,336 @@ +//---------------------------------------------------------------------------// +// 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_DETAIL_PREPARED_COMPOSITION_COEFFICIENTS_HPP +#define CRYPTO3_MATH_DETAIL_PREPARED_COMPOSITION_COEFFICIENTS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace nil::crypto3::math::detail { + + template + class prepared_composition_coefficients { + public: + static constexpr bool supported = false; + }; + + template + class prepared_composition_coefficients> { + using params = typename Value::field_type::extension_policy::fast_params; + using kernel = algebra::fields::detail::fp12_fast::fp12_fast; + using prepared_type = typename kernel::prepared_product_operand; + using simd_kernel = typename kernel::rank23; + using limb = algebra::fields::detail::fp12_fast::limb; + + public: + static constexpr bool supported = true; + static constexpr std::size_t block_tile_size = 4; + static constexpr std::size_t coefficient_tile_size = 4; + + // Scratch belongs to a composition call, so concurrent readers never mutate the shared cache. + struct workspace { + std::array, block_tile_size> scalars; + std::array, block_tile_size> exponents; + std::array, block_tile_size> simd_scalars; + }; + + prepared_composition_coefficients() = default; + + template + prepared_composition_coefficients(const std::vector &powers, std::size_t byte_limit, + std::size_t simd_lanes = 8) { + if (simd_lanes != 0 && simd_lanes != 4 && simd_lanes != 8) { + throw std::invalid_argument("composition SIMD lanes must be 0, 4 or 8"); + } + if (byte_limit == 0 || powers.size() < 2) { + return; + } + const std::size_t terms = powers.size() - 1; // The constant power is handled by direct addition. + std::size_t width = 0; + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + width = std::max(width, powers[exponent].size()); + } +#if defined(CRYPTO3_HAS_AVX512IFMA) + if constexpr (simd_kernel::supported) { + // Cache dimensions rule out short dots and narrow polynomials. use_for_block also checks + // the actual nonzero outer terms before selecting either prepared kernel. + if (simd_lanes != 0 && width >= 4 && terms >= minimum_prepared_terms && + prepare_simd(powers, width, terms, byte_limit, simd_lanes)) { + return; + } + } +#endif + const std::size_t capacity = std::min(byte_limit / sizeof(prepared_type), coefficients_.max_size()); + // Division before multiplication checks both the budget and size_t overflow without allocating. + if (width == 0 || width > capacity / terms) { + return; + } + coefficients_.resize(width * terms); + for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + if (coefficient < powers[exponent].size()) { + coefficients_[coefficient * terms + exponent - 1] = + kernel::prepare_product_operand(powers[exponent][coefficient]); + } + // Default-constructed backends are zero, padding shorter powers without borrowed views. + } + } + terms_ = terms; + width_ = width; + } + + bool available() const { + return !coefficients_.empty() || !simd_coefficients_.empty(); + } + + std::size_t bytes() const { + return coefficients_.size() * sizeof(prepared_type) + simd_coefficients_.size() * sizeof(limb); + } + + std::size_t simd_lanes() const { + return simd_lanes_; + } + + // Short or sparse sums reuse ordinary Fp12 multiplication instead of paying for operand preparation + // and recombination. The constant coefficient needs only an addition and does not count as a term. + template + bool use_for_block(const Polynomial &outer, std::size_t block_index) const { + if (!available()) { + return false; + } + const std::size_t block_size = terms_ + 1; + const std::size_t first = block_index * block_size; + assert(first < outer.size()); + const std::size_t count = std::min(block_size, outer.size() - first); + if (count <= minimum_prepared_terms) { + return false; + } + std::size_t active_terms = 0; + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero() && ++active_terms == minimum_prepared_terms) { + return true; + } + } + return false; + } + + // Evaluate up to four adjacent outer blocks. Each scalar is prepared once and reused across all + // output coefficients. The coefficient-major cache keeps each dot contiguous; tiles reuse the same + // cached columns across adjacent blocks without storing a matrix of every block's output. + template + void evaluate_blocks(std::array &output, const Polynomial &outer, + std::size_t first_block, std::size_t block_count, workspace &scratch) const { + assert(available() && block_count > 0 && block_count <= block_tile_size); +#if defined(CRYPTO3_HAS_AVX512IFMA) + if constexpr (simd_kernel::supported) { + if (simd_lanes_ == 4) { + evaluate_simd_blocks<4>(output, outer, first_block, block_count, scratch); + return; + } + if (simd_lanes_ == 8) { + evaluate_simd_blocks<8>(output, outer, first_block, block_count, scratch); + return; + } + } +#endif + const std::size_t block_size = terms_ + 1; + for (std::size_t block = 0; block < block_count; ++block) { + const std::size_t first = (first_block + block) * block_size; + assert(first < outer.size()); + const std::size_t count = std::min(block_size, outer.size() - first); + auto &scalars = scratch.scalars[block]; + auto &exponents = scratch.exponents[block]; + scalars.clear(); + exponents.clear(); + scalars.reserve(count - 1); + exponents.reserve(count - 1); + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero()) { + scalars.push_back(kernel::prepare_product_operand(outer[first + exponent])); + exponents.push_back(exponent - 1); + } + } + output[block].assign(scalars.empty() ? 1 : width_, Value::zero()); + output[block][0] = outer[first]; + } + for (std::size_t first_coefficient = 0; first_coefficient < width_;) { + const std::size_t end = first_coefficient + std::min(coefficient_tile_size, width_ - first_coefficient); + for (std::size_t coefficient = first_coefficient; coefficient < end; ++coefficient) { + const auto select_power = [this, coefficient](std::size_t exponent) -> const prepared_type & { + return coefficients_[coefficient * terms_ + exponent]; + }; + for (std::size_t block = 0; block < block_count; ++block) { + const auto &scalars = scratch.scalars[block]; + if (!scalars.empty()) { + const auto powers = scratch.exponents[block] | std::views::transform(select_power); + output[block][coefficient] += + kernel::template inner_product(scalars.begin(), scalars.end(), powers.begin()); + } + } + } + first_coefficient = end; + } + for (std::size_t block = 0; block < block_count; ++block) { + condense(output[block]); + } + } + + private: + // Conservative dispatch heuristic shared with the SIMD cache cutoff, not an arithmetic bound. + static constexpr std::size_t minimum_prepared_terms = 8; +#if defined(CRYPTO3_HAS_AVX512IFMA) + template + bool prepare_simd(const std::vector &powers, std::size_t width, std::size_t terms, + std::size_t byte_limit, std::size_t lanes) { + using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; + constexpr std::size_t coordinates = simd_kernel::product_coordinate_count; + const std::size_t panel_count = 1 + (width - 1) / lanes; + const std::size_t stride = coordinates * 5 * lanes; + const std::size_t capacity = std::min(byte_limit / sizeof(limb), simd_coefficients_.max_size()); + if (terms > capacity / stride || panel_count > capacity / stride / terms) { + return false; + } + simd_coefficients_.resize(panel_count * stride * terms); + // Owned [panel][coordinate][term][digit][lane] storage. Only this representation is retained; + // 40 bytes per form replace the scalar backend object, and padding is included in the budget. + for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { + const std::size_t panel = coefficient / lanes, lane = coefficient % lanes; + for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { + if (coefficient >= powers[exponent].size()) { + continue; + } + const auto prepared = simd_kernel::prepare_product_operand(powers[exponent][coefficient]); + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const auto packed = radix::from_montgomery64(prepared[coordinate]); + const std::size_t offset = ((panel * coordinates + coordinate) * terms + exponent - 1) * 5; + for (std::size_t digit = 0; digit < 5; ++digit) { + simd_coefficients_[(offset + digit) * lanes + lane] = packed[digit]; + } + } + } + } + terms_ = terms; + width_ = width; + simd_lanes_ = lanes; + return true; + } + + template + void evaluate_simd_blocks(std::array &output, const Polynomial &outer, + std::size_t first_block, std::size_t block_count, workspace &scratch) const { + namespace fast = algebra::fields::detail::fp12_fast; + using radix = fast::montgomery_radix52; + constexpr std::size_t coordinates = simd_kernel::product_coordinate_count; + bool any_active = false; + for (std::size_t block = 0; block < block_count; ++block) { + const std::size_t first = (first_block + block) * (terms_ + 1); + assert(first < outer.size()); + const std::size_t count = std::min(terms_ + 1, outer.size() - first); + auto &indices = scratch.exponents[block]; + indices.clear(); + for (std::size_t exponent = 1; exponent < count; ++exponent) { + if (!outer[first + exponent].is_zero()) { + indices.push_back(exponent - 1); + } + } + auto &scalars = scratch.simd_scalars[block]; + if (indices.size() > scalars.max_size() / (coordinates * 5)) { + throw std::length_error("composition SIMD scratch exceeds maximum size"); + } + scalars.resize(indices.size() * coordinates * 5); + for (std::size_t term = 0; term < indices.size(); ++term) { + const auto prepared = simd_kernel::prepare_product_operand(outer[first + indices[term] + 1]); + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const auto packed = radix::from_montgomery64(prepared[coordinate]); + std::copy(packed.begin(), packed.end(), + scalars.begin() + (coordinate * indices.size() + term) * 5); + } + } + any_active = any_active || !indices.empty(); + output[block].assign(indices.empty() ? 1 : width_, Value::zero()); + output[block][0] = outer[first]; + } + for (std::size_t first = 0; any_active && first < width_;) { + const std::size_t active_lanes = std::min(Lanes, width_ - first); + // Keep the dot products in canonical radix-52 form until tower recombination. + std::array, block_tile_size> products; + for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { + const limb *right = + simd_coefficients_.data() + ((first / Lanes * coordinates + coordinate) * terms_) * 5 * Lanes; + // Reuse this coordinate panel across neighboring blocks before moving to the next one. + for (std::size_t block = 0; block < block_count; ++block) { + const auto &indices = scratch.exponents[block]; + if (indices.empty()) { + continue; + } + fast::montgomery_inner_product52>( + products[block].data() + coordinate * 5 * Lanes, + scratch.simd_scalars[block].data() + coordinate * indices.size() * 5, right, indices.data(), + indices.size()); + } + } + for (std::size_t block = 0; block < block_count; ++block) { + if (!scratch.exponents[block].empty()) { + std::array combined; + simd_kernel::template recombine_products52>( + combined.data(), products[block].data()); + for (std::size_t lane = 0; lane < active_lanes; ++lane) { + Value value; + for (std::size_t coordinate = 0; coordinate < 12; ++coordinate) { + typename radix::value_type packed; + for (std::size_t digit = 0; digit < 5; ++digit) { + packed[digit] = combined[(coordinate * 5 + digit) * Lanes + lane]; + } + radix::to_montgomery64(value.coordinate(coordinate).data.backend().base_data(), packed); + } + output[block][first + lane] += value; + } + } + } + first += active_lanes; + } + for (std::size_t block = 0; block < block_count; ++block) { + condense(output[block]); + } + } +#endif + std::size_t terms_ = 0, width_ = 0, simd_lanes_ = 0; + std::vector coefficients_; + std::vector simd_coefficients_; + }; + +} // namespace nil::crypto3::math::detail + +#endif // CRYPTO3_MATH_DETAIL_PREPARED_COMPOSITION_COEFFICIENTS_HPP diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp index 957a5898c..ca110b07f 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_composition.hpp @@ -27,316 +27,17 @@ #include #include -#include #include -#include #include -#include #include #include -#include #include +#include #include namespace nil::crypto3::math { - namespace detail { - template - class prepared_composition_coefficients { - public: - static constexpr bool supported = false; - }; - - template - class prepared_composition_coefficients< - Value, std::void_t> { - using params = typename Value::field_type::extension_policy::fast_params; - using kernel = algebra::fields::detail::fp12_fast::fp12_fast; - using prepared_type = typename kernel::prepared_product_operand; - using simd_kernel = typename kernel::rank23; - using limb = algebra::fields::detail::fp12_fast::limb; - - public: - static constexpr bool supported = true; - static constexpr std::size_t block_tile_size = 4; - static constexpr std::size_t coefficient_tile_size = 4; - - // Scratch belongs to a composition call, so concurrent readers never mutate the shared cache. - struct workspace { - std::array, block_tile_size> scalars; - std::array, block_tile_size> exponents; - std::array, block_tile_size> simd_scalars; - }; - - prepared_composition_coefficients() = default; - - template - prepared_composition_coefficients(const std::vector &powers, std::size_t byte_limit, - std::size_t simd_lanes = 8) { - if (simd_lanes != 0 && simd_lanes != 4 && simd_lanes != 8) { - throw std::invalid_argument("composition SIMD lanes must be 0, 4 or 8"); - } - if (byte_limit == 0 || powers.size() < 2) { - return; - } - const std::size_t terms = powers.size() - 1; // The constant power is handled by direct addition. - std::size_t width = 0; - for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { - width = std::max(width, powers[exponent].size()); - } -#if defined(CRYPTO3_HAS_AVX512IFMA) - if constexpr (simd_kernel::supported) { - // Cache dimensions rule out short dots and narrow polynomials. use_for_block also checks - // the actual nonzero outer terms before selecting either prepared kernel. - if (simd_lanes != 0 && width >= 4 && terms >= minimum_prepared_terms && - prepare_simd(powers, width, terms, byte_limit, simd_lanes)) { - return; - } - } -#endif - const std::size_t capacity = std::min(byte_limit / sizeof(prepared_type), coefficients_.max_size()); - // Division before multiplication checks both the budget and size_t overflow without allocating. - if (width == 0 || width > capacity / terms) { - return; - } - coefficients_.resize(width * terms); - for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { - for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { - if (coefficient < powers[exponent].size()) { - coefficients_[coefficient * terms + exponent - 1] = - kernel::prepare_product_operand(powers[exponent][coefficient]); - } - // Default-constructed backends are zero, padding shorter powers without borrowed views. - } - } - terms_ = terms; - width_ = width; - } - - bool available() const { - return !coefficients_.empty() || !simd_coefficients_.empty(); - } - - std::size_t bytes() const { - return coefficients_.size() * sizeof(prepared_type) + simd_coefficients_.size() * sizeof(limb); - } - - std::size_t simd_lanes() const { - return simd_lanes_; - } - - // Short or sparse sums reuse ordinary Fp12 multiplication instead of paying for operand preparation - // and recombination. The constant coefficient needs only an addition and does not count as a term. - template - bool use_for_block(const Polynomial &outer, std::size_t block_index) const { - if (!available()) { - return false; - } - const std::size_t block_size = terms_ + 1; - const std::size_t first = block_index * block_size; - assert(first < outer.size()); - const std::size_t count = std::min(block_size, outer.size() - first); - if (count <= minimum_prepared_terms) { - return false; - } - std::size_t active_terms = 0; - for (std::size_t exponent = 1; exponent < count; ++exponent) { - if (!outer[first + exponent].is_zero() && ++active_terms == minimum_prepared_terms) { - return true; - } - } - return false; - } - - // Evaluate up to four adjacent outer blocks. Each scalar is prepared once and reused across all - // output coefficients. The coefficient-major cache keeps each dot contiguous; tiles reuse the same - // cached columns across adjacent blocks without storing a matrix of every block's output. - template - void evaluate_blocks(std::array &output, const Polynomial &outer, - std::size_t first_block, std::size_t block_count, workspace &scratch) const { - assert(available() && block_count > 0 && block_count <= block_tile_size); -#if defined(CRYPTO3_HAS_AVX512IFMA) - if constexpr (simd_kernel::supported) { - if (simd_lanes_ == 4) { - evaluate_simd_blocks<4>(output, outer, first_block, block_count, scratch); - return; - } - if (simd_lanes_ == 8) { - evaluate_simd_blocks<8>(output, outer, first_block, block_count, scratch); - return; - } - } -#endif - const std::size_t block_size = terms_ + 1; - for (std::size_t block = 0; block < block_count; ++block) { - const std::size_t first = (first_block + block) * block_size; - assert(first < outer.size()); - const std::size_t count = std::min(block_size, outer.size() - first); - auto &scalars = scratch.scalars[block]; - auto &exponents = scratch.exponents[block]; - scalars.clear(); - exponents.clear(); - scalars.reserve(count - 1); - exponents.reserve(count - 1); - for (std::size_t exponent = 1; exponent < count; ++exponent) { - if (!outer[first + exponent].is_zero()) { - scalars.push_back(kernel::prepare_product_operand(outer[first + exponent])); - exponents.push_back(exponent - 1); - } - } - output[block].assign(scalars.empty() ? 1 : width_, Value::zero()); - output[block][0] = outer[first]; - } - for (std::size_t first_coefficient = 0; first_coefficient < width_;) { - const std::size_t end = - first_coefficient + std::min(coefficient_tile_size, width_ - first_coefficient); - for (std::size_t coefficient = first_coefficient; coefficient < end; ++coefficient) { - const auto select_power = [this, coefficient](std::size_t exponent) -> const prepared_type & { - return coefficients_[coefficient * terms_ + exponent]; - }; - for (std::size_t block = 0; block < block_count; ++block) { - const auto &scalars = scratch.scalars[block]; - if (!scalars.empty()) { - const auto powers = scratch.exponents[block] | std::views::transform(select_power); - output[block][coefficient] += kernel::template inner_product( - scalars.begin(), scalars.end(), powers.begin()); - } - } - } - first_coefficient = end; - } - for (std::size_t block = 0; block < block_count; ++block) { - condense(output[block]); - } - } - - private: - // Conservative dispatch heuristic shared with the SIMD cache cutoff, not an arithmetic bound. - static constexpr std::size_t minimum_prepared_terms = 8; -#if defined(CRYPTO3_HAS_AVX512IFMA) - template - bool prepare_simd(const std::vector &powers, std::size_t width, std::size_t terms, - std::size_t byte_limit, std::size_t lanes) { - using radix = algebra::fields::detail::fp12_fast::montgomery_radix52; - constexpr std::size_t coordinates = simd_kernel::product_coordinate_count; - const std::size_t panel_count = 1 + (width - 1) / lanes; - const std::size_t stride = coordinates * 5 * lanes; - const std::size_t capacity = std::min(byte_limit / sizeof(limb), simd_coefficients_.max_size()); - if (terms > capacity / stride || panel_count > capacity / stride / terms) { - return false; - } - simd_coefficients_.resize(panel_count * stride * terms); - // Owned [panel][coordinate][term][digit][lane] storage. Only this representation is retained; - // 40 bytes per form replace the scalar backend object, and padding is included in the budget. - for (std::size_t coefficient = 0; coefficient < width; ++coefficient) { - const std::size_t panel = coefficient / lanes, lane = coefficient % lanes; - for (std::size_t exponent = 1; exponent < powers.size(); ++exponent) { - if (coefficient >= powers[exponent].size()) { - continue; - } - const auto prepared = simd_kernel::prepare_product_operand(powers[exponent][coefficient]); - for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { - const auto packed = radix::from_montgomery64(prepared[coordinate]); - const std::size_t offset = ((panel * coordinates + coordinate) * terms + exponent - 1) * 5; - for (std::size_t digit = 0; digit < 5; ++digit) { - simd_coefficients_[(offset + digit) * lanes + lane] = packed[digit]; - } - } - } - } - terms_ = terms; - width_ = width; - simd_lanes_ = lanes; - return true; - } - - template - void evaluate_simd_blocks(std::array &output, const Polynomial &outer, - std::size_t first_block, std::size_t block_count, workspace &scratch) const { - namespace fast = algebra::fields::detail::fp12_fast; - using radix = fast::montgomery_radix52; - constexpr std::size_t coordinates = simd_kernel::product_coordinate_count; - bool any_active = false; - for (std::size_t block = 0; block < block_count; ++block) { - const std::size_t first = (first_block + block) * (terms_ + 1); - assert(first < outer.size()); - const std::size_t count = std::min(terms_ + 1, outer.size() - first); - auto &indices = scratch.exponents[block]; - indices.clear(); - for (std::size_t exponent = 1; exponent < count; ++exponent) { - if (!outer[first + exponent].is_zero()) { - indices.push_back(exponent - 1); - } - } - auto &scalars = scratch.simd_scalars[block]; - if (indices.size() > scalars.max_size() / (coordinates * 5)) { - throw std::length_error("composition SIMD scratch exceeds maximum size"); - } - scalars.resize(indices.size() * coordinates * 5); - for (std::size_t term = 0; term < indices.size(); ++term) { - const auto prepared = simd_kernel::prepare_product_operand(outer[first + indices[term] + 1]); - for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { - const auto packed = radix::from_montgomery64(prepared[coordinate]); - std::copy(packed.begin(), packed.end(), - scalars.begin() + (coordinate * indices.size() + term) * 5); - } - } - any_active = any_active || !indices.empty(); - output[block].assign(indices.empty() ? 1 : width_, Value::zero()); - output[block][0] = outer[first]; - } - for (std::size_t first = 0; any_active && first < width_;) { - const std::size_t active_lanes = std::min(Lanes, width_ - first); - // Keep the dot products in canonical radix-52 form until tower recombination. - std::array, block_tile_size> products; - for (std::size_t coordinate = 0; coordinate < coordinates; ++coordinate) { - const limb *right = simd_coefficients_.data() + - ((first / Lanes * coordinates + coordinate) * terms_) * 5 * Lanes; - // Reuse this coordinate panel across neighboring blocks before moving to the next one. - for (std::size_t block = 0; block < block_count; ++block) { - const auto &indices = scratch.exponents[block]; - if (indices.empty()) { - continue; - } - fast::montgomery_inner_product52>( - products[block].data() + coordinate * 5 * Lanes, - scratch.simd_scalars[block].data() + coordinate * indices.size() * 5, right, - indices.data(), indices.size()); - } - } - for (std::size_t block = 0; block < block_count; ++block) { - if (!scratch.exponents[block].empty()) { - std::array combined; - simd_kernel::template recombine_products52>( - combined.data(), products[block].data()); - for (std::size_t lane = 0; lane < active_lanes; ++lane) { - Value value; - for (std::size_t coordinate = 0; coordinate < 12; ++coordinate) { - typename radix::value_type packed; - for (std::size_t digit = 0; digit < 5; ++digit) { - packed[digit] = combined[(coordinate * 5 + digit) * Lanes + lane]; - } - radix::to_montgomery64(value.coordinate(coordinate).data.backend().base_data(), - packed); - } - output[block][first + lane] += value; - } - } - } - first += active_lanes; - } - for (std::size_t block = 0; block < block_count; ++block) { - condense(output[block]); - } - } -#endif - std::size_t terms_ = 0, width_ = 0, simd_lanes_ = 0; - std::vector coefficients_; - std::vector simd_coefficients_; - }; - } // namespace detail - /** * Reusable baby-step precomputation for Brent-Kung composition with a fixed inner polynomial and divisor. For a * block size k, the object stores the reduced powers diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_exponentiation.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_exponentiation.hpp index 017e992ad..99752ff19 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_exponentiation.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_exponentiation.hpp @@ -32,7 +32,6 @@ #include #include -#include #include namespace nil::crypto3::math { @@ -167,20 +166,7 @@ namespace nil::crypto3::math { // X^(2*k+b): square result, then multiply by X only when b is one. squaremod(result, result, divisor_context, arithmetic_context); if (bits.back()) { - if (result.size() == divisor_degree) { - // With d = divisor_degree, X*result has degree d and leading coefficient result[d-1]. - // Subtract c*H once, with c = result[d-1]/lead(H). Store its lower d coefficients directly, - // descending to preserve the old result[i-1]; the canceled X^d coefficient needs no storage. - const value_type c = result[divisor_degree - 1] * inverse_lead; - for (std::size_t i = divisor_degree - 1; i > 0; --i) { - result[i] = result[i - 1] - c * modulus[i]; - } - result[0] = -c * modulus[0]; - condense(result); - } else { - // The shift stays below degree(H); this helper also preserves canonical zero. - shift_left(result, result, 1); - } + detail::multiply_by_x_mod(result, divisor_context); } bits.pop_back(); } diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp index 90d0bce30..e760b5525 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_frobenius.hpp @@ -163,10 +163,7 @@ namespace nil::crypto3::math { using value_type = typename polynomial_type::value_type; using cache_type = polynomial_composition_precomputation; const std::size_t degree = divisor_context.degree(); - const auto canonical = [](const polynomial_type &value) { - return !value.empty() && (value.size() == 1 || value.back() != value_type::zero()); - }; - if (!canonical(input) || !canonical(x_to_frobenius_power) || + if (!is_canonical_polynomial(input) || !is_canonical_polynomial(x_to_frobenius_power) || x_to_frobenius_power.size() > std::max(1, degree)) { throw std::invalid_argument( "Frobenius products require canonical inputs and a reduced Frobenius image"); @@ -290,10 +287,8 @@ namespace nil::crypto3::math { using value_type = typename polynomial_type::value_type; using cache_type = polynomial_composition_precomputation; const std::size_t degree = divisor_context.degree(); - const auto canonical = [](const polynomial_type &value) { - return !value.empty() && (value.size() == 1 || value.back() != value_type::zero()); - }; - if (!canonical(input) || !canonical(weight) || !canonical(x_to_frobenius_power) || + if (!detail::is_canonical_polynomial(input) || !detail::is_canonical_polynomial(weight) || + !detail::is_canonical_polynomial(x_to_frobenius_power) || x_to_frobenius_power.size() > std::max(1, degree)) { throw std::invalid_argument("Frobenius sums require canonical inputs and a reduced Frobenius image"); } diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_modular_arithmetic.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_modular_arithmetic.hpp index 3f29a5976..300e6bc12 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_modular_arithmetic.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_modular_arithmetic.hpp @@ -25,10 +25,42 @@ #ifndef CRYPTO3_MATH_POLYNOMIAL_MODULAR_ARITHMETIC_HPP #define CRYPTO3_MATH_POLYNOMIAL_MODULAR_ARITHMETIC_HPP +#include +#include + #include +#include namespace nil::crypto3::math { + namespace detail { + + // Multiply a canonical reduced representative by X modulo a nonconstant divisor B, in place. + // Reuse 1/lead(B) from the divisor context; no polynomial product or division is needed. + template + void multiply_by_x_mod(typename Backend::polynomial_type &value, + const polynomial_divisor_context &divisor_context) { + using value_type = typename Backend::polynomial_type::value_type; + const std::size_t degree = divisor_context.degree(); + assert(degree > 0 && value.size() <= degree && is_canonical_polynomial(value)); + if (value.size() == degree) { + const auto &modulus = divisor_context.divisor(); + const value_type c = value[degree - 1] * divisor_context.reversed_divisor_inverse()[0]; + // Subtract c*B from X*value. Descending writes preserve the old value[i-1], and the + // canceled X^degree coefficient needs no storage. + for (std::size_t i = degree - 1; i > 0; --i) { + value[i] = value[i - 1] - c * modulus[i]; + } + value[0] = -c * modulus[0]; + condense(value); + } else { + // The shift stays below degree(B), including canonical zero. + shift_left(value, value, 1); + } + } + + } // namespace detail + /** * Compute output = (left * right) mod B, where B is the nonzero polynomial stored in divisor_context. B need not * be monic or irreducible. The result is canonical, has degree less than degree(B) when B is nonconstant, and may diff --git a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp index aaf066db1..c8d3942c8 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/quotient_ring/polynomial_square_root.hpp @@ -36,7 +36,6 @@ #include #include -#include #include #include @@ -141,8 +140,7 @@ namespace nil::crypto3::math { } void require_canonical_reduced(const polynomial_type &candidate) const { - if (candidate.empty() || candidate.size() > divisor_context().degree() || - (candidate.size() > 1 && candidate.back() == value_type::zero())) { + if (candidate.size() > divisor_context().degree() || !detail::is_canonical_polynomial(candidate)) { throw std::invalid_argument( "a polynomial square-root nonresidue candidate must be a nonempty canonical reduced " "representative"); @@ -160,30 +158,6 @@ namespace nil::crypto3::math { polynomial_type non_residue_to_odd_order_; }; - /** - * Return whether input is a square in K[X]/(B), where K is a finite field and B is the irreducible polynomial - * stored in divisor_context. If Q is the order of K and d = degree(B), the quotient is a field of order Q^d. - * For a nonzero residue a in odd characteristic, Euler's criterion gives - * - * a^((Q^d - 1) / 2) = 1 - * - * exactly when a is a square. Zero is a square. In characteristic two, squaring is an automorphism of every - * finite field, so every residue is a square. - * - * The canonical indeterminate X uses the closed-form field norm - * - * Norm(X) = (-1)^d B(0) / leading_coefficient(B). - * - * An element of a finite extension of an odd-order field is a square exactly when its norm is a square in the - * coefficient field. When the coefficient value type provides is_square(), this avoids quotient-ring - * exponentiation for X. Other inputs use Euler's criterion as a reference implementation. - * - * The irreducibility precondition is essential: for a reducible B the quotient has zero divisors and Euler's - * criterion does not characterize its squares. The function does not repeat an irreducibility test. - * - * @throws std::invalid_argument if B is constant or input is not a reduced quotient-field representative. - * @pre B is irreducible and input is a nonempty canonical coefficient polynomial. - */ namespace detail { template @@ -243,6 +217,30 @@ namespace nil::crypto3::math { } // namespace detail + /** + * Return whether input is a square in K[X]/(B), where K is a finite field and B is the irreducible polynomial + * stored in divisor_context. If Q is the order of K and d = degree(B), the quotient is a field of order Q^d. + * For a nonzero residue a in odd characteristic, Euler's criterion gives + * + * a^((Q^d - 1) / 2) = 1 + * + * exactly when a is a square. Zero is a square. In characteristic two, squaring is an automorphism of every + * finite field, so every residue is a square. + * + * The canonical indeterminate X uses the closed-form field norm + * + * Norm(X) = (-1)^d B(0) / leading_coefficient(B). + * + * An element of a finite extension of an odd-order field is a square exactly when its norm is a square in the + * coefficient field. When the coefficient value type provides is_square(), this avoids quotient-ring + * exponentiation for X. Other inputs use Euler's criterion as a reference implementation. + * + * The irreducibility precondition is essential: for a reducible B the quotient has zero divisors and Euler's + * criterion does not characterize its squares. The function does not repeat an irreducibility test. + * + * @throws std::invalid_argument if B is constant or input is not a reduced quotient-field representative. + * @pre B is irreducible and input is a nonempty canonical coefficient polynomial. + */ template requires algebra::FieldValue bool is_square_mod(const typename Backend::polynomial_type &input, @@ -253,23 +251,6 @@ namespace nil::crypto3::math { namespace detail { - // Multiply a canonical reduced representative by X, reducing the single possible leading term. - template - void square_root_multiply_by_x(typename Backend::polynomial_type &value, - const polynomial_divisor_context &divisor_context) { - using value_type = typename Backend::polynomial_type::value_type; - const std::size_t degree = divisor_context.degree(); - shift_left(value, value, 1); - if (value.size() > degree) { - const value_type leading = value[degree] * divisor_context.reversed_divisor_inverse()[0]; - value.resize(degree); - for (std::size_t i = 0; i < degree; ++i) { - value[i] -= leading * divisor_context.divisor()[i]; - } - condense(value); - } - } - // Return {X^((Q-1)/2), X^Q} modulo B. Called after the root checks, with d > 1 and B(0) != 0. template std::pair @@ -290,7 +271,7 @@ namespace nil::crypto3::math { polynomial_type characteristic_seed, x_to_p; powmod_x(characteristic_seed, exponent, divisor_context, arithmetic_context); squaremod(x_to_p, characteristic_seed, divisor_context, arithmetic_context); - square_root_multiply_by_x(x_to_p, divisor_context); + multiply_by_x_mod(x_to_p, divisor_context); frobenius_product_impl(seed, characteristic_seed, field_type::arity, x_to_p, divisor_context, arithmetic_context); } else { @@ -299,7 +280,7 @@ namespace nil::crypto3::math { powmod_x(seed, exponent, divisor_context, arithmetic_context); } squaremod(x_to_q, seed, divisor_context, arithmetic_context); - square_root_multiply_by_x(x_to_q, divisor_context); + multiply_by_x_mod(x_to_q, divisor_context); return {std::move(seed), std::move(x_to_q)}; } @@ -362,7 +343,7 @@ namespace nil::crypto3::math { } auto [seed, x_to_q] = detail::square_root_x_seed(divisor_context, arithmetic_context); - detail::square_root_multiply_by_x(seed, divisor_context); + detail::multiply_by_x_mod(seed, divisor_context); const std::size_t term_count = (degree - 1) / 2; polynomial_type product, x_to_q_squared; @@ -379,7 +360,7 @@ namespace nil::crypto3::math { frobenius_product(product, product, term_count, x_to_q_squared, divisor_context, arithmetic_context); } - detail::square_root_multiply_by_x(product, divisor_context); + detail::multiply_by_x_mod(product, divisor_context); scalar_multiplication(output, product, norm_root.inversed()); return true; } @@ -447,20 +428,19 @@ namespace nil::crypto3::math { break; } candidate = generator(); - if (candidate.empty() || candidate.size() > degree || - (candidate.size() > 1 && candidate.back() == value_type::zero())) { + if (candidate.size() > degree || !detail::is_canonical_polynomial(candidate)) { throw std::invalid_argument("a trace square-root candidate must be a canonical reduced representative"); } } polynomial_type constant; squaremod(constant, sum, divisor_context, arithmetic_context); - detail::square_root_multiply_by_x(constant, divisor_context); + detail::multiply_by_x_mod(constant, divisor_context); if (constant.size() != 1 || constant[0] == value_type::zero() || !constant[0].is_square()) { throw std::logic_error("a nonzero weighted trace must yield a nonzero coefficient-field square"); } const value_type scalar = algebra::fields::sqrt_known_square(constant[0]).inversed(); - detail::square_root_multiply_by_x(sum, divisor_context); + detail::multiply_by_x_mod(sum, divisor_context); scalar_multiplication(output, sum, scalar); return true; } diff --git a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_rational_reconstruction.hpp b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_rational_reconstruction.hpp index 7f102fd1e..6c725a009 100644 --- a/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_rational_reconstruction.hpp +++ b/libs/math/include/nil/crypto3/math/polynomial/reconstruction/polynomial_rational_reconstruction.hpp @@ -85,11 +85,7 @@ namespace nil::crypto3::math { if (std::addressof(numerator) == std::addressof(denominator)) { throw std::invalid_argument("rational-reconstruction outputs must be distinct objects"); } - const auto is_canonical = [](const polynomial_type &polynomial) { - return !polynomial.empty() && - (polynomial.size() == 1 || polynomial[polynomial.size() - 1] != value_type {}); - }; - if (!is_canonical(residue) || !is_canonical(modulus)) { + if (!detail::is_canonical_polynomial(residue) || !detail::is_canonical_polynomial(modulus)) { throw std::invalid_argument("rational reconstruction requires canonical nonempty inputs"); } if (modulus.size() == 1) { 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 9ece608c2..0ffc989dc 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 @@ -55,12 +55,7 @@ namespace nil::crypto3::math { 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); + return is_canonical_polynomial(representation.p) && is_canonical_polynomial(representation.q); } } // namespace detail @@ -183,7 +178,7 @@ namespace nil::crypto3::math { 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 {})) { + if (!detail::is_canonical_polynomial(g)) { throw std::invalid_argument("polynomial X-norm recovery requires a canonical nonempty polynomial"); } if (g.size() == 1) { @@ -440,7 +435,7 @@ namespace nil::crypto3::math { using value_type = typename polynomial_type::value_type; using representation_type = polynomial_x_norm_representation; - if (h.empty() || (h.size() > 1 && h[h.size() - 1] == value_type::zero())) { + if (!detail::is_canonical_polynomial(h)) { throw std::invalid_argument("polynomial X-norm recovery requires a canonical nonempty polynomial"); } @@ -478,8 +473,10 @@ namespace nil::crypto3::math { std::vector factor_representations; bool factor_recovery_failed = false; - const auto factorization = complete_factorization( - h, arithmetic_context, coefficient_generator, [&](const polynomial_factor &factor) { + value_type leading_coefficient; + const auto control = detail::complete_factorization_impl( + leading_coefficient, h, arithmetic_context, coefficient_generator, + [&](const polynomial_factor &factor) { auto representation = detail::recover_polynomial_x_norm_factor_power( factor, arithmetic_context, coefficient_generator); if (!representation) { @@ -493,13 +490,12 @@ namespace nil::crypto3::math { if (factor_recovery_failed) { return std::nullopt; } - if (!factorization.complete || factor_representations.empty()) { + if (control != factorization_control::continue_factorization || factor_representations.empty()) { throw std::logic_error("complete factorization did not produce all nonconstant factors"); } representation_type result = detail::combine_polynomial_x_norm_representations_balanced( std::move(factor_representations), arithmetic_context); - const value_type leading_coefficient = factorization.leading_coefficient; if (leading_coefficient.is_zero()) { throw std::logic_error("complete factorization produced a zero leading coefficient"); } diff --git a/libs/math/test/complete_factorization.cpp b/libs/math/test/complete_factorization.cpp index 73901245d..eadce8d9b 100644 --- a/libs/math/test/complete_factorization.cpp +++ b/libs/math/test/complete_factorization.cpp @@ -191,7 +191,7 @@ BOOST_AUTO_TEST_CASE(irreducible_callback_runs_before_later_degree_group_work) { BOOST_REQUIRE_EQUAL(square_free.factors.size(), 1); std::vector> groups; const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( - groups, std::move(square_free.factors.front().polynomial), + &groups, std::move(square_free.factors.front().polynomial), math::detail::kaltofen_shoup_block_size(input.degree()), prefix_context, [&](const auto &group) { BOOST_CHECK(group.polynomial == linear); return math::factorization_control::stop_factorization; @@ -202,34 +202,46 @@ BOOST_AUTO_TEST_CASE(irreducible_callback_runs_before_later_degree_group_work) { // This first degree group is already irreducible: delivering its factor requires no further backend calls // or random coefficients. Returning or throwing from the callback must both avoid all later degree-group work. - for (bool throw_from_callback : {false, true}) { - std::size_t calls = 0, callback_count = 0; - polynomial_arithmetic::polynomial_context context( - counting_schoolbook_backend(calls), options); - auto generator = []() -> value_type { throw std::logic_error("unexpected random coefficient request"); }; - auto callback = [&](const math::polynomial_factor &factor) { - ++callback_count; - BOOST_CHECK(factor.polynomial == linear); - BOOST_CHECK_EQUAL(factor.multiplicity, 1); + for (bool collect : {false, true}) { + for (bool throw_from_callback : {false, true}) { + std::size_t calls = 0, callback_count = 0; + polynomial_arithmetic::polynomial_context context( + counting_schoolbook_backend(calls), options); + auto generator = []() -> value_type { throw std::logic_error("unexpected random coefficient request"); }; + auto callback = [&](const math::polynomial_factor &factor) { + ++callback_count; + BOOST_CHECK(factor.polynomial == linear); + BOOST_CHECK_EQUAL(factor.multiplicity, 1); + BOOST_CHECK_EQUAL(calls, prefix_calls); + if (throw_from_callback) + throw std::runtime_error("factor callback failure"); + return math::factorization_control::stop_factorization; + }; + auto run = [&] { + if (collect) { + const auto result = + math::complete_factorization(input, context, generator, callback); + BOOST_CHECK(!result.complete); + BOOST_REQUIRE_EQUAL(result.factors.size(), 1); + BOOST_CHECK(result.factors.front().polynomial == linear); + } else { + value_type leading_coefficient; + BOOST_CHECK(math::detail::complete_factorization_impl( + leading_coefficient, input, context, generator, callback) == + math::factorization_control::stop_factorization); + BOOST_CHECK(leading_coefficient == value_type::one()); + } + }; + if (throw_from_callback) { + BOOST_CHECK_EXCEPTION(run(), std::runtime_error, [](const auto &error) { + return std::string(error.what()) == "factor callback failure"; + }); + } else { + run(); + } + BOOST_CHECK_EQUAL(callback_count, 1); BOOST_CHECK_EQUAL(calls, prefix_calls); - if (throw_from_callback) - throw std::runtime_error("factor callback failure"); - return math::factorization_control::stop_factorization; - }; - if (throw_from_callback) { - BOOST_CHECK_EXCEPTION( - math::complete_factorization(input, context, generator, callback), - std::runtime_error, - [](const auto &error) { return std::string(error.what()) == "factor callback failure"; }); - } else { - const auto result = - math::complete_factorization(input, context, generator, callback); - BOOST_CHECK(!result.complete); - BOOST_REQUIRE_EQUAL(result.factors.size(), 1); - BOOST_CHECK(result.factors.front().polynomial == linear); } - BOOST_CHECK_EQUAL(callback_count, 1); - BOOST_CHECK_EQUAL(calls, prefix_calls); } } @@ -261,42 +273,64 @@ BOOST_AUTO_TEST_CASE(incremental_groups_preserve_factor_order_multiplicities_and polynomial_arithmetic::polynomial_context_options options; options.distinct_degree_screening_limit = screening_limit; for (std::size_t stop_after = 0; stop_after <= expected.size(); ++stop_after) { - polynomial_arithmetic::polynomial_context context(backend_type {}, options); - std::size_t next_coefficient = 0; - auto generator = [&] { return coefficients.at(next_coefficient++); }; - std::vector> reported; - const auto result = - math::complete_factorization(input, context, generator, [&](const auto &factor) { + for (bool collect : {false, true}) { + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + std::size_t next_coefficient = 0; + auto generator = [&] { return coefficients.at(next_coefficient++); }; + std::vector> reported; + auto callback = [&](const auto &factor) { reported.push_back(factor); BOOST_CHECK(generator() == value_type(100 + reported.size())); return reported.size() == stop_after ? math::factorization_control::stop_factorization : math::factorization_control::continue_factorization; - }); - const auto count = stop_after ? stop_after : expected.size(); - BOOST_REQUIRE_EQUAL(result.factors.size(), count); - BOOST_CHECK(result.factors == reported); - for (std::size_t i = 0; i < count; ++i) - BOOST_CHECK(result.factors[i] == expected[i]); - BOOST_CHECK(result.leading_coefficient == expected_result.leading_coefficient); - BOOST_CHECK_EQUAL(result.complete, stop_after == 0); - BOOST_CHECK_EQUAL(next_coefficient, consumed_at_factor[count - 1]); - if (result.complete) - BOOST_CHECK(reconstruct(backend, result) == input); + }; + value_type leading_coefficient; + bool complete; + if (collect) { + const auto result = math::complete_factorization(input, context, generator, callback); + BOOST_CHECK(result.factors == reported); + leading_coefficient = result.leading_coefficient; + complete = result.complete; + if (complete) + BOOST_CHECK(reconstruct(backend, result) == input); + } else { + complete = math::detail::complete_factorization_impl(leading_coefficient, input, + context, generator, callback) == + math::factorization_control::continue_factorization; + } + const auto count = stop_after ? stop_after : expected.size(); + BOOST_REQUIRE_EQUAL(reported.size(), count); + for (std::size_t i = 0; i < count; ++i) + BOOST_CHECK(reported[i] == expected[i]); + BOOST_CHECK(leading_coefficient == expected_result.leading_coefficient); + BOOST_CHECK_EQUAL(complete, stop_after == 0); + BOOST_CHECK_EQUAL(next_coefficient, consumed_at_factor[count - 1]); + } } } - polynomial_arithmetic::polynomial_context context; - auto throwing_generator = []() -> value_type { throw std::runtime_error("coefficient generator failure"); }; - std::size_t callback_count = 0; - BOOST_CHECK_EXCEPTION( - math::complete_factorization(input, context, throwing_generator, - [&](const auto &) { - ++callback_count; - return math::factorization_control::continue_factorization; - }), - std::runtime_error, - [](const auto &error) { return std::string(error.what()) == "coefficient generator failure"; }); - BOOST_CHECK_EQUAL(callback_count, 0); + for (bool collect : {false, true}) { + polynomial_arithmetic::polynomial_context context; + auto throwing_generator = []() -> value_type { throw std::runtime_error("coefficient generator failure"); }; + std::size_t callback_count = 0; + auto callback = [&](const auto &) { + ++callback_count; + return math::factorization_control::continue_factorization; + }; + auto run = [&] { + if (collect) { + math::complete_factorization(input, context, throwing_generator, callback); + } else { + value_type leading_coefficient; + math::detail::complete_factorization_impl(leading_coefficient, input, context, + throwing_generator, callback); + } + }; + BOOST_CHECK_EXCEPTION(run(), std::runtime_error, [](const auto &error) { + return std::string(error.what()) == "coefficient generator failure"; + }); + BOOST_CHECK_EQUAL(callback_count, 0); + } } BOOST_AUTO_TEST_CASE(screening_delivers_a_linear_factor_before_building_the_baby_table) { diff --git a/libs/math/test/distinct_degree_factorization.cpp b/libs/math/test/distinct_degree_factorization.cpp index 2c48600b1..39e702913 100644 --- a/libs/math/test/distinct_degree_factorization.cpp +++ b/libs/math/test/distinct_degree_factorization.cpp @@ -413,13 +413,13 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_fine_splitting_separates_exact_degrees) { }; polynomial_type giant_step = precomputation.baby_step(block_size); - BOOST_CHECK(math::detail::kaltofen_shoup_split_coarse_block(factors, first_block, giant_step, 1, block_size, + BOOST_CHECK(math::detail::kaltofen_shoup_split_coarse_block(&factors, first_block, giant_step, 1, block_size, precomputation, arithmetic_context, record_factor) == math::factorization_control::continue_factorization); precomputation.apply_block_frobenius(giant_step, giant_step, frobenius_context, arithmetic_context); - BOOST_CHECK(math::detail::kaltofen_shoup_split_coarse_block(factors, second_block, giant_step, 3, 1, precomputation, - arithmetic_context, record_factor) == + BOOST_CHECK(math::detail::kaltofen_shoup_split_coarse_block(&factors, second_block, giant_step, 3, 1, + precomputation, arithmetic_context, record_factor) == math::factorization_control::continue_factorization); BOOST_CHECK_EQUAL(callback_count, 3); @@ -448,7 +448,7 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_fine_splitting_honors_an_early_stop) { std::vector factors; const auto control = math::detail::kaltofen_shoup_split_coarse_block( - factors, input, giant_step, 1, block_size, precomputation, arithmetic_context, + &factors, input, giant_step, 1, block_size, precomputation, arithmetic_context, [](const factor_type &) { return math::factorization_control::stop_factorization; }); BOOST_CHECK(control == math::factorization_control::stop_factorization); @@ -477,18 +477,18 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_fine_splitting_rejects_invalid_degree_interv return math::factorization_control::continue_factorization; }; - BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(factors, input, giant_step, 0, 1, precomputation, + BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(&factors, input, giant_step, 0, 1, precomputation, arithmetic_context, continue_factorization), std::invalid_argument); - BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(factors, input, giant_step, 1, 0, precomputation, + BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(&factors, input, giant_step, 1, 0, precomputation, arithmetic_context, continue_factorization), std::invalid_argument); - BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(factors, input, giant_step, 1, 3, precomputation, + BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(&factors, input, giant_step, 1, 3, precomputation, arithmetic_context, continue_factorization), std::invalid_argument); // The stated interval contains only degree one, so the unclassified quadratic must be rejected. - BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(factors, input, giant_step, 1, 1, precomputation, + BOOST_CHECK_THROW(math::detail::kaltofen_shoup_split_coarse_block(&factors, input, giant_step, 1, 1, precomputation, arithmetic_context, continue_factorization), std::invalid_argument); } @@ -514,7 +514,7 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_internal_driver_processes_multiple_blocks) { std::size_t callback_count = 0; // A block size of two forces a full degree-1-to-2 block followed by a partial block containing degree three. const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( - factors, input, 2, arithmetic_context, [&](const factor_type &) { + &factors, input, 2, arithmetic_context, [&](const factor_type &) { ++callback_count; return math::factorization_control::continue_factorization; }); @@ -537,7 +537,7 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_internal_driver_emits_the_final_irreducible_ polynomial_arithmetic::polynomial_context arithmetic_context; std::vector factors; const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( - factors, cubic_factor, 2, arithmetic_context, + &factors, cubic_factor, 2, arithmetic_context, [](const factor_type &) { return math::factorization_control::continue_factorization; }); BOOST_CHECK(control == math::factorization_control::continue_factorization); @@ -555,7 +555,7 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_internal_driver_rejects_a_zero_block_size) { std::vector factors; BOOST_CHECK_THROW(math::detail::kaltofen_shoup_factor_monic_square_free( - factors, input, 0, arithmetic_context, + &factors, input, 0, arithmetic_context, [](const factor_type &) { return math::factorization_control::continue_factorization; }), std::invalid_argument); } @@ -576,7 +576,7 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_internal_driver_propagates_an_early_stop) { polynomial_arithmetic::polynomial_context arithmetic_context; std::vector factors; const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( - factors, input, 2, arithmetic_context, [](const factor_type &factor) { + &factors, input, 2, arithmetic_context, [](const factor_type &factor) { return factor.irreducible_factor_degree == 2 ? math::factorization_control::stop_factorization : math::factorization_control::continue_factorization; }); @@ -587,6 +587,53 @@ BOOST_AUTO_TEST_CASE(kaltofen_shoup_internal_driver_propagates_an_early_stop) { BOOST_CHECK(factors[1] == factor_type({quadratic_factor, 2})); } +BOOST_AUTO_TEST_CASE(kaltofen_shoup_streaming_matches_collected_groups_and_stop_prefixes) { + using backend_type = polynomial_arithmetic::schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + using factor_type = math::distinct_degree_factor; + + backend_type backend; + const polynomial_type linear = {fq_value_type::one(), fq_value_type::one()}; + const polynomial_type quadratic = {-fq_value_type(3), fq_value_type::zero(), fq_value_type::one()}; + const polynomial_type cubic = {-fq_value_type(3), fq_value_type::zero(), fq_value_type::zero(), + fq_value_type::one()}; + const polynomial_type second_cubic = {-fq_value_type(9), fq_value_type::zero(), fq_value_type::zero(), + fq_value_type::one()}; + const auto mixed = multiply(backend, multiply(backend, linear, quadratic), multiply(backend, cubic, second_cubic)); + for (const auto &input : {linear, cubic, mixed}) { + for (std::size_t screening_limit : {0, 1, 100}) { + for (std::size_t block_size : {1, 2}) { + polynomial_arithmetic::polynomial_context_options options; + options.distinct_degree_screening_limit = screening_limit; + polynomial_arithmetic::polynomial_context context(backend_type {}, options); + std::vector expected; + BOOST_REQUIRE(math::detail::kaltofen_shoup_factor_monic_square_free( + &expected, input, block_size, context, [&](const factor_type &factor) { + // The public collecting path must append before reporting the group. + BOOST_REQUIRE(!expected.empty()); + BOOST_CHECK(&factor == &expected.back()); + return math::factorization_control::continue_factorization; + }) == math::factorization_control::continue_factorization); + + for (std::size_t stop_after = 0; stop_after <= expected.size(); ++stop_after) { + polynomial_arithmetic::polynomial_context streaming_context(backend_type {}, options); + std::size_t reported = 0; + const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( + nullptr, input, block_size, streaming_context, [&](const factor_type &factor) { + BOOST_REQUIRE_LT(reported, expected.size()); + BOOST_CHECK(factor == expected[reported++]); + return reported == stop_after ? math::factorization_control::stop_factorization : + math::factorization_control::continue_factorization; + }); + BOOST_CHECK_EQUAL(reported, stop_after ? stop_after : expected.size()); + BOOST_CHECK(control == (stop_after ? math::factorization_control::stop_factorization : + math::factorization_control::continue_factorization)); + } + } + } + } +} + BOOST_AUTO_TEST_CASE(screening_continuation_matches_reference_and_unscreened_blocks) { using backend_type = polynomial_arithmetic::schoolbook_backend; using mixed_backend = polynomial_arithmetic::mixed_radix_backend; @@ -634,7 +681,7 @@ BOOST_AUTO_TEST_CASE(screening_continuation_matches_reference_and_unscreened_blo const auto reference = math::distinct_degree_factorization_reference(mixed_degrees, context); std::vector> groups; const auto control = math::detail::kaltofen_shoup_factor_monic_square_free( - groups, mixed_degrees, block_size, context, + &groups, mixed_degrees, block_size, context, [](const auto &) { return math::factorization_control::continue_factorization; }); BOOST_CHECK(control == math::factorization_control::continue_factorization); BOOST_CHECK(groups == reference.factors); diff --git a/libs/math/test/polynomial_algorithm_constraints.cpp b/libs/math/test/polynomial_algorithm_constraints.cpp index 9bb0e09ef..2a9d5d70d 100644 --- a/libs/math/test/polynomial_algorithm_constraints.cpp +++ b/libs/math/test/polynomial_algorithm_constraints.cpp @@ -44,6 +44,10 @@ using coefficient_polynomial = math::polynomial; using evaluation_polynomial = math::polynomial_dfs; using interpolation_points = std::vector>; +template +concept SupportsCanonicalPolynomialCheck = + requires(const Polynomial &input) { math::detail::is_canonical_polynomial(input); }; + template concept SupportsPolynomialMultiplication = requires(Range &output, const Range &input) { math::multiplication(output, input, input); }; @@ -103,6 +107,10 @@ static_assert(math::detail::PolynomialCoefficientRange static_assert(!math::detail::MutablePolynomialCoefficientRange>); static_assert(!math::detail::PolynomialCoefficientRange>); +static_assert(SupportsCanonicalPolynomialCheck); +static_assert(!SupportsCanonicalPolynomialCheck); +static_assert(!SupportsCanonicalPolynomialCheck); + static_assert(SupportsPolynomialMultiplication); static_assert(!SupportsPolynomialMultiplication>); diff --git a/libs/math/test/polynomial_exponentiation.cpp b/libs/math/test/polynomial_exponentiation.cpp index dd8642fce..0a6528516 100644 --- a/libs/math/test/polynomial_exponentiation.cpp +++ b/libs/math/test/polynomial_exponentiation.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -245,6 +246,51 @@ BOOST_AUTO_TEST_CASE(powmod_supports_extension_coefficients_and_base_field_roots BOOST_CHECK(result == expected); } +BOOST_AUTO_TEST_CASE(multiply_by_x_mod_matches_general_products_for_reduced_inputs) { + using backend_type = counting_schoolbook_backend; + using polynomial_type = typename backend_type::polynomial_type; + const auto zero = fq_value_type::zero(); + const auto one = fq_value_type::one(); + const polynomial_type x = {zero, one}; + const polynomial_type divisors[] = {{zero, one}, + {fq_value_type(3), fq_value_type(2)}, + {one, zero, one}, + {fq_value_type(3), fq_value_type(2), fq_value_type(4)}, + {fq_value_type(3), fq_value_type(2), zero, fq_value_type(4)}, + {zero, zero, zero, one}}; + std::size_t squares = 0, products = 0; + polynomial_arithmetic::polynomial_context context(backend_type {squares, products}); + for (const auto &modulus : divisors) { + const std::size_t degree = modulus.size() - 1; + const math::polynomial_divisor_context divisor(modulus, degree, context); + std::vector inputs = {{zero}}; + for (std::size_t width = 1; width <= degree; ++width) { + polynomial_type input(width); + for (std::size_t i = 0; i < width; ++i) { + input[i] = fq_value_type(3 * i + 2); + } + inputs.push_back(std::move(input)); + } + // X * ((B-B(0))/X) reduces to a constant, or zero when B(0)=0. + polynomial_type collapsing(degree); + for (std::size_t i = 0; i < degree; ++i) { + collapsing[i] = modulus[i + 1]; + } + inputs.push_back(std::move(collapsing)); + for (const auto &input : inputs) { + polynomial_type expected; + math::mulmod(expected, input, x, divisor, context); + const auto saved_squares = squares, saved_products = products; + auto actual = input; + math::detail::multiply_by_x_mod(actual, divisor); + BOOST_CHECK(actual == expected); + BOOST_CHECK_LE(actual.size(), degree); + BOOST_CHECK_EQUAL(squares, saved_squares); + BOOST_CHECK_EQUAL(products, saved_products); + } + } +} + BOOST_AUTO_TEST_CASE(powmod_x_matches_generic_exponentiation_across_moduli_and_reduction_paths) { using backend_type = polynomial_arithmetic::schoolbook_backend; using polynomial_type = typename backend_type::polynomial_type; diff --git a/libs/math/test/small_prime_fft_arithmetic.cpp b/libs/math/test/small_prime_fft_arithmetic.cpp index 0b32f3864..7f6ab99e5 100644 --- a/libs/math/test/small_prime_fft_arithmetic.cpp +++ b/libs/math/test/small_prime_fft_arithmetic.cpp @@ -24,8 +24,8 @@ #define BOOST_TEST_MODULE small_prime_fft_arithmetic_test -#include -#include +#include +#include #include #include diff --git a/libs/math/test/small_prime_fft_conversion.cpp b/libs/math/test/small_prime_fft_conversion.cpp index 91a81c82a..22011ee6f 100644 --- a/libs/math/test/small_prime_fft_conversion.cpp +++ b/libs/math/test/small_prime_fft_conversion.cpp @@ -24,7 +24,7 @@ #define BOOST_TEST_MODULE small_prime_fft_conversion_test -#include +#include #include #include diff --git a/libs/math/test/small_prime_fft_conversion_32bit.cpp b/libs/math/test/small_prime_fft_conversion_32bit.cpp index d04383eb1..386f95d84 100644 --- a/libs/math/test/small_prime_fft_conversion_32bit.cpp +++ b/libs/math/test/small_prime_fft_conversion_32bit.cpp @@ -27,7 +27,7 @@ // Boost first so subsequent includes cannot re-enable its 128-bit integer support. #undef BOOST_HAS_INT128 -#include +#include #include diff --git a/libs/math/test/small_prime_fft_crt.cpp b/libs/math/test/small_prime_fft_crt.cpp index 25d46a055..30b03c775 100644 --- a/libs/math/test/small_prime_fft_crt.cpp +++ b/libs/math/test/small_prime_fft_crt.cpp @@ -24,7 +24,8 @@ #define BOOST_TEST_MODULE small_prime_fft_crt_test -#include +#include +#include #include #include diff --git a/libs/math/test/small_prime_fft_params.cpp b/libs/math/test/small_prime_fft_params.cpp index 9c6901341..8ffc28565 100644 --- a/libs/math/test/small_prime_fft_params.cpp +++ b/libs/math/test/small_prime_fft_params.cpp @@ -35,7 +35,7 @@ #include #include -#include +#include namespace { namespace detail = nil::crypto3::math::detail; diff --git a/libs/math/test/small_prime_fft_transform.cpp b/libs/math/test/small_prime_fft_transform.cpp index 39d60aa1d..30f8c4527 100644 --- a/libs/math/test/small_prime_fft_transform.cpp +++ b/libs/math/test/small_prime_fft_transform.cpp @@ -24,7 +24,7 @@ #define BOOST_TEST_MODULE small_prime_fft_transform_test -#include +#include #include #include