Lab 09 — Matrix Exponentiation & Linear Recurrences

Goal

Compute the n-th term of any fixed-order linear recurrence mod p for n up to 10^18 in O(k³ log n) via matrix power, where k is the recurrence order.

Background

A linear recurrence f(n) = c₁f(n−1) + … + c_k f(n−k) is one multiplication of the state vector [f(n−1), …, f(n−k)] by its companion matrix:

[c₁ c₂ … c_{k−1} c_k]
[ 1  0 …    0     0 ]
[ 0  1 …    0     0 ]
[ 0  0 …    1     0 ]

So f(n) reads off M^(n−k+1) applied to the initial state. Binary exponentiation (phase-07 lab-02) needs only associativity — lift it from scalars to matrices unchanged: O(log n) multiplications, O(k³) each.

Fibonacci is the canonical k = 2 case: [[1,1],[1,0]]^n = [[F(n+1), F(n)], [F(n), F(n−1)]].

Graph interpretation: for adjacency matrix A, (A^k)[i][j] counts length-k walks i → j. Every “count sequences of length n under local rules” problem is a walk count on an automaton → matrix power.

The pattern: DP with a huge time dimension. If one DP step is linear (state vector × fixed transition matrix), 10^18 steps collapse to ~60 matrix squarings.

Beyond k³: Kitamasa — exponentiate x^n mod the characteristic polynomial — gives O(k² log n), or O(k log k log n) with FFT multiplication (phase-12 lab-07). Name-drop it; implement only if k ≥ ~10³.

Interview Context

  • Codeforces / ICPC: regular — Div 2 D/E staple
  • Quant: common (recurrences, n-step Markov transitions)
  • FAANG: rare but real — shows up as the “now n = 10^18” follow-up to climbing-stairs / tiling / string-counting DPs
  • Google-style “count binary strings of length 10^15 avoiding pattern X”: exactly this
  • Standard FAANG loop: low probability, but the cheapest hard topic in this phase

When to Skip This Topic

Skip if any of these are true:

  • Scalar binary exponentiation (phase-07 lab-02) isn’t automatic yet
  • You can’t write bottom-up Fibonacci DP cold
  • You’re on a frontend/mobile interview track

Otherwise don’t skip — this is the highest ROI-per-hour lab in the phase.

Problem Statement

Given a, b, c and n ≤ 10^18, compute

f(n) = a·f(n−1) + b·f(n−2) + c (mod 10^9+7)

with f(0), f(1) given. The constant c makes the recurrence affine: a 2×2 matrix no longer suffices — carry a constant “1” lane, giving a 3×3 matrix.

Constraints

  • 0 ≤ n ≤ 10^18
  • Coefficients and initial values in [0, 10^9+7)
  • Modulus 10^9+7, prime
  • Wall-clock: sub-millisecond per query — ~60 squarings of a 3×3

Clarifying Questions

  1. Answer mod what — and is it prime? (10^9+7, prime; primality only matters if division shows up.)
  2. Can coefficients be negative? (In general yes — reduce into [0, mod) before building the matrix.)
  3. One n or many queries? (One here; for many, see follow-ups.)
  4. Is k fixed and small? (k ≤ ~50 for the k³ method; larger k needs Kitamasa.)

Examples

Fibonacci identity at n = 5 (F: 0, 1, 1, 2, 3, 5, 8):

[[1,1],[1,0]]^5 = [[8,5],[5,3]] = [[F(6),F(5)],[F(5),F(4)]]  ✓

Affine 3×3 for f(n) = a·f(n−1) + b·f(n−2) + c:

[f(n)  ]   [a b c]   [f(n−1)]
[f(n−1)] = [1 0 0] × [f(n−2)]
[  1   ]   [0 0 1]   [  1   ]

a=2, b=1, c=3, f(0)=0, f(1)=1: f(2) = 2+0+3 = 5, f(3) = 10+1+3 = 14, f(4) = 28+5+3 = 36.

Brute Force

Iterate the recurrence keeping the last k values. O(n) time. At n = 10^18 and 10^9 steps/sec: ~30 years. Dead.

Brute Force Complexity

  • Time: O(n·k)
  • Space: O(k)

