i. Matrices and Matrix Power

Matrix arithmetic for contests: a fast multiplication loop, binary exponentiation, many power queries, geometric sums, semiring products ($(\min ,+)$, $(\max ,\min )$, boolean), and Strassen. How to build the transition matrix of a DP is in Matrix Exponentiation DP; recurrences of order $k$ in $O(k\log k\log n)$ are in Linear Recurrences. mint is from Modular Arithmetic.

$(AB{)}_{ij}={\sum }_k{A}_{ik}{B}_{kj}$ for $A$ of size $n\times m$ and $B$ of size $m\times p$. The product is associative but not commutative; $(AB{)}^T=B^T{A}^T$, $(AB{)}^{-1}=B^{-1}{A}^{-1}$, $\mathrm{tr}(AB)=\mathrm{tr}(BA)$.

Task	Method	Cost	Note
$A^e$	binary exponentiation	$O(n^3\log e)$	this note
$A^{e_q}v$ for many queries	precomputed $A^{2^j}$, vector products	$O(n^3\log E+q{n}^2\log E)$	this note
${\sum }_{i<e}{A}^i$	block matrix ${\left(\begin{array}{cc}A& I\\ 0& I\end{array}\right)}^e$	$O(n^3\log e)$	this note
solve $Ax=b$, rank, inverse	Gaussian elimination	$O(n^3)$	Gaussian Elimination
determinant (field, composite modulus, exact integers)	elimination, Euclid-style rows, Bareiss	$O(n^3)$	Determinants
characteristic polynomial	Hessenberg reduction	$O(n^3)$	Characteristic Polynomial
sparse system, sparse determinant	Wiedemann / Lanczos	$O(n\cdot \mathrm{nnz})$	Sparse Matrices
spanning trees	Matrix-Tree	$O(n^3)$	Matrix-Tree Theorem
perfect matching exists	Tutte matrix rank (random)	$O(n^3)$	General Matching
verify $AB=C$	Freivalds	$O(n^2)$ per round	Hashing and Verification
boolean product	bitset rows	$O(n^3/64)$	Bitset

Everything works over ${\mathbb{F}}_p$ with division by inverses; over a composite modulus division may fail, so use Euclid-style row reduction (repeatedly subtract multiples, as in the gcd) instead of pivoting by inverses.

using Mat = vector<vector<mint>>;
Mat eye(int n) { Mat r(n, vector<mint>(n, 0)); for (int i = 0; i < n; i++) r[i][i] = 1; return r; }
Mat operator*(const Mat& a, const Mat& b) {                // i, k, j order: rows of b and c are walked contiguously
    int n = sz(a), K = sz(b), m = sz(b[0]); Mat c(n, vector<mint>(m, 0));
    for (int i = 0; i < n; i++) for (int k = 0; k < K; k++) {
        if (a[i][k] == mint(0)) continue;
        mint x = a[i][k]; for (int j = 0; j < m; j++) c[i][j] += x * b[k][j];
    }
    return c;
}
Mat matPow(Mat a, ll e) { Mat r = eye(sz(a)); for (; e; e >>= 1, a = a * a) if (e & 1) r = r * a; return r; }
vector<mint> mulVec(const Mat& a, const vector<mint>& v) {
    vector<mint> r(sz(a), 0); for (int i = 0; i < sz(a); i++) for (int j = 0; j < sz(v); j++) r[i] += a[i][j] * v[j]; return r;
}
 
struct PowerQueries {                                      // A^e v in O(n^2 log e) per query after O(n^3 log E) setup
    vector<Mat> p;
    PowerQueries(const Mat& a, int LOG = 60) : p(LOG) { p[0] = a; for (int j = 1; j < LOG; j++) p[j] = p[j - 1] * p[j - 1]; }
    vector<mint> query(vector<mint> v, ull e) const { for (int j = 0; e; j++, e >>= 1) if (e & 1) v = mulVec(p[j], v); return v; }
};
 
Mat geometricSum(const Mat& a, ll e) {                     // sum_{i < e} a^i
    int n = sz(a); Mat b(2 * n, vector<mint>(2 * n, 0));
    for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) b[i][j] = a[i][j]; b[i][n + i] = 1; b[n + i][n + i] = 1; }
    Mat r = matPow(b, e), s(n, vector<mint>(n));
    for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) s[i][j] = r[i][n + j];
    return s;
}
 
