In high school algebra, it is demonstrated that a rational function (a quotient of two polynomial functions) can be expressed as the sum of partial fractions, which are terms of the form

A(xx0)k \frac{A}{(x-x_0)^k}

The importance of the partial fraction decomposition lies in the fact that it provides algorithms for various computations with rational functions, including the explicit computation of antiderivatives (discussed later in the blog post), Taylor series expansions and inverse Laplace transforms.

In this blog post, I assume that degPdegQ\deg P \le \deg Q, i.e. that the denominator has a lower degree than the numerator. Otherwise, one can use Euclidean division and seek partial fractions of R(x)/Q(x)R(x) / Q(x), where R(x)R(x) is the remainder of the division guaranteed to have a lower degree of QQ.

The high school method

Consider the following function:

s22s+1s(s2)2(s2+1)=As+Bs2+C(s2)2+Ds+Es2+1 \frac{s^2-2s+1}{s(s-2)^2(s^2+1)} = \frac{A}{s} + \frac{B}{s-2} + \frac{C}{(s-2)^2} + \frac{Ds + E}{s^2+1}

We clear the fractions, i.e. multiply both sides by the denominator:

s22s+1=A(s2)2(s2+1)+Bs(s2)(s2+1)+Cs(s2+1)+(Ds+E)s(s2)2 s^2-2s+1=A(s-2)^2(s^2+1) + Bs(s-2)(s^2+1) + Cs(s^2+1) + (Ds + E)s(s-2)^2

Clearly, A=0.25A = 0.25 (found by substituting 00). Substituting 22 tells us that C=0.1C=0.1. Consequently, B=0.07B=0.07. Finally, D=0.32D = 0.32 and E=0.24E = 0.24. Hence:

s22s+1s(s2)2(s2+1)=0.25s+0.07s2+0.1(s2)2+0.32s+0.24s2+1 \frac{s^2-2s+1}{s(s-2)^2(s^2+1)} = \frac{0.25}{s} + \frac{0.07}{s-2} + \frac{0.1}{(s-2)^2} + \frac{0.32s + 0.24}{s^2+1}

Clearly, this method is very tedious and difficult to automate (the algorithm needs the ability to simplify and extensively manipulate rational functions).

Heaviside method

The Heaviside method in the base case can deal only with rational functions, the denominators of which do not contain any repeated roots. To illustrate the intuition behind this method, consider the following function:

2s+1s(s1)(s+1)=As+Bs1+Cs+1 \frac{2s+1}{s(s-1)(s+1)} = \frac{A}{s} + \frac{B}{s-1} + \frac{C}{s+1}

Olivier Heaviside proposes to find the value of CC by covering the s+1s+1 term and substituting s=1s=-1:

2s+1s(s1)(s+1)1=Cs+1 \frac{2s+1}{s(s-1)\color{transparent}{(s+1)}} \Bigg|_{-1} = \frac{C}{\color{transparent}{s+1}}

Clearly, C=0.5C = 0.5. To justify this maneuver, multiply both sides by s+1s+1 in the original expression:

(2s+1)(s+1)s(s1)(s+1)=A(s+1)s+B(s+1)s1+C(s+1)s+1 \frac{(2s+1)(s+1)}{s(s-1)(s+1)} = \frac{A(s+1)}{s} + \frac{B(s+1)}{s-1} + \frac{C(s+1)}{s+1}

After a chain of simplifications the problem reduces to the following:

(2s+1)s(s1)=A(s+1)s+B(s+1)s1+C \frac{(2s+1)}{s(s-1)} = \frac{A(s+1)}{s} + \frac{B(s+1)}{s-1} + C

Setting s=1s = -1 annihilates the first two terms on the right-hand side, and the result follows:

(2s+1)s(s1)1=C \frac{(2s+1)}{s(s-1)} \Bigg|_{-1} = C

Most textbooks, courses and online resources stop here, but the method can also be extended to rational functions with repeated roots.

Domain extension

