Cantor–Zassenhaus algorithm

The Cantor–Zassenhaus algorithm factors polynomials over a finite field ${\mathbb{F}}_p$ by a randomized, divide-and-conquer method^[1]. It works in three stages. Square-free factorization removes repeated factors. Distinct-degree factorization separates the remaining irreducible factors into groups of equal degree, by taking gcds with $x^{p^d}-x$. Equal-degree factorization splits each group: for a random polynomial $a$, the gcd of the group with $a^{(p^d-1)/2}-1$ collects a random subset of its irreducible factors, and repeating on the pieces separates all of them in an expected logarithmic number of rounds. David Cantor and Hans Zassenhaus published the method in 1981^[1], building on Berlekamp's randomized algorithm and Rabin's work on probabilistic algorithms in finite fields^[2][7]. Factoring over finite fields is the first step of factoring over the integers by Hensel lifting.

Over ${\mathbb{F}}_p$, the polynomial $x^{p^d}-x$ is the product of all monic irreducible polynomials whose degree divides $d$, because its roots are exactly the elements of ${\mathbb{F}}_{p^d}$. For a square-free $f$, the gcd of $f$ with $x^p-x$ is therefore the product of its linear factors; removing them and taking the gcd with $x^{p^2}-x$ gives the product of its quadratic factors, and so on^[3]. The powers $x^{p^d}\mathrm{mod}f$ are computed by repeated $p$-th powering modulo $f$, never expanding $x^{p^d}$ itself.

Let $f=f_1\cdots f_r$ be a product of distinct monic irreducibles of degree $d$, with $p$ odd. By the Chinese remainder theorem,

In a finite field of odd order $q=p^d$, the power $b^{(q-1)/2}$ is 1 for the nonzero squares, $-1$ for the non-squares and 0 for 0, and half of the nonzero elements are squares. For a random $a$ of degree less than $\mathrm{deg}f$, the components of $a^{(q-1)/2}$ are independent and each is $\pm 1$ with probability close to $1/2$. The gcd of $f$ with $a^{(q-1)/2}-1$ is the product of the $f_i$ whose component is 1^[1][3].

The split is nontrivial unless all components agree, which happens with probability about $2^{1-r}$. The pieces are split recursively with fresh random elements until each is irreducible. For $p=2$, the power is replaced by the trace map $a+a^2+a^4+\cdots +a^{2^{d-1}}$, which takes each value in ${\mathbb{F}}_2$ half the time^[3].

(* Factoring a square-free polynomial over F_p, p an odd prime:
   distinct-degree factorization, then Cantor–Zassenhaus equal-degree
   splitting. *)
open Modp
 
(* Split f into products of the irreducible factors of each degree d,
   using gcd(f, x^(p^d) - x). *)
let distinct_degree p f =
  let rec go f xp d acc =
    if deg f < 2 * d then (if deg f > 0 then (deg f, monic p f) :: acc else acc)
    else begin
      let xp = powmod p xp p f in                         (* x^(p^d) mod f *)
      let g = gcd p f (sub p xp [| 0; 1 |]) in
      if deg g > 0 then go (fst (divmod p f g)) (rem p xp (fst (divmod p f g))) (d + 1) ((d, g) :: acc)
      else go f xp (d + 1) acc
    end
  in
  List.rev (go (monic p f) [| 0; 1 |] 1 [])
 
let attempts = ref 0
 
(* f is a product of irreducibles of degree d. A random a gives
   a^((p^d - 1)/2) = +1 or -1 in each factor's field independently, so
   gcd(f, a^((p^d-1)/2) - 1) splits f with probability about 1/2. *)
let rec equal_degree p d f =
  if deg f <= d then [ monic p f ]
  else begin
    incr attempts;
    let a = norm (Array.init (deg f) (fun _ -> Random.int p)) in
    let pd = Z.pow (Z.of_int p) d in
    let exponent = Z.div (Z.sub pd Z.one) (Z.of_int 2) in
    let rec powz acc b n = if Z.equal n Z.zero then acc
      else powz (if Z.is_odd n then rem p (mul p acc b) f else acc) (rem p (mul p b b) f) (Z.shift_right n 1) in
    let b = powz [| 1 |] (rem p a f) exponent in
    let g = gcd p f (sub p b [| 1 |]) in
    if deg g > 0 && deg g < deg f then equal_degree p d g @ equal_degree p d (fst (divmod p f g))
    else equal_degree p d f
  end
 
