Rational functions are functions of the form $f(x) = \frac{P(x)}{Q(x)}$, where $P(x)$ and $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.

Warmup⌗

We know that for an infinite sum of a rational function to be convergent, $\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

$$\sum_{n=1}^{\infty} \frac{1}{n^2}$$

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

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

is divergent, since $\deg 1 - \deg n = 0 - 1 = -1 \nless -1$.

Consider the sum

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

where $a$ is a constant. The following equalities hold:

$$\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 $z \in \mathbb{C} \setminus \mathbb{Z}$:

$$\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 $b_n$ as the residues:

$$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 $C_N$ is a circle enclosing first $N$ poles of $f$:

$$I_N = \oint_{C_N} \frac{f(\omega)d\omega}{\omega(\omega - z)}$$

Using the residue theorem we find:

$$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 $\cot z - z^{-1}$ to remove a singularity. $b_n$ is found using the residue theorem as follows:

$$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$$

Hence:

$$\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:

$$\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:

$$\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:

$$\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:

$$\psi (1-x)-\psi (x)=\pi \cot \pi x$$

The problem reduces to:

$$\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.

Conjecture⌗

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 \cdot (\psi\circ g)$ where $f$ and $g$ are functions of $x$.

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)$ as $R_0 \dots R_m$. The following equality holds if $P(x)=1$:

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

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

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

Demonstration⌗

Consider the following sum:

$$\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) = 12x^3 + 24x^2 + 12x + 8$$

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

• $R_0 = -0.2906762263 - 1.176486157i$
• $R_1 = -0.2906762263 + 1.176486157i$
• $R_2 = 0.2870228254$
• $R_3 = -2.372337039$

Recall the formula:

$$\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.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.

Implementation⌗

A KamilaLisp implementation follows. Define a few helper methods: one to substitute $x+1$ into some polynomial $P(x)$, compute $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)
  (let-seq
    (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)
  (let-seq
    (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)))))