Berlekamp's algorithm

Berlekamp's algorithm factors a square-free polynomial over a finite field $F_{p}$ using linear algebra. The residues $v$ modulo $f$ with $v^{p} \equiv v$ form a vector space over $F_{p}$ , the Berlekamp subalgebra, whose dimension equals 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 \mapsto v^{p}$ , so $r$ comes from one Gaussian elimination. Any element of the space that is not a constant splits $f$ through greatest common divisors^[1]. The algorithm, published by Berlekamp in 1967, was the first polynomial-time factorization method for small fields, and together with Hensel lifting it made factorization of integer polynomials practical^[7].

§ 01

Overview

Factoring a polynomial over $Z$ or $Q$ is usually done by factoring it modulo a prime and lifting the factors. Over $F_{p}$ there are only finitely many polynomials of each degree, so factoring by trial division is possible in principle but takes time exponential in the degree. Berlekamp's algorithm reduces the problem to linear algebra over $F_{p}$ , with a cost polynomial in the degree $n$ and in $p$ . The input is assumed square-free, which square-free decomposition arranges beforehand without any factoring.

§ 02

The Berlekamp subalgebra

Let $f = f_{1} f_{2} \dots f_{r}$ be square-free with distinct monic irreducible factors $f_{i}$ of degrees $d_{i}$ . By the Chinese remainder theorem^[10],

F_{p} [x] / (f) ≅ i = 1 \prod r F_{p} [x] / (f_{i}) ≅ i = 1 \prod r F_{p^{d_{i}}}

The Frobenius map $v \mapsto v^{p}$ is $F_{p}$ -linear on the ring on the left, since $(a + b)^{p} = a^{p} + b^{p}$ in characteristic $p$ and $c^{p} = c$ for $c \in F_{p}$ . In each field $F_{p^{d_{i}}}$ its fixed points are exactly the elements of $F_{p}$ . So the fixed points in the whole ring are the residues that are constant modulo every $f_{i}$ :

B = {v \in F_{p} [x] / (f) : v^{p} \equiv v (mod f)} ≅ F_{p}^{r}, dim B = r

The number of irreducible factors is therefore a dimension, and it can be computed without finding a single factor.

The matrix Q

On the basis $1, x, \dots, x^{n - 1}$ , the Frobenius map has the matrix $Q$ whose row $i$ holds the coefficients of $x^{i p} mod f$ . A residue $v$ , written as a row vector of coefficients, is fixed exactly when $v Q = v$ , so

r = n - rank (Q - I)

The first row of $Q - I$ is always zero, since the constant 1 is fixed, which accounts for the trivial solution. The matrix was considered by Petr in 1937^[6], and is sometimes called the Petr–Berlekamp matrix.

§ 03

Splitting

Every $v \in B$ is a vector of constants $(s_{1}, \dots, s_{r})$ in the product decomposition, with $v \equiv s_{i} (mod f_{i})$ . Since $x^{p} - x = \prod_{s \in F_{p}} (x - s)$ , $v^{p} - v$ is divisible by $f$ and^[1]^[3]:

f = s \in F_{p} \prod g cd (f, v - s)

where the factor for $s$ is the product of the $f_{i}$ with $s_{i} = s$ . If $v$ is not a constant, two of its components differ, and at least two of the gcds are nontrivial. Running through a basis of $B$ and splitting the current factors with each basis vector separates all $r$ factors, because for any two factors some basis vector has different components on them.

§ 04

Example

Polynomials over F_p, as in the square-free decomposition article.

(* Polynomials over the prime field F_p, as arrays of coefficients,
   constant term first, with no trailing zeros. *)