let factor p f =
  List.concat_map (fun (d, g) -> equal_degree p d g) (distinct_degree p f)

(* Factor x^15 - 1 over F_31 and over F_101, showing the degree groups
   produced by distinct-degree factorization. *)

x^15 - 1 over F_31
  degree-1 part: x^15 + 30
  15 factors, 21 random trials: (x + 26) (x + 13) (x + 15) (x + 6) (x + 11) (x + 21) (x + 24) (x + 12) (x + 3) (x + 27) (x + 17) (x + 30) (x + 29) (x + 22) (x + 23)
x^15 - 1 over F_101
  degree-1 part: x^5 + 100
  degree-2 part: x^10 + x^5 + 1
  10 factors, 12 random trials: (x + 100) (x + 17) (x + 14) (x + 6) (x + 65) (x^2 + 87x + 95) (x^2 + x + 1) (x^2 + 36x + 84) (x^2 + 95x + 36) (x^2 + 84x + 87)

Both results agree with SymPy. The degrees of the factors are determined by arithmetic: a root of unity of order $n$ lies in ${\mathbb{F}}_{p^d}$ exactly when $p^d\equiv 1(\mathrm{mod}n)$. Since $31\equiv 1(\mathrm{mod}15)$, all fifteen roots are in ${\mathbb{F}}_{31}$; since $101\equiv 1(\mathrm{mod}5)$ but $101\equiv 2(\mathrm{mod}3)$, over ${\mathbb{F}}_{101}$ the fifth roots are linear and the other ten roots pair into quadratics.

(* How often does one random element split a product of r irreducible
   quadratics over F_10007? A split fails when all r components agree. *)
let p = 10007
 
let random_irreducible d =
  let rec try_ () =
    let f = Modp.norm (Array.init (d + 1) (fun i -> if i = d then 1 else Random.int p)) in
    if List.length (Cz.factor p f) = 1 then f else try_ () in
  try_ ()

r   splits per 1000 trials   predicted 1 - 2^(1-r)
2             517                      500
3             777                      750
4             865                      875
6             959                      969
8             996                      992

The dominant cost is computing powers $a^{(q-1)/2}\mathrm{mod}f$ and $x^{p^d}\mathrm{mod}f$, each needing $O(d\log p)$ multiplications of polynomials of degree less than $n=\mathrm{deg}f$. With classical arithmetic the whole algorithm takes $O(n^3\log p)$ expected operations in ${\mathbb{F}}_p$ up to logarithmic factors. Von zur Gathen and Shoup reduced the cost of the Frobenius powers with modular composition^[5], Kaltofen and Shoup obtained subquadratic exponents in $n$^[4], and Kedlaya and Umans a nearly optimal exponent of $1.5$^[6]. Shoup's NTL implementation made factoring polynomials of degree in the thousands routine^[8].

Berlekamp's algorithm solves a linear system to find the subalgebra of ${\mathbb{F}}_p[x]/(f)$ fixed by the Frobenius map, whose dimension is the number of factors, and splits $f$ with gcds of its elements^[2]. Its linear algebra is cubic in the degree. Cantor–Zassenhaus avoids the matrix and replaces it with powering, which is simpler and parallels the probabilistic version of Berlekamp's own method for large fields. For small $p$, Berlekamp's deterministic algorithm remains competitive; for large $p$, no deterministic polynomial-time algorithm is known, and randomization is essential.

Berlekamp gave the first practical factoring algorithm for finite fields in 1967 and a randomized variant for large fields in 1970^[2]. Rabin analyzed probabilistic algorithms for finding roots and irreducible factors in 1980^[7]. Cantor and Zassenhaus's paper of 1981 combined distinct-degree factorization with the equal-degree splitting now named after them^[1]. The asymptotically faster algorithms of von zur Gathen, Kaltofen, Shoup, Kedlaya and Umans followed between 1992 and 2011^[5][4][6].

