Factoring polynomials over finite fields

Let \mathbb{F}_q be a finite field with q = p^k elements and let f \in \mathbb{F}_q[x] have degree n. Every factorisation algorithm used by FLINT reduces the problem to three steps, which go back to Berlekamp [Ber1967] [Ber1970] and to Cantor and Zassenhaus [CanZas1981]:

1. Squarefree factorisation writes f = c \prod_i g_i^{e_i} with squarefree, pairwise coprime g_i and distinct e_i (Squarefree factorisation).
2. Distinct-degree factorisation (DDF) splits a squarefree g into g = \prod_d g_{[d]}, where g_{[d]} is the product of the irreducible factors of g of degree d.
3. Equal-degree factorisation (EDF) splits each g_{[d]} into its irreducible factors, by a randomised algorithm.

The algorithms differ in how they perform steps 2 and 3. The Cantor–Zassenhaus algorithm computes the DDF with one exponentiation by q per degree; the Kaltofen–Shoup algorithm [KalSho1998] [Sho1995] uses a baby-step giant-step method in which the exponentiations are replaced by modular compositions; Berlekamp's algorithm skips the DDF and splits g with the kernel of a linear map. Textbook accounts are [GG2003, Ch. 14] and [Sho2009, Ch. 20].

All of them are implemented once, generically, in gr_poly (gr_poly_factor_finite_field and the functions it calls), for any gr context that is a finite field: nmod, fmpz_mod, fq, fq_nmod and fq_zech. The modules nmod_poly_factor, fmpz_mod_poly_factor and the fq_*_poly_factor modules are wrappers. Throughout, M(n) denotes the cost of multiplying polynomials of degree n over \mathbb{F}_q, counted in field operations, and C(n) the cost of one modular composition a(b) \bmod f for a, b of degree less than n (Modular composition).

Finite fields and the Frobenius map#

The gcd is computed with x^{q^d} \bmod g, never with the huge polynomial x^{q^d} - x itself. The powers x^{q^d} \bmod g are linked by the Frobenius endomorphism a \mapsto a^q of \mathbb{F}_q[x]:

So there are two ways to step from h_i to h_{i+1}: an exponentiation h_{i+1} = h_i^q \bmod g, which costs about \tfrac{3}{2} \log_2 q multiplications modulo g, and a composition h_{i+1} = h_i(h_1) \bmod g, which costs C(n) independently of q. With the Brent–Kung algorithm and the matrix of powers of h_1 precomputed, a composition costs about \sqrt{n} multiplications modulo g plus a product of an \sqrt{n} \times \sqrt{n} and an \sqrt{n} \times n matrix, so composition wins when \log q is large compared to \sqrt{n}. (The Frobenius map on elements of \mathbb{F}_q itself is discussed in Frobenius map, norm, trace and roots.) The generic code uses composition when q has more than \tfrac{3}{4}(\lfloor \sqrt{n} \rfloor + 1) bits.

Squarefree factorisation and deflation#

Squarefree factorisation needs only gcds, divisions, derivatives and, in characteristic p, p-th roots of coefficients; it costs O(M(n) \log n) operations per gcd and is described in Squarefree factorisation. It reduces the problem to monic squarefree polynomials, and all steps below assume this.

Before the squarefree factorisation, _gr_poly_factor_finite_field deflates the input: if F(x) = G(x^m) for some m \ge 2 (all exponents of nonzero terms are multiples of m), it factors G, and then factors g_i(x^m) for each factor g_i of G (without deflating again). Each factor of F divides exactly one g_i(x^m), so this splits the problem into independent problems of smaller size; polynomials like x^m - a or x^{2m} + a x^m + b benefit most.

Distinct-degree factorisation#

The classical algorithm#

By induction on i, when u is computed v has no factors of degree less than i, so by Corollary 2 u is the product of the factors of degree exactly i. When \deg v < 2(i + 1), the remaining v has no factor of degree at most i and therefore cannot be a product of two or more factors; it is irreducible. The loop runs at most n/2 times, and each iteration costs an exponentiation, O(M(n) \log q), and a gcd, for a total of O(n M(n) \log q).

Baby steps and giant steps#

Shoup [Sho1995] and Kaltofen and Shoup [KalSho1998] remove the factor \log q by grouping the degrees into intervals. Choose \ell \approx \sqrt{n/2} and m = \lceil n / (2\ell) \rceil, so that the m intervals (\ell j, \ell (j+1)], 0 \le j < m, cover all degrees up to n/2.

