Hensel lifting

Hensel lifting is the process of refining a factorization modulo a prime $p$ into factorizations modulo higher and higher powers of $p$ ^[2]. If $f \equiv g h (mod p)$ and $g$ and $h$ are coprime modulo $p$ , then for every $k$ there are unique factors $g_{k} \equiv g$ and $h_{k} \equiv h$ with $f \equiv g_{k} h_{k} (mod p^{k})$ . This is Hensel's lemma, from Kurt Hensel's theory of $p$ -adic numbers of 1908^[3]. Each lifting step solves a linear congruence, and the quadratic version doubles the exponent each time, as Newton's method doubles the number of correct digits. Hans Zassenhaus made it the basis of polynomial factorization over the integers in 1969^[1]: a factorization modulo a small prime, which is easy to find, is lifted until the modulus is large enough for the integer coefficients of the true factors, bounded by the Landau–Mignotte bound, to be read off.

§ 01

Roots

The simplest case is a root. If $f (r) \equiv 0 (mod p^{k})$ and $f^{'} (r)$ is invertible modulo $p$ , then the Newton step

r^{'} = r - \frac{f ( r )}{f ^{'} ( r )} (mod p^{2 k})

is a root modulo $p^{2 k}$ , and the only one congruent to $r$ . The Taylor expansion $f (r + ε) = f (r) + f^{'} (r) ε + O (ε^{2})$ shows why: with $ε$ divisible by $p^{k}$ , the error term vanishes modulo $p^{2 k}$ .

Lifting the square root of 2 modulo 7, where 3² ≡ 2.

(* Hensel's lemma for a root: if f(r) = 0 mod p and f'(r) is a unit mod p,
   Newton's step r - f(r)/f'(r) doubles the number of correct p-adic
   digits. Here f = x^2 - 2 and p = 7, since 3^2 = 9 = 2 mod 7. *)
let base_digits p n r = (* the base-p digits of r, lowest first *)
  let rec go r k = if k = 0 then [] else Z.to_int (Z.erem r (Z.of_int p)) :: go (Z.div r (Z.of_int p)) (k - 1) in
  go r n

Running it; the digits are base 7, the lowest on the right.

mod 7^1   r = 3                              digits 3
mod 7^2   r = 10                             digits 13
mod 7^4   r = 2166                           digits 6213
mod 7^8   r = 1802916                        digits 21216213
mod 7^16  r = 9611769806236                  digits 2011266421216213
mod 7^32  r = 891834753593197726755637001    digits 54365536623164112011266421216213

The correct base-7 digits of √2 in the 7-adic integers after each quadratic step. The number of correct digits doubles each time; linear lifting, which gains one digit per step, would need fifteen steps for the same sixteen digits.

§ 02

Factors

For a factorization $f \equiv g h (mod m)$ with $h$ monic and $s g + t h \equiv 1 (mod m)$ , the error $e = f - g h$ is divisible by $m$ . Corrections $g + Δ g$ and $h + Δ h$ change the product by $g Δ h + h Δ g$ modulo $m^{2}$ , so it suffices to solve $g Δ h + h Δ g \equiv e$ , which the Bézout relation does: $Δ g = t e$ and $Δ h = se$ , reduced so that the degrees do not grow. The Bézout coefficients $s$ and $t$ are lifted along with the factors^[2].

One quadratic step, after von zur Gathen and Gerhard, Algorithm 15.10, together with linear lifting for comparison.

(* One quadratic step: from f = g h and s g + t h = 1 modulo m, with h
   monic, to the same modulo m^2 (von zur Gathen and Gerhard, 15.10). *)
let hensel_step m f g h s t =
  incr steps;
  let m = Z.mul m m in
  let e = red m (sub f (mul g h)) in
  let q, r = divmod_monic m (mulm m s e) h in
  let g' = red m (add g (add (mul t e) (mul q g))) and h' = red m (add h r) in
  let b = red m (sub (add (mul s g') (mul t h')) one) in
  let c, d = divmod_monic m (mulm m s b) h' in
  let s' = red m (sub s d) and t' = red m (sub t (add (mul t b) (mul c g'))) in
  (m, g', h', s', t')
 
(* Lift f = g h mod p to modulo at least target, quadratically. *)
let lift_two p target f g h =
  let pz = Z.of_int p in
  let _, s, t = Modp.ext_gcd p (Modp.of_z p g) (Modp.of_z p h) in
  let s = Array.map Z.of_int s and t = Array.map Z.of_int t in
  let rec go m g h s t = if Z.geq m target then (m, g, h) else
      let m, g, h, s, t = hensel_step m f g h s t in go m g h s t in
  go pz (red pz g) (red pz h) s t
 
(* One linear step: from modulo p^k to p^(k+1), with s and t fixed mod p. *)
let lift_linear p target f g h =
  let pz = Z.of_int p in
  let _, s, t = Modp.ext_gcd p (Modp.of_z p g) (Modp.of_z p h) in
  let rec go m g h =
    if Z.geq m target then (m, g, h)
    else begin
      incr steps;
      let e = Modp.of_z p (Array.map (fun c -> Z.divexact c m) (sub (red (Z.mul m pz) f) (red (Z.mul m pz) (mul g h)))) in
      let q, r = Modp.divmod p (Modp.mul p s e) (Modp.of_z p h) in
      let dg = Modp.add p (Modp.mul p t e) (Modp.mul p q (Modp.of_z p g)) in
      let up d = Array.map (fun c -> Z.mul m (Z.of_int c)) d in
      let m' = Z.mul m pz in
      go m' (red m' (add g (up dg))) (red m' (add h (up r)))
    end in
  go pz (red pz g) (red pz h)

Lifting a factorization modulo 3.

open Zpoly
 
(* Lifting the factorization of f = (x^2 - 2)(x^3 + 5x + 11) modulo 3. *)

Running it; coefficients are shown in the symmetric range.

f = x^5 + 3x^3 + 11x^2 - 10x - 22
mod 3: (x^2 + 1) · (x^3 + 2x + 2)
mod 3            x^3 - x - 1                        x^2 + 1
mod 9            x^3 - 4x + 2                       x^2 - 2
mod 81           x^3 + 5x + 11                      x^2 - 2
mod 6561         x^3 + 5x + 11                      x^2 - 2
mod 43046721     x^3 + 5x + 11                      x^2 - 2

Modulo 3 the factors are not the true ones, but they are their images. Modulo 81 the symmetric representatives already equal the integer factors, and further steps leave them unchanged.

§ 03

Linear and quadratic lifting

Linear lifting goes from $p^{k}$ to $p^{k + 1}$ with $s$ and $t$ fixed modulo $p$ , so all its arithmetic on corrections is modulo the small prime; quadratic lifting goes from $p^{k}$ to $p^{2 k}$ but must also lift $s$ and $t$ and works with large numbers. Quadratic lifting needs logarithmically many steps.

Steps needed to lift a two-factor split modulo 101 past twice the coefficient bound of the product.

open Zpoly
 
(* Steps needed to lift a two-factor split of a product of two random
   polynomials past twice its coefficient bound: linear against quadratic. *)
let random n = normalize (Array.init (n + 1) (fun i -> if i = n then Z.one else Z.of_int (Random.int 19 - 9)))

Running it.

degree  target bits  linear steps  quadratic steps
    20           27             4                3
    40           48             7                3
    80           88            13                4
   160          169            25                5
   320          329            49                6

Since the work of a quadratic step grows with the size of the modulus, the total cost of quadratic lifting is dominated by its last step and is comparable to a few multiplications at full precision, which is asymptotically better than linear lifting^[2]. Implementations often combine the two.

§ 04

Factoring over the integers

Zassenhaus's algorithm factors a square-free primitive polynomial $f \in Z [x]$ as follows^[1]^[4]:

1. Choose a prime $p$ that does not divide the leading coefficient and for which $f mod p$ is square-free. 2. Factor $f mod p$ into irreducibles, by Berlekamp's algorithm^[7] or Cantor–Zassenhaus. 3. Lift the factorization to modulus $p^{k} > 2 B$ , where $B$ bounds the coefficients of any factor of $f$ times its leading coefficient. 4. Recombine: for subsets of the lifted factors of increasing size, multiply them together with the leading coefficient, take symmetric representatives, and test whether the result divides $f$ .

The whole algorithm, with lifting of several factors by splitting off one at a time.

(* Lift all monic factors u1..ur of f (mod p) at once, splitting off one
   factor at a time. *)
let lift_all p target f us =
  let lcf = Z.erem (lc f) target in
  let rec go f = function
    | [] | [ _ ] -> [ ]
    | u :: rest ->
      let cof = List.fold_left (fun acc v -> red (Z.of_int p) (mul acc (Array.map Z.of_int v))) (const lcf) rest in
      let m, g, h = lift_two p target f cof (Array.map Z.of_int u) in
      let g = red target g and h = red target h in
      ignore m;
      h :: (if List.length rest = 1 then
              (* the last cofactor, made monic *)
              let inv = Z.invert (lc g) target in [ red target (scale g inv) ]
            else go g rest) in
  go (red target f) us
 
(* Zassenhaus: factor a square-free primitive polynomial over Z. *)
let trials = ref 0
let factor_squarefree f =
  let rec pick p = if Modp.is_prime p && not (Z.equal (Z.erem (lc f) (Z.of_int p)) Z.zero)
                      && Modp.deg (Modp.gcd p (Modp.of_z p f) (Modp.deriv p (Modp.of_z p f))) = 0
    then p else pick (p + 2) in
  let p = pick 3 in
  let us = Cz.factor p (Modp.of_z p f) in
  let bound = Z.mul (Z.abs (lc f)) (Z.mul (Z.bin (Z.of_int (deg f)) (deg f / 2)) (norm2_ceil f)) in
  let rec target m = if Z.gt m (Z.shift_left bound 1) then m else target (Z.mul m (Z.of_int p)) in
  let m = target (Z.of_int p) in
  let lifted = if List.length us = 1 then [ Array.map Z.of_int (List.hd us) ] else lift_all p m f us in
  (* try subsets of increasing size *)
  let rec subsets k l = if k = 0 then [ [] ] else match l with
      | [] -> [] | x :: rest -> List.map (fun s -> x :: s) (subsets (k - 1) rest) @ subsets k rest in
  let rec recombine f us k acc =
    if 2 * k > List.length us then (if deg f > 0 then f :: acc else acc)
    else
      match List.find_opt (fun s ->
          incr trials;
          let g = primitive (symmetric m (red m (List.fold_left (mulm m) (const (lc f)) s))) in
          divides g f) (subsets k us) with
      | Some s ->
        let g = primitive (symmetric m (red m (List.fold_left (mulm m) (const (lc f)) s))) in
        recombine (fst (divmod_exact f g)) (List.filter (fun u -> not (List.memq u s)) us) k (g :: acc)
      | None -> recombine f us (k + 1) acc
  in
  (p, List.length us, recombine f lifted 1 [])

The factorizer above agrees with SymPy on more than a hundred random products tested. Each true factor over the integers is the product of some subset of the modular factors, and the recombination step can be exponential in their number; the Swinnerton-Dyer polynomials are the extreme case, irreducible over the integers but splitting into many small factors modulo every prime. Van Hoeij's algorithm of 2002 replaced subset enumeration with lattice reduction and made recombination polynomial in practice^[9].

§ 05

Multivariate lifting

The same idea works with a polynomial ideal in place of the prime: to factor $f (x, y)$ , substitute a value for $y$ , factor the univariate image, and lift the factors modulo powers of $y - a$ . Musser and Wang and Rothschild developed this for multivariate factorization over the integers in 1975^[6]^[5], and it is also the basis of the EEZ-GCD algorithm for multivariate gcds.

§ 06

History

Hensel introduced the $p$ -adic numbers around 1897 and presented the lemma that bears his name in his book of 1908^[3]. Zassenhaus proposed in 1969 to factor integer polynomials by lifting a modular factorization^[1], combining it with Berlekamp's algorithm for factoring over finite fields^[7]; the combination, the Berlekamp–Zassenhaus algorithm, became the standard method. Musser, Wang and Rothschild extended lifting to several variables in 1975^[6]^[5], and van Hoeij's lattice-based recombination of 2002 removed the exponential worst case from practice^[9].