Matrix Exponentiation

Start from brute force

The obvious approach starts from the identity and multiplies by mat once per unit of power. It is correct, since M^n is by definition n copies of M multiplied together. It costs O(k³ · power) time for a k×k matrix, because each schoolbook multiply is k³ operations, and O(k²) space for the running product. The multiply count grows linearly with the exponent: power = 10⁹ means a billion matrix multiplies, while repeated squaring needs about 60. That linear growth is exactly what the problems built on this trick, such as F(10¹⁸) mod a prime, are designed to rule out.

The approach

Write the exponent in binary: M^10 = M^8 · M^2 because 10 is 1010. Repeated squaring yields M, M², M⁴, M⁸, … at one multiply each, and the answer is the product of exactly the powers whose bit is set. The reference keeps result, starting at the identity (the matrix version of 1), and base, a copy of mat. Each iteration tests the low bit with power & 1 and, if it is set, does result = matmul(result, base); then base = matmul(base, base) and power >>= 1.

The invariant: result · base^power always equals the original matrix to the original power. It holds at the start, since the identity times mat^power is the answer, and both branches preserve it, so when power reaches 0, result is the answer. matmul builds entry [i][j] as the dot product of row i of A with column j of B.

Walk through an example

For the first example, power = 10 is 1010 in binary, so result picks up M² at step 2 and M⁸ at step 4, and the last squaring is computed but never used.

Complexity

The loop runs ⌊log₂ power⌋ + 1 times, with one squaring each plus one multiply per set bit, so time is O(k³ log power) for a k×k matrix and space is O(k²) for result, base and one product. Repeated multiplication is O(k³ · power) time with the same O(k²) space. Big-integer entries make each multiply slower as numbers grow.

Edge cases to test

• power = 0 returns the identity, even for the zero matrix: matrix_power([[0,0],[0,0]], 0) gives [[1, 0], [0, 1]], by the convention M⁰ = I.
• power = 1 returns the input unchanged: matrix_power([[1,2],[3,4]], 1) gives [[1, 2], [3, 4]].
• A 1×1 matrix is ordinary fast exponentiation: matrix_power([[3]], 4) returns [[81]].
• A negative power is not rejected: the loop never runs, so matrix_power([[1,1],[1,0]], -1) silently returns the identity. Validate the input; a true inverse needs fractions or a modulus.
• Large entries: Python integers are exact at any size, so matrix_power([[1,1],[1,0]], 90)[0][1] is exactly F(90) = 2880067194370816120. TypeScript numbers are 64-bit floats, exact only up to 2⁵³ − 1, so at power 78 the top-left entry, F(79) = 14472334024676221, already comes back as 14472334024676220.

Common mistakes

• Naive repeated multiplication O(n) instead of the O(log n) squaring.
• Wrong base case — the accumulator must start as the identity matrix.
• Losing exactness on large entries. Python integers never overflow, but every TypeScript number is a 64-bit float, exact only up to 2⁵³ − 1; use BigInt, or reduce by a modulus after every multiply when the problem asks for one.
• Raising each entry to the power instead of multiplying matrices: for [[1,2],[3,4]] cubed that gives [[1, 8], [27, 64]] instead of [[37, 54], [81, 118]].

Follow-ups an interviewer asks next

Compute the n-th Fibonacci number mod 10⁹ + 7 for n up to 10¹⁸.

Raise [[1,1],[1,0]] to the n-th power, reducing every entry mod 1_000_000_007 after each multiply, and read [0][1], which is F(n). Since log₂(10¹⁸) is about 60, that is at most about 120 two-by-two multiplies. In TypeScript use BigInt for the products: two numbers near 10⁹ multiply to about 10¹⁸, past 2⁵³.

How does this solve a general linear recurrence?

Put the last k terms in a state vector and write one step as a k×k companion matrix. For a(n) = 2·a(n−1) + 3·a(n−2) the matrix is [[2,3],[1,0]], and [a(n), a(n−1)] = M^(n−1) · [a(1), a(0)]. Any order-k recurrence with constant coefficients then costs O(k³ log n).

What does a power of an adjacency matrix count?

Entry [i][j] of A^n counts the walks of exactly n edges from i to j, because each multiply extends every walk by one edge. For the graph with edges 0→1, 0→2, 1→2 and 2→0, A³ has a 1 at [0][0]: the single closed walk 0 → 1 → 2 → 0.

Can each multiply be faster than k³?

Strassen’s algorithm multiplies in about O(k^2.81) by trading multiplications for additions, and it pays off only for large k. For the small matrices in recurrence problems, the practical wins are a modulus, an unrolled 2×2 multiply, and skipping the final squaring, which the reference computes and never uses.

Where this shows up

Matrix exponentiation evaluates linear recurrences in O(log n) — the classic [[1,1],[1,0]]^n computes Fibonacci fast. It also counts length-n paths in a graph (powers of the adjacency matrix) and advances Markov chains, making it a go-to for “huge n” recurrence and counting problems.

Reference solution

def matrix_power(mat, power):
    n = len(mat)

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

    # start at the identity matrix
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    base = [row[:] for row in mat]
    while power > 0:
        if power & 1:
            result = matmul(result, base)   # multiply in when the bit is set
        base = matmul(base, base)           # square the base
        power >>= 1                         # halve the exponent
    return result


print(matrix_power([[1, 1], [1, 0]], 10))   # -> [[89, 55], [55, 34]]
print(matrix_power([[2, 0], [0, 2]], 3))    # -> [[8, 0], [0, 8]]

function matrixPower(mat: number[][], power: number): number[][] {
  const n = mat.length;
  const matmul = (A: number[][], B: number[][]): number[][] =>
    Array.from({ length: n }, (_, i) =>
      Array.from({ length: n }, (_, j) => {
        let s = 0;
        for (let k = 0; k < n; k++) s += A[i][k] * B[k][j];
        return s;
      }));

  // identity matrix
  let result = Array.from({ length: n }, (_, i) =>
    Array.from({ length: n }, (_, j) => (i === j ? 1 : 0)));
  let base = mat.map((row) => row.slice());
  while (power > 0) {
    if (power & 1) result = matmul(result, base);  // multiply in on a set bit
    base = matmul(base, base);                     // square the base
    power = Math.floor(power / 2);                 // halve the exponent
  }
  return result;
}

console.log(JSON.stringify(matrixPower([[1,1],[1,0]], 10))); // -> [[89,55],[55,34]]
console.log(JSON.stringify(matrixPower([[2,0],[0,2]], 3)));  // -> [[8,0],[0,8]]