vector<vector<ll>> mulFast(const vector<vector<ll>>& a, const vector<vector<ll>>& b, ll mod) {   // entries < mod < 2^30: reduce every 16 terms
    int n = sz(a), K = sz(b), m = sz(b[0]); vector<vector<ll>> c(n, vector<ll>(m));
    vector<ull> acc(m);
    for (int i = 0; i < n; i++) {
        fill(all(acc), 0);
        for (int k = 0; k < K; k++) {
            ull x = a[i][k]; for (int j = 0; j < m; j++) acc[j] += x * b[k][j];
            if ((k & 15) == 15) for (int j = 0; j < m; j++) acc[j] %= mod;
        }
        for (int j = 0; j < m; j++) c[i][j] = acc[j] % mod;
    }
    return c;
}

Proof: exponentiation, power queries, geometric sums and the delayed reduction

Binary exponentiation is correct in any monoid: with $e={\sum }_j{b}_j{2}^j$, the loop keeps $a=A^{2^j}$ and multiplies $r$ by it when $b_j=1$, so at the end $r={\prod }_j{A}^{b_j{2}^j}=A^e$ (all factors are powers of $A$, so they commute, and associativity makes the grouping irrelevant). Power queries: the same decomposition applied to a vector, $A^{e}v={\prod }_j(A^{2^j}{)}^{b_j}v$, with each factor applied as an $O(n^2)$ matrix-vector product. Geometric sum: $B^e=\left(\begin{array}{cc}{A}^e& S_e\\ 0& I\end{array}\right)$ with $S_e={\sum }_{i<e}{A}^i$, by induction ($B^{e+1}=B^{e}B$ has top-right block $A^e+S_e$). Delayed reduction: after a reduction every accumulator is $<mod$, and 16 more products add at most $16(mod-1{)}^2$, so the accumulator stays below $16mo{d}^2+mod<2^{64}$ whenever $mod<2^{30}$ (true for $10^9+7$ and 998244353). Checked against the plain product.

Semirings

Matrix multiplication and fast exponentiation work over any semiring $(\oplus ,\otimes )$: $\oplus$ associative and commutative with identity $\bar{0}$, $\otimes$ associative with identity $\bar{1}$, $\otimes$ distributing over $\oplus$, and $\bar{0}$ absorbing. $(A^{\otimes n}{)}_{ij}$ then aggregates (with $\oplus$) the $\otimes$-products of edge weights over all walks of exactly $n$ edges.

Semiring $(\oplus ,\otimes )$	$\bar{0}$, $\bar{1}$	$(A^n{)}_{ij}$ computes
$(+,\times )$	0, 1	number (or total weight) of walks of length $n$
$(\min ,+)$	$+\infty$, 0	shortest walk with exactly $n$ edges
$(\max ,+)$	$-\infty$, 0	longest walk with exactly $n$ edges
$(\vee ,\wedge )$	false, true	reachability in exactly $n$ steps
$(\max ,\min )$	$-\infty$, $+\infty$	widest (bottleneck) walk
$({\oplus }_{\text{xor}},\wedge )$	0, 1	parity of the number of walks

“At most $n$ edges” with $(\min ,+)$: add a 0-weight self-loop at every vertex. Max-plus: negate the weights. Cheapest route using exactly $k$ flights for $k\le 10^9$ is $O(V^3\log k)$ with this, a problem with no other reasonable approach. In $(\min ,+)$ there is no subtraction, so no Strassen and no $O(n^{3-\epsilon })$ algorithm is known (the APSP conjecture, Theoretical Algorithms).

struct MinPlusMatrix {
    static constexpr ll INF = numeric_limits<ll>::max() / 4;
    int n; vector<vector<ll>> d;
    explicit MinPlusMatrix(int size) : n(size), d(size, vector<ll>(size, INF)) {}
    static MinPlusMatrix identity(int size) { MinPlusMatrix r(size); for (int i = 0; i < size; i++) r.d[i][i] = 0; return r; }
    friend MinPlusMatrix operator*(const MinPlusMatrix& a, const MinPlusMatrix& b) {
        MinPlusMatrix c(a.n);
        for (int i = 0; i < a.n; i++) for (int k = 0; k < a.n; k++) {
            if (a.d[i][k] >= INF) continue;
            for (int j = 0; j < a.n; j++) if (b.d[k][j] < INF) c.d[i][j] = min(c.d[i][j], a.d[i][k] + b.d[k][j]);
        }
        return c;
    }
    MinPlusMatrix pow(ll k) const { MinPlusMatrix res = identity(n), b = *this; for (; k > 0; k >>= 1, b = b * b) if (k & 1) res = res * b; return res; }
};
 