The baby steps are computed in _gr_poly_iterated_frobenius_preinv, by repeated exponentiation by q if q is small (see Finite fields and the Frobenius map) and otherwise by modular composition using a doubling scheme: given h_1, \ldots, h_s, the vector composition h_{s+i} = h_i(h_s), 1 \le i \le s, is a single matrix product with the powers of h_s (gr_poly_preinv_compose_mod_brent_kung_vec). The giant steps always use composition with the precomputed matrix of powers of h_\ell.

With C(n) = O(\sqrt{n} M(n) + n^{(\omega+1)/2}) the bound of Theorem 5 is O(n M(n) + n^{(\omega+2)/2}) and no longer depends on q except through h_1; this is what makes the algorithm fast for large q. Kaltofen and Shoup also show how to make the interval products subquadratic, which together with fast matrix multiplication gives a factorisation algorithm with O(n^{1.815} \log q) operations [KalSho1998]; the asymptotically fastest known method, based on the modular composition algorithm of Kedlaya and Umans, uses O(n^{1.5 + o(1)} + n^{1 + o(1)} \log q) bit operations per field operation [KU2011]. FLINT uses neither of these refinements.

Equal-degree factorisation#

Let now g be monic, squarefree, of degree n = rd, the product of r \ge 2 irreducible factors u_1, \ldots, u_r of degree d. By the Chinese remainder theorem

The Berlekamp subalgebra B = \{ b \in R : b^q = b \} consists of the elements whose components all lie in \mathbb{F}_q, so B \cong \mathbb{F}_q^r. Every b \in B whose components are not all equal gives a proper factor \gcd(g, b - c) for some c \in \mathbb{F}_q; the algorithms below produce such b at random and then find a c, or a polynomial in b playing the same role.

By Lemma 7, each attempt succeeds with probability at least 4/9, and because each factor goes to either side independently, the expected depth of the recursion is O(\log r). The cost of an attempt is that of the trace, an exponentiation, O(M(n) \log q), and a gcd. FLINT computes the trace in one of three ways (_gr_poly_trace):

• for small q (the same test as in Finite fields and the Frobenius map), as a + a^q + \cdots with d - 1 exponentiations by q;
• for d < 16, with d - 1 compositions a^{q^{i+1}} = a^{q^i}(h_1) using the precomputed matrix of powers of h_1;
• otherwise by the doubling method of von zur Gathen and Shoup [GatSho1992]: with T_s(a) = a + a^q + \cdots + a^{q^{s-1}}, the identities T_{2s}(a) = T_s(a) + T_s(a)(x^{q^s}), T_{s+1}(a) = a + T_s(a)(h_1) and x^{q^{2s}} = x^{q^s}(x^{q^s}) compute T_d(a) along the binary expansion of d with O(\log d) compositions.

Roots (d = 1). Every element of R is then in B, so no trace is needed, and b = x + c with random c suffices (for even q, b = \alpha x + c with random \alpha \ne 0): the components of b are the values \rho + c at the roots \rho of g, so \gcd(g, (x + c)^{(q-1)/2} - 1) collects the roots for which \rho + c is a nonzero square. For q of more than 16 bits, the power is computed as (x + c)^e \bmod g = R(x + c) where R = x^e \bmod g(x - c), because powers of x modulo a polynomial are cheaper (the multiplications by x in the exponentiation are shifts); this costs two Taylor shifts. See Roots of polynomials over finite fields.

Splitting by the minimal polynomial. When q is large and the number r of factors is moderate, it is cheaper to split g into all its factors at once [Mar2014]. If b \in B has components c_1, \ldots, c_r, its minimal polynomial over \mathbb{F}_q is \mu(y) = \prod_{c \in C} (y - c), where C is the set of distinct values c_i, and

FLINT computes the sequence s_i = w(b^i \bmod g), 0 \le i < 2r, for a random linear form w, obtains its minimal polynomial with the Berlekamp–Massey algorithm [Mas1969] (which equals \mu unless w is degenerate, in which case the attempt is repeated), finds the roots C of \mu by the d = 1 case above, and splits g by divide and conquer: with \mu_1 = \prod_{c \in C_1} (y - c) for the first half C_1 of the roots, \gcd(g, \mu_1(b) \bmod g) is the product of the u_i with c_i \in C_1, and \mu_1(b) \bmod g is a modular composition. This uses 2r multiplications modulo g and no exponentiation by (q-1)/2 for the splitting element. A random b has r distinct components with probability at least 1 - r^2/(2q), in which case the factorisation is complete after one round. The generic code uses this method when d > 1, q \ge 2^7 and r \le 8 \lambda(d + 1) + \lambda(q), where \lambda denotes the bit length.

Berlekamp's algorithm#