• Berlekamp's algorithmBerlekamp's algorithm factors a square-free polynomial f over a finite field F_p by linear algebra. The polynomials v with v^p = v modulo f form a vector space whose dimension is exactly the number r of irreducible factors of f; it is the null space of Q - I, where Q is the matrix of the Frobenius map v -> v^p. Any nonconstant v in it splits f as the product of gcd(f, v - s) over s in F_p. For large p a random element of the space and one exponentiation split f with probability about one half.
• Hensel liftingHensel lifting turns a factorization of a polynomial modulo a prime p into a factorization modulo p², p⁴ and so on, as long as the factors modulo p are coprime. Each step solves a linear equation using the Bézout relation between the factors, and in quadratic form it doubles the precision, the p-adic analogue of Newton's method. It is the core of Zassenhaus's algorithm for factoring integer polynomials: factor modulo a small prime, lift until the modulus exceeds twice the Landau–Mignotte bound, and recombine the lifted factors into true factors over the integers.
• Square-free decompositionThe square-free decomposition of a polynomial writes it as f = f_1 f_2^2 f_3^3 ..., where the f_i are square-free and pairwise coprime, so f_i collects the irreducible factors of multiplicity i. It needs only greatest common divisors with the derivative, not factorization, because gcd(f, f') contains every repeated factor with its multiplicity reduced by one. Yun's algorithm computes it at the cost of about one gcd. In characteristic p a p-th power has zero derivative, and the algorithm must take p-th roots.
• Swinnerton-Dyer polynomialThe Swinnerton-Dyer polynomials are the integer polynomials SD_k(x) = ∏(x ± √p₁ ± √p₂ ± … ± √p_k), taken over all 2^k sign choices, for the first k primes. They are irreducible over the integers but factor into polynomials of degree at most two modulo every prime, because their Galois group (ℤ/2)^k has no element of order greater than two. They are the worst case for the Berlekamp–Zassenhaus factoring algorithm: every prime yields at least 2^(k−1) modular factors, and before concluding that the polynomial is irreducible the algorithm must try exponentially many subsets of them. They motivated the lattice-based factoring algorithms of Lenstra, Lenstra and Lovász and of van Hoeij.
• Landau–Mignotte boundThe Landau–Mignotte bound limits the coefficients of any factor of an integer polynomial in terms of the polynomial itself: if g of degree k divides f, then each coefficient of g is at most binom(k, j) times the Euclidean norm of f, up to the ratio of leading coefficients, so all are below 2^k ‖f‖₂. It follows from Landau's inequality for the Mahler measure. Computer algebra algorithms use it to decide how large a modulus is needed to recover a factor or a gcd from its image modulo a prime power or a product of primes, as in Hensel lifting and the modular gcd. Factors can have larger coefficients than the polynomial they divide, so some such bound is necessary.

1. [1]D. G. Cantor, H. Zassenhaus, “A new algorithm for factoring polynomials over finite fields”, Mathematics of Computation 36 (1981).
2. [2]E. R. Berlekamp, “Factoring polynomials over large finite fields”, Mathematics of Computation 24 (1970).
3. [3]J. von zur Gathen, J. Gerhard, Modern Computer Algebra, ch. 14, Cambridge University Press (3rd ed., 2013).
4. [4]E. Kaltofen, V. Shoup, “Subquadratic-time factoring of polynomials over finite fields”, Mathematics of Computation 67 (1998).
5. [5]J. von zur Gathen, V. Shoup, “Computing Frobenius maps and factoring polynomials”, Computational Complexity 2 (1992).
6. [6]K. S. Kedlaya, C. Umans, “Fast polynomial factorization and modular composition”, SIAM Journal on Computing 40 (2011).
7. [7]M. O. Rabin, “Probabilistic algorithms in finite fields”, SIAM Journal on Computing 9 (1980).
8. [8]V. Shoup, “A new polynomial factorization algorithm and its implementation”, Journal of Symbolic Computation 20 (1995).

last updated 28 September 2026