SpaceComplexity

Matrix Exponentiation: The O(log n) Trick for Any Linear Recurrence

June 5, 202615 min read
dsaalgorithmsinterview-prepdata-structures
Matrix Exponentiation: The O(log n) Trick for Any Linear Recurrence
TL;DR
  • Matrix exponentiation computes M^n in O(k³ log n) by applying binary exponentiation to matrix multiplication instead of scalar numbers.
  • The Fibonacci transition matrix [[1,1],[1,0]] raised to the nth power gives F(n) at position [0][1].
  • Any k-th order linear recurrence maps to a k×k companion matrix, achieving the same O(k³ log n) time.
  • Start the accumulator with the identity matrix, not zeros; it plays the same role as 1 in scalar exponentiation.
  • Apply % MOD after every matrix multiply to prevent overflow, not once at the end.
  • Reach for this when the recurrence is linear, n is large (10^9 or 10^18), and you need a modular answer.

You know the Fibonacci DP solution. Two variables, one loop, O(n) time. Perfectly respectable. Solid citizen of the algorithm world.

What about n = 10^18?

At a billion operations per second, that's a billion seconds of runtime. Roughly 31 years. Your interviewer will not wait. Matrix exponentiation solves it in 60 operations. Yes, 60.

Alternative Big O notation, O(log n) = O(nice)

O(n) for n = 10^18 is technically correct. O(log n) is correct and also your career.

You Can Already Do This With Numbers

Binary exponentiation computes x^n in O(log n) multiplications instead of O(n). The trick is repeated squaring:

x^8 = ((x^2)^2)^2     # 3 multiplications, not 8
x^13 = x^8 * x^4 * x  # read the bits of 13 (1101₂)

Every time the exponent halves, you save half the work. Matrix exponentiation applies the exact same logic to square matrices. The algorithm is identical. You are literally just doing the same thing, except your "numbers" are tiny grids of numbers. Computers love tiny grids.

The Fibonacci Matrix

The Fibonacci recurrence F(n) = F(n-1) + F(n-2) rewrites as a matrix-vector product:

|F(n+1)|   |1  1|   |F(n)  |
|F(n)  | = |1  0| * |F(n-1)|

Apply this matrix once and you advance one step. Apply it n times and you jump from the base case straight to F(n). No intermediate steps needed. No loop counting to a billion.

[[1,1],[1,0]]^n = [[F(n+1), F(n)],
                  [F(n),   F(n-1)]]

To compute F(n), raise the transition matrix to the nth power and read position [0][1]. You can verify this by hand in about two minutes, which I recommend because it's the kind of thing that feels like magic until you do it once:

[[1,1],[1,0]]^2 = [[2,1],[1,1]]   → [0][1] = 1 = F(2) ✓
[[1,1],[1,0]]^3 = [[3,2],[2,1]]   → [0][1] = 2 = F(3) ✓
[[1,1],[1,0]]^5 = [[8,5],[5,3]]   → [0][1] = 5 = F(5) ✓

Classic distracted boyfriend meme, interviewer turning away from your GitHub profile toward "algorithms 101 questions"

They don't care about your five years of production code. They want Fibonacci for n = 10^18.

Binary Exponentiation, but for Matrices

The iterative approach maintains one invariant: result * M^remaining_exponent = M^original_exponent. Say that five times fast. Or just look at the loop:

result   = identity matrix I  (like 1 for numbers)
base     = M
exponent = n

while exponent > 0:
    if exponent is odd:
        result = result * base     # absorb one factor
    base = base * base             # square the base
    exponent = exponent >> 1       # halve the exponent

Start with the identity (the matrix equivalent of 1, the thing that does nothing when multiplied). Each time the exponent is odd, multiply the accumulator by the current base. Square the base, halve the exponent. After O(log n) iterations, result holds M^n.

This is identical to binary exponentiation for scalars. Same bits, same loop structure, different multiplication operation. The only mental leap is accepting that matrices can be "multiplied" just like scalars. They can. Linear algebra said so.

Time and Space

One k×k matrix multiply costs O(k³). Binary exponentiation does O(log n) of them.

Total: O(k³ log n).

For Fibonacci, k = 2, so O(8 log n) = O(log n). For tribonacci (F(n) = F(n-1) + F(n-2) + F(n-3)), k = 3, so O(27 log n). Still O(log n). Space is O(k²) for the two tiny matrices you're keeping in memory.

