Hi everyone!

There are already dozens of blogs on linear recurrences, why not make another one? In this article, my main goal is to highlight the possible approaches to solving linear recurrence relations, their applications and implications. I will try to derive the results with different approaches independently from each other, but will also highlight similarities between them after they're derived.

Definitions

Def. 1. An order $$$d$$$ homogeneous linear recurrence with constant coefficients (or linear recurrence) is an equation of the form

$$$ F_n = \sum\limits_{k=1}^d a_k F_{n-k}. $$$

Def. 2. In the equation above, the coefficients $$$a_1, \dots, a_d$$$ are called the recurrence parameters,

Def. 3. and a sequence $$$F_0, F_1, \dots$$$ is called an order $$$d$$$ linear recurrence sequence.

---

The most common task with linear recurrences is, given initial coefficients $$$F_0, F_1, \dots, F_{d-1}$$$, to find the value of $$$F_n$$$.

---

Example 1. A famous Fibonacci sequence $$$F_n = F_{n-1} + F_{n-2}$$$ is an order 2 linear recurrence sequence.

Example 2. Let $$$F_n = n^2$$$. One can prove that $$$F_n = 3 F_{n-1} - 3 F_{n-2} + F_{n-3}$$$.

Example 3. Moreover, for $$$F_n = P(n)$$$, where $$$P(n)$$$ is a degree $$$d$$$ polynomial, it holds that

$$$ F_n = \sum\limits_{k=1}^{d+1} (-1)^{k+1}\binom{d+1}{k} F_{n-k} = 0. $$$

If this fact is not obvious to you, do not worry as it will be explained further below.

---

Finally, before proceeding to next sections, we'll need one more definition.

---

Def. 4. A polynomial $$$ A(x) = x^d - \sum\limits_{k=1}^d a_k x^{d-k} $$$

is called the characteristic polynomial of the linear recurrence defined by $$$a_1, \dots, a_d$$$.

Example 4. For Fibonacci sequence, the characteristic polynomial is $$$A(x) = x^2-x-1$$$.

Matrix approach

For linear recurrence sequences it holds that

$$$ \mathbf{F}_{n} = \begin{bmatrix}F_{n+d-1} \\ F_{n+d-2} \\ \dots \\ F_{n}\end{bmatrix} = \begin{bmatrix} a_1 & a_2 & \dots & a_{d-1} & a_d \\ 1 & 0 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} \begin{bmatrix}F_{n+d-2} \\ F_{n+d-3} \\ \dots \\ F_{n-1}\end{bmatrix} =

\left(\begin{bmatrix} a_1 & a_2 & \dots & a_{d-1} & a_d \\ 1 & 0 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} \right)^{n} \begin{bmatrix}F_{d-1} \\ F_{d-2} \\ \dots \\ F_{0}\end{bmatrix} = A^n \mathbf{F}_0. $$$

This representation allows to find $$$F_n$$$ in $$$O(d^3 \log n)$$$ with binary matrix exponentiation.

Improving the complexity of this solution is possible with the following theorem:

---

Theorem 1 (Cayley–Hamilton theorem). Every square matrix over a commutative ring satisfies its own characteristic polynomial.

---

For the matrix $$$A$$$ used above, it can be proven that its characteristic polynomial is exactly $$$A(x)$$$, thus $$$ A^d = \sum\limits_{k=1}^d a_k A^{d-k}. $$$

For the purpose of calculating $$$F_n$$$ it means that we may find $$$b_0, b_1, \dots, b_{d-1}$$$ such that

$$$ x^n \equiv b_0 + b_1 x + \dots + b_{d-1} x^{d-1} \mod{A(x)}, $$$

and this equivalence would also hold if we substitute $$$x$$$ with $$$A$$$. This also implies that

$$$ F_n = (\mathbf F_n)_d = (A^n \mathbf F_0)_d= \sum\limits_{k=0}^{d-1} b_k \left(A^k \mathbf{F}_0\right)_d = \sum\limits_{k=0}^{d-1} b_k F_k, $$$