template <class T, class Add, class Mul>                   // generic semiring product (zero = identity of add)
vector<vector<T>> semiringMul(const vector<vector<T>>& a, const vector<vector<T>>& b, T zero, Add add, Mul mul) {
    int n = sz(a), K = sz(b), m = sz(b[0]); vector<vector<T>> c(n, vector<T>(m, zero));
    for (int i = 0; i < n; i++) for (int k = 0; k < K; k++) for (int j = 0; j < m; j++) c[i][j] = add(c[i][j], mul(a[i][k], b[k][j]));
    return c;
}
 
ll graphPathsI(int n, const vector<pair<int, int>>& edges, ll k) {   // CSES 1723: walks 1 -> n with exactly k edges, mod p
    Mat a(n, vector<mint>(n, 0)); for (auto [u, v] : edges) a[u - 1][v - 1] += 1;
    return matPow(a, k)[0][n - 1].v;
}
ll graphPathsII(int n, const vector<array<ll, 3>>& edges, ll k) {   // CSES 1724: min weight walk 1 -> n with exactly k edges, -1 if none
    MinPlusMatrix a(n); for (auto [u, v, w] : edges) a.d[u - 1][v - 1] = min(a.d[u - 1][v - 1], w);
    ll r = a.pow(k).d[0][n - 1];
    return r >= MinPlusMatrix::INF ? -1 : r;
}

Graph Paths I and II (CSES). The DP $dp[t][v]$ = ways (resp. minimum cost) to be at $v$ after exactly $t$ edges has the same transition at every step, $dp[t+1][v]={\sum }_{u\to v}dp[t][u]$ (resp. ${\min }_{u\to v}(dp[t][u]+w)$), and $t$ goes up to $10^9$: the state is the vector $dp[t][\cdot ]$, the transition is the adjacency matrix in the $(+,\times )$ (resp. $(\min ,+)$) semiring, and the answer is entry $(1,n)$ of its $k$-th power. On the samples (edges $1\to 2,2\to 3,3\to 1,3\to 2$, $k=8$) the answers are 2 walks and minimum cost 27 (weights 5, 4, 1, 2).

Proof: semiring powers aggregate walks

Associativity of the semiring matrix product: $((AB)C{)}_{ij}={⨁}_l({⨁}_k{A}_{ik}\otimes B_{kl})\otimes C_{lj}={⨁}_{k,l}{A}_{ik}\otimes B_{kl}\otimes C_{lj}=(A(BC){)}_{ij}$, using distributivity on both sides, associativity of $\otimes$, and associativity and commutativity of $\oplus$ to regroup the finite sum. The identity matrix has $\bar{1}$ on the diagonal and $\bar{0}$ elsewhere ($\bar{0}$ absorbing makes the off-diagonal terms vanish). So binary exponentiation applies. Walks: by induction on $n$, $(A^n{)}_{ij}=⨁$ over walks $i=v_0,v_1,\dots ,v_n=j$ of $A_{v_0{v}_1}\otimes \cdots \otimes A_{v_{n-1}{v}_n}$: a walk of length $n+1$ is a walk of length $n$ to some $k$ followed by the edge $k\to j$, and distributivity gives $(A^{n+1}{)}_{ij}={⨁}_k(A^n{)}_{ik}\otimes A_{kj}$. Each row of the table is an instance: $(\min ,+)$ minimises total weight, $(\max ,\min )$ maximises the minimum edge, $(\vee ,\wedge )$ asks for existence. A 0-weight self-loop turns “exactly $n$” into “at most $n$” because padding a shorter walk with loops does not change its weight. The two CSES functions were checked on the samples and against step-by-step DP on random graphs.

Fast matrix multiplication

Algorithm	Exponent	In practice
schoolbook	$n^3$	the default; tuned loops reach several GFLOP/s
Strassen (1969)	${\log }_{2}7\approx 2.807$	crossover in the hundreds; numerically less stable for doubles; fine modulo $p$
Coppersmith–Winograd and successors (laser method)	$\omega <2.372$	galactic
boolean, word-RAM	$n^3/w$ with bitsets	best for every practical boolean size
lower bound	$n^2$ (reading the input)	$\omega =2$ is open

