Matrix Exponentiation for Linear Recurrences

From Fibonacci to finite fields, the trick of lifting a linear recurrence to a matrix and powering via repeated squaring turns $\Theta(n)$ iteration into $O(k^3 \log n)$ arithmetic operations. But Cayley–Hamilton says we are doing too much work: the characteristic polynomial knows everything, and Fiduccia's algorithm rides that observation down to $O(k^2 \log n)$, or $O(k \log k \log n)$ with FFT. Here is the theory, tight and honest about what is proved.

The Lift

Given a linear recurrence of order $k$ with constant coefficients over a commutative ring $R$, $$x_{n+k} = c_1 x_{n+k-1} + c_2 x_{n+k-2} + \cdots + c_k x_n,$$ stack $k$ consecutive terms into a state vector $v_n = (x_{n+k-1}, \ldots, x_n)^T$. Then $v_{n+1} = A v_n$ with $A$ the companion matrix $$A = \begin{pmatrix} c_1 & c_2 & \cdots & c_{k-1} & c_k \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{pmatrix}.$$ So $v_n = A^n v_0$, and the entire computational question becomes: how fast can we form $A^n$?

Repeated Squaring

Naive iteration is $\Theta(nk)$. The classical win uses the binary expansion $n = \sum_{i} b_i 2^i$: $$A^n = \prod_{i\,:\,b_i=1} A^{2^i}.$$ Compute $A, A^2, A^4, \ldots, A^{2^{\lfloor \log_2 n\rfloor}}$ by squaring, then multiply the selected factors. That is $O(\log n)$ matrix multiplications, each $O(k^\omega)$ where $\omega$ is the matrix multiplication exponent. Schoolbook gives $O(k^3 \log n)$; Strassen gives $O(k^{\log_2 7} \log n) \approx O(k^{2.807} \log n)$. Correctness is associativity plus binary expansion — no subtlety. The only failure mode is floating point with $|A| > 1$, where error compounds; in the integer or modular setting (where this technique actually earns its keep) there is no issue.

Go, Briefly

For the modular setting — which is the regime where this dominates Fiduccia in constants for small $k$:

const MOD = 1_000_000_007

type Mat [][]int64

func mul(A, B Mat) Mat {
    k := len(A)
    C := make(Mat, k)
    for i := range C {
        C[i] = make([]int64, k)
        for l := 0; l < k; l++ {
            if A[i][l] == 0 { continue }
            a := A[i][l]
            for j := 0; j < k; j++ {
                C[i][j] = (C[i][j] + a*B[l][j]) % MOD
            }
        }
    }
    return C
}

func pow(A Mat, n uint64) Mat {
    k := len(A)
    R := make(Mat, k)
    for i := range R { R[i] = make([]int64, k); R[i][i] = 1 }
    for n > 0 {
        if n&1 == 1 { R = mul(R, A) }
        A = mul(A, A)
        n >>= 1
    }
    return R
}

The loop order i, l, j matters: it keeps access to B row-wise and is worth a factor of 2–5 over i, j, l once $k$ is out of L1. That is not pedantry; it is the reason the $O(k^3)$ constant is tolerable up to $k \approx 200$.

Cayley–Hamilton Says You Are Wasting Work

Now the real point. Let $p(\lambda) = \det(\lambda I - A)$ be the characteristic polynomial of $A$. Because $A$ is a companion matrix, $$p(\lambda) = \lambda^k - c_1 \lambda^{k-1} - c_2 \lambda^{k-2} - \cdots - c_k,$$ which you can read directly off the recurrence without computing a determinant.

Cayley–Hamilton theorem. $p(A) = 0$.

Load-bearing consequence: $A^k$ is an $R$-linear combination of $I, A, \ldots, A^{k-1}$. By induction, every power satisfies $$A^n = \sum_{i=0}^{k-1} r_i(n)\, A^i$$ for some scalars $r_i(n)$. We never need the full $k \times k$ matrix $A^n$; we need $k$ scalars.

How to get them without ever forming $A$. Let $q_n(\lambda) = \lambda^n \bmod p(\lambda)$ — the remainder when dividing $\lambda^n$ by $p$ in $R[\lambda]$. The evaluation map $R[\lambda] \to R[A]$, $\lambda \mapsto A$, is a ring homomorphism whose kernel contains $p$ (by Cayley–Hamilton). Therefore $\lambda^n$ and $q_n(\lambda)$ have the same image, i.e. $$A^n = q_n(A).$$ That is the step doing the work. Every other part of the argument is bookkeeping.

