Matrix Exponentiation for Linear Recurrences

Linear recurrences are secretly linear maps in disguise, and once you see the disguise you get an $O(\log n)$ algorithm for the $n$-th term almost for free. This is one of those beautiful cases where switching representations — from "next term = combination of previous terms" to "next state = matrix times current state" — turns a naive $O(n)$ walk into repeated squaring. We'll build the trick from scratch on Fibonacci, generalize it, count operations honestly, and finally connect the whole thing to Cayley–Hamilton and Binet.

The idea, on Fibonacci

Fibonacci says $F_0 = 0$, $F_1 = 1$, and each subsequent term is the sum of its two predecessors. The naive definition is a recipe for computing $F_n$ by taking $n$ additions. That's fine when $n$ is a hundred; it is dreadful when $n$ is $10^{18}$.

The trick starts with a shift of perspective. Instead of tracking a single number $F_n$, we track a state vector — the current term together with the previous one:

Here $v_n$ is a column vector in $\mathbb{R}^2$; its top entry is the "next" Fibonacci number and its bottom entry is the "current" one. The recurrence $F_{n+2} = F_{n+1} + F_n$ says: the new top entry is the sum of both old entries, and the new bottom entry is the old top entry. That is a linear map, and any linear map on $\mathbb{R}^2$ is a $2 \times 2$ matrix. Let

Read: the top row 1 1 says "new top = 1·(old top) + 1·(old bottom)", i.e. $F_{n+2} = F_{n+1} + F_n$. The bottom row 1 0 says "new bottom = 1·(old top) + 0·(old bottom)", i.e. the new bottom is just the old top. So $v_{n+1} = M v_n$, and iterating,

In words: to find the $n$-th Fibonacci number, raise a fixed $2 \times 2$ matrix to the $n$-th power and read off an entry.

Why this is a win: fast exponentiation

Naively, computing $M^n$ takes $n-1$ matrix multiplications, and each matrix multiplication is a fixed amount of work (eight multiplies and four adds for a $2\times 2$). We saved nothing.

The win comes from repeated squaring. Any positive integer $n$ has a binary expansion, so we can write, say, $n=13 = 8+4+1 = 1101_2$, and use

The powers $M^2, M^4, M^8$ are obtained by successive squaring: $M^2 = M \cdot M$, then $M^4 = M^2 \cdot M^2$, then $M^8 = M^4 \cdot M^4$. Three squarings plus two multiplications: five matrix products to reach $M^{13}$, instead of twelve.

In general, if $n$ has $k$ bits ($k = \lceil \log_2(n+1) \rceil$) and $h$ of them are set, this scheme costs $(k-1) + (h-1)$ matrix multiplications. The worst case is $2k - 2$, which is $O(\log n)$. Here "cost" is counted in matrix multiplications, not in scalar operations — we'll come back to what that means in bits when the numbers get huge.

Small example, hand-cranked

Let's confirm. With $M = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$:

Applying $M^4$ to $v_0 = (1, 0)^\top$ gives $(5, 3)^\top$, which is $(F_5, F_4)$. Correct. Squaring again, $M^8 = \begin{pmatrix} 34 & 21 \\ 21 & 13 \end{pmatrix}$, and $M^8 v_0 = (34, 21)^\top = (F_9, F_8)$. Notice the entries of $M^n$ are themselves Fibonacci numbers:

That's not a coincidence — it drops out of induction on the recurrence. It also gives you the beautiful Cassini identity $F_{n+1} F_{n-1} - F_n^2 = (-1)^n$ for free, by taking the determinant of both sides: $\det M = -1$, so $\det(M^n) = (-1)^n$.

The general trick

Any linear recurrence of order $d$ over a ring $R$,

with constants $c_1, \dots, c_d \in R$, becomes a matrix–vector recursion. Let

be the $d$-dimensional state (the last $d$ terms). Define the companion matrix

The top row encodes the recurrence itself, and each subrow just shifts a coordinate down. Then $s_{n+1} = C s_n$, so $s_n = C^n s_0$. Same trick, wider matrix.

The algorithm, in Python

Here is a small, honest implementation. It computes the $n$-th term of the recurrence $x_{n+2} = a x_{n+1} + b x_n$ modulo some prime p, which is the standard competitive-programming setting.

def matmul(A, B, p):
    n = len(A)
    return [[sum(A[i][k] * B[k][j] for k in range(n)) % p
             for j in range(n)] for i in range(n)]

def matpow(M, k, p):
    # returns M^k mod p, assuming k >= 0
    n = len(M)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]  # identity
    base = [row[:] for row in M]
    while k > 0:
        if k & 1:
            result = matmul(result, base, p)
        base = matmul(base, base, p)
        k >>= 1
    return result

def linrec(a, b, x0, x1, n, p):
    # x_{k+2} = a * x_{k+1} + b * x_k, returns x_n mod p
    if n == 0: return x0 % p
    if n == 1: return x1 % p
    M = [[a % p, b % p], [1, 0]]
    Mn = matpow(M, n - 1, p)
    return (Mn[0][0] * x1 + Mn[0][1] * x0) % p

