diff --git a/libs/math/README.md b/libs/math/README.md index d50ac677b..2dd6a2e1a 100644 --- a/libs/math/README.md +++ b/libs/math/README.md @@ -1,7 +1,8 @@ # Crypto3.Math -Matrix and polynomial arithmetic, evaluation domains, and Fast Fourier -Transform algorithms for the Crypto3 suite. +Matrix and polynomial arithmetic, evaluation domains, exact geometric-domain +Lagrange weights, Fast Fourier Transform algorithms, finite-field polynomial +factorization, and polynomial recovery primitives for the Crypto3 suite. This header-only component is maintained as part of the Crypto3 monorepo. See the [root build instructions](../../README.md#clone-and-build) to configure the diff --git a/libs/math/docs/concepts.md b/libs/math/docs/concepts.md index 313bb9acc..900c39b55 100644 --- a/libs/math/docs/concepts.md +++ b/libs/math/docs/concepts.md @@ -1,3 +1,117 @@ -# Concepts # {#fft_concepts} +# C++ concepts {#fft_concepts} -@tableofcontents \ No newline at end of file +@tableofcontents + +This page lists the C++20 concepts exposed by Crypto3.Math. A concept is a compile-time predicate used to constrain a +template argument. It describes the interface an argument must provide; it does not construct an object or generate +runtime code. + +Crypto3.Math currently exposes concepts for matrix and vector expressions, polynomial representations, and polynomial +multiplication backends. Concepts declared in a `detail` namespace are implementation constraints and are summarized +separately at the end of this page. + +## Matrix and vector concepts + +The matrix concepts are declared in: + + + +### Readable expressions + +`VectorExpression` describes a readable vector-like expression. The type must provide `value_type`, `size_type`, +`size()`, and indexed access through `value(i)`. + +`MatrixExpression` describes a readable matrix-like expression. The type must provide `value_type`, `size_type`, +indexed access through `value(row, column)`, and dimensions through either: + +* `rows()` and `columns()`; or +* `size1()` and `size2()`. + +The free `rows(expression)` and `columns(expression)` functions provide one interface for both dimension naming +conventions. + +Expression concepts are deliberately weak. In particular, an expression-template result may be readable without +owning storage or permitting mutation. + +### Writable backends + +`VectorBackend` extends `VectorExpression`. It additionally requires a semiregular type and writable indexed +access through `value(i)`. + +`MatrixBackend` similarly extends `MatrixExpression` with semiregular value semantics and writable access +through `value(row, column)`. + +The resizable variants add their corresponding resize operation: + +| Concept | Additional operation | +|---|---| +| `ResizableVectorBackend` | `value.resize(size)` | +| `ResizableMatrixBackend` | `value.resize(rows, columns)` | + +Algorithms that only read their operands should accept expression types. Algorithms that create or overwrite results +use backend or resizable-backend constraints as appropriate. + +## Polynomial representation concepts + +The polynomial representation concepts are declared in: + + + +Both concepts require `value_type`, `size_type`, `size()`, `degree()`, indexed read access, and a representation tag. +They accept const and reference-qualified types. + +### CoefficientPolynomial + +`CoefficientPolynomial` identifies a polynomial stored as coefficients in ascending degree order: + + [a0, a1, ..., an] represents a0 + a1 X + ... + an X^n. + +Its `representation_type` must be `coefficient_representation`. Crypto3's `polynomial`, `polynomial_view`, and +`polymorphic_polynomial` satisfy this concept. + +```cpp +template +void consume_coefficients(const Polynomial &polynomial); +``` + +### EvaluationPolynomial + +`EvaluationPolynomial` identifies a polynomial stored as evaluations over a domain. Indexed entries are samples, +not coefficients, and the logical degree may differ from `size() - 1`. + +Its `representation_type` must be `evaluation_representation`. Crypto3's `polynomial_dfs`, `polynomial_dfs_view`, and +`polymorphic_polynomial_dfs` satisfy this concept. + +```cpp +template +void consume_evaluations(const Polynomial &polynomial); +``` + +The representation tag prevents a coefficient algorithm from accidentally accepting an evaluation polynomial merely +because both types provide similar container operations. A bare `std::vector` satisfies neither concept because it +does not identify what its entries represent. + +These concepts describe readable representation only. They do not require field-valued entries, mutability, +canonical storage, arithmetic operators, or alias-safe operations. + +## PolynomialBackend + +The polynomial multiplication backend concept is declared in: + + + +`polynomial_arithmetic::PolynomialBackend` requires an associated `polynomial_type` that satisfies +`CoefficientPolynomial` and three operations: + +```cpp +backend.multiply(output, left, right); +backend.square(output, input); +backend.multiply_low(output, left, right, coefficient_count); +``` + +The operations produce canonical coefficient polynomials and permit output to alias an input. `multiply_low` computes +the product modulo `X^coefficient_count`. + +Backend operations receive a mutable backend object so implementations can reuse plans and scratch storage. The +[polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) documents the supplied schoolbook and +mixed-radix backends and the `polynomial_context` that owns them. diff --git a/libs/math/docs/geometric_lagrange.md b/libs/math/docs/geometric_lagrange.md new file mode 100644 index 000000000..338dc9212 --- /dev/null +++ b/libs/math/docs/geometric_lagrange.md @@ -0,0 +1,97 @@ +# Exact geometric-domain Lagrange weights {#math_geometric_lagrange} + +@tableofcontents + +`geometric_sequence_domain` represents the exact-size domain + + x_i = r^i, 0 <= i < m, + +where `r` is the field's configured geometric generator. The domain provides a linear-time field-element API for +evaluating every Lagrange basis polynomial at an arbitrary point. + +The relevant headers are: + +| Facility | Header | +|---|---| +| Montgomery batch inversion | `` | +| Evaluation-domain interface | `` | +| Exact geometric domain | `` | + +## Domain precomputation + +Construction requires `m > 1`, a nonzero geometric generator, and distinct points `1, r, ..., r^(m-1)`. It rejects a +size that would repeat a point. + +For the Lagrange path, the constructor precomputes: + +* every domain point `x_i`; +* the barycentric weights `w_i = 1 / Z'(x_i)`; and +* the coefficients of the vanishing polynomial + + Z(X) = product(X - x_i, i = 0 .. m - 1). + +For general points, constructing every `w_i` from pairwise differences would take quadratic work. Geometric points +instead satisfy the recurrence + + w_0 = product((1 - x_j)^-1, j = 1 .. m - 1), + w_i = -w_(i-1) * x_(m-1-i)^-1 * (1 - x_(m-i)) / (1 - x_i). + +All required inverses of `1 - r^i`, together with `r^-1`, are obtained through one Montgomery batch inversion. For +`m` nonzero inputs, `batch_inverse_nonzero` performs one field inversion and exactly `3 * (m - 1)` field +multiplications. + +Writing + + Z(X) = sum(c_j * X^(m-j), j = 0 .. m), c_0 = 1, + +the vanishing-polynomial coefficients use the recurrence + + c_j = -c_(j-1) * r^(j-1) * (1 - r^(m-j+1)) / (1 - r^j). + +If `r` has exact order `m`, the final denominator vanishes and the implementation uses the resulting special case +`Z(X) = X^m - 1`. + +Construction therefore takes `O(m)` field operations, one field inversion, and `O(m)` stored field elements. + +## Evaluating all weights + +The combined overload returns both the weights and `Z(t)` without repeating the product: + +```cpp +#include + +namespace math = nil::crypto3::math; + +math::geometric_sequence_domain domain(domain_size); + +value_type vanishing_at_t; +const std::vector weights = + domain.evaluate_all_lagrange_polynomials(t, vanishing_at_t); +``` + +Away from the domain, it uses the barycentric formula + + L_i(t) = Z(t) * w_i / (t - x_i). + +The implementation forms all denominators `t - x_i`, computes `Z(t)` as their product, and batch-inverts the +denominators with one field inversion. If `t = x_i`, it instead returns the corresponding unit vector and sets +`Z(t) = 0` without attempting an inversion. + +For evaluations `f(x_i)`, the returned weights satisfy + + f(t) = sum(f(x_i) * L_i(t), i = 0 .. m - 1) + +for every polynomial of degree below `m`. The weights belong to `FieldType::value_type`; they may also scale values in +a compatible extension field. + +## Rough current costs + +| Operation | Rough current cost | +|---|---| +| Construct a size-`m` domain | `O(m)` field operations and one inversion | +| Evaluate all weights at one off-domain point | `O(m)` field operations and one inversion | +| Evaluate at a domain point | `O(m)` work and no inversion | +| Evaluate all weights at `k` independent points | `O(k * m)` field operations and at most `k` inversions | + +Returning `m` weights already requires linear output work. All scratch storage used by an evaluation is local, while +the domain precomputation is immutable, so concurrent weight evaluations on one domain do not share mutable scratch. diff --git a/libs/math/docs/index.md b/libs/math/docs/index.md index dd81bf2c6..a014df3ee 100644 --- a/libs/math/docs/index.md +++ b/libs/math/docs/index.md @@ -1,3 +1,15 @@ # Crypto3.Math {#fft_index} -@subpage fft_introduction @subpage fft_manual @subpage fft_concepts +@subpage fft_introduction + +@subpage fft_manual + +@subpage fft_concepts + +@subpage math_geometric_lagrange + +@subpage math_polynomial_arithmetic + +@subpage math_polynomial_factorization + +@subpage math_polynomial_recovery diff --git a/libs/math/docs/introduction.md b/libs/math/docs/introduction.md index d84e52150..a8f54d0e5 100644 --- a/libs/math/docs/introduction.md +++ b/libs/math/docs/introduction.md @@ -2,16 +2,20 @@ @tableofcontents -Crypto3.Math extends the Crypto3 algebra library and provides Fast Fourier Transform -evaluation algorithms implemented in way C++ standard library implies: concepts, algorithms, predictable behavior, -latest standard features support and clean architecture without compromising security and performance. +Crypto3.Math extends the Crypto3 algebra library with polynomial and matrix arithmetic, evaluation domains, FFTs over +finite fields, polynomial factorization, and quotient-field recovery algorithms. Its interfaces use modern C++ +concepts, generic algorithms, explicit arithmetic contexts, and interchangeable multiplication backends. Crypto3.Math consists of several parts to review: * [Manual](@ref fft_manual). * [Concepts](@ref fft_concepts). +* [Exact geometric-domain Lagrange weights](@ref math_geometric_lagrange). +* [Polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic). +* [Polynomial factorization](@ref math_polynomial_factorization). +* [Polynomial recovery](@ref math_polynomial_recovery). -## Background +## FFT background There is currently a variety of algorithms for computing the Fast Fourier Transform (FFT) over the field of complex numbers. For this situation, there exists many libraries, such as [FFTW](http://www.fftw.org/), that have been @@ -84,6 +88,9 @@ conversion algorithm between the monomial and the Newton bases. The domain takes Newton evaluation and interpolation by choosing sample points that form a geometric progression, _a\_n = r^(n-1)_, see \[BS05\]. +The field-element path also provides [exact geometric-domain Lagrange weights](@ref math_geometric_lagrange) with +linear domain precomputation and linear work per evaluation point. + ## Dependencies ## {#fft_dependencies} Internal dependencies: diff --git a/libs/math/docs/manual.md b/libs/math/docs/manual.md index 9e63803e2..a4e025101 100644 --- a/libs/math/docs/manual.md +++ b/libs/math/docs/manual.md @@ -1,3 +1,20 @@ # Manual # {#fft_manual} -@tableofcontents \ No newline at end of file +@tableofcontents + +Crypto3.Math provides coefficient-form and evaluation-form polynomials, evaluation domains, FFT algorithms, and +backend-aware polynomial arithmetic. The higher-level polynomial algorithms share a +`polynomial_arithmetic::polynomial_context`. Reusing a context lets its multiplication backend retain plans, +configuration, and scratch storage across an operation. + +The following pages describe the optimized domain and polynomial facilities: + +* [Exact geometric-domain Lagrange weights](@ref math_geometric_lagrange) covers exact-size geometric points, + barycentric precomputation, batch inversion, vanishing polynomials, and linear-time weight evaluation. +* [Polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) describes multiplication backends, + reusable contexts, division, GCD, quotient-ring arithmetic, modular composition, and Frobenius maps shared by the + higher-level algorithms. +* [Polynomial factorization](@ref math_polynomial_factorization) covers square-free, distinct-degree, equal-degree, + complete, and staged factorization. +* [Polynomial recovery](@ref math_polynomial_recovery) covers square testing and square roots in polynomial quotient + fields, bounded rational reconstruction, and the relation between these operations and polynomial norms. diff --git a/libs/math/docs/polynomial_arithmetic.md b/libs/math/docs/polynomial_arithmetic.md new file mode 100644 index 000000000..044ecec91 --- /dev/null +++ b/libs/math/docs/polynomial_arithmetic.md @@ -0,0 +1,281 @@ +# Polynomial arithmetic infrastructure {#math_polynomial_arithmetic} + +@tableofcontents + +The higher-level polynomial algorithms in Crypto3.Math share a backend-aware arithmetic layer. It separates the +mathematical algorithms from the multiplication implementation and gives callers explicit ownership of expensive +precomputation. + +The main dependency chain is: + +| Layer | Main types and operations | +|---|---| +| Multiplication | `PolynomialBackend`, `schoolbook_backend`, `mixed_radix_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` | +| Modular composition | `compose_mod_reference`, `polynomial_composition_precomputation`, `compose_mod` | +| Frobenius | `polynomial_frobenius_context`, `frobenius_map` | + +### Header map + +| Facility | Header | +|---|---| +| Representation concepts | `` | +| Backend concept and arithmetic context | `` | +| Schoolbook and mixed-radix backends | ``, `` | +| Power-series inversion | `` | +| Divisor context and division | `` | +| GCD | `` | +| Quotient-ring multiplication | `` | +| Quotient-ring exponentiation | `` | +| Modular composition | `` | +| Frobenius maps | `` | + +## Multiplication backends + +`PolynomialBackend` is the compile-time interface used by the arithmetic layer. A backend supplies an associated +coefficient polynomial type and three alias-safe operations: + +```cpp +backend.multiply(output, left, right); +backend.square(output, input); +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: + +* `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. + +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. + +## 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: + +```cpp +namespace pa = nil::crypto3::math::polynomial_arithmetic; + +pa::polynomial_context_options options; +options.gcd_half_gcd_cutoff = 0; + +pa::polynomial_context arithmetic_context( + backend_type{}, options); +``` + +The context forwards `multiply`, `square`, and `multiply_low` to its backend. It is intended for sequential reuse. +Separate contexts are required for concurrent calls when the backend owns mutable scratch storage. + +### Algorithm-selection options + +| Option | Default | Meaning | +|---|---:|---| +| `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. | +| `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. | + +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 +coefficient type and multiplication backend; callers should choose a nonzero cutoff from representative benchmarks. + +## Newton power-series inversion + +`inverse_series(output, input, coefficient_count, arithmetic_context)` computes + + output * input = 1 mod X^coefficient_count. + +The constant coefficient of `input` must be nonzero. Starting from its scalar inverse, Newton iteration doubles the +known precision on each step. The implementation preserves the already correct low block and computes only the new +high coefficients using `multiply_low`. Output may alias input, and a zero requested precision returns the canonical +zero polynomial. + +Power-series inversion is the precomputation behind fast polynomial division. + +## Polynomial divisor contexts + +`polynomial_divisor_context` stores a canonical nonzero divisor `B` and a truncated inverse of its reversed +polynomial. If `d = degree(B)`, define + + reverse(B) = X^d * B(X^-1). + +Construction precomputes + + reverse(B)^-1 mod X^inverse_precision. + +The chosen precision is the maximum number of quotient coefficients that later fast divisions may require. Reusing the +context avoids repeating this Newton inversion for every division or modular reduction by the same `B`. + +```cpp +const std::size_t quotient_coefficient_bound = + dividend.size() >= divisor.size() + ? dividend.size() - divisor.size() + 1 + : 1; +math::polynomial_divisor_context divisor_context( + divisor, quotient_coefficient_bound, arithmetic_context); +``` + +The context is immutable after construction. Operations receiving it assume it represents the intended divisor; when +another precomputation also depends on `B`, callers must keep those contexts paired with the same polynomial. + +## Polynomial division + +`divrem` computes canonical `quotient` and `remainder` satisfying + + dividend = quotient * B + remainder, + degree(remainder) < degree(B). + +The quotient and remainder must be distinct output objects, but either may alias the dividend. For small divisors or +quotients, the arithmetic options select quadratic long division. Otherwise, for + + k = degree(dividend) - degree(B) + 1, + +reversal turns division into the truncated product + + reverse(quotient) = reverse(dividend) * reverse(B)^-1 mod X^k. + +The divisor context must have inverse precision of at least `k` on this fast path. Long division does not use the +precomputed inverse and therefore does not require that precision. + +Two convenience operations share the same dispatch: + +* `remainder` returns only `dividend mod B`. +* `exact_division` returns the quotient and throws `std::invalid_argument` if the remainder is nonzero. + +## GCD and half-GCD + +`gcd(output, left, right, arithmetic_context)` returns the canonical monic greatest common divisor. Output may alias an +input, `gcd(0, B)` is the monic form of `B`, and `gcd(0, 0)` is zero. + +The Euclidean algorithm repeatedly replaces `(A, B)` by `(B, A mod B)`. Each step preserves the common divisors and +strictly lowers the second degree until the final nonzero remainder is the GCD up to a scalar. + +For sufficiently large operands, half-GCD batches several Euclidean steps. One quotient step acts on the polynomial +pair through + + [ 0 1 ] [ A ] [ B ] + [ 1 -q ] [ B ] = [ A - q * B ]. + +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. + +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 +field. + +## Quotient-ring arithmetic + +The modular operations work in `K[X]/(B)` for any nonzero `B`. Irreducibility is not required unless a caller needs the +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. + +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. + +## Brent-Kung modular composition + +Modular composition computes + + outer(inner(X)) mod B. + +`compose_mod_reference` uses Horner's rule and one modular multiplication per nonleading coefficient of `outer`. It is +the simple correctness reference. + +The faster `compose_mod` overload uses blocked Brent-Kung composition. For block size `k`, write + + outer(Y) = F0(Y) + F1(Y) * Y^k + F2(Y) * Y^(2k) + ... . + +`polynomial_composition_precomputation` caches the baby powers + + 1, inner, ..., inner^(k - 1) mod B + +and the giant step `inner^k mod B`. Each `Fi(inner)` is a linear combination of the baby powers, and the block values +are combined by Horner's rule in the giant step. The default `k = ceil(sqrt(L))`, where `L` is the maximum outer +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. + +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. + +## Iterated Frobenius maps + +Let the finite coefficient field `K` contain `Q` elements. Every coefficient satisfies `a^Q = a`, so for a polynomial +`A`: + + 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. + +```cpp +math::polynomial_frobenius_context frobenius_context( + divisor, arithmetic_context); + +math::frobenius_map(output, input, frobenius_context, arithmetic_context); +math::frobenius_map(output, input, iteration_count, + frobenius_context, arithmetic_context); +``` + +One map raises a quotient-ring element to the `Q`-th power. The iterated overload applies that map +`iteration_count` times, producing the `Q^iteration_count` power; zero iterations only reduce the input. The divisor +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. + +## Context ownership and reuse + +The arithmetic context normally has the longest lifetime. Divisor, composition, Frobenius, and square-root contexts +are immutable precomputations tied to particular polynomials and may reference or logically depend on one another. +Construct them once per fixed input and reuse them sequentially. Do not combine precomputations built for different +divisors merely because their degrees match. + +## Cost overview + +Let `n` be the polynomial size, `M(n)` the active backend's multiplication cost, `e` an exponent, and `Q` the +coefficient-field order. These are rough bounds for the current implementation: + +| Operation | Rough time bound | +|---|---| +| Schoolbook multiplication or squaring | `O(n^2)` | +| Mixed-radix multiplication | `O(N * sum(radices))`, where `N` is the selected transform size | +| Newton series inversion | `O(M(n))` | +| Divisor-context construction | `O(M(n))` | +| Long division | `O(n^2)` | +| Newton division with a reused divisor context | `O(M(n))` | +| Euclidean GCD | `O(n^2)` with the default small-quotient path | +| 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 | +| 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)` | + +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. diff --git a/libs/math/docs/polynomial_factorization.md b/libs/math/docs/polynomial_factorization.md new file mode 100644 index 000000000..ce701113e --- /dev/null +++ b/libs/math/docs/polynomial_factorization.md @@ -0,0 +1,159 @@ +# Polynomial factorization {#math_polynomial_factorization} + +@tableofcontents + +Crypto3.Math factors univariate coefficient polynomials over finite fields. The algorithms operate on the polynomial +type supplied by a multiplication backend, so the same factorization code can use the reference schoolbook backend or +a faster backend suitable for the coefficient field and operand sizes. + +## Result representation + +A complete factorization has the form + + input = leading_coefficient * product(factor.polynomial ^ factor.multiplicity). + +`polynomial_factorization_result` stores the original leading coefficient, the monic irreducible factors +and their positive multiplicities, and a `complete` flag. Nonconstant factors are canonical and monic. The canonical +zero polynomial and nonzero constants produce no polynomial factors; their scalar value is returned as +`leading_coefficient`. + +Distinct-degree factorization uses a related result type. Each `distinct_degree_factor` contains the product of all +irreducible factors having one degree. Such a group is generally not itself irreducible. For example, a group reported +with `irreducible_factor_degree == 2` may be the product of several distinct irreducible quadratics. + +## Factorization pipeline + +Complete factorization composes three stages: + +| Stage | Public API | Purpose | +|---|---|---| +| Square-free factorization | `square_free_factorization` | Separate factors by their multiplicity using Yun's algorithm. | +| Distinct-degree factorization | `distinct_degree_factorization_kaltofen_shoup` | Group the square-free factors by irreducible degree using blocked Frobenius steps and GCDs. | +| Equal-degree factorization | `equal_degree_factorization` | Split one degree group into individual irreducible factors using Cantor-Zassenhaus. | + +`distinct_degree_factorization_reference` provides the unblocked distinct-degree algorithm as a correctness reference. +Production callers normally use the Kaltofen-Shoup implementation selected by `complete_factorization`. + +The [polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) provides the reusable divisor contexts, +Newton division, GCD and half-GCD, quotient-ring operations, Brent-Kung modular composition, and iterated Frobenius +maps used by these stages. + +### Header map + +| Facility | Header | +|---|---| +| Result and callback types | `` | +| Square-free factorization | `` | +| Reference distinct-degree factorization | `` | +| Kaltofen-Shoup distinct-degree factorization | `` | +| Equal-degree factorization | `` | +| Complete factorization | `` | + +## Complete factorization + +The simplest entry point takes a polynomial arithmetic context and a caller-owned random generator: + +```cpp +#include +#include +#include +#include + +namespace math = nil::crypto3::math; +namespace pa = math::polynomial_arithmetic; +namespace fields = nil::crypto3::algebra::fields; + +using field_type = fields::babybear; +using value_type = field_type::value_type; +using backend_type = pa::schoolbook_backend; +using polynomial_type = backend_type::polynomial_type; + +// (X + 1) * (X + 2)^2 = X^3 + 5 X^2 + 8 X + 4. +polynomial_type input = { + value_type(4), value_type(8), value_type(5), value_type::one() +}; + +pa::polynomial_context arithmetic_context; +nil::crypto3::random::algebraic_engine generator(42); +auto result = math::complete_factorization( + input, arithmetic_context, generator); + +for (const auto &factor : result.factors) { + // factor.polynomial is monic and irreducible. + // factor.multiplicity is its multiplicity in input. +} +``` + +The generator is injected rather than created internally. It must return independent, uniformly distributed elements +of the coefficient field. Supplying a seeded generator makes a run reproducible. + +Cantor-Zassenhaus repeats randomized split trials until a nontrivial factor is found; the implementation does not set +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. + +## Staged factorization + +The factorization APIs have callback overloads for callers that can decide after each emitted factor whether more work +is necessary: + +```cpp +auto result = math::complete_factorization( + input, arithmetic_context, generator, + [](const math::polynomial_factor &factor) { + return should_stop(factor) + ? math::factorization_control::stop_factorization + : math::factorization_control::continue_factorization; + }); +``` + +The factor that triggers the stop is included in `result.factors`, and `result.complete` is `false`. A complete run has +`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. + +## Preconditions and current restrictions + +* Inputs are nonempty coefficient polynomials. Algorithms canonicalize trailing zero coefficients and return monic + nonconstant factors. +* Square-free factorization currently requires the coefficient-field characteristic to exceed the input degree. This + avoids the separate polynomial p-th-root path required when differentiation erases p-th powers. +* Distinct-degree and equal-degree entry points require square-free input. Complete factorization establishes this + precondition through its first stage. +* Equal-degree and complete factorization currently implement the odd-characteristic Cantor-Zassenhaus algorithm and + reject characteristic two. +* An equal-degree input must be a group whose irreducible factors all have the supplied degree. The API checks whether + the factor degree divides the total degree but does not repeat distinct-degree factorization to verify the full + precondition. + +## Verifying a result + +A complete result can be checked by starting with the constant polynomial containing `leading_coefficient`, multiplying +each monic factor into it `multiplicity` times through the same backend, and comparing the canonical result with the +input. This reconstruction identity is the primary contract of `polynomial_factorization_result`. + +## Cost overview + +Let `n` be the input degree, `d` the common irreducible-factor degree in an equal-degree group, and `Q` the +coefficient-field order. The current implementations perform the rough work shown below. + +The table counts expensive polynomial operations. Their individual bounds—particularly GCD, exact division, modular +multiplication, and cached Frobenius maps—are given in the +[polynomial arithmetic cost overview](@ref math_polynomial_arithmetic) and depend on the selected multiplication +backend. + +| Stage | Rough current work | +|---|---| +| Yun square-free factorization | Up to `n` iterations, each with one GCD and two exact divisions; divisor inverses are rebuilt | +| Reference distinct-degree factorization | Up to `n / 2` cached Frobenius maps and `n / 2` GCDs | +| Kaltofen-Shoup coarse phase | `O(sqrt(n))` Frobenius maps, `O(n)` modular multiplications, and `O(sqrt(n))` coarse GCDs | +| Kaltofen-Shoup fine splitting | Up to `O(n)` additional GCDs and exact divisions | +| One Cantor-Zassenhaus trial | `O(d * log(Q))` modular multiplications plus at most two polynomial GCDs | +| Equal-degree factorization | An expected linear number of Cantor-Zassenhaus trials in the number of output factors | +| 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. diff --git a/libs/math/docs/polynomial_recovery.md b/libs/math/docs/polynomial_recovery.md new file mode 100644 index 000000000..13ef8f0cc --- /dev/null +++ b/libs/math/docs/polynomial_recovery.md @@ -0,0 +1,241 @@ +# Polynomial recovery {#math_polynomial_recovery} + +@tableofcontents + +Crypto3.Math provides two complementary recovery primitives: + +* square testing and square roots in a finite polynomial quotient field; and +* bounded rational reconstruction from a residue modulo a polynomial. + +Together they can recover representations of irreducible factors by the polynomial norm form + + P(X)^2 - X * Q(X)^2. + +The library exposes the generic arithmetic operations. Applications remain responsible for selecting factors, degree +bounds, random sources, and any policy for combining or rejecting recovered representations. + +### Header map + +| Facility | Header | +|---|---| +| Field orders and multiplicative-group decomposition | `` | +| Coefficient-field square-root helpers | `` | +| Quotient-field square testing and roots | `` | +| Bounded rational reconstruction | `` | + +## Field orders and coefficient square roots + +The quotient-field algorithms derive their exponents from the public field-order utilities: + +* `field_characteristic()` returns the prime characteristic; +* `field_order()` returns the number of elements, including the extension degree described by `Field`; +* `extension_field_order(d)` returns `field_order()^d`; and +* the corresponding multiplicative-group decomposition functions write an order minus one as + `odd_order * 2^two_adicity`. + +The templates accept either a Crypto3 field type or its value type and return `boost::multiprecision::cpp_int`, since +an extension-field order can exceed the fixed-width field representation. These integers are used as public loop +bounds and exponents; polynomial coefficients remain native field values. + +BN254 Fq12 values provide `is_square()` and `sqrt()`. `is_square()` handles zero and applies the finite-field square +criterion. `sqrt()` requires a square input and asserts that precondition; callers handling arbitrary values should +test first. The implementation uses the field's multiplicative-group decomposition and Tonelli-Shanks rather than +converting the value to a generic integer representation. + +## 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 +`order(K)^d` elements. A polynomial with degree below `d` is the canonical representative of one quotient-field +element. + +`is_square_mod(input, divisor_context, arithmetic_context)` reports whether that element is a square. It treats zero +as a square and uses Euler's criterion for a general nonzero representative. For the canonical indeterminate `X`, the +implementation can use + + Norm(X) = (-1)^d * B(0) / leading_coefficient(B) + +when the coefficient value type provides `is_square()`. This replaces an expensive quotient-ring exponentiation with +one coefficient-field square test. + +Irreducibility of `B` is a caller precondition. The function does not factor `B`, and Euler's criterion does not +characterize squares if the quotient has zero divisors. + +## Quotient-field square roots + +`square_root_mod` uses Tonelli-Shanks in `K[X]/(B)`. Repeated calls should share two immutable precomputations: + +1. `polynomial_divisor_context` caches the divisor and the inverse needed for modular reduction. +2. `polynomial_square_root_context` caches the decomposition + `order(K)^d - 1 = odd_order * 2^two_adicity` and a suitable quadratic nonresidue raised to `odd_order`. + +```cpp +#include +#include + +// B must be irreducible. Its degree is d. +const std::size_t d = B.size() - 1; +const std::size_t inverse_precision = std::max(1, d - 1); + +pa::polynomial_context arithmetic_context; +math::polynomial_divisor_context divisor_context( + B, inverse_precision, arithmetic_context); + +// known_nonresidue is canonical, reduced modulo B, and nonsquare in K[X]/(B). +math::polynomial_square_root_context square_root_context( + known_nonresidue, divisor_context, arithmetic_context); + +polynomial_type root; +if (math::square_root_mod(root, value, square_root_context, + arithmetic_context)) { + // root^2 == value mod B. +} +``` + +Instead of supplying a known nonresidue, a second context constructor accepts a caller-owned generator. The generator +returns canonical reduced polynomial representatives until the context finds a nonsquare. The divisor context must +outlive the square-root context that references it. + +The search has no fixed retry limit. The generator must be capable of producing a quotient-field nonsquare; a +pathological generator can make context construction fail to progress. Applications that need a work limit should +enforce one in the supplied generator. + +The square-root operation supports odd characteristic. It returns `false` and stores zero when the input is not a +square; zero is returned as its own root. Output may alias input. The coefficient field may itself be an extension +field, including BN254 Fq12. + +## Bounded rational reconstruction + +Given a reduced residue `R` modulo a nonconstant polynomial `B`, `rational_reconstruct` searches for polynomials `P` +and `Q` such that + + P = R * Q mod B, + +subject to caller-supplied degree bounds. It follows the Euclidean remainder sequence until the remainder reaches the +numerator bound, while tracking the corresponding coefficient of `R`. The returned denominator is monic, and the +numerator is scaled by the same field element. + +```cpp +#include + +polynomial_type numerator; +polynomial_type denominator; +const bool recovered = math::rational_reconstruct( + numerator, denominator, residue, modulus, + maximum_numerator_degree, maximum_denominator_degree, + arithmetic_context); +``` + +The strict condition + + maximum_numerator_degree + maximum_denominator_degree < degree(modulus) + +ensures uniqueness. The function returns `false` if the Euclidean candidate exceeds the denominator bound and leaves +both outputs unchanged. It rejects noncanonical inputs, a constant modulus, an unreduced residue, nonunique bounds, or +using the same object for both outputs. Either output may otherwise alias an input. + +For example, over any field in which the displayed small integers are distinct, the residue + + R = 2 + X + X^2 + X^3 + +modulo + + B = 3 + 3X + 3X^2 + 2X^3 + X^4 + +reconstructs with bounds `degree(P) <= 0` and `degree(Q) <= 2` as + + P = 1, + Q = 2 + 2X + X^2. + +## Recovering a polynomial norm representation + +The square-root and reconstruction APIs fit together as follows. For one irreducible factor `B` of degree `d`: + +1. Represent the indeterminate by `x = {0, 1}` and test whether it is a square modulo `B`. +2. If it is square, compute `R` such that `R^2 = X mod B`. +3. Rationally reconstruct `R` with + + maximum_numerator_degree = floor(d / 2), + maximum_denominator_degree = floor((d - 1) / 2). + + This produces `P = R * Q mod B`. +4. Squaring the congruence gives + + P^2 - X * Q^2 = 0 mod B. + + The degree bounds make the left side have degree at most `d`, so it is a scalar multiple of `B`. + +This is the local representation of `B` by the norm from adjoining a square root of `X`. Representations compose +multiplicatively: + + (P1^2 - X Q1^2) * (P2^2 - X Q2^2) + = (P1 P2 + X Q1 Q2)^2 - X * (P1 Q2 + P2 Q1)^2. + +Consequently, a caller can factor a target polynomial, recover eligible irreducible factors independently, account for +multiplicities and the scalar leading coefficient, and combine the local representations. Crypto3.Math deliberately +keeps that application-level policy separate from the generic factorization, quotient-field square-root, and bounded +rational-reconstruction primitives. + +### Worked example over F7 + +Consider + + H = X^2 + 1 + +over the field with seven elements. Its only possible roots would square to `-1 = 6`, but the squares in F7 are +`0, 1, 2, 4`, so `H` is irreducible and its complete factorization contains just `H`. + +In the quotient field `F7[X]/(H)`, the polynomial + + R = 2 + 2X + +is a square root of `X`: + + R^2 = 4 + 8X + 4X^2 = X mod H. + +For `d = 2`, reconstruction uses numerator degree at most one and denominator degree at most zero. It returns + + P = 2 + 2X, + Q = 1, + +so that + + P^2 - X Q^2 = 4 + 7X + 4X^2 = 4H. + +The reconstruction determines a representation up to this nonzero scalar. Since `4^-1 = 2` and `3^2 = 2` in F7, +scale both recovered polynomials by `3`: + + P' = 6 + 6X, + Q' = 3. + +Then the scalar is corrected exactly: + + P'^2 - X Q'^2 = H. + +This example performs the same local steps required for a higher-degree factor: factor, test whether `X` is square in +the quotient, compute its root, reconstruct bounded `P` and `Q`, and normalize the remaining coefficient-field +scalar. Combining several eligible factors then uses the multiplicative identity above. + +## Reuse and performance + +The [polynomial arithmetic infrastructure](@ref math_polynomial_arithmetic) describes these contexts in detail. Reuse +one polynomial arithmetic context throughout a recovery, one divisor context for every operation modulo the same `B`, +and one square-root context for repeated roots modulo that divisor. This avoids rebuilding divisor inverses, +multiplicative-group decompositions, nonresidue powers, backend plans, and scratch storage. + +Let `Q` be the coefficient-field order, `d` the irreducible divisor degree, and `s` the two-adicity of `Q^d - 1`. +The current implementations have the rough bounds shown below. + +The table counts quotient-ring multiplications and Euclidean steps. The cost of each underlying polynomial +multiplication, reduction, and division 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 | +| 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.