thus allowing to compute $$$F_n$$$ in $$$O(d \log d \log n)$$$ using the binary lifting modulo $$$A(x)$$$ (see on CP-Algorithms).

---

Approach above is intuitive and good if you want just to find the $$$n$$$-th element of the recurrence. However it might be hard to analyze sequences with it (e.g. estimate the asymptotic behavior or find a recurrence for a sum of sequences).

Problem example: CodeChef — Random Number Generator.

Generating function approach

Consider the generating function $$$G(x) = F_0 + F_1 x + F_2 x^2 + \dots$$$ of the sequence. For $$$G(x)$$$ it holds that

$$$ G(x) = P(x) + \sum\limits_{k=1}^d a_k x^k G(x), $$$

where $$$P(x)$$$ is some polynomial of degree less than $$$d$$$ used to calibrate the fact that $$$F_0, \dots, F_{d-1}$$$ are pretty much arbitrary and might not adhere to the recurrence with lower terms. From this equation we see that

$$$ G(x) = \frac{P(x)}{Q(x)}, $$$

where $$$Q(x)$$$ is a polynomial obtained by reversing the coefficients of $$$A(x)$$$, that is

$$$ Q(x) = 1 - \sum\limits_{k=1}^d a_k x^k = x^d A(x^{-1}). $$$

Now, if we allow negative powers near $$$x$$$, the following equation will hold:

$$$ F_n = [x^0] x^n G(x^{-1}), $$$

where $$$[x^k] P(x)$$$ is the coefficient near $$$x^k$$$ in the series $$$P(x)$$$. We should notice that

$$$ A(x) G(x^{-1}) = x^d P(x^{-1}) $$$

has only positive coefficients, so is $$$x^k A(x) G(x^{-1})$$$ and so is all their linear combinations. In other words,

$$$ [x^0] X(x) G(x^{-1}) = [x^0] Y(x) G(x^{-1}) $$$

whenever $$$X(x)-Y(x)$$$ is divisible by $$$A(x)$$$. From this, we conclude that

$$$ F_n = [x^0] (x^n \bmod A)G(x^{-1}) = \sum\limits_{k=0}^{d-1} b_k [x^0] x^k G(x^{-1}) = \sum\limits_{k=0}^{d-1} b_k F_k, $$$

which is essentially the same result as in the matrix approach.

Closed-form solution to the recurrence

While it could be hard to analyze the asymptotic behavior of $$$F_n$$$ with matrices (you probably can, but you'd need to work with eigenvalues of $$$A$$$), it is way simpler with the generating function representation, thanks to the existence of the partial fraction decomposition:

$$$ \frac{P(x)}{Q(x)} = \sum\limits_{i=1}^m \frac{C_i(x)}{(1-\lambda_i x)^{d_i}}, $$$

where $$$\lambda_1, \dots, \lambda_m$$$ are the roots of $$$A(x)$$$ with multiplicities $$$d_1,\dots,d_m$$$, and $$$C_i(x)$$$ is a polynomial of degree less than $$$d_i$$$. In other words,

$$$ \begin{gather} Q(x) = \prod\limits_{i=1}^m (1-\lambda_i x)^{d_i},\\ A(x) = \prod\limits_{i=1}^m (x-\lambda_i)^{d_i}.\\ \end{gather} $$$

It is generally known that

$$$ \frac{1}{(1-\lambda x)^d} = \sum\limits_{k=0}^\infty \binom{k+d-1}{d-1} \lambda^k x^k, $$$

from which we may conclude that the general solution to the linear recurrence is given by

$$$ F_n = \sum\limits_{i=1}^m p_i(n)\lambda_i^n, $$$

where $$$p_i(n)$$$ is a polynomial of degree less than $$$d_i$$$, determined from $$$F_0, \dots, F_{d-1}$$$.

Umbral calculus approach