Plain iterative DP runs in O(n) time with O(1) space. Matrix exponentiation wins only when n is large enough that O(n) is too slow, roughly n > 10^7 in practice. Below that, DP's smaller constants win and the whole matrix machinery is overkill. Don't be a hero.

Beyond Fibonacci: Any Linear Recurrence

Here's the actual point, and it's bigger than Fibonacci. Any k-th order linear recurrence can be computed in O(k³ log n). The Fibonacci matrix is just the most photogenic example.

A k-th order recurrence looks like this:

a(n) = c₁·a(n-1) + c₂·a(n-2) + ... + cₖ·a(n-k)

The transition matrix is the companion matrix. For tribonacci (k = 3, all coefficients 1):

|1  1  1|
|1  0  0|
|0  1  0|

Climbing stairs with 1, 2, or 3 steps is tribonacci. Count Vowels Permutation (LeetCode 1220) is a 5-state linear recurrence. Any time the next value is a fixed linear combination of the previous k values, you get a k×k matrix and O(k³ log n) time. Non-homogeneous recurrences (ones with a constant term) fit too, by augmenting with one extra row and column.

When Does This Actually Show Up?

Most of the time it doesn't. If n ≤ 50, just write the DP loop and move on.

The signal that you need matrix exponentiation is three things together:

  • The recurrence is linear
  • n is enormous: 10^9, 10^12, 10^18
  • The answer needs a modulus (just apply % MOD after every matrix multiply, not once at the end, or intermediate values will overflow before you can say "wrong answer")

In standard FAANG-style interviews this technique almost never appears. In quant trading rounds, competitive programming, and at firms like Jane Street or HRT, it shows up in harder problem sets.

Competitive programmers vs LeetCode grinders: "look what they need to mimic a fraction of our power"

The guy who solved 600 LeetCode mediums encountering a Codeforces grandmaster at a quant shop.

The value in knowing this isn't memorizing the implementation. It's pattern recognition: when you see a linear recurrence plus a suspiciously large n, you immediately know an O(log n) path exists. Naming that approach earns signal even if you don't fully implement it under pressure.

That recognition is exactly what SpaceComplexity surfaces in mock interviews. The first time you hear "compute this for n up to 10^18" should not be live in front of an interviewer.

One Structure, Every Language

Each implementation below has the same three pieces: mat_mul multiplies two square matrices, mat_pow applies binary exponentiation iteratively, and fib wraps both. The n = 0 base case is handled separately since M^n[0][1] gives F(n) only for n ≥ 1.

Python 3

from typing import List

def mat_mul(A: List[List[int]], B: List[List[int]]) -> List[List[int]]:
    n = len(A)
    C = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                C[i][j] += A[i][k] * B[k][j]
    return C

def mat_pow(M: List[List[int]], p: int) -> List[List[int]]:
    n = len(M)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    while p > 0:
        if p & 1:
            result = mat_mul(result, M)
        M = mat_mul(M, M)
        p >>= 1
    return result

def fib(n: int) -> int:
    if n == 0:
        return 0
    return mat_pow([[1, 1], [1, 0]], n)[0][1]

Python 2

def mat_mul(A, B):
    n = len(A)
    C = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                C[i][j] += A[i][k] * B[k][j]
    return C

def mat_pow(M, p):
    n = len(M)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    while p > 0:
        if p & 1:
            result = mat_mul(result, M)
        M = mat_mul(M, M)
        p >>= 1
    return result

def fib(n):
    if n == 0:
        return 0
    return mat_pow([[1, 1], [1, 0]], n)[0][1]

JavaScript

function matMul(A, B) {
    const n = A.length;
    const C = Array.from({ length: n }, () => new Array(n).fill(0));
    for (let i = 0; i < n; i++)
        for (let j = 0; j < n; j++)
            for (let k = 0; k < n; k++)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

function matPow(M, p) {
    const n = M.length;
    let result = Array.from({ length: n }, (_, i) =>
        Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))
    );
    while (p > 0) {
        if (p & 1) result = matMul(result, M);
        M = matMul(M, M);
        p >>= 1;
    }
    return result;
}

function fib(n) {
    if (n === 0) return 0;
    return matPow([[1, 1], [1, 0]], n)[0][1];
}

TypeScript

function matMul(A: number[][], B: number[][]): number[][] {
    const n = A.length;
    const C: number[][] = Array.from({ length: n }, () => new Array(n).fill(0));
    for (let i = 0; i < n; i++)
        for (let j = 0; j < n; j++)
            for (let k = 0; k < n; k++)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

function matPow(M: number[][], p: number): number[][] {
    const n = M.length;
    let result: number[][] = Array.from({ length: n }, (_, i) =>
        Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))
    );
    while (p > 0) {
        if (p & 1) result = matMul(result, M);
        M = matMul(M, M);
        p >>= 1;
    }
    return result;
}