The basic idea is factoring out the repeated roots and recursing (with memoisation) as follows:

f(s)=1(s1)2(s2)=1s1(1(s1)(s2))=1s1(1s21s1)by Heaviside=1(s1)2+1(s1)(s2)=1(s1)2+1s21s1by Heaviside \begin{aligned} f(s) &= \frac{1}{(s-1)^2(s-2)}\\ &= \frac{1}{s-1} \left(\frac{1}{(s-1)(s-2)}\right)\\ &= \frac{1}{s-1} \left(\frac{1}{s-2} - \frac{1}{s-1}\right) && \text{by Heaviside}\\ &= \frac{-1}{(s-1)^2} + \frac{1}{(s-1)(s-2)}\\ &= \frac{-1}{(s-1)^2} + \frac{1}{s-2} - \frac{1}{s-1} && \text{by Heaviside} \end{aligned}

Before we allow repeated roots, we must assume that the polynomial Q(x)F[x]Q(x) \in F[x] does not have factors that are irreducible over FF. This is not an unreasonable demand, as the polynomial can often be solved numerically.

Over the complex field, the rational function is by definition factored as follows:

f(x)=i(ai1xxi+ai2(xxi)2++aiki(xxi)ki) f(x)=\sum_{i}\left({\frac {a_{i1}}{x-x_{i}}}+{\frac {a_{i2}}{(x-x_{i})^{2}}}+\cdots +{\frac {a_{ik_{i}}}{(x-x_{i})^{k_{i}}}}\right)

Let gij(x)=(xxi)j1f(x){\displaystyle g_{ij}(x)=(x-x_{i})^{j-1}f(x)}. Then, according to the uniqueness of Laurent series, aija_{ij} is the coefficient of the term (xxi)1(x − x_i)^{−1} in the Laurent expansion of gij(x)g_{ij}(x) about the point xix_i, i.e., its residue aij=Res(gij,xi){\displaystyle a_{ij}=\operatorname {Res} (g_{ij},x_{i})}.

This is given by the formula:

aij=1(kij)!limxxidkijdxkij((xxi)kiP(x)Q(x)) {\displaystyle a_{ij}={\frac {1}{(k_{i}-j)!}}\lim_{x\to x_{i}}{\frac {d^{k_{i}-j}}{dx^{k_{i}-j}}}\left((x-x_{i})^{k_{i}}\frac{P(x)}{Q(x)}\right)}

When xix_i is a root with multiplicity of one, the formula reduces to:

ai1=P(xi)Q(xi) {\displaystyle a_{i1}={\frac {P(x_{i})}{Q’(x_{i})}}}

Consistently with Lagrange interpolation.

Improvements to the formula

Unfortunately, Heaviside’s original approach requires successive differentiation. This can be avoided in the following way:

Consider the following derivation for the partial fraction’s last term of a given root (given kik_i as the multiplicity of the ii-th unique root xix_i):

aiki=limxxif(x)(xxi)ki a_{i k_{i}}=\lim_{x \to x_i} f(x) (x-x_i)^{k_i}

By multiplying both sides of the definition by (xxi)ki(x - x_i)^{k_i}, we cancel out the (xxi)(x - x_i) factor in Q(x)Q(x) and include (xxi)(x - x_i) or its higher power in each partial fraction, except the constant term aikia_{i k_i}.

Consider the following function:

Pi1(x)Qi1(x)=P(x)Q(x)aiki(xxi)ki \frac{P_{i1}(x)}{Q_{i1}(x)} = \frac{P(x)}{Q(x)} - \frac{a_{i k_i}}{(x - x_i)^{k_i}}

Due to the uniqueness of the partial fraction expansion, the highest power of (xxi)(x - x_i) in Qi1(x)Q_{i1}(x) is ki1k_i-1, further implying that (xxi)(x - x_i) must be cancelled out when simplifying Pi1(x)Qi1(x)\frac{P_{i1}(x)}{Q_{i1}(x)}. Thus, we can reduce the multiplicity of the pole at xix_i by one via a single polynomial division. By applying the technique again, we can determine aiki1a_{i k_{i-1}} readily:

aiki1=limxxiPi1(x)Qi1(x)(xxi)ki1 a_{i k_{i-1}} = \lim_{x \to x_i} \frac{P_{i1}(x)}{Q_{i1}(x)} (x-x_i)^{k_i-1}

As such, we consider:

Pij(x)Qij(x)=P(x)Q(x)m=0j1aikim(xxi)kim \frac{P_{ij}(x)}{Q_{ij}(x)} = \frac{P(x)}{Q(x)} - \sum_{m=0}^{j-1} \frac{a_{i k_i-m}}{(x - x_i)^{k_i-m}}

To determine the formula for aikija_{i k_i-j} for i[1,s]i \in [1, s] (ss being the amount of unique roots) and j[1,ki1]j \in [1, k_i-1], proceed as:

aikij=limxxiPij(x)Qij(x)(xxi)kij=limxxi(P(x)Q(x)m=0j1aikim(xxi)kim)(xxi)kij a_{i k_i-j} = \lim_{x \to x_i} \frac{P_{ij}(x)}{Q_{ij}(x)} (x-x_i)^{k_i-j} = \lim_{x \to x_i} \left(\frac{P(x)}{Q(x)} - \sum_{m=0}^{j-1} \frac{a_{i k_i-m}}{(x - x_i)^{k_i-m}}\right) (x-x_i)^{k_i-j}

An example

We want to find the partial fractions of the following function:

2x+1x(x1)(x2)2=Ax+Bx1+Cx2+D(x2)2 \frac{2x+1}{x(x-1)(x-2)^2} = \frac{A}{x} + \frac{B}{x-1} + \frac{C}{x-2} + \frac{D}{(x-2)^2}

Applying the first step to AA, BB and DD yields:

A=limxxif(x)(xxi)kiby definition=limx0f(x)xby values of x1 and k1=limx0x(2x+1)x(x1)(x2)2by definition of f=limx0(2x+1)(x1)(x2)2cancel out x=1(01)(02)2substitute into the limit=14B=limxxif(x)(xxi)ki=limx1f(x)(x1)=limx1(x1)(2x+1)x(x1)(x2)2=limx1(2x+1)x(x2)2=2+1(12)2=3D=limxxif(x)(xxi)ki=limx2f(x)(x2)2=limx2(x2)2(2x+1)x(x1)(x2)2=limx2(2x+1)x(x1)=(4+1)2=2.5 \begin{aligned} A &= \lim_{x\to x_i} f(x) (x-x_i)^{k_i} && \text{by definition} \\ &= \lim_{x\to 0} f(x) x && \text{by values of x1x_1 and k1k_1} \\ &= \lim_{x\to 0} \frac{x(2x+1)}{x(x-1)(x-2)^2} && \text{by definition of ff} \\ &= \lim_{x\to 0} \frac{(2x+1)}{(x-1)(x-2)^2} && \text{cancel out xx} \\ &= \frac{1}{(0-1)(0-2)^2} && \text{substitute into the limit} \\ &= -\frac{1}{4} \\ B &= \lim_{x\to x_i} f(x) (x-x_i)^{k_i} = \lim_{x\to 1} f(x) (x-1) \\ &= \lim_{x\to 1} \frac{(x-1)(2x+1)}{x(x-1)(x-2)^2} = \lim_{x\to 1} \frac{(2x+1)}{x(x-2)^2} \\ &= \frac{2+1}{(1-2)^2} = 3 \\ D &= \lim_{x\to x_i} f(x) (x-x_i)^{k_i} = \lim_{x\to 2} f(x) (x-2)^2 \\ &= \lim_{x\to 2} \frac{(x-2)^2(2x+1)}{x(x-1)(x-2)^2} = \lim_{x\to 2} \frac{(2x+1)}{x(x-1)} \\ &= \frac{(4+1)}{2} = 2.5 \end{aligned}