Two-line summary: matpow is exponentiation-by-squaring on matrices; linrec plugs the recurrence coefficients into the companion, raises it to the right power, and reads off the answer. This matpow is the standard, slightly-suboptimal variant: it multiplies result by the identity on the lowest set bit and squares base one last time after the top bit is consumed, so it does $k$ squarings and $h$ multiplications rather than the $(k-1)+(h-1)$ optimum — a couple of wasted products, harmless in practice. Under p = 10**9 + 7, linrec(1, 1, 0, 1, 10**18, p) returns in about a millisecond.

What is the cost, really?

I owe an honest accounting. For a companion matrix of size $d$:

• Each matrix multiplication costs $O(d^3)$ scalar operations with the schoolbook algorithm, or $O(d^\omega)$ with fast matrix multiplication (currently $\omega < 2.372$). We do $O(\log n)$ of them, so $O(d^3 \log n)$ scalar operations.
• If we work modulo a fixed-size integer, each scalar op is $O(1)$ and we're done: $O(d^3 \log n)$ bit operations.
• If we work over $\mathbb{Z}$, entries of $C^n$ grow like $\lambda^n$ where $\lambda$ is the dominant eigenvalue. Their bit length is $\Theta(n \log \lambda)$, so a "scalar multiplication" is really a big-integer multiplication and the true bit cost is closer to $O(d^3 \mu(n \log \lambda))$, where $\mu(b)$ is the cost of multiplying $b$-bit integers — the entry bit-lengths roughly double each squaring, so the total is a geometric series dominated by the final squaring. Fast exponentiation still beats naive iteration for exact integer $F_n$, but the win is more modest than the $\log n$ suggests — the last few squarings dominate. This is where fast-doubling identities (like $F_{2k} = F_k(2F_{k+1} - F_k)$) can shave constants.

There is a stronger complexity claim available. By the Cayley–Hamilton theorem, $C$ satisfies its own characteristic polynomial $\chi(t) = t^d - c_1 t^{d-1} - \cdots - c_d$. That means the vector $s_n = C^n s_0$ is a linear combination of $s_0, C s_0, \dots, C^{d-1} s_0$ with coefficients depending on $n$, and those coefficients are exactly $t^n \bmod \chi(t)$. So instead of doing matrix arithmetic in $d^2$ entries, we can do polynomial arithmetic modulo $\chi$ — $d$ coefficients. This is Kitamasa's algorithm (Fiduccia 1985 in its full generality), and it runs in $O(d^2 \log n)$ modular operations, or $O(d \log d \log n)$ using FFT-based polynomial multiplication. For large $d$ it is a real improvement; for $d = 2$ (Fibonacci) it just recovers the fast-doubling identities.

Where does $\phi$ come from?

Diagonalize $M$. Its characteristic polynomial is

with roots $\phi = (1+\sqrt{5})/2$ and $\psi = (1-\sqrt{5})/2$. Because the eigenvalues are distinct, $M = P D P^{-1}$ with $D = \mathrm{diag}(\phi, \psi)$, and therefore $M^n = P D^n P^{-1}$. Chasing the algebra through gives Binet's formula, though de Moivre had it a century earlier:

In words: the $n$-th Fibonacci number is a difference of two pure exponentials. Since $|\psi| < 1$, the second term is tiny for large $n$, and $F_n \approx \phi^n / \sqrt{5}$; in fact $F_n$ is the nearest integer to $\phi^n / \sqrt{5}$ for all $n \geq 0$. The load-bearing step here is diagonalizability, which fails precisely when $\chi$ has repeated roots — in that case you get a Jordan block and the closed form picks up polynomial-in-$n$ factors like $n \lambda^n$.

What is proved, sketched, and asserted

Proved above: the reformulation $v_n = M^n v_0$ (routine linear algebra); the closed form of $M^n$ for Fibonacci (induction); the operation count of fast exponentiation ($O(\log n)$ multiplications from binary expansion); Binet's formula from diagonalization.

Sketched: the Cayley–Hamilton reduction to polynomial-mod-$\chi$ arithmetic. The theorem itself — every square matrix satisfies its own characteristic polynomial — has a slick proof over $\mathbb{C}$ by density of diagonalizable matrices, and a coordinate-free proof over any commutative ring using the adjugate identity $\chi(C) = 0$; I did neither.

Asserted, not proved: that fast matrix multiplication achieves $O(d^\omega)$ with $\omega < 2.372$. The bound $\omega < 2.3729$ was already crossed before 2021 — by Vassilevska Williams 2012 ($\approx 2.37287$) and Le Gall 2014 ($\approx 2.37286$); Alman–Williams 2021 refined it to $\omega < 2.372860$, and later work (Duan–Wu–Zhou 2022; Williams–Xu–Xu–Zhou 2024) pushed the bound below 2.372. Whether $\omega = 2$ is a famous open problem — the current world record is a hair below 2.372, and the "galactic" nature of the constants means schoolbook $O(d^3)$ is what you actually run.

Where the trick genuinely fails: if the coefficients of the recurrence themselves depend on $n$ (non-constant coefficients), $M$ is no longer fixed and squaring doesn't apply — you're back to $O(n)$. That includes the Catalan-like recurrences and anything driven by a variable-length combinatorial sum.

The moral: whenever the state of an integer sequence can be captured in a bounded-size vector and the update rule is linear in that vector, you can leap ahead exponentially. The disguise is small, but it turns "walk $n$ steps" into "square $\log_2 n$ times", and that is essentially always worth the change of variables.