Compute $q_n$ by repeated squaring in the quotient ring $R[\lambda]/(p)$: $$q_{2m} \equiv q_m^2 \pmod{p}, \qquad q_{2m+1} \equiv \lambda \cdot q_m^2 \pmod{p}.$$ Each square-and-reduce step multiplies two degree-$(k-1)$ polynomials and reduces a degree-$(2k-2)$ polynomial modulo $p$. Schoolbook: $O(k^2)$ per step. With $O(\log n)$ steps, $$T(n) = O(k^2 \log n).$$ This is **Fiduccia's algorithm** (Fiduccia, 1985), often called *Kitamasa's method* in competitive-programming circles. For $k$ in the hundreds and $n$ astronomical, the gap from $O(k^3 \log n)$ is several orders of magnitude.

Once $q_n(\lambda) = \sum_{i=0}^{k-1} r_i \lambda^i$ is computed, extracting $x_n$ is one dot product against the initial conditions, $$x_n = \sum_{i=0}^{k-1} r_i\, x_i,$$ because $A^n v_0 = q_n(A) v_0$ and the last coordinate of $A^i v_0$ is $x_i$ by construction of the companion form. You touch only $x_0, \ldots, x_{k-1}$.

What the FFT Buys

If $R$ supports suitable roots of unity (prime fields with $2^s \mid (p-1)$, or $\mathbb{C}$), polynomial multiplication drops to $O(k \log k)$. Division by $p$ uses the Newton-iteration reciprocal trick — compute the reciprocal of the reverse of $p$ modulo $\lambda^k$, then multiply — also at $O(k \log k)$. The final bound is $$T(n) = O(k \log k \cdot \log n).$$ For $k = 10^5$ and $n = 10^{18}$, this is the difference between seconds and never. Over $\mathbb{Z}$ with no friendly roots of unity, Schönhage–Strassen gives $O(k \log k \log \log k)$ bit-complexity per multiplication and the asymptotic still holds modulo a log-log factor.

Asymptotics vs. Diagonalization

Why not diagonalize? If $A = PDP^{-1}$ with $D = \mathrm{diag}(\lambda_1, \ldots, \lambda_k)$, then $A^n = P D^n P^{-1}$ looks free. Two honest objections:

1. Exactness. Over $\mathbb{Q}$ or $\mathbb{F}_q$, eigenvalues live in an extension of degree up to $k$; arithmetic there costs. Fiduccia stays in the base ring throughout and requires no inverses — which matters over $\mathbb{Z}/m$ with $m$ composite, where most elements are not invertible but the algorithm still runs verbatim.
2. Repeated roots. If $p$ has repeated roots, $A$ is only similar to a Jordan form; powers acquire $\binom{n}{j} \lambda^{n-j}$ terms. Cayley–Hamilton is oblivious to Jordan structure and always correct.

For analytic asymptotics, diagonalization is the right tool. If $\lambda_1$ is the unique eigenvalue of largest modulus, the Perron–Frobenius theorem (when $A \geq 0$ is irreducible and aperiodic) gives $x_n \sim c \lambda_1^n$ with $c$ read off the dominant eigenvector's pairing with $v_0$. Binet's formula $F_n = (\phi^n - \psi^n)/\sqrt{5}$ is the $k=2$ instance with $\phi, \psi$ roots of $\lambda^2 - \lambda - 1$.

Scope and Pitfalls

The technique is a black box for constant-coefficient linear recurrences over any commutative ring. It degrades in three directions:

• Variable coefficients. $x_{n+1} = a_n x_n + b_n x_{n-1}$ gives a product $A_{n-1} \cdots A_0$ with no squaring structure; you pay $\Theta(n)$ unless the $A_i$ have exploitable algebraic structure.
• D-finite (holonomic) recurrences. Coefficients are rational functions of $n$. Bostan–Mori (2021, building on Chudnovsky–Chudnovsky) achieve $\tilde O(\sqrt{n})$ for one term via baby-step/giant-step polynomial evaluation — a genuinely different algorithm, sharper only for enormous $n$.
• Nonlinear recurrences. Nothing here applies; the orbit structure is generically intractable.

What Is Actually Proved Here

• Matrix exponentiation correctness: trivial, by associativity. Proved.
• Fiduccia/Kitamasa correctness: the load-bearing step is Cayley–Hamilton plus the ring-homomorphism identity $A^n = q_n(A)$ where $q_n = \lambda^n \bmod p$. Proved (sketched above).
• Complexity $O(k^2 \log n)$ and $O(k \log k \log n)$: asserted using standard polynomial-arithmetic bounds (FFT is rigorously $O(k \log k)$ over fields with roots of unity; Schönhage–Strassen over $\mathbb{Z}$). Asserted, not reproved.
• Perron–Frobenius asymptotics: asserted, not proved here.
• Matching lower bounds: open. No known $\Omega(k \log k \log n)$ lower bound for the $n$-th term problem rules out subquadratic-in-$k$ improvements beyond FFT.