Berlekamp's algorithm [Ber1967] [Ber1970] does not need the DDF. Let g be monic squarefree of degree n with r irreducible factors (of arbitrary degrees). The map \phi(a) = a^q - a on R = \mathbb{F}_q[x]/(g) is \mathbb{F}_q-linear, and its matrix Q - I in the basis 1, x, \ldots, x^{n-1} has as column i the coefficients of x^{qi} \bmod g minus the unit vector e_i.

gr_poly_factor_berlekamp computes x^q \bmod g, the n columns x^{qi} \bmod g by repeated multiplication, and a basis of the kernel with gr_mat_nullspace. If the nullity is 1, g is irreducible. Otherwise it takes random linear combinations b of the basis vectors (uniform in B), tries \gcd(g, b) and then \gcd(g, b^{(q-1)/2} - 1) (or \gcd(g, T(b)) for even q), and recurses on the two parts of the first proper split, recomputing the matrix for each part. The kernel costs O(n^\omega), the matrix O(n M(n)), and each split attempt O(M(n) \log q); since random splits are balanced on average, the expected total is O(n^\omega + n M(n) \log r + M(n) \log q \log r). The cubic linear algebra makes the algorithm uncompetitive for large n in FLINT, where it serves mainly as an independent check.

The complete algorithms and FLINT's choice#

All three algorithms begin with gr_poly_factor_squarefree and treat each squarefree factor separately.

• gr_poly_factor_cantor_zassenhaus runs Algorithm 1 on each squarefree factor, computing h \gets h^q \bmod v by binary exponentiation, and calls the equal-degree factorisation (_gr_poly_factor_equal_deg_with_frob) on each g_{[d]} with d < \deg g_{[d]}, passing x^q \bmod g_{[d]} obtained from the first exponentiation.
• gr_poly_factor_kaltofen_shoup computes x^q \bmod g once, runs Algorithm 2 (_gr_poly_factor_distinct_deg_with_frob), and calls the equal-degree factorisation on each piece with x^q reduced modulo the piece.
• gr_poly_factor_berlekamp is described above.

gr_poly_factor_finite_field (after deflation) uses Cantor–Zassenhaus for degree

and Kaltofen–Shoup otherwise. The rule fits measured crossovers for dense random polynomials over nmod: about n = 125 for q = 2, 70 for q = 3, 32 for q = 7, 16 for q = 31, 8 for q \approx 2^8 and 4 for q \approx 2^{16}. For small q the exponentiation by q is a few squarings and the classical DDF has the smaller constant; for large q the \log q factor dominates and the modular compositions win.

Parallelism. When threads are available (flint_set_num_threads) and the context is thread-safe, polynomials of degree above 200 (GR_POLY_FACTOR_THREADED_CUTOFF) use _gr_poly_factor_distinct_deg_with_frob_threaded: the giant steps are computed in blocks of B (the number of threads) independent compositions H_{j+B} = H_j(H_{B-1}) with a precomputed matrix for H_{B-1}, the interval polynomials of a block in parallel, and the fine DDF in parallel over the intervals; v is only reduced after each block, which costs a little extra work when there are factors of small degree. The equal-degree factorisation processes the pieces on the top of its work stack in parallel, each with its own random state. Matrix products use the threads of gr_mat_mul.

Randomness. Each top-level function initialises its own flint_rand_t with the default seed, so results, including the order of the factors, are reproducible from run to run (but not across FLINT versions). The factors are returned in no particular order.

Use in FLINT#

The generic implementations are gr_poly_factor_finite_field (with _gr_poly_factor_finite_field for a fixed algorithm), gr_poly_factor_cantor_zassenhaus, gr_poly_factor_kaltofen_shoup, gr_poly_factor_berlekamp, gr_poly_factor_squarefree, gr_poly_factor_distinct_deg, gr_poly_factor_equal_deg and gr_poly_factor_equal_deg_prob, with the helpers gr_poly_iterated_frobenius_preinv, _gr_poly_factor_distinct_deg_with_frob and _gr_poly_factor_equal_deg_with_frob. The gr method GR_METHOD_POLY_FACTOR of the finite field contexts points to gr_poly_factor_finite_field. The type-specific modules wrap these functions without copying the polynomials:

The fq_poly_factor, fq_nmod_poly_factor and fq_zech_poly_factor modules are generated from the same template and wrap the same functions; the arithmetic of the finite fields themselves is described in Arithmetic in finite fields. The modular factorisations used by the Zassenhaus algorithm over \mathbb{Z} (Factoring polynomials over the integers) and by the multivariate factorisation code (Factoring over finite fields) call nmod_poly_factor and fmpz_mod_poly_factor. Irreducibility tests are described in Irreducibility testing over finite fields, and root finding, which is the case d = 1 with a cheaper first split, on Roots of polynomials over finite fields.