Factoring polynomials over the integers

Factoring a polynomial f \in \mathbb{Z}[x] means writing it as f = c \prod_i f_i^{e_i} with c \in \mathbb{Z} and distinct primitive irreducible f_i of positive degree. By Gauss' lemma this is the same as factoring over \mathbb{Q}, up to the content. All practical methods follow the plan of Zassenhaus [Zas1969]: factor f modulo a small prime p, which is easy (Factoring polynomials over finite fields), lift the modular factorisation to one modulo a large power p^a (Hensel lifting), and find out which products of modular factors are true factors. The last step, the recombination, is the difficult one: a polynomial with r modular factors has 2^r candidate subsets. Lenstra, Lenstra and Lovász showed that factoring is possible in polynomial time [LLL1982], but their algorithm was not competitive. Van Hoeij [vHo2002] turned the recombination into a knapsack problem solved by lattice reduction, which is fast in practice and, with the refinements of Belabas, van Hoeij, Klüners and Steel [BHKS2009], van Hoeij and Novocin [vHN2010], and Hart, van Hoeij and Novocin [HvHN2011], provably polynomial. The implementation in FLINT is due to Novocin and Hart and follows [HvHN2011] and [Nov2008].

Throughout, f has degree n and leading coefficient \ell = \operatorname{lc}(f) > 0, and \|f\|_2 is the Euclidean norm of its coefficient vector.

Reductions#

fmpz_poly_factor first reduces to the factorisation of squarefree primitive polynomials with nonzero constant term:

1. The content is removed, with the sign of \ell. Polynomials of degree at most 3 are handled directly: a quadratic splits if and only if its discriminant is a square, and a cubic is reducible if and only if it has a rational root, which _fmpz_poly_factor_cubic locates from real and 2-adic approximations of the roots (Cubics).
2. The factor x^k is split off.
3. If f = g(x^d) with d > 1, the polynomial g of degree n/d is factored, and the inflation h(x^d) of each irreducible factor h of g is handled separately (Deflation and Capelli's theorem).
4. Otherwise fmpz_poly_factor_squarefree computes the squarefree factorisation f = \prod_i g_i^i by Yun's algorithm [Yun1976] (Squarefree factorisation), with the GCDs of GCDs of polynomials over Z, and each g_i is factored by _fmpz_poly_factor_zassenhaus.

Deflation and Capelli's theorem#

If h is irreducible, h(x^d) need not be: x - 4 is irreducible but x^2 - 4 is not. For prime exponents, Capelli's theorem decides the question [Sch2000].

FLINT applies the prime case along a chain h(x^d) = U(x^{d/p}) with U = h(x^p). A witness that \theta is not a p-th power is a prime q and a residue field \mathbb{F}_{q^f} of K above q in which the image of \theta is not a p-th power, which by Euler's criterion means \theta^{(q^f - 1)/p} \ne 1 when p \mid q^f - 1. _fmpz_poly_factor_inflation_is_irreducible_capelli searches for a witness with a budget linear in \deg h; if it finds one, h(x^p) is irreducible, and otherwise h(x^p) is factored by the general algorithm at degree p \deg h, which is still smaller than factoring f from scratch.

Cubics#

_fmpz_poly_factor_cubic looks for rational roots instead of factoring modulo a prime, by a method from Schultz's notes [Sch-mpoly]. Let f = a x^3 + b x^2 + c x + d with a > 0. With X = 3a x + b,

and a rational root of the monic g = X^3 - 3AX - B \in \mathbb{Z}[X] is an integer. Each integer root r of g gives the factor \operatorname{pp}(3a x + b - r) of f. The integer roots of g are found as follows (the cases B = 0 and B = \pm 1 are immediate).

1. Small primes. If A and B are odd, g \equiv X^3 + X + 1 \pmod 2 has no root, so f is irreducible. If g has no root modulo one of the primes 5, 7, \ldots, 29, f is irreducible; this is tested at once by computing \prod_{x=0}^{28} g(x) modulo the product of these primes, which fits in a word.
2. Discriminant. Let D = B^2 - 4A^3, so that the discriminant of g is -27D. If D = 0, the roots are B/A and the double root -B/(2A). If D > 0, g has one real root, \bigl(\sqrt[3]{4B + 4\sqrt D} + \sqrt[3]{4B - 4\sqrt D}\bigr)/2 by Cardano's formula; it is computed with integer square and cube roots, and the three integers around the approximation are tested.
3. Three real roots (D < 0, hence A > 0). With X = 2\sqrt{A} \cos \theta the equation becomes 2 A^{3/2} \cos 3\theta = B, so all roots satisfy |r| \le 2\sqrt A. After a test modulo the primes 31, \ldots, 47, the roots are computed in the 2-adic integers \mathbb{Z}_2 to a precision 2^k with 2^{k - 1} > 2\sqrt A, reduced to the symmetric range, and tested exactly.

For the 2-adic roots write -3A = 2^{\alpha} a' and -B = 2^{\beta} b' with a', b' odd, so that g = X^3 + 2^\alpha a' X + 2^\beta b'. The Newton polygon of g at 2 gives three cases.

• 2\beta = 3\alpha: with X = 2^{\beta/3} Y, g = 2^\beta (Y^3 + a' Y + b'), which has no root modulo 2; g has no root in \mathbb{Z}_2, and f is irreducible.
• 2\beta < 3\alpha: all roots have valuation \beta/3, so if 3 \nmid \beta there is none in \mathbb{Z}_2. Otherwise X = 2^{\beta/3} Y gives 2^\beta (Y^3 + 2^{\alpha - 2\beta/3} a' Y + b'), which is (Y + 1)(Y^2 + Y + 1) modulo 2 with an irreducible quadratic factor, so there is exactly one root in \mathbb{Z}_2, found by Hensel lifting.

• 2\beta > 3\alpha: with X = 2^{\beta - \alpha} Y and e = 2\beta - 3\alpha > 0, g = 2^\beta (2^e Y^3 + a' Y + b'), which has the simple root Y \equiv 1 modulo 2. Newton iteration, doubling the precision at each step, lifts the factorisation

  2^e Y^3 + a' Y + b' = (Y + r)(2^e Y^2 - 2^e r Y + s), \qquad r, s \in \mathbb{Z}_2 \text{ odd},

  that is, g = (X + 2^{\beta - \alpha} r)(X^2 - 2^{\beta - \alpha} r X + 2^\alpha s). The roots of the quadratic have valuation \alpha/2. If \alpha is odd, -2^{\beta - \alpha} r is the only root in \mathbb{Z}_2. If \alpha is even, the other roots are

  2^{\beta - \alpha - 1} r \pm 2^{\alpha/2} \sqrt{2^{2\beta - 3\alpha - 2} r^2 - s},

  with a 2-adic square root (an odd 2-adic integer is a square if and only if it is 1 modulo 8). The square root loses precision, so r and s are lifted further when the result is not precise enough.

Once an integer root r_0 has been found, the other roots are those of X^2 + r_0 X + r_0^2 - 3A, rational if and only if 3(4A - r_0^2) is a square. The method costs a bounded number of integer square and cube roots and 2-adic Newton iterations on numbers of the size of the coefficients, so it is much faster than the general algorithm for cubics with large coefficients.

The modular factorisation#

Let f be squarefree and primitive with f(0) \ne 0. A prime p is admissible if p \nmid \ell, p \nmid f(0) and f \bmod p is squarefree; all but finitely many primes are, since f \bmod p is squarefree when p divides neither \ell nor the discriminant. For an admissible prime,

with distinct monic irreducible g_i \in \mathbb{F}_p[x], and every irreducible factor g of f in \mathbb{Z}[x] reduces to \operatorname{lc}(g) \prod_{i \in S} g_i for a subset S, the subsets of distinct factors being disjoint. In particular f has at most r irreducible factors, and it is irreducible if r = 1.

_fmpz_poly_factor_zassenhaus factors f modulo the first three admissible primes, starting from 2, and keeps the one with the fewest factors, since the cost of recombination grows with r. The factor degrees modulo all three primes are combined: a true factor of degree k is a product of modular factors of total degree k modulo each prime, so k must be a subset sum of the modular degrees for all three. The zassenhaus_prune functions maintain the set of such degrees; if it contains no degree strictly between 0 and n, the polynomial is irreducible, and otherwise subsets of an impossible degree are skipped during recombination. By Chebotarev's density theorem, the factor degrees modulo a random prime are the cycle lengths of a random element of the Galois group, so for a polynomial with Galois group S_n, the typical case, r is about \log n on average.

Bounds and precision#

Recombination recovers a factor from its image modulo p^a, which requires a bound for its coefficients.

fmpz_poly_factor_mignotte computes B, the maximum of this bound over 1 \le j < n with \|f\|_2 rounded up, of |f(0)| and of \ell, so that every proper factor of f has coefficients bounded by B (the extreme coefficients of a factor divide those of f).

The Zassenhaus algorithm therefore lifts to the least a with p^a \ge 2 \ell B + 1.

Zassenhaus recombination#

The number of subsets tried is at most \sum_{k \le r/2} \binom{r}{k} \le 2^{r-1}, and each test is a trial division, made cheap by fmpz_poly_divides, which first compares the constant coefficients and the values at 1. This is fast when r is small, the typical case, and hopeless when r is large. FLINT uses Zassenhaus recombination for r \le 8.

Van Hoeij's algorithm#

Van Hoeij's idea [vHo2002] is to represent subsets by indicator vectors v \in \{0, 1\}^r and to find the vectors of the true factors as short vectors of a lattice. Let g_1, \ldots, g_r be monic with f \equiv \ell \prod g_i \pmod{p^a}, and for S \subseteq \{1, \ldots, r\} let g_S = \prod_{i \in S} g_i.

Logarithmic derivatives#

The map g \mapsto f g'/g turns products into sums:

The coefficients of the logarithmic derivative (CLD) are the coefficients of the polynomials f g_i' / g_i, which have degree less than n and are well defined modulo p^a, since g_i divides f there. If g is a true factor of f, then f g'/g = (f/g) g' \in \mathbb{Z}[x], and its coefficients are much smaller than p^a: on |x| = \rho, f g'/g = f \sum_\alpha 1/(x - \alpha) over the roots \alpha of g, and splitting the coefficients of f at index j gives the bound computed by fmpz_poly_CLD_bound [HvHN2011]: with f = \sum a_i x^i and any \rho > 0, the coefficient of x^j in f g'/g is at most

which is minimised numerically over \rho. For the k lowest and k highest coefficients j, let c_{ij} be the coefficient of x^j in f g_i'/g_i reduced into (-p^a/2, p^a/2] (_fmpz_poly_factor_CLD_mat). If v is the indicator vector of a true factor g, then \sum_i v_i c_{ij} \equiv [x^j]\, f g'/g \pmod{p^a} is small. The low coefficients are computed with fmpz_poly_mullow and fmpz_poly_divlow_smodp, the high ones with fmpz_poly_mulhigh_n and fmpz_poly_divhigh_smodp, so that only O(k) coefficients of each product are needed.

The knapsack lattice#

For a set J of CLD columns, let L be the lattice spanned by the rows of

with each column of C scaled down so that the values produced by true factors become small. The vector (2^U v, \ldots) of a true factor then has length at most about 2^U \sqrt{r + 1}, while a random 0-1 combination produces entries of the size of p^a. The vectors of the true factors span a sublattice W whose rank is the number of irreducible factors, and the reduced echelon basis of W consists of the indicator vectors of the factors, a 0-1 basis. Lattice reduction finds a basis of a lattice L' \supseteq W of decreasing rank, by the following lemma.

After each LLL reduction, fmpz_lll_with_removal_ulll drops the final basis vectors whose Gram–Schmidt norm exceeds the bound \beta^2 = (r + 1) 2^{2U} (LLL with removals, with the truncation strategy Truncation: ULLL); by the lemma, W survives. Each new column of CLD data can only lower the rank of L', and enough columns at enough precision make L' = W [vHo2002], [HvHN2011].

The test in step 2 is fmpz_poly_factor_van_hoeij_check_if_solved: it reads a partition from the echelon form of the first r columns (fmpz_mat_col_partition), forms the products of the lifted factors in each part, multiplied by \ell and reduced modulo p^a, and checks that their primitive parts divide f, first modulo 2 and then over \mathbb{Z}. Since the answer is verified by division, the output is always correct; the bounds only ensure that the algorithm makes progress and terminates.

Gradual feeding, adding one column at a time, keeps each LLL call on a matrix whose new column has only a few more bits than the reduced basis, which is what makes the total work polynomial [vHN2010], [HvHN2011]. The starting precision is about the least a with

where b is the bit size of the smaller of the CLD bounds of the lowest and the highest coefficient, but not more than the Zassenhaus precision. This is usually far below the Mignotte bound and enough for the first few columns to determine the factorisation; otherwise the precision is doubled and the Hensel lift continued with _fmpz_poly_hensel_continue_lift. The first rounds use k = 30 (or 50 for r > 200) coefficients from each end instead of 10 when 3r > n + 1.

Complexity#

For a polynomial of degree n with b-bit coefficients, the Mignotte precision is a \log p = O(n + b) bits, and the Hensel lift to that precision costs O(M(n) \log r) operations on integers of that size (Hensel lifting). Zassenhaus recombination adds up to 2^{r-1} trial divisions. Van Hoeij's algorithm with gradual feeding runs in time polynomial in n and b [HvHN2011], and usually needs much less than the Mignotte precision; the LLL work is small unless r is large and there are many factors of similar degree. The classical worst cases of the Zassenhaus algorithm, such as the Swinnerton-Dyer polynomials of degree n = 2^k, which have n/2 factors modulo every prime, are easy: proving S_7(x + 1), of degree 128, irreducible takes about 0.1 seconds on the x86-64 machine used for this documentation.

Use in FLINT#

• fmpz_poly_factor performs the reductions of Reductions and calls _fmpz_poly_factor_zassenhaus with cutoff 8 and van Hoeij's algorithm enabled; fmpz_poly_factor_zassenhaus uses only the Zassenhaus recombination.
• _fmpz_poly_factor_zassenhaus chooses the prime, lifts with fmpz_poly_hensel_lift_once to the Mignotte precision and calls fmpz_poly_factor_zassenhaus_recombination_with_prune, or calls fmpz_poly_factor_van_hoeij when r > 8.
• fmpz_poly_factor_van_hoeij uses _fmpz_poly_hensel_start_lift, _fmpz_poly_hensel_continue_lift, _fmpz_poly_factor_CLD_mat, fmpz_mat_next_col_van_hoeij, fmpz_lll_with_removal_ulll (with \delta = 0.75 and \eta = 0.81; see Lattice reduction (LLL)) and fmpz_poly_factor_van_hoeij_check_if_solved.
• fmpz_poly_factor_squarefree, fmpz_poly_factor_mignotte, fmpz_poly_CLD_bound, _fmpz_poly_factor_inflation_is_irreducible_capelli, _fmpz_poly_factor_quadratic and _fmpz_poly_factor_cubic provide the remaining steps.
• Multivariate polynomials over \mathbb{Z} are factored by fmpz_mpoly_factor, which reduces to univariate factorisations of this kind (Multivariate polynomials: representations and algorithms); polynomials over number fields and finite fields have their own modules.