function fib(n: number): number {
    if (n === 0) return 0;
    return matPow([[1, 1], [1, 0]], n)[0][1];
}

Java

static long[][] matMul(long[][] A, long[][] B) {
    int n = A.length;
    long[][] C = new long[n][n];
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

static long[][] matPow(long[][] M, long p) {
    int n = M.length;
    long[][] result = new long[n][n];
    for (int i = 0; i < n; i++) result[i][i] = 1;
    while (p > 0) {
        if ((p & 1) == 1) result = matMul(result, M);
        M = matMul(M, M);
        p >>= 1;
    }
    return result;
}

static long fib(int n) {
    if (n == 0) return 0;
    long[][] base = {{1, 1}, {1, 0}};
    return matPow(base, n)[0][1];
}

C++

#include <vector>
using Matrix = std::vector<std::vector<long long>>;

Matrix matMul(const Matrix& A, const Matrix& B) {
    int n = A.size();
    Matrix C(n, std::vector<long long>(n, 0));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Matrix matPow(Matrix M, long long p) {
    int n = M.size();
    Matrix result(n, std::vector<long long>(n, 0));
    for (int i = 0; i < n; i++) result[i][i] = 1;
    while (p > 0) {
        if (p & 1) result = matMul(result, M);
        M = matMul(M, M);
        p >>= 1;
    }
    return result;
}

long long fib(int n) {
    if (n == 0) return 0;
    return matPow({{1, 1}, {1, 0}}, n)[0][1];
}

C

#include <string.h>
typedef long long ll;

typedef struct { ll v[2][2]; } Mat;

static Mat mat_mul(Mat A, Mat B) {
    Mat C = {{{0, 0}, {0, 0}}};
    for (int i = 0; i < 2; i++)
        for (int j = 0; j < 2; j++)
            for (int k = 0; k < 2; k++)
                C.v[i][j] += A.v[i][k] * B.v[k][j];
    return C;
}

static Mat mat_pow(Mat M, long long p) {
    Mat result = {{{1, 0}, {0, 1}}};
    while (p > 0) {
        if (p & 1) result = mat_mul(result, M);
        M = mat_mul(M, M);
        p >>= 1;
    }
    return result;
}

ll fib(int n) {
    if (n == 0) return 0;
    Mat base = {{{1, 1}, {1, 0}}};
    return mat_pow(base, n).v[0][1];
}

Go

type Matrix [2][2]int64

func matMul(A, B Matrix) Matrix {
    var C Matrix
    for i := 0; i < 2; i++ {
        for j := 0; j < 2; j++ {
            for k := 0; k < 2; k++ {
                C[i][j] += A[i][k] * B[k][j]
            }
        }
    }
    return C
}

func matPow(M Matrix, p int64) Matrix {
    result := Matrix{{1, 0}, {0, 1}}
    for p > 0 {
        if p&1 == 1 {
            result = matMul(result, M)
        }
        M = matMul(M, M)
        p >>= 1
    }
    return result
}

func fib(n int) int64 {
    if n == 0 {
        return 0
    }
    base := Matrix{{1, 1}, {1, 0}}
    return matPow(base, int64(n))[0][1]
}

Rust

type Matrix = [[i64; 2]; 2];

fn mat_mul(a: &Matrix, b: &Matrix) -> Matrix {
    let mut c = [[0i64; 2]; 2];
    for i in 0..2 {
        for j in 0..2 {
            for k in 0..2 {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
    c
}

fn mat_pow(mut m: Matrix, mut p: u64) -> Matrix {
    let mut result: Matrix = [[1, 0], [0, 1]];
    while p > 0 {
        if p & 1 == 1 {
            result = mat_mul(&result, &m);
        }
        m = mat_mul(&m, &m);
        p >>= 1;
    }
    result
}

fn fib(n: u64) -> i64 {
    if n == 0 {
        return 0;
    }
    mat_pow([[1, 1], [1, 0]], n)[0][1]
}

Swift

typealias Matrix = [[Int]]

func matMul(_ A: Matrix, _ B: Matrix) -> Matrix {
    let n = A.count
    var C = Array(repeating: Array(repeating: 0, count: n), count: n)
    for i in 0..<n {
        for j in 0..<n {
            for k in 0..<n {
                C[i][j] += A[i][k] * B[k][j]
            }
        }
    }
    return C
}

func matPow(_ M: Matrix, _ p: Int) -> Matrix {
    let n = M.count
    var result = (0..<n).map { i in (0..<n).map { j in i == j ? 1 : 0 } }
    var base = M
    var exp = p
    while exp > 0 {
        if exp & 1 == 1 { result = matMul(result, base) }
        base = matMul(base, base)
        exp >>= 1
    }
    return result
}

func fib(_ n: Int) -> Int {
    if n == 0 { return 0 }
    return matPow([[1, 1], [1, 0]], n)[0][1]
}

Kotlin

typealias Matrix = Array<LongArray>

fun matMul(A: Matrix, B: Matrix): Matrix {
    val n = A.size
    val C = Array(n) { LongArray(n) }
    for (i in 0 until n)
        for (j in 0 until n)
            for (k in 0 until n)
                C[i][j] += A[i][k] * B[k][j]
    return C
}

fun matPow(M: Matrix, p: Long): Matrix {
    val n = M.size
    var result = Array(n) { i -> LongArray(n) { j -> if (i == j) 1L else 0L } }
    var base = M
    var exp = p
    while (exp > 0) {
        if (exp and 1L == 1L) result = matMul(result, base)
        base = matMul(base, base)
        exp = exp shr 1
    }
    return result
}

fun fib(n: Int): Long {
    if (n == 0) return 0L
    val base = arrayOf(longArrayOf(1, 1), longArrayOf(1, 0))
    return matPow(base, n.toLong())[0][1]
}

C#

static long[,] MatMul(long[,] A, long[,] B) {
    int n = A.GetLength(0);
    var C = new long[n, n];
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                C[i, j] += A[i, k] * B[k, j];
    return C;
}

static long[,] MatPow(long[,] M, long p) {
    int n = M.GetLength(0);
    var result = new long[n, n];
    for (int i = 0; i < n; i++) result[i, i] = 1;
    while (p > 0) {
        if ((p & 1) == 1) result = MatMul(result, M);
        M = MatMul(M, M);
        p >>= 1;
    }
    return result;
}

static long Fib(int n) {
    if (n == 0) return 0;
    var baseMatrix = new long[,] { { 1, 1 }, { 1, 0 } };
    return MatPow(baseMatrix, n)[0, 1];
}

Ruby

def mat_mul(a, b)
  n = a.size
  c = Array.new(n) { Array.new(n, 0) }
  n.times do |i|
    n.times do |j|
      n.times do |k|
        c[i][j] += a[i][k] * b[k][j]
      end
    end
  end
  c
end

def mat_pow(m, p)
  n = m.size
  result = Array.new(n) { |i| Array.new(n) { |j| i == j ? 1 : 0 } }
  while p > 0
    result = mat_mul(result, m) if p.odd?
    m = mat_mul(m, m)
    p >>= 1
  end
  result
end

def fib(n)
  return 0 if n == 0
  mat_pow([[1, 1], [1, 0]], n)[0][1]
end

PHP

function matMul(array $A, array $B): array {
    $n = count($A);
    $C = array_fill(0, $n, array_fill(0, $n, 0));
    for ($i = 0; $i < $n; $i++)
        for ($j = 0; $j < $n; $j++)
            for ($k = 0; $k < $n; $k++)
                $C[$i][$j] += $A[$i][$k] * $B[$k][$j];
    return $C;
}

function matPow(array $M, int $p): array {
    $n = count($M);
    $result = array_fill(0, $n, array_fill(0, $n, 0));
    for ($i = 0; $i < $n; $i++) $result[$i][$i] = 1;
    while ($p > 0) {
        if ($p & 1) $result = matMul($result, $M);
        $M = matMul($M, $M);
        $p >>= 1;
    }
    return $result;
}

function fib(int $n): int {
    if ($n === 0) return 0;
    return matPow([[1, 1], [1, 0]], $n)[0][1];
}

Key Takeaways

  • Matrix exponentiation computes M^n in O(k³ log n) by applying binary exponentiation to matrix multiplication, where k is the matrix size.
  • The Fibonacci transition matrix is [[1,1],[1,0]]. Position [0][1] of M^n gives F(n).
  • Any k-th order linear recurrence maps to a k×k companion matrix, giving the same O(k³ log n) time.
  • Start the accumulator with the identity matrix, not zeros. The identity is the matrix equivalent of 1.
  • Apply % MOD after each multiply, not once at the end, or intermediate values overflow.
  • Related: recursion time complexity, Master Theorem.

Further Reading