Partial fractions

Partial fractions express a quotient of polynomials as a sum of fractions whose denominators are powers of the factors of the original denominator^[1]^[2]. Each term is simpler to integrate, to expand as a power series or to invert under a transform than the whole, which makes the decomposition a standard step in calculus, in the analysis of linear systems, and in computer algebra. Its existence follows from Bézout's identity for coprime polynomials, and computing it costs little beyond a factorization of the denominator; the level of factorization chosen, square-free, over the rationals or over the complex numbers, determines both the cost and the shape of the result^[1]^[3]. Leibniz and Johann Bernoulli used partial fractions in 1702 and 1703 to integrate rational functions^[6]^[8].

§ 01

Example

\frac{x ^{5} + 2 x + 1}{( x - 1 ) ^{2} ( x ^{2} + 1 ) ( x + 2 )} = 1 + \frac{5/18}{x - 1} + \frac{2/3}{( x - 1 ) ^{2}} + \frac{( x - 1 ) /2}{x ^{2} + 1} - \frac{7/9}{x + 2}

The polynomial part comes from dividing the numerator by the denominator; each factor of the denominator contributes one term for each power up to its multiplicity, with a numerator of lower degree than the factor. The existence of such a decomposition follows from Bézout's identity: for coprime $a$ and $b$ there are $s$ and $t$ with $s a + t b = 1$ , so $n / (ab) = n t / a + n s / b$ .

§ 02

Statement

Let $k$ be a field, $n, d \in k [x]$ with $d \neq = 0$ , and $d = \prod_{i = 1}^{r} f_{i}^{e_{i}}$ with the $f_{i}$ pairwise coprime. Then there are unique polynomials $q$ and $c_{ij}$ such that

\frac{n}{d} = q + i = 1 \sum r j = 1 \sum e_{i} \frac{c _{ij}}{f _{i}^{j}}, de g c_{ij} < de g f_{i}

The factors need not be irreducible: any pairwise coprime factorization gives a decomposition, and irreducible factors give the finest one^[2].

§ 03

Algorithm

Divide first: $n = q d + r$ with $de g r < de g d$ . Then split the proper fraction $r / d$ across one coprime factor at a time: with $d = ab$ , solve $s b + t a = r$ with $de g s < de g a$ by the extended Euclidean algorithm, so that $r / (ab) = s / a + t / b$ . Finally, for a numerator $s$ over a power $f^{e}$ , divide repeatedly by $f$ : the successive remainders are the numerators of $1/ f^{e}, 1/ f^{e - 1}, \dots$ , the digits of $s$ in base $f$ .

The decomposition of the example. Division gives the polynomial part; Bézout relations split the remainder among the coprime factors; the numerator over (x − 1)² is expanded in powers of x − 1 to give one term for each power.

Partial fractions for a given coprime factorization of the denominator.

(* Partial fractions over Q for a given factorization of the denominator
   into pairwise coprime powers d = prod f_i^e_i. *)
open Qpoly
 
(* n / (a b) with a, b coprime: find s, t with s b + t a = n, so that
   n / (a b) = s / a + t / b. *)
let split n a b = let t, s = diophantine a b n in (s, t)
 
(* r / f^e as sum over j of c_j / f^j, by repeated division by f. *)
let expand_power r f e =
  let rec go r j acc = if j = 0 then acc else
      let q, c = divmod r f in go q (j - 1) ((c, j) :: acc) in
  List.filter (fun (c, _) -> not (is_zero c)) (go r e [])
 
(* Returns the polynomial part and a list of (numerator, factor, power). *)
let apart n (factors : (t * int) list) =
  let d = List.fold_left (fun acc (f, e) -> mul acc (pow f e)) one factors in
  let q, r = divmod n d in
  let rec go r = function
    | [] -> []
    | [ (f, e) ] -> List.map (fun (c, j) -> (c, f, j)) (expand_power r f e)
    | (f, e) :: rest ->
      let rest_d = List.fold_left (fun acc (g, k) -> mul acc (pow g k)) one rest in
      let a, b = split r (pow f e) rest_d in
      List.map (fun (c, j) -> (c, f, j)) (expand_power a f e) @ go b rest in
  (q, go r factors)
 
(* Recombine, to check the result. *)
let recombine (q, terms) =
  let d = List.fold_left (fun acc (_, f, j) -> let g = pow f j in mul acc (quo g (gcd acc g))) one terms in
  let n = List.fold_left (fun acc (c, f, j) -> add acc (mul c (quo d (pow f j)))) (mul q d) terms in
  (n, d)
 
(* One term c / f^j, written with the factor and its power. *)
let show_term ?(ascending = false) (c, f, j) =
  let par s = if String.contains s ' ' || String.contains s '/' then "(" ^ s ^ ")" else s in
  let fs = "(" ^ to_string ~ascending f ^ ")" ^ (if j > 1 then "^" ^ string_of_int j else "") in
  par (to_string c) ^ " / " ^ fs

The example above, checked by recombining the terms.

open Qpoly

Running it.

polynomial part: 1
  (5/18) / (x - 1)
  (2/3) / (x - 1)^2
  ((1/2)x - 1/2) / (x^2 + 1)
  (-7/9) / (x + 2)
recombines to the original: true

The result agrees with SymPy's apart. The whole computation costs $O (n^{2})$ field operations with classical arithmetic for a denominator of degree $n$ , and Kung and Tong showed how to do it in $O (n lo g^{2} n)$ with fast arithmetic^[4]^[5].

§ 04

Levels of factorization

The irreducible factorization over $Q$ requires polynomial factoring. The square-free factorization $d = \prod d_{i}^{i}$ , in which each $d_{i}$ collects the factors of multiplicity $i$ , requires only gcds, and it is already a coprime factorization, so it gives a coarser decomposition cheaply. For integration, this coarser level is all that is needed: Hermite reduction removes the repeated factors, and the logarithmic part can be computed without factoring^[1]. Over the complex numbers every irreducible factor is linear, and the numerator of a simple pole $a$ is the residue $n (a) / d^{'} (a)$ ; Bronstein and Salvy showed how to express this full decomposition exactly, as a sum over the roots of a polynomial, without computing the roots^[3].

One denominator at both levels.

open Qpoly
 
(* The square-free factorization of the denominator needs only gcds; the
   full factorization over Q needs factoring. Partial fractions can be
   taken at either level. *)

Running it.

1 / (x^9 - 3x^8 - 6x^7 + 8x^6 + 26x^5 + 24x^4 - 15x^3 - 47x^2 - 42x - 18)
square-free partial fractions:
  (-(33/784)x^2 + 41/392) / (x^3 + x^2 - 2x - 2)
  ((33/784)x^2 - (99/784)x + 17/784) / (x^3 - 2x^2 - 2x - 3)
  (-(3/28)x^2 + (9/28)x + 1/28) / (x^3 - 2x^2 - 2x - 3)^2
over the irreducible factors:
  ((1/49)x - 1/49) / (x^2 - 2)
  (-1/16) / (x + 1)
  (-795/1722448) / (x - 3)
  (1/4732) / (x - 3)^2
  ((4581/107653)x + 176/107653) / (x^2 + x + 1)
  ((31/1183)x - 6/1183) / (x^2 + x + 1)^2

§ 05

Generating functions

A rational generating function has coefficients given by a linear recurrence, and partial fractions give their closed form: $c / (1 - a x)$ contributes $c a^{n}$ to the $n$ -th coefficient and $c / (1 - a x)^{2}$ contributes $c (n + 1) a^{n}$ ^[7].

The coefficients of 1 / ((1 − x)²(1 − 2x)) from its partial fractions, against direct series division.

open Qpoly
 
(* Coefficients of the generating function 1 / ((1 - x)^2 (1 - 2x)):
   partial fractions turn it into sums of geometric series. *)

Running it.

  -2 / (1 - x)
  -1 / (1 - x)^2
  4 / (1 - 2x)
from partial fractions: 1 4 11 26 57 120 247 502 1013 2036
by series division:     1 4 11 26 57 120 247 502 1013 2036

The terms give the closed form $4 \cdot 2^{n} - (n + 1) - 2 = 2^{n + 2} - n - 3$ .

§ 06

History

Leibniz and Johann Bernoulli introduced the decomposition in 1702 and 1703 to reduce integrals of rational functions to known forms, and argued over whether every real polynomial splits into real linear and quadratic factors, which the fundamental theorem of algebra later confirmed^[6]^[8]. Euler used it systematically. In computer algebra, Horowitz analyzed partial fraction algorithms for rational function integration in 1971^[5], Kung and Tong gave fast algorithms in 1977^[4], and Bronstein and Salvy the full decomposition over the algebraic closure in 1993^[3].