Rational functions are functions of the form f(x)=P(x)Q(x)f(x) = \frac{P(x)}{Q(x)}, where P(x)P(x) and Q(x)Q(x) are polynomials. Their sums often yield curious results and are very commonly used in mathematics. Today, I will be explaining what was essentially a product of my curiosity turning into a very elegant solution for evaluating them.


We know that for an infinite sum of a rational function to be convergent, degP(x)degQ(x)<1\deg P(x) - \deg Q(x) \lt -1 must hold. This is trivially proven using a combination of the p-test and the comparison test. For example, the sum

n=11n2\sum_{n=1}^{\infty} \frac{1}{n^2}

is clearly convergent, since deg1degn2=02=2<1\deg 1 - \deg n^2 = 0 - 2 = -2 \lt -1. On the other hand,

n=11n\sum_{n=1}^{\infty} \frac{1}{n}

is divergent, since deg1degn=01=11\deg 1 - \deg n = 0 - 1 = -1 \nless -1.

Consider the sum

k=11k2a\sum_{k=1}^\infty \frac{1}{k^2 - a}

where aa is a constant. The following equalities hold:

k=11k2a=1aπcot(aπ)2a=12a(ψ(1+a)ψ(1a))\sum_{k=1}^\infty \frac{1}{k^2 - a} = \frac{1-\sqrt{a}\pi\cot(\sqrt{a}\pi)}{2a} = \frac{1}{2\sqrt{a}}\left(\psi\left(1+\sqrt{a}\right)-\psi\left(1-\sqrt{a}\right)\right)

To attempt at understanding why is this happening, recall the following identity due to Euler (1770) for zCZz \in \mathbb{C} \setminus \mathbb{Z}:

πcot(πz)=z1+v=12zz2v2 \pi \cot(\pi z) = z^{-1} + \sum_{v=1}^\infty \frac{2z}{z^2 - v^2}

A concise proof using Mittag-Leffler’s theorem follows. Given bnb_n as the residues:

f(z)=f(0)+n=1bn(1zzn+1zn)=f(0)+n=1zbnzn(znz) f(z) = f(0) + \sum_{n=1}^\infty b_n \left(\frac{1}{z-z_n} + \frac{1}{z_n}\right) = f(0) + \sum_{n=1}^\infty \frac{zb_n}{z_n(z_n-z)}

Using the contour integral where CNC_N is a circle enclosing first NN poles of ff:

IN=CNf(ω)dωω(ωz) I_N = \oint_{C_N} \frac{f(\omega)d\omega}{\omega(\omega - z)}

Using the residue theorem we find:

IN=2πif(0)f(z)+2πif(z)z+n=1N2πibnzn(znz) I_N = -2\pi i \frac{f(0)}{f(z)} + 2\pi i \frac{f(z)}{z} + \sum_{n=1}^N \frac{2\pi i b_n}{z_n(z_n - z)}

Consider cotzz1\cot z - z^{-1} to remove a singularity. bnb_n is found using the residue theorem as follows:

bn=limznπ(znπ)cotz=limznπ(znπ)zcoszsinzzsinz=1 b_n = \lim_{z \to n\pi} (z-n\pi)\cot z = \lim_{z \to n\pi} (z-n\pi) \frac{z \cos z - \sin z}{z \sin z} = 1


cotzz1=n=1N1znπ+1nπ+1z+nπ1nπ \cot z - z^{-1} = \sum_{n=1}^N \frac{1}{z-n\pi} + \frac{1}{n\pi} + \frac{1}{z+n\pi} - \frac{1}{n\pi}

Substitute and rearrange:

πcot(πz)=z1+n=1N(1zn+1z+n)=z1+n=1N2zz2n2 \pi \cot(\pi z) = z^{-1} + \sum_{n=1}^N \left(\frac{1}{z-n} + \frac{1}{z+n}\right) = z^{-1} + \sum_{n=1}^N \frac{2z}{z^2-n^2}

Knowing where the first part of the equality comes from, the next logical step is reasoning about the second equality:

1aπcot(aπ)2a=12a(ψ(1+a)ψ(1a)) \frac{1-\sqrt{a}\pi\cot(\sqrt{a}\pi)}{2a} = \frac{1}{2\sqrt{a}}\left(\psi\left(1+\sqrt{a}\right)-\psi\left(1-\sqrt{a}\right)\right)

Use the recurrence relation of the digamma function once:

1aπcot(aπ)2a=12a(1a+ψ(a)ψ(1a)) \frac{1-\sqrt{a}\pi\cot(\sqrt{a}\pi)}{2a} = \frac{1}{2\sqrt{a}}\left(\frac{1}{\sqrt{a}}+\psi\left(\sqrt{a}\right)-\psi\left(1-\sqrt{a}\right)\right)