Optimization Path

  1. Express one step as state_new = M × state_old (companion matrix; add the “1” lane for +c).
  2. n steps = M^n × state₀ — the loop becomes a matrix power.
  3. Compute M^n by repeated squaring: log₂(10^18) ≈ 60 levels.
  4. Read f(n) out of one row of the result.
  5. If k grows past ~50: Kitamasa, O(k² log n).

Final Expected Approach

MOD = 10**9 + 7

def mat_mul(A, B, mod=MOD):
    n, m, r = len(A), len(B[0]), len(B)
    C = [[0] * m for _ in range(n)]
    for i in range(n):
        Ci = C[i]
        for k in range(r):
            a = A[i][k]
            if a:
                Bk = B[k]
                for j in range(m):
                    Ci[j] = (Ci[j] + a * Bk[j]) % mod
    return C

def mat_pow(A, e, mod=MOD):
    n = len(A)
    R = [[int(i == j) for j in range(n)] for i in range(n)]
    while e:
        if e & 1:
            R = mat_mul(R, A, mod)
        A = mat_mul(A, A, mod)
        e >>= 1
    return R

def kth_term(coeffs, initial, n, mod=MOD):
    k = len(coeffs)
    if n < k:
        return initial[n] % mod
    M = [[0] * k for _ in range(k)]
    M[0] = [c % mod for c in coeffs]
    for i in range(1, k):
        M[i][i - 1] = 1
    P = mat_pow(M, n - k + 1, mod)
    state = [initial[k - 1 - j] % mod for j in range(k)]
    return sum(P[0][j] * state[j] for j in range(k)) % mod

def affine_term(a, b, c, f0, f1, n, mod=MOD):
    if n == 0:
        return f0 % mod
    M = [[a % mod, b % mod, c % mod], [1, 0, 0], [0, 0, 1]]
    P = mat_pow(M, n - 1, mod)
    return (P[0][0] * f1 + P[0][1] * f0 + P[0][2]) % mod

def fib(n, mod=MOD):
    def pair(n):
        if n == 0:
            return 0, 1
        a, b = pair(n >> 1)
        c = a * (2 * b - a) % mod
        d = (a * a + b * b) % mod
        return (d, (c + d) % mod) if n & 1 else (c, d)
    return pair(n)[0]

Data Structures

  • k×k list-of-lists matrices; no numpy (int64 overflows at this modulus)
  • State vector ordered newest-first: [f(t), f(t−1), …, f(t−k+1)]

Correctness Argument

  • One recurrence step is exactly one companion-matrix multiply: the first row computes Σ cᵢ f(t−i+1); the shifted-identity rows slide each older value down one slot.
  • Matrix multiplication is associative, so binary exponentiation is valid: M^(2e) = (M^e)², M^(e+1) = M·M^e.
  • Induction: s_t = [f(t), …, f(t−k+1)] satisfies s_{t+1} = M·s_t, hence f(n) = (M^(n−k+1)·s_{k−1})[0].
  • The [0 0 1] row preserves the constant lane, so +c is added exactly once per step.

Complexity

  • Time: O(k³ log n) — k = 3, n = 10^18: 27 × ~120 ≈ 3×10³ scalar mulmods
  • Space: O(k²)
  • Fast doubling (Fibonacci only): O(log n), ~3 mulmods per level

Implementation Requirements

  • Reduce mod inside mat_mul, per accumulation — mandatory for overflow in C++/Java, mandatory for speed in Python
  • Seed exponentiation with the identity matrix; e = 0 must return I
  • Answer n < k directly from initial values — no matrix work
  • Exponent is n − k + 1, not n; verify at n = k before trusting anything
  • Reduce negative coefficients into [0, mod) before building M

Tests

  • fib(n) == kth_term([1,1], [0,1], n) == iterative Fibonacci for n = 0..30
  • Affine a=2, b=1, c=3, f0=0, f1=1: f(2)=5, f(3)=14, f(4)=36 (hand values)
  • Edge n = 0, 1, k−1 (pure initial-value reads)
  • k = 1 degenerate case: f(n) = c₁^n · f(0) — matches pow(c1, n, mod) route
  • fib(10^18) returns in milliseconds; matrix and fast-doubling versions agree
  • Randomized k = 4 recurrence vs O(n) iteration for n ≤ 10^4

