xiii. Matrix Exponentiation DP

If the state vector evolves by a fixed linear map, $v_{i+1}=M{v}_i$, then $v_n=M^n{v}_0$ and $n$ up to $10^{18}$ costs $O(k^3\log n)$ for $k$ states. The same works in any semiring: $(\min ,+)$ gives shortest walks with exactly $n$ edges, $(\max ,\min )$ bottleneck walks, $(\vee ,\wedge )$ reachability in exactly $n$ steps. Matrix arithmetic, determinants and fast recurrence methods are in Matrices and Linear Recurrences.

using Mat = vector<vector<ll>>; const ll MOD = 1e9 + 7;
Mat mul(const Mat& a, const Mat& b) {
    int n = sz(a), m = sz(b[0]), K = sz(b); Mat c(n, vector<ll>(m));
    for (int i = 0; i < n; i++) for (int k = 0; k < K; k++) if (a[i][k]) 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, ll e) { Mat r(sz(a), vector<ll>(sz(a))); for (int i = 0; i < sz(a); i++) r[i][i] = 1; for (; e; e >>= 1, a = mul(a, a)) if (e & 1) r = mul(r, a); return r; }

(Loop order $i,k,j$ walks rows contiguously; for $k\le 100$ with a $10^9$ modulus, accumulate in unsigned __int128 or reduce every 16 terms with ull to save the inner %.)

Building the matrix

Recurrence	State vector
$f_n={\sum }_{i=1}^k{c}_i{f}_{n-i}$	$(f_n,\dots ,f_{n-k+1})$, companion matrix
$+$ constant $c$	append a constant $1$
$+$ polynomial $n^d$	append $n^d,n^{d-1},\dots ,1$; update with binomial coefficients
$+$ geometric $r^n$	append $r^n$ (multiplied by $r$ each step)
prefix sums $S_n={\sum }_{i\le n}{f}_i$	append $S$
several interleaved sequences	stack them
periodic transitions $M_1,\dots ,M_p$	exponentiate $P=M_p\cdots M_1$, then apply the remainder
walks of length $n$ in a graph	adjacency matrix
strings of length $n$ avoiding patterns	automaton (KMP / Aho–Corasick) transition counts
tilings of a $k\times n$ board	profile-to-profile transfer matrix (Profile DP)
Markov chain after $n$ steps	transition matrix (probabilities mod $p$ or doubles)

Proof: matrix powers count walks

By induction on $n$: $(A^n{)}_{ij}$ = number of walks of length $n$ from $i$ to $j$. For $n=1$ it’s the adjacency matrix. A walk of length $n+1$ from $i$ to $j$ is a walk of length $n$ to some $k$ followed by an edge $k\to j$, so the count is ${\sum }_k(A^n{)}_{ik}{A}_{kj}=(A^{n+1}{)}_{ij}$. Any linear DP is a weighted walk count in its state graph, and in a semiring the same argument holds with $(\oplus ,\otimes )$ in place of $(+,\times )$, which is exactly why $(\min ,+)$ powers give shortest walks with exactly $n$ edges. Associativity of the semiring matrix product justifies fast exponentiation.

Sums of powers and many queries

• ${\sum }_{i=0}^{n-1}{M}^i$: exponentiate the block matrix $\left(\begin{array}{cc}M& I\\ 0& I\end{array}\right)$, whose $n$-th power has ${\sum }_{i<n}{M}^i$ in the top-right block. Or divide and conquer: $S_{2m}=S_m(I+M^m)$.
• Many queries $M^{n_q}v$: precompute $M^{2^j}$ once ($O(k^3\log n)$), then each query multiplies the vector by $O(\log n)$ matrices: $O(k^2\log n)$ per query instead of $O(k^3\log n)$.
• Min-plus powers with “at most $n$ edges”: add a self-loop of weight 0 at every vertex.
• Large $k$ but a linear recurrence: use Kitamasa ($O(k^2\log n)$) or Bostan–Mori ($O(k\log k\log n)$); if only the sequence is known, recover the recurrence with Berlekamp–Massey from $2k$ terms. Concretely, for $k=100$, $n=10^{18}$: matrix exponentiation costs $\sim 6\cdot 10^7$ operations, Kitamasa $\sim 6\cdot 10^5$ — a real difference once $k$ is in the hundreds.
• $M$ with special structure: circulant matrices multiply as cyclic convolutions ($O(k\log k)$ or $O(k^2)$ per product); upper-triangular with equal diagonals is a polynomial in a shift.

Proof: the block matrix gives the geometric sum

Let $B=\left(\begin{array}{cc}M& I\\ 0& I\end{array}\right)$. Claim $B^n=\left(\begin{array}{cc}{M}^n& S_n\\ 0& I\end{array}\right)$ with $S_n={\sum }_{i<n}{M}^i$. True for $n=1$ ($S_1=I$). Then $B^{n+1}=B^{n}B$ has top-left block $M^{n}M=M^{n+1}$ and top-right block $M^n\cdot I+S_n\cdot I=M^n+S_n=S_{n+1}$, while the bottom row stays $(0,I)$. So the claim holds by induction.

Pitfalls

• The state must be closed under the transition: every quantity used on the right-hand side must be in the vector.
• Check that $M$ does not depend on $i$; if it depends on $i\mathrm{mod}p$ use the period product; if it depends on $i$ polynomially, you usually need extra state for powers of $i$.
• $(\min ,+)$ matrices need $+\infty$ handling: skip infinite entries rather than adding them.