Now it is possible to apply the digamma reflection formula given as follows:

ψ(1x)ψ(x)=πcotπx \psi (1-x)-\psi (x)=\pi \cot \pi x

The problem reduces to:

12a(1aπcotπa)=12a(1aπcotπaa)=1aπcotπa2a \frac{1}{2\sqrt{a}}\left(\frac{1}{\sqrt{a}}-\pi\cot\pi\sqrt{a}\right) = \frac{1}{2\sqrt{a}}\left(\frac{1-\sqrt{a}\pi\cot\pi\sqrt{a}}{\sqrt{a}}\right) = \frac{1-\sqrt{a}\pi\cot\pi\sqrt{a}}{2a}

Which concludes the proof.


I state the following conjecture:

Every infinite sum of a rational function, the denominator of which does not contain roots of multiplicity \ge 2, can be meaningfully expressed as a finite sum of the terms in form f(ψg)f \cdot (\psi\circ g) where ff and gg are functions of xx.

The following part of the post focuses on sheding light on this theorem from multiple angles and using it to derive a method for efficiently evaluating infinite sums of rational functions. While it is possible to special-case the general result for certain functions (cubic and quadratic), it will not be done in the example code for the sake of simplicity.

The algorithm for convergent series is as follows. Denote the roots of the polynomial Q(x+1)Q(x+1) as R0RmR_0 \dots R_m. The following equality holds if P(x)=1P(x)=1:

n=11Q(n)=ωRψ(ω)Q(ω+1) \sum_{n=1}^\infty \frac{1}{Q(n)} = -\sum_{\omega \in R} \frac{\psi(-\omega)}{Q’(\omega + 1)}

If P(x)1P(x) \ne 1 (i.e. the rational function is not an inverse of some polynomial), the formula is as follows:

n=1P(n)Q(n)=ωRP(ω+1)ψ(ω)Q(ω+1) \sum_{n=1}^\infty \frac{P(n)}{Q(n)} = -\sum_{\omega \in R} \frac{P(\omega + 1) \psi(-\omega)}{Q’(\omega + 1)}


Consider the following sum:

k=1k2+4k23k44k3+8k10 \sum_{k=1}^\infty \frac{k^2+4k-2}{3k^4-4k^3+8k-10}

Compute the derivative of Q(x+1)Q’(x+1):

Q(x+1)=12x3+24x2+12x+8 Q’(x + 1) = 12x^3 + 24x^2 + 12x + 8

The approximate roots of Q(x+1)Q(x + 1) are found using numerical methods as follows:

  • R0=0.29067622631.176486157iR_0 = -0.2906762263 - 1.176486157i
  • R1=0.2906762263+1.176486157iR_1 = -0.2906762263 + 1.176486157i
  • R2=0.2870228254R_2 = 0.2870228254
  • R3=2.372337039R_3 = -2.372337039

Recall the formula:

n=1k2+4k23k44k3+8k10=ωRP(ω+1)ψ(ω)Q(ω+1) \sum_{n=1}^\infty \frac{k^2+4k-2}{3k^4-4k^3+8k-10} = -\sum_{\omega \in R} \frac{P(\omega + 1) \psi(-\omega)}{Q’(\omega + 1)}

Computing the partial sums yields the sum 0.323477+0.392877i0.3234770.392877i+0.806386+0.0784787=0.2379107-0.323477 +0.392877i -0.323477 - 0.392877 i + 0.806386 + 0.0784787 = 0.2379107, and this value negated is the approximate value of the sum when computed by iteration until convergence.


A KamilaLisp implementation follows. Define a few helper methods: one to substitute x+1x+1 into some polynomial P(x)P(x), compute P(x)P’(x) and evaluate a polynomial at a given point.

; Utility functions.
(def polyevl
  [(inner-product cmplx64:+ cmplx64:*) #0 ^cmplx64:**&[#0 range@tally]])

(def polyderv cdr@[* #0 range@tally])

(defun poly+1 (c)
    (def bin-mat (:[:^binomial #0 \range $(+ 1)]@range@tally c))
    (def coeff-mat (:$(take (tally c)) (* c bin-mat)))
    (:$(foldl1 +)@transpose coeff-mat)))

Further, the direct implementation of the algorithm is given as follows:

(defun rsum (p q)
    (def q+1 (poly+1 q))
    (def p+1 $(polyevl (poly+1 p)))
    (def q+1p $(polyevl (polyderv q+1)))
    (def R (cmplx64:solve q+1))
    (def V [/ [* p+1 cmplx64:digamma@-] q+1p])
    (cmplx64:neg (foldl + 0 (:V R)))))