Follow-up Questions

  • For Fibonacci specifically, can you beat the 2×2 power? → Fast doubling: F(2k) = F(k)·(2F(k+1) − F(k)), F(2k+1) = F(k)² + F(k+1)². Same O(log n), but ~3 mulmods per level vs ~8+ for naive 2×2 squaring — roughly an 8× constant-factor win over unoptimized matrix code.
  • Count walks of length 10^12 between two nodes in a 50-vertex graph. → (A^k)[i][j] by matrix power mod p: 50³ × ~80 multiplications. Same machinery counts length-n strings avoiding a pattern: build the KMP/Aho-Corasick automaton, power its transition matrix.
  • Why must the mod live inside mat_mul? → In C++/Java, one product a·b already fills 64 bits and summing k of them overflows without intermediate reduction — wrong answers, not exceptions. Python can’t overflow; it just degrades into multi-word big-int arithmetic and slows several-fold.
  • k = 10^5 — now what? → k³ log n ≈ 10^15·60: dead. Kitamasa computes x^n mod the characteristic polynomial in O(k² log n), or O(k log k log n) with NTT (phase-12 lab-07). If only the sequence is given, Berlekamp-Massey first recovers the minimal recurrence from 2k terms.
  • All of f(l..r), not a single n? → Matrix-power to f(l), then iterate linearly: O(k³ log l + k·(r−l)). For many scattered queries, precompute the ~60 squarings M^(2^i) once and answer each query with ≤ 60 multiplications, or sort queries and reuse prefix products.
  • Does the recurrence have to be homogeneous? → No. Constants get the “1” lane; polynomial terms like +n or +n² get extra lanes updated Pascal-style ((n+1)² = n² + 2n + 1). Anything linear in the augmented state folds in.

Product Extension

  • Markov-chain systems projected far forward (population models, failure states)
  • State machines in rate limiters: probability of reaching a state within n events
  • Graphics pipelines: composing affine transforms is exactly 3×3/4×4 matrix products
  • Financial recurrences (annuity-style balance updates) at huge horizons

Language/Runtime Follow-ups

  • Python: big ints make it correct by default; keep list-of-lists — numpy int64 silently overflows at 10^9+7 (object dtype works but is slower than plain lists)
  • C++: accumulate in unsigned 64-bit with periodic reduction, or __int128; hand-unroll k = 2, 3
  • Java: long with reduction inside the loop; Math.floorMod for negatives
  • Rust: u128 intermediates; const-generic fixed-size matrices optimize well

Common Bugs

  1. Off-by-one exponent: M^n vs M^(n−1) vs M^(n−k+1) depending on the seeded state — verify at n = k first.
  2. State vector reversed: companion matrix assumes newest-first order; feeding [f(0), …, f(k−1)] gives plausible garbage.
  3. Mod applied too late: overflow in C++/Java; silent multi-fold slowdown in Python.
  4. Missing constant lane: modeling +c in a 2×2 — the error grows by c per step and small tests may still pass.
  5. Aliased squaring: mat_mul writing into an input matrix corrupts A = A·A.
  6. n < k unhandled: matrix path assumes at least one transition; n = 0 must short-circuit.

Debugging Strategy

  • Diff against the O(n) iterative recurrence for n = 0..1000 — catches every off-by-one
  • Check M^0 = I and M^1 = M explicitly
  • Verify [[1,1],[1,0]]^5 = [[8,5],[5,3]] by hand
  • Print the state vector after 1, 2, 3 matrix applications and match the direct recurrence

Mastery Criteria

  • Derive the companion matrix for any dictated recurrence in ≤ 3 min
  • Implement mat_pow + kth_term bug-free in ≤ 25 min
  • Add the affine “+c” lane without hesitation
  • Write Fibonacci fast doubling from memory in ≤ 10 min
  • Spot “linear DP + astronomical n” and say “matrix power” within 2 min of reading a problem
  • Name the large-k escape hatch (Kitamasa / Berlekamp-Massey) with its complexity