module Make (P : sig val p : int end) = struct
  let p = P.p
  let md x = ((x mod p) + p) mod p
  let rec pw b e = if e = 0 then 1 else let h = pw (b * b mod p) (e / 2) in if e land 1 = 1 then h * b mod p else h
  let inv a = pw (md a) (p - 2)
 
  let trim a = let n = ref (Array.length a) in while !n > 0 && a.(!n - 1) = 0 do decr n done; Array.sub a 0 !n
  let make l = trim (Array.of_list (List.map md l))
  let deg a = Array.length a - 1
  let zero = [||] and one = [| 1 |]
  let x_plus c = make [ c; 1 ]
  let lead a = a.(deg a)
 
  let add a b = trim (Array.init (max (Array.length a) (Array.length b)) (fun i ->
    md ((if i < Array.length a then a.(i) else 0) + if i < Array.length b then b.(i) else 0)))
  let neg a = Array.map (fun c -> md (-c)) a
  let sub a b = add a (neg b)
  let scale k a = trim (Array.map (fun c -> md (k * c)) a)
  let mul a b =
    if a = zero || b = zero then zero else begin
      let r = Array.make (Array.length a + Array.length b - 1) 0 in
      Array.iteri (fun i x -> Array.iteri (fun j y -> r.(i + j) <- (r.(i + j) + x * y) mod p) b) a;
      trim r end
  let rec pow a e = if e = 0 then one else mul a (pow a (e - 1))
 
  let divmod a b =
    let q = Array.make (max 0 (deg a - deg b + 1)) 0 and r = Array.copy a in
    let il = inv (lead b) in
    for i = deg a - deg b downto 0 do
      let c = r.(i + deg b) * il mod p in
      q.(i) <- c;
      Array.iteri (fun j y -> r.(i + j) <- md (r.(i + j) - c * y)) b
    done;
    (trim q, trim r)
  let div a b = fst (divmod a b)
  let rem a b = snd (divmod a b)
 
  let monic a = if a = zero then a else scale (inv (lead a)) a
  let rec gcd a b = if b = zero then monic a else gcd b (rem a b)
 
  let deriv a = trim (Array.init (max 0 (Array.length a - 1)) (fun i -> md ((i + 1) * a.(i + 1))))
 
  (* If a' = 0, every exponent is a multiple of p, and since c^p = c in
     F_p, a is the p-th power of the polynomial with those coefficients. *)
  let pth_root a = trim (Array.init (deg a / p + 1) (fun i -> a.(i * p)))
 
  (* Printed with coefficients in (-p/2, p/2], highest degree first. *)
  let to_string a =
    if a = zero then "0" else
    let signed c = if c > p / 2 then c - p else c in
    let terms = List.rev (List.filter (fun (_, c) -> c <> 0) (List.mapi (fun i c -> (i, signed c)) (Array.to_list a))) in
    let mono i c =
      let k = abs c in
      let cs = if k = 1 && i > 0 then "" else string_of_int k in
      match i with 0 -> string_of_int k | 1 -> cs ^ "x" | _ -> cs ^ "x^" ^ string_of_int i in
    String.concat "" (List.mapi (fun n (i, c) ->
      (if n = 0 then (if c < 0 then "-" else "") else if c < 0 then " - " else " + ") ^ mono i c) terms)
end

The Q matrix, its null space, and splitting with gcd(g, v - s).