Doubling DP: when the transition isn’t linear in a fixed-size vector

Not every “huge $n$, small extra parameter $K$” DP fits the fixed-matrix mould above. When a quantity over $\{1,\dots ,n\}$ supports both “append one more element” and “combine two halves directly” without going through a transition matrix at all, exponentiation-by-doubling on $n$ itself can beat even Kitamasa.

Example: $e_k(1,\dots ,n)={\sum }_{1\le i_1<\cdots <i_k\le n}{i}_1{i}_2\cdots i_k$ (the elementary symmetric polynomial of $\{1,\dots ,n\}$) for every $k\le K$ simultaneously, $n$ up to $10^{18}$, in $O(K^2\log n)$ — a transition-matrix approach would need $K\times K$ state and cost $O(K^3\log n)$.

template <class Z>
vector<Z> elementarySymmetricUpTo(i64 n, int K) {
    // f(n) = e_k over {0, 1, ..., n}; since 0 contributes nothing, this is e_k over {1..n}
    vector<Z> inv(K + 2, Z(1));
    for (int i = 1; i <= K + 1; i++) inv[i] = Z(1) / Z(i);
    auto rec = [&](auto&& self, i64 m) -> vector<Z> {
        vector<Z> ans(K + 1, Z(0));
        if (m == 0) {
            ans[0] = 1;
            return ans;
        }
        if (m & 1) {   // {0..m} = {0..m-1} + {m}
            vector<Z> l = self(self, m - 1);
            for (int i = 0; i <= K; i++) {
                ans[i] = l[i];
                if (i) ans[i] += l[i - 1] * Z(m);
            }
            return ans;
        }
        // {0..m} = {0..p} + {p + 1 .. 2p} with p = m / 2: elements p + j for j = 1..p
        i64 p = m / 2;
        vector<Z> l = self(self, p);
        // r = e over {p + 1, ..., 2p} = sum_j e_j(1..p) * p^{i-j} * C(p - j, i - j)
        vector<Z> r(K + 1, Z(0)), pw(K + 1);
        pw[0] = 1;
        for (int i = 1; i <= K; i++) pw[i] = pw[i - 1] * Z(p);
        for (int j = 0; j <= K; j++) {
            if (p - j < 0) break;
            Z c = 1;   // C(p - j, t) built incrementally in t
            for (int t = 0; j + t <= K; t++) {
                if (p - j - t < 0) break;
                r[j + t] += l[j] * pw[t] * c;
                c *= Z(p - j - t);
                c *= inv[t + 1];
            }
        }
        for (int i = 0; i <= K; i++) {
            for (int j = 0; i + j <= K; j++) ans[i + j] += l[i] * r[j];
        }
        return ans;
    };
    return rec(rec, n);
}

Proof: the odd and even steps of the doubling recurrence

Work with generating functions: $[x^k]{\prod }_{i=1}^m(1+ix)=e_k(1,\dots ,m)$, since choosing the $ix$ term from exactly $k$ of the $m$ factors and $1$ from the rest sums every $k$-subset product exactly once.
Odd step ($m\to m-1,m$): ${\prod }_{i=1}^m(1+ix)=({\prod }_{i=1}^{m-1}(1+ix))(1+mx)$, so $e_k(1,\dots ,m)=e_k(1,\dots ,m-1)+m\cdot e_{k-1}(1,\dots ,m-1)$ — exactly ans[i] = l[i] + m*l[i-1].
Even step ($m=2p$, splitting into $\{1,\dots ,p\}$ and $\{p+1,\dots ,2p\}$): the second half is the first half’s generating function with every element shifted by $c=p$. For a shift by a constant $c$ applied to $p$ values with elementary symmetric polynomials $e_t$, write $1+(a_j+c)x=(1+cx)(1+\frac{a_{j}x}{1+cx})$, so ${\prod }_j(1+(a_j+c)x)=(1+cx{)}^{p}F(\frac{x}{1+cx})$ where $F(x)={\prod }_j(1+a_{j}x)$. Substituting $F(x)={\sum }_t{e}_t{x}^t$ gives ${\sum }_t{e}_t{x}^t(1+cx{)}^{p-t}$, and extracting $[x^i]$: $e_i(a+c)={\sum }_{t\le i}{e}_t\left(\genfrac{}{}{0ex}{}{p-t}{i-t}\right)c^{i-t}$ — exactly the code’s r[j+t] += l[j] * pw[t] * C(p-j, t) with $a=\{1,\dots ,p\}$, $c=p$. The final combination ans[i+j] += l[i]*r[j] is just multiplying the two halves’ generating functions. $■$

Tested against a direct $O(nK)$ DP (append elements $1,\dots ,n$ one at a time, truncating at degree $K$) for random $n\le 60$, $K\le 8$.

This is the general pattern for a DP over a huge range $n$ that supports both “append one element” ($n\to n+1$) and “double the range” ($n\to 2n$) — not every such DP needs a transition matrix. The $e_k(1,\dots ,n)$ here are, up to reindexing, unsigned Stirling numbers of the first kind; for $K$ close to $n$ the direct $O(nK)$ or $O(n\log n)$ polynomial methods (Stirling Numbers) are the better choice instead.