Square-free decomposition

The square-free decomposition of a polynomial writes it as a product of powers of polynomials that are square-free and pairwise coprime, $f = f_{1} f_{2}^{2} f_{3}^{3} \dots$ , so that $f_{i}$ is the product of the irreducible factors that occur exactly $i$ times. Computing it does not require factoring: the greatest common divisor of $f$ and its derivative contains every repeated factor with its multiplicity reduced by one, so a sequence of gcds and exact divisions separates the multiplicities^[3]. Yun's algorithm does this at a cost comparable to a single gcd^[1]. The decomposition is the first step of polynomial factorization, of partial fraction decomposition and of symbolic integration of rational functions. Over fields of characteristic $p$ , a $p$ -th power has zero derivative and the algorithm must also take $p$ -th roots^[5].

§ 01

Definition

A polynomial is square-free if no irreducible polynomial divides it twice. Over a field, every nonzero polynomial $f$ has a decomposition

f = c i \geq 1 \prod f_{i}^{i}, f_{i} monic and square-free, g cd (f_{i}, f_{j}) = 1 (i \neq = j)

with $c$ the leading coefficient, and it is unique. If $f = c \prod_{k} g_{k}^{e_{k}}$ is the factorization into distinct monic irreducibles, then $f_{i} = \prod_{e_{k} = i} g_{k}$ . Some $f_{i}$ may be 1. The square-free part of $f$ is $f_{1} f_{2} f_{3} \dots$ , the product of the distinct irreducible factors.

§ 02

Derivatives and repeated factors

The method rests on the derivative. If $f = g^{e} h$ then $f^{'} = e g^{e - 1} g^{'} h + g^{e} h^{'}$ , so $g^{e - 1}$ divides $f^{'}$ . In characteristic zero $g$ does not divide $e g^{'} h$ when $g$ and $h$ are coprime, so exactly $g^{e - 1}$ divides the gcd^[3]:

g cd (f, f^{'}) = i \geq 2 \prod f_{i}^{i - 1} \frac{f}{g cd ( f , f ^{'} )} = f_{1} f_{2} f_{3} \dots