Strassen computes a $2\times 2$ block product with 7 multiplications instead of 8 (the “trade multiplications for additions” idea of Karatsuba):
$M_1=(A_{11}+A_{22})(B_{11}+B_{22})$, $M_2=(A_{21}+A_{22})B_{11}$, $M_3=A_{11}(B_{12}-B_{22})$, $M_4=A_{22}(B_{21}-B_{11})$, $M_5=(A_{11}+A_{12})B_{22}$, $M_6=(A_{21}-A_{11})(B_{11}+B_{12})$, $M_7=(A_{12}-A_{22})(B_{21}+B_{22})$;
$C_{11}=M_1+M_4-M_5+M_7$, $C_{12}=M_3+M_5$, $C_{21}=M_2+M_4$, $C_{22}=M_1-M_2+M_3+M_6$.

Mat strassen(const Mat& A, const Mat& B) {                 // square, size a power of two
    int n = sz(A);
    if (n <= 64) return A * B;
    int h = n / 2;
    auto blk = [&](const Mat& X, int r, int c) { Mat Y(h, vector<mint>(h)); for (int i = 0; i < h; i++) for (int j = 0; j < h; j++) Y[i][j] = X[r + i][c + j]; return Y; };
    auto add = [&](const Mat& X, const Mat& Y, int s) { Mat Z(h, vector<mint>(h)); for (int i = 0; i < h; i++) for (int j = 0; j < h; j++) Z[i][j] = s > 0 ? X[i][j] + Y[i][j] : X[i][j] - Y[i][j]; return Z; };
    Mat a11 = blk(A, 0, 0), a12 = blk(A, 0, h), a21 = blk(A, h, 0), a22 = blk(A, h, h);
    Mat b11 = blk(B, 0, 0), b12 = blk(B, 0, h), b21 = blk(B, h, 0), b22 = blk(B, h, h);
    Mat m1 = strassen(add(a11, a22, 1), add(b11, b22, 1)), m2 = strassen(add(a21, a22, 1), b11), m3 = strassen(a11, add(b12, b22, -1));
    Mat m4 = strassen(a22, add(b21, b11, -1)), m5 = strassen(add(a11, a12, 1), b22), m6 = strassen(add(a21, a11, -1), add(b11, b12, 1));
    Mat m7 = strassen(add(a12, a22, -1), add(b21, b22, 1));
    Mat C(n, vector<mint>(n));
    for (int i = 0; i < h; i++) for (int j = 0; j < h; j++) {
        C[i][j] = m1[i][j] + m4[i][j] - m5[i][j] + m7[i][j]; C[i][j + h] = m3[i][j] + m5[i][j];
        C[i + h][j] = m2[i][j] + m4[i][j]; C[i + h][j + h] = m1[i][j] - m2[i][j] + m3[i][j] + m6[i][j];
    }
    return C;
}

Proof: Strassen's identities and running time

Expand: $M_1+M_4-M_5+M_7=A_{11}{B}_{11}+A_{11}{B}_{22}+A_{22}{B}_{11}+A_{22}{B}_{22}+A_{22}{B}_{21}-A_{22}{B}_{11}-A_{11}{B}_{22}-A_{12}{B}_{22}+A_{12}{B}_{21}+A_{12}{B}_{22}-A_{22}{B}_{21}-A_{22}{B}_{22}=A_{11}{B}_{11}+A_{12}{B}_{21}$. $M_3+M_5=A_{11}{B}_{12}-A_{11}{B}_{22}+A_{11}{B}_{22}+A_{12}{B}_{22}=A_{11}{B}_{12}+A_{12}{B}_{22}$. $M_2+M_4=A_{21}{B}_{11}+A_{22}{B}_{11}+A_{22}{B}_{21}-A_{22}{B}_{11}=A_{21}{B}_{11}+A_{22}{B}_{21}$. $M_1-M_2+M_3+M_6=A_{11}{B}_{11}+A_{11}{B}_{22}+A_{22}{B}_{11}+A_{22}{B}_{22}-A_{21}{B}_{11}-A_{22}{B}_{11}+A_{11}{B}_{12}-A_{11}{B}_{22}+A_{21}{B}_{11}+A_{21}{B}_{12}-A_{11}{B}_{11}-A_{11}{B}_{12}=A_{21}{B}_{12}+A_{22}{B}_{22}$. The identities use only ring axioms (no commutativity), so they hold for blocks. $T(n)=7T(n/2)+O(n^2)$ gives $O(n^{{\log }_{2}7})$ by the master theorem. Checked against the cubic product for sizes up to 256.

Related notes

• Matrix Exponentiation DP (building transition matrices)
• Linear Recurrences (Kitamasa, Bostan–Mori, Berlekamp–Massey)
• Bitset (boolean products, transitive closure)