(* Berlekamp's algorithm over F_p, for monic square-free f. *)
module Make (P : sig val p : int end) = struct
  include Fpoly.Make (P)
 
  let rec powmod a e f = if e = 0 then one else
    let h = powmod (rem (mul a a) f) (e / 2) f in
    if e land 1 = 1 then rem (mul h a) f else h
 
  (* Row i of Q holds the coefficients of x^(ip) mod f. *)
  let q_matrix f =
    let n = deg f in
    let xp = powmod (make [ 0; 1 ]) p f in
    let row = ref one in
    Array.init n (fun _ ->
      let r = Array.init n (fun j -> if j < Array.length !row then !row.(j) else 0) in
      row := rem (mul !row xp) f; r)
 
  (* A basis of { v : v (Q - I) = 0 }, by Gaussian elimination on the
     transpose of Q - I. *)
  let kernel q =
    let n = Array.length q in
    let a = Array.init n (fun i -> Array.init n (fun j -> md (q.(j).(i) - if i = j then 1 else 0))) in
    let pivot_col = Array.make n (-1) in
    let row = ref 0 in
    for col = 0 to n - 1 do
      match List.find_opt (fun r -> a.(r).(col) <> 0) (List.init (n - !row) (fun k -> !row + k)) with
      | None -> ()
      | Some r ->
          let t = a.(r) in a.(r) <- a.(!row); a.(!row) <- t;
          let iv = inv a.(!row).(col) in
          a.(!row) <- Array.map (fun x -> x * iv mod p) a.(!row);
          for r' = 0 to n - 1 do
            if r' <> !row && a.(r').(col) <> 0 then begin
              let k = a.(r').(col) in
              a.(r') <- Array.mapi (fun j x -> md (x - k * a.(!row).(j))) a.(r')
            end
          done;
          pivot_col.(!row) <- col; incr row
    done;
    let pivots = Array.to_list (Array.sub pivot_col 0 !row) in
    let free = List.filter (fun c -> not (List.mem c pivots)) (List.init n Fun.id) in
    List.map (fun fc ->
      let v = Array.make n 0 in
      v.(fc) <- 1;
      List.iteri (fun r pc -> v.(pc) <- md (- a.(r).(fc))) pivots;
      trim v) free
 
  (* Split with gcd(g, v - s) for every s in F_p, until there are r factors. *)
  let factor f =
    let basis = kernel (q_matrix f) in
    let r = List.length basis in
    let factors = ref [ f ] in
    List.iter (fun v ->
      if deg v > 0 then
        factors := List.concat_map (fun g ->
          if List.length !factors >= r || deg g <= 1 then [ g ] else
          let parts = List.filter (fun h -> deg h > 0)
            (List.init p (fun s -> gcd g (sub v (make [ s ])))) in
          if List.length parts > 1 then parts else [ g ]) !factors) basis;
    (r, List.sort compare !factors)
end

Knuth's example over F_13.

module B = Berl.Make (struct let p = 13 end)
open B

Running it.

f = x^8 + x^6 - 3x^4 - 3x^3 - 5x^2 + 2x - 5,  gcd(f, f') = 1
Q - I:
   0  0  0  0  0  0  0  0
   2  0  7 11 10 12  5 11
   3  6  3  3  0  4  7  2
   4  3  6  4  1  6  2  3
   2 11  8  8  2  1  3 11
   6 11  8  6  2  6 10  9
   5 11  7 10  0 11  6 12
   3  3 12  5  0 11  9 11
kernel basis (3 vectors):
  1
  x^6 + 5x^5 - 4x^4 + 5x^2 + 5x
  x^7 - x^5 - 3x^4 - 4x^3 - 2x^2 - 4x
r = 3 factors:
  x + 3   (its own r = 1)
  x^3 - 5x^2 + 4x - 1   (its own r = 1)
  x^4 + 2x^3 + 3x^2 + 4x + 6   (its own r = 1)
product equals f: true

This is the example Knuth works through in Seminumerical Algorithms^[3]. The rank of $Q - I$ is 5, so the null space has dimension 3 and the polynomial has three irreducible factors, of degrees 1, 3 and 4. Each factor on its own gives $r = 1$ , which confirms its irreducibility.

Checking the count

Over $F_{3}$ all monic polynomials of degree 6 can be checked. For every square-free one, the dimension of the null space is compared with the number of irreducible factors found by trial division:

An exhaustive check over F_3.

(* For every monic polynomial of degree 6 over F_3: when it is square-free,
   compare dim ker(Q - I) with the number of irreducible factors found by
   trial division. *)
module B = Berl.Make (struct let p = 3 end)
open B
 
let monics d = (* all monic polynomials of degree d *)
  List.init (int_of_float (3. ** float d)) (fun k ->
    let c = Array.make (d + 1) 1 in
    let k = ref k in
    for i = 0 to d - 1 do c.(i) <- !k mod 3; k := !k / 3 done;
    c)
 
let irreducibles = (* monic irreducibles of degree 1 to 3, by trial division *)
  List.fold_left (fun acc d ->
    acc @ List.filter (fun f -> List.for_all (fun g -> rem f g <> zero) acc) (monics d)) [] [ 1; 2; 3 ]
 
let count_factors f =
  let rec go f = function
    | [] -> if deg f > 0 then 1 else 0
    | g :: gs -> if deg g > deg f then go f [] else if rem f g = zero then 1 + go (div f g) (g :: gs) else go f gs
  in
  go f irreducibles

Running it.

729 monic polynomials of degree 6 over F_3, 486 square-free
dim ker(Q - I) equals the number of irreducible factors for all of them: true
  1 factor : 116 polynomials
  2 factors: 226 polynomials
  3 factors: 127 polynomials
  4 factors:  17 polynomials

That the count needs only a rank makes the algorithm useful even when factors are found by other means. The Cantor–Zassenhaus algorithm shares none of its code, so comparing the number of factors it finds with $r$ is an independent check on an implementation.

§ 05

Large fields

Trying every $s \in F_{p}$ costs $p$ gcds per split, which is impractical for large $p$ . Berlekamp's 1970 paper gave methods for large fields^[2], and the randomized approach now standard takes a random $v \in B$ and computes one gcd. For odd $p$ , each component $s_{i}$ of a random $v$ is uniform on $F_{p}$ , and $s_{i}^{(p - 1) /2}$ is $+ 1$ for about half of the nonzero values and $- 1$ for the other half:

g cd (f, v^{(p - 1) /2} - 1) = i : s_{i}^{(p - 1) /2} = 1 \prod f_{i}

which is a proper factor with probability at least about one half^[4]. This is the same splitting step as in the Cantor–Zassenhaus algorithm, applied to an element of $B$ instead of a random polynomial of the right degree^[5].

The randomized split over F_1000003.

(* For large p, trying every s in F_p is too slow. A random element v of
   the Berlekamp subalgebra has v^((p-1)/2) equal to +1 or -1 in each
   component, independently, so gcd(g, v^((p-1)/2) - 1) splits g with
   probability about 1/2. *)
module B = Berl.Make (struct let p = 1_000_003 end)
open B
 
let () = Random.init 5
 
let split_random basis g =
  let v = List.fold_left (fun acc b -> add acc (scale (Random.int p) b)) zero basis in
  let w = sub (powmod (rem v g) ((p - 1) / 2) g) one in
  gcd g w
 
let rec factor_rand basis f =
  if List.length (kernel (q_matrix f)) = 1 then [ f ]
  else
    let rec try_ n = let h = split_random basis f in
      if deg h > 0 && deg h < deg f then (h, n) else try_ (n + 1) in
    let h, _ = try_ 1 in
    factor_rand basis h @ factor_rand basis (div f h)

Running it.

f = x^10 + 7x^9 + 5x^8 + 41x^7 + 49x^6 + 64x^5 + 115x^4 + 81x^3 + 93x^2 + 122x + 70 over F_1000003, square-free: true
r = 5
  x + 7
  x + 410588
  x - 410588
  x^3 + 3x + 5
  x^4 + x + 1
product equals f: true

§ 06

Complexity

Building $Q$ takes the power $x^{p} mod f$ , computed by repeated squaring in $O (lo g p)$ multiplications modulo $f$ , followed by $n$ multiplications by it. The null space costs $O (n^{3})$ field operations by Gaussian elimination. Deterministic splitting costs $O (p r n^{2})$ and randomized splitting an expected $O (r n^{2} lo g p)$ , so the randomized algorithm runs in time polynomial in $n$ and $lo g p$ ^[4]. The cubic linear algebra dominates for large degrees. Kaltofen and Lobo replaced it with Wiedemann's black-box linear algebra, which only multiplies vectors by $Q$ ^[8], and Kaltofen and Shoup gave algorithms subquadratic in the degree^[9].

§ 07

Uses

Factoring polynomials over the integers reduces modulo a prime $p$ for which the polynomial remains square-free, factors there, and lifts the factors to a larger modulus by Hensel lifting, as proposed by Zassenhaus^[7]; the lifted factors are then recombined by trial or, in modern systems, by lattice reduction. The Swinnerton-Dyer polynomials, which split into many factors modulo every prime, are the standard hard case for the recombination step. Finite-field factorization is also used in coding theory, where $x^{n} - 1$ is factored to construct cyclic codes, and in cryptography.

§ 08

History

Petr studied the matrix of the Frobenius map in 1937^[6]. Berlekamp, working on the decoding of error-correcting codes at Bell Labs, published the factorization algorithm in 1967^[1] and methods for large fields in 1970^[2]. Zassenhaus combined it with Hensel lifting to factor integer polynomials in 1969^[7], Knuth presented it in Seminumerical Algorithms^[3], and Cantor and Zassenhaus's 1981 algorithm replaced the linear algebra with exponentiation for most purposes^[5].