In particular $f$ is square-free exactly when $g cd (f, f^{'}) = 1$ , and dividing by the gcd gives the square-free part. Repeating the idea separates the multiplicities.

Musser's algorithm

Let $c = g cd (f, f^{'})$ and $w = f / c$ , the square-free part. Then $g cd (w, c)$ is the product of the factors of multiplicity at least 2, and $w / g cd (w, c)$ those of multiplicity exactly 1, which is $f_{1}$ . Replacing $w$ by $g cd (w, c)$ and $c$ by $c / g cd (w, c)$ shifts every multiplicity down by one, and the next round yields $f_{2}$ ^[2].

Yun's algorithm

Yun rearranged the computation so that every gcd after the first directly produces one $f_{i}$ ^[1]:

a_{0} f_{i} = g cd (f, f^{'}), = g cd (b_{i}, d_{i}), b_{1} b_{i + 1} = f / a_{0}, = b_{i} / f_{i}, c_{1} c_{i + 1} = f^{'} / a_{0}, = d_{i} / f_{i}, d_{1} d_{i + 1} = c_{1} - b_{1}^{'} = c_{i + 1} - b_{i + 1}^{'}

stopping when $b_{i} = 1$ . The invariant is that $b_{i} = f_{i} f_{i + 1} \dots$ and $d_{i} = f_{i} \cdot (something coprime to b_{i})$ , so the gcd is exactly $f_{i}$ . The degrees of $b_{i}$ and $d_{i}$ shrink as factors are removed, and Yun showed that the total cost is of the same order as computing $g cd (f, f^{'})$ once^[1]^[3].

§ 03

Implementation

Polynomials over the prime field F_p.

(* 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

Musser's and Yun's algorithms over F_p with p = 1000003, larger than any degree used, so the characteristic-zero reasoning applies.

module F = Fpoly.Make (struct let p = 1_000_003 end)
open F
 
(* Cost measure: the total degree of the arguments of every gcd. *)
let gcd_work = ref 0
let gcd a b = gcd_work := !gcd_work + deg a + deg b; gcd a b
 
(* Musser: c = gcd(f, f') holds every repeated factor once less; each
   round takes the part of w that divides c, and what does not is the
   factor of the current multiplicity. *)
let musser f =
  let c = ref (gcd f (deriv f)) and w = ref (div f (gcd f (deriv f))) and i = ref 1 and out = ref [] in
  gcd_work := !gcd_work - (deg f + deg (deriv f));
  while deg !w > 0 do
    let y = gcd !w !c in
    let z = div !w y in
    if deg z > 0 then out := (z, !i) :: !out;
    w := y; c := div !c y; incr i
  done;
  List.rev !out
 
(* Yun: after one gcd with f', each round's gcd has the next factor as
   its answer, and the arguments shrink as fast as the factors do. *)
let yun f =
  let a0 = gcd f (deriv f) in
  let b = ref (div f a0) and c = ref (div (deriv f) a0) in
  let d = ref (sub !c (deriv !b)) and i = ref 1 and out = ref [] in
  while deg !b > 0 do
    let a = gcd !b !d in
    if deg a > 0 then out := (a, !i) :: !out;
    b := div !b a;
    c := div !d a;
    d := sub !c (deriv !b);
    incr i
  done;
  List.rev !out
 
let expand fs = List.fold_left (fun acc (g, i) -> mul acc (pow g i)) one fs
let show fs = String.concat " * " (List.map (fun (g, i) -> Printf.sprintf "(%s)^%d" (to_string g) i) fs)

Running it.

f has degree 17
Yun:     (x + 1)^1 * (x + 2)^2 * (x^3 - 3x^2 + x - 3)^3 * (x + 5)^5
Musser:  (x + 1)^1 * (x + 2)^2 * (x^3 - 3x^2 + x - 3)^3 * (x + 5)^5
both multiply back to f: true true
total degree of gcd arguments: Musser 70, Yun 61

The polynomial is $(x + 1) (x + 2)^{2} (x^{2} + 1)^{3} (x - 3)^{3} (x + 5)^{5}$ . The two factors of multiplicity 3 come out together as $x^{3} - 3 x^{2} + x - 3 = (x^{2} + 1) (x - 3)$ : square-free decomposition groups factors by multiplicity and does not split them further, which is the job of factorization. Yun's version does less gcd work on the same input.

§ 04

Positive characteristic

Over a field of characteristic $p$ the derivative of $x^{p}$ is $p x^{p - 1} = 0$ . Any polynomial in $x^{p}$ has zero derivative, and the argument above fails for every factor whose multiplicity is divisible by $p$ : such a factor divides $f^{'}$ with its full multiplicity, and $g cd (f, f^{'})$ can be $f$ itself. Over a finite field every polynomial with zero derivative is a $p$ -th power, because the Frobenius map $a \mapsto a^{p}$ is a bijection on the coefficients^[8]:

f^{'} = 0 ⟹ f = k \sum a_{k} x^{k p} = (k \sum a_{k}^{1/ p} x^{k})^{p}, a^{1/ p} = a in F_{p}

The algorithm for finite fields therefore runs the loop above, which finds the factors whose multiplicity is not divisible by $p$ , and recurses on the $p$ -th root of what remains, multiplying the multiplicities it finds by $p$ ^[3]^[4]:

Yun's algorithm used unchanged over F_3, and the version that takes p-th roots.

module F = Fpoly.Make (struct let p = 3 end)
open F
 
let yun f =
  let a0 = gcd f (deriv f) in
  let b = ref (div f a0) and c = ref (div (deriv f) a0) in
  let d = ref (sub !c (deriv !b)) and i = ref 1 and out = ref [] in
  while deg !b > 0 do
    let a = gcd !b !d in
    if deg a > 0 then out := (a, !i) :: !out;
    b := div !b a; c := div !d a; d := sub !c (deriv !b); incr i
  done;
  List.rev !out
 
(* Characteristic p: the loop finds the factors whose multiplicity is not
   divisible by p; what is left over, c, has zero derivative, so it is a
   p-th power, and the algorithm recurses on its p-th root. *)
let rec sqf f =
  let c = ref (gcd f (deriv f)) and w = ref (div f (gcd f (deriv f))) and i = ref 1 and out = ref [] in
  while deg !w > 0 do
    let y = gcd !w !c in
    let z = div !w y in
    if deg z > 0 then out := (z, !i) :: !out;
    w := y; c := div !c y; incr i
  done;
  let rest = if deg !c > 0 then List.map (fun (g, k) -> (g, k * p)) (sqf (pth_root !c)) else [] in
  List.sort (fun (_, a) (_, b) -> compare a b) (List.rev !out @ rest)
 
let expand fs = List.fold_left (fun acc (g, i) -> mul acc (pow g i)) one fs
let show fs = if fs = [] then "(nothing)" else String.concat " * " (List.map (fun (g, i) -> Printf.sprintf "(%s)^%d" (to_string g) i) fs)

Running it.

(x-1)^9 = x^9 - 1
  f' = 0
  Yun, as in characteristic 0:  (nothing)   (multiplies back: false)
  with p-th roots:              (x - 1)^9   (multiplies back: true)
(x^2+1)^3 (x-1)^4 (x+1) = x^11 - x^9 - x^8 + x^6 + x^5 - x^3 - x^2 + 1
  f' = -x^10 + x^7 - x^4 + x
  Yun, as in characteristic 0:  (x^2 - 1)^1   (multiplies back: false)
  with p-th roots:              (x + 1)^1 * (x^2 + 1)^3 * (x - 1)^4   (multiplies back: true)

For $(x - 1)^{9} = x^{9} - 1$ the derivative vanishes, and the characteristic-zero algorithm returns nothing. In the second example $(x^{2} + 1)^{3}$ is invisible to the derivative, and the unmodified algorithm finds only a wrong factor. Over fields of characteristic $p$ that are not perfect, such as $F_{p} (t)$ , $p$ -th roots need not exist, and Gianni and Trager gave algorithms for that case^[5].

§ 05

Over the integers and rationals

Over $Q$ the algorithms apply directly, but the gcds suffer from coefficient growth in the intermediate remainders. Practical systems work with primitive polynomials over $Z$ , compute the gcds with subresultant sequences or modular methods, and divide exactly. A polynomial over the integers that is square-free stays square-free modulo all but finitely many primes, those dividing its discriminant, which is what lets a factorizer reduce modulo a suitable prime before Berlekamp's algorithm and Hensel lifting.

§ 06

Uses

Polynomial factorization over finite fields and over the integers begins with square-free decomposition, since the algorithms that follow, distinct-degree factorization, Berlekamp and Cantor–Zassenhaus, require square-free input^[6]. In symbolic integration, the denominator of a rational function is decomposed first, so that Hermite reduction can remove the parts of multiplicity above one without factoring, and partial fraction decomposition can be carried out on the square-free factors^[7]. Real root isolation also works on the square-free part, since a repeated root is a common root of the polynomial and its derivative.

§ 07

History

Using the gcd with the derivative to find repeated factors is classical. Musser gave the algorithm above in his 1971 thesis on polynomial factorization^[2], and Horowitz used square-free decomposition in rational function integration in the same year^[7]. Yun's 1976 paper gave the improved algorithm and the analysis of its cost^[1]. Gianni and Trager treated fields of positive characteristic that are not perfect in 1996^[5].