Matrix Exponentiation for Linear Recurrences

This is an extension topic beyond most interviews; it pays off when a simple DP or recurrence has a tiny state but the number of steps is astronomical (n up to 10^18), which shows up in hard OAs and as a follow-up to Fibonacci-style questions.

Intuition

If one step of a process transforms a small state vector linearly, that step is a matrix M, and n steps are M^n. Fast exponentiation computes M^n with O(log n) matrix multiplications, so a recurrence that would need 10^18 iterations finishes in about 60 squarings.

Fibonacci is the classic example. The state (F(i+1), F(i)) becomes (F(i+1) + F(i), F(i+1)):

| F(i+2) |   | 1 1 |   | F(i+1) |
| F(i+1) | = | 1 0 | x | F(i)   |

so [[1, 1], [1, 0]]^n = [[F(n+1), F(n)], [F(n), F(n-1)]]. Any linear recurrence f(i) = c1*f(i-1) + ... + ck*f(i-k) works the same way with a k x k "companion" matrix: the first row holds the coefficients and the rows below shift the window by one.

Two more ideas extend the reach:

• Constants and extra terms: f(i) = f(i-1) + f(i-2) + 1 is not linear in (f(i-1), f(i-2)), but it is linear in (f(i-1), f(i-2), 1). Add the constant (or i, or a running sum) to the state.
• Walks in a graph: if A is the adjacency matrix, (A^k)[i][j] counts walks with exactly k edges from i to j. Every DP of the form "state = current node, one step per round" is a matrix power in disguise.

When to use it in interviews / OAs

• A recurrence or step-by-step DP with at most a few dozen states, but n up to 10^9-10^18.
• Counting sequences or strings of huge length whose next character depends only on the previous one (or previous few).
• Counting walks of an exact length in a small graph.
• The follow-up to "climbing stairs / Fibonacci / tiling": "now n is 10^18".

Template

def mat_mul(A: list[list[int]], B: list[list[int]], mod: int) -> list[list[int]]:
    cols = list(zip(*B))
    return [[sum(a * b for a, b in zip(row, col)) % mod for col in cols] for row in A]


def mat_pow(M: list[list[int]], n: int, mod: int) -> list[list[int]]:
    k = len(M)
    result = [[int(i == j) for j in range(k)] for i in range(k)]  # identity matrix
    while n:
        if n & 1:
            result = mat_mul(result, M, mod)
        M = mat_mul(M, M, mod)
        n >>= 1
    return result


def linear_recurrence(coeffs: list[int], initial: list[int], n: int, mod: int) -> int:
    # f(i) = coeffs[0]*f(i-1) + coeffs[1]*f(i-2) + ...; initial = [f(0), ..., f(k-1)]
    k = len(coeffs)
    if n < k:
        return initial[n] % mod
    # companion matrix: first row = coefficients, then shift the window down by one
    companion = [coeffs[:]] + [[int(c == r - 1) for c in range(k)] for r in range(1, k)]
    P = mat_pow(companion, n - k + 1, mod)
    state = initial[::-1]  # [f(k-1), f(k-2), ..., f(0)]
    return sum(P[0][j] * state[j] for j in range(k)) % mod

For Fibonacci specifically, "fast doubling" squeezes the 2 x 2 matrix power into two identities, F(2k) = F(k) * (2F(k+1) - F(k)) and F(2k+1) = F(k)^2 + F(k+1)^2:

def fib_pair(n: int, mod: int) -> tuple[int, int]:
    # returns (F(n), F(n+1))
    if n == 0:
        return 0, 1
    a, b = fib_pair(n >> 1, mod)  # a = F(k), b = F(k+1) with k = n // 2
    c = a * (2 * b - a) % mod  # F(2k)
    d = (a * a + b * b) % mod  # F(2k+1)
    return (d, (c + d) % mod) if n & 1 else (c, d)

Counting walks

def count_walks(adj: list[list[int]], k: int, mod: int) -> list[list[int]]:
    # entry [i][j] = number of walks with exactly k edges from i to j
    return mat_pow(adj, k, mod)


def knight_dialer(n: int, mod: int = 10**9 + 7) -> int:
    # phone keypad: digit -> digits a chess knight can jump to
    jumps = {0: (4, 6), 1: (6, 8), 2: (7, 9), 3: (4, 8), 4: (0, 3, 9),
             5: (), 6: (0, 1, 7), 7: (2, 6), 8: (1, 3), 9: (2, 4)}
    adj = [[int(v in jumps[u]) for v in range(10)] for u in range(10)]
    walks = count_walks(adj, n - 1, mod)  # a number with n digits is a walk with n - 1 jumps
    return sum(map(sum, walks)) % mod

The usual O(n * 10) DP is fine for n up to 10^5 or so; the matrix version answers n = 10^18 in a few hundred tiny multiplications.

Complexity

O(k^3 log n) time and O(k^2) space for a k x k matrix. Fast doubling for Fibonacci is O(log n). Because of the k^3 factor, this is practical in Python for k up to roughly 30-50; beyond that, use the plain DP if n allows it.

Common pitfalls

• Exponent off-by-one: write out the first two or three terms by hand and check that M^(n-k+1) (or whatever you derived) lands on f(n).
• Row vs column vectors: decide whether the state is a column vector (v' = M v, coefficients in the first row) or a row vector (v' = v M, coefficients in the first column) and stay consistent.
• Forgetting the modulus inside multiplication: entries grow exponentially; reduce after every product. In C++/Java make sure a sum of k products fits in 64 bits.
• Negative intermediate values: fast doubling computes 2*b - a, which is negative before the modulo in C++/Java; add mod before taking %.
• Non-linear steps: max, min or conditions on the value cannot be expressed as an ordinary matrix product. (The "min-plus" product, where multiplication becomes + and addition becomes min, handles shortest paths with exactly k edges.)

Interview variations

• Fibonacci, climbing stairs, tribonacci for huge n.
• Tilings of a 2 x n or 3 x n board for huge n: the DP over column profiles is a matrix.
• Knight dialer, count vowel strings with follow rules: walks in a small transition graph.
• Painting an n x 3 grid: transitions between the 12 valid row colorings.
• Recurrences with a constant or polynomial term: augment the state with 1, i, or a running sum.