Matrix Exponentiation on DP

When a DP’s transition is linear and the same at every step, the whole DP is a matrix, and $n$ steps become one matrix power — $O(k^3\log n)$ instead of $O(nk)$. This is what makes $n=10^{18}$ tractable.

The condition

Your recurrence must be expressible as

${\vec{v}}_i=M{\vec{v}}_{i-1}$

with a constant matrix $M$. Then ${\vec{v}}_n=M^n{\vec{v}}_0$, computable by binary exponentiation on matrices.

If $M$ depends on $i$, this fails — unless the dependence is periodic, in which case multiply one period’s matrices together and exponentiate that.

Building the matrix

Fibonacci.
$\left[\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}{F}_{n-1}\\ F_{n-2}\end{array}\right]$

General linear recurrence $f_n={\sum }_{i=1}^k{c}_i{f}_{n-i}$: the companion matrix

$M=\left[\begin{array}{ccccc}{c}_1& c_2& \cdots & c_{k-1}& c_k\\ 1& 0& \cdots & 0& 0\\ 0& 1& \cdots & 0& 0\\ ⋮& & \ddots & & ⋮\\ 0& 0& \cdots & 1& 0\end{array}\right]$

Adding a constant term $f_n=f_{n-1}+c$: extend the vector with a permanent 1.

$\left[\begin{array}{c}{f}_n\\ 1\end{array}\right]=\left[\begin{array}{cc}1& c\\ 0& 1\end{array}\right]\left[\begin{array}{c}{f}_{n-1}\\ 1\end{array}\right]$

Prefix sums $S_n=S_{n-1}+f_n$: add a row for $S$.

Polynomial terms $f_n=f_{n-1}+n^2$: extend with $n^2,n,1$ and encode the binomial updates $(n+1{)}^2=n^2+2n+1$.

Code

using Mat = vector<vector<long long>>;
const long long MOD = 1e9 + 7;
 
Mat mul(const Mat& a, const Mat& b) {
    int n = a.size(), m = b[0].size(), K = b.size();
    Mat c(n, vector<long long>(m, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < K; k++) {
            if (!a[i][k]) continue;
            for (int j = 0; j < m; j++)
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % MOD;
        }
    return c;
}
 
Mat power(Mat a, long long e) {
    int n = a.size();
    Mat r(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++) r[i][i] = 1;
    while (e) { if (e & 1) r = mul(r, a); a = mul(a, a); e >>= 1; }
    return r;
}

The if (!a[i][k]) continue; and the $i,k,j$ loop order together give a 3-5× speedup over the naive $i,j,k$ version — the inner loop then walks both b and c sequentially.

Where it applies

Problem	Matrix size
Linear recurrence of order $k$	$k\times k$
Count walks of length $n$ in a graph	$V\times V$ — the adjacency matrix
Count strings of length $n$ avoiding a pattern	states of the Aho-Corasick automaton
Count strings with a forbidden substring	KMP automaton states
Tilings of a $k\times n$ board	$2^k\times 2^k$ over profile masks
Markov chain after $n$ steps	state count
DP over a small state set, huge $n$	number of states
Shortest path with exactly $n$ edges	$V\times V$ in the min-plus semiring

The min-plus semiring

Replace $(+,\times )$ with $(\min ,+)$ in the matrix multiplication and the same exponentiation computes shortest paths using exactly $n$ edges:

c[i][j] = min(c[i][j], a[i][k] + b[k][j]);

$O(V^3\log n)$. This is how “cheapest route with exactly $k$ flights” problems are solved for $k$ up to $10^9$.

Beating $O(k^3\log n)$

For a linear recurrence of order $k$, you do not need the full matrix:

Method	Time
Matrix exponentiation	$O(k^3\log n)$
Kitamasa	$O(k^2\log n)$
Bostan-Mori	$O(k\log k\log n)$

All three compute $x^n\mathrm{mod}f(x)$ where $f$ is the characteristic polynomial; they differ only in how fast the polynomial multiplication is. For $k=100$ and $n=10^{18}$, matrix exponentiation is $6\cdot 10^7$ operations and Kitamasa is $6\cdot 10^5$.

If you have the sequence but not the recurrence, Berlekamp-Massey recovers it from the first $2k$ terms — a powerful combination: compute a few terms by brute force, find the recurrence automatically, then jump to term $10^{18}$.

See also: Matrix Exponentiation · Kitamasa · Berlekamp-Massey