The value of C=a31C = a_{31} is calculated by considering the formula for aikija_{i k_i-j} with i=3i = 3, ki=2k_i = 2, j=1j = 1 and xi=2x_i = 2:

C=a31=limx2P31(x)Q31(x)(x2)=limx2(P(x)Q(x)a32(x2)2)(x2)=limx2(2x+1x(x1)(x2)22.5(x2)2)(x2)=limx22x+1x(x1)(x2)2.5(x2)=limx22.5x0.5x(x1)=2.75 \begin{aligned} C = a_{31} &= \lim_{x \to 2} \frac{P_{31}(x)}{Q_{31}(x)} (x-2) \\ &= \lim_{x \to 2} \left(\frac{P(x)}{Q(x)} - \frac{a_{32}}{(x - 2)^{2}}\right) (x-2) \\ &= \lim_{x \to 2} \left(\frac{2x+1}{x(x-1)(x-2)^2} - \frac{2.5}{(x - 2)^{2}}\right) (x-2) \\ &= \lim_{x \to 2} \frac{2x+1}{x(x-1)(x-2)} - \frac{2.5}{(x - 2)} \\ &= \lim_{x \to 2} \frac{-2.5x-0.5}{x(x-1)} \\ &= -2.75 \end{aligned}

Hence, the partial fraction decomposition follows as:

2x+1x(x1)(x2)2=0.25x+3x1+2.75x2+2.5(x2)2 \frac{2x+1}{x(x-1)(x-2)^2} = \frac{-0.25}{x} + \frac{3}{x-1} + \frac{-2.75}{x-2} + \frac{2.5}{(x-2)^2}

Application to antiderivatives

Partial fraction decomposition, especially when performed using this simple algorithm can be used to compute antiderivatives and infinite sums of rational functions.

Consider this motivating example:

3x4+4x3+3x2(x+6)(x7)(x3)2dx=4361104(x7)212(x3)233524(x3)11639(x+6)+3dx \int \dfrac{3x^4+ 4x^3 + 3x^2}{(x+6)(x-7)(x-3)^2} \mathrm dx = \int \frac{4361}{104 (x-7)} - \frac{21}{2(x-3)^2} - \frac{335}{24 (x-3)} - \frac{116}{39 (x+6)} + 3 \mathrm dx

Which is trivially computed using a elementary integration techniques, starting with a rule to break up the expression into parts:

f(x)+g(x)dx=f(x)dx+g(x)dx \displaystyle \int f(x)+g(x) \mathrm dx = \int f(x) \mathrm dx + \int g(x) \mathrm dx

Then, each expression following the same rule is integrated as follows:

a(xb)ndx=a(xb)1n1n+C \displaystyle \int \frac{a}{(x-b)^n} \mathrm dx = \frac{a(x-b)^{1-n}}{1-n} + C

An alternative way to approach problems like these is using the Ostrogradsky-Hermite’s method. Consider the following integral:

I=27x8+18x723x6+15x599x4+3x3+2x2+10x55x763x6+267x5349x4324x3+912x2448xdx I=\int \dfrac{27x^8+18x^7-23x^6+15x^5-99x^4+3x^3+2x^2+10x-5}{5x^7-63x^6+267x^5-349x^4-324x^3+912x^2-448x} \mathrm dx

We begin by finding Q(x)Q’(x) as follows:

Q(x)=(x4)2(35x498x39x2+100x28) Q’(x) = (x-4)^2 (35x^4-98x^3-9x^2+100x-28)

Then, it becomes apparent that

gcd(Q,Q)=Q1=(x4)2(x1) \text{gcd}(Q, Q’) = Q_1 = (x-4)^2(x-1)

We find Q2Q_2 as follows:

Q2=QQ1=5x418x315x2+28x Q_2 = \frac{Q}{Q_1} = 5 x^4 - 18 x^3 - 15 x^2 + 28 x

P1(x)P_1(x) and P2(x)P_2(x) are currently unknown, but they will have degrees one less than those of Q1(x)Q_1(x) and Q2(x)Q_2(x), respectively. Thus:

I=ax3+bx2+cx+d5x418x315x2+28x+ex3+fx2+gx+h5x418x315x2+28xdx I=\frac{ax^3+bx^2+cx+d}{5 x^4 - 18 x^3 - 15 x^2 + 28 x}+\int \frac{ex^3+fx^2+gx+h}{5 x^4 - 18 x^3 - 15 x^2 + 28 x} \mathrm dx

The merit of this method is obvious: The integral of a rational function which has roots with non-one multiplicity was reduced to an integral of a rational function with simple roots, akin to the domain extension to the Heaviside method. The coefficients must be equated now to find the values from aa to hh, which is a rather tedious and simple task left to a determined reader.

Application to infinite sums

The following identity involving the polygamma function can be used to compute infinite sums of rational functions:

n=0P(n)Q(n)=n=0k=1mak(n+bk)rk=k=1m(1)rk(rk1)!akψ(rk1)(bk) \sum_{n=0}^{\infty}\frac{P(n)}{Q(n)}=\sum_{n=0}^{\infty}\sum_{k=1}^{m}{\frac{a_{k}}{(n+b_{k})^{r_{k}}}}=\sum_{k=1}^{m}{\frac{(-1)^{r_{k}}}{(r_{k}-1)!}}a_{k}\psi ^{(r_{k}-1)}(b_{k})

The second term is a partial fraction decomposition of the rational function, which can be computed using methods described earlier. The notation ψ(a)(x)\psi ^{(a)}(x) denotes the polygamma function of the aa-th order, defined as the aa-th derivative of the digamma function, the logarithmic derivative of Γ(x)\Gamma(x). While in this scenario a>0a>0 and aNa\in \mathbb{N} is almost taken for granted, the following formula can be used to compute the polygamma function also given that z>0\Re z > 0:

ψ(m)(z)=(1)m+10tmezt1etdt=01tz11t(lnt)mdt=(1)m+1m!ζ(m+1,z) {\begin{aligned}\psi ^{(m)}(z)&=(-1)^{m+1}\int _{0}^{\infty }{\frac {t^{m}e^{-zt}}{1-e^{-t}}}\mathrm {d} t\\&=-\int _{0}^{1}{\frac {t^{z-1}}{1-t}}(\ln t)^{m}\mathrm {d} t\\&=(-1)^{m+1}m!\zeta (m+1,z)\end{aligned}}

If a=0a = 0 and z>0\Re z > 0, then:

ψ(z)=0ettezt1etdt \psi(z) = \int_0^\infty {\frac {e^{-t} } t - \frac {e^{-z t} } {1 - e^{-t} } } \mathrm dt

If zz is outside of the domain, this is easily fixed using the digamma reflection formula: ψ(1x)ψ(x)=πcotπx \psi (1-x)-\psi (x)=\pi \cot \pi x

Meanwhile the same defect with the general formula is fixed using the polygamma reflection formula: (1)mψ(m)(1z)ψ(m)(z)=πdmdzmcotπz=πm+1Pm(cosπz)sinm+1(πz) (-1)^{m}\psi ^{(m)}(1-z)-\psi ^{(m)}(z)=\pi {\frac {\mathrm {d} ^{m}}{\mathrm {d} z^{m}}}\cot {\pi z}=\pi ^{m+1}{\frac {P_{m}(\cos {\pi z})}{\sin ^{m+1}(\pi z)}}

Given that:

P0(x)=xPm+1(x)=((m+1)xPm(x)+(1x2)Pm(x)) {\begin{aligned}P_{0}(x)&=x\\P_{m+1}(x)&=-\left((m+1)xP_{m}(x)+\left(1-x^{2}\right)P’_{m}(x)\right)\end{aligned}}


The Heaviside method is a simple and elegant way to compute partial fraction decompositions of rational functions. Researching this topic has been a very interesting experience, and I hope that this blog post will be useful to someone who is also dealing with real and complex analysis in computer algebra.