iv. Characteristic Polynomial

${\chi }_A(x)=\det (xI-A)=x^n+c_{n-1}{x}^{n-1}+\cdots +c_0$, with $c_{n-1}=-\mathrm{tr}A$ and $c_0=(-1{)}^n\det A$; its roots are the eigenvalues, so $\mathrm{tr}{A}^k=\sum {\lambda }_i^k$ (closed walks of length $k$ in a graph). Cayley–Hamilton: ${\chi }_A(A)=0$, so every power $A^k$ is a polynomial of degree $<n$ in $A$: this is the bridge from matrix powers to linear recurrences (any coordinate of $A^{k}v$ satisfies the recurrence with characteristic polynomial ${\chi }_A$). mint from Modular Arithmetic.

Frobenius (rational canonical) form is the formal name for this fact as a change of basis: $A$ is similar to a block-diagonal matrix of companion matrices, one per invariant factor of the minimal polynomial (a single block, the companion matrix of ${\chi }_A$ itself, whenever the minimal and characteristic polynomials agree — the generic case for a random matrix). A companion matrix’s powers acting on $e_1$ are exactly the coefficient-shift recurrence above, so computing $A^{k}v$ via “reduce $x^k\mathrm{mod}{\chi }_A(x)$ then apply to $v$” (Kitamasa / Bostan–Mori) is fast exponentiation in the Frobenius form, without ever constructing the change-of-basis matrix explicitly.

Method	Cost	Field requirement
Hessenberg reduction + recurrence	$O(n^3)$ ($n=500$ in about 0.2 s)	any field
Faddeev–LeVerrier	$O(n^4)$ (matrix products and traces)	divides by $\mathrm{1..}n$: characteristic 0 or $p>n$
interpolation of $\det (xI-A)$ at $n+1$ points	$O(n^4)$	$p>n$
Berlekamp–Massey on $u^T{A}^{k}v$ (random $u,v$)	$O(n^3)$ for $2n$ terms	finds the minimal polynomial w.h.p. (Wiedemann)

Use	How
$A^{k}v$ for huge $k$ and $n\approx 100$–$500$	$r(x)=x^k\mathrm{mod}{\chi }_A$ in $O(n^2\log k)$ (or $O(n\log n\log k)$ with NTT), then $A^{k}v=\sum r_i{A}^{i}v$ in $O(n^3)$: beats $O(n^3\log k)$
$\det (A+xB)$ as a polynomial	$B$ invertible: $\det B\cdot {\chi }_{-B^{-1}A}(x)$; in general interpolate $n+1$ values
inverse	Faddeev–LeVerrier byproduct $A^{-1}=-M_{n-1}/c_0$
closed walks of each length	$\mathrm{tr}{A}^k$ = Newton power sums of the roots

template <class T> vector<T> characteristicPolynomial(vector<vector<T>> a) {   // p[0..n], p[n] = 1, p[0] = (-1)^n det A
    int n = sz(a);
    for (int j = 0; j + 2 < n; j++) {                      // similarity transforms to upper Hessenberg form
        int piv = j + 1; while (piv < n && a[piv][j] == T(0)) piv++;
        if (piv == n) continue;
        if (piv != j + 1) { swap(a[piv], a[j + 1]); for (int i = 0; i < n; i++) swap(a[i][piv], a[i][j + 1]); }
        T inv = T(1) / a[j + 1][j];
        for (int i = j + 2; i < n; i++) {
            T c = a[i][j] * inv; if (c == T(0)) continue;
            for (int k = j; k < n; k++) a[i][k] -= c * a[j + 1][k];   // row i -= c row j+1
            for (int k = 0; k < n; k++) a[k][j + 1] += c * a[k][i];   // column j+1 += c column i
        }
    }
    vector<vector<T>> p(n + 1); p[0] = {T(1)};
    for (int i = 0; i < n; i++) {                          // p_{i+1} = (x - h_ii) p_i - sum_t h_ti (h_{t+1,t}..h_{i,i-1}) p_t
        p[i + 1].assign(i + 2, T(0));
        for (int k = 0; k <= i; k++) p[i + 1][k + 1] += p[i][k], p[i + 1][k] -= a[i][i] * p[i][k];
        T prod = 1;
        for (int t = i - 1; t >= 0; t--) {
            prod *= a[t + 1][t]; if (prod == T(0)) break;
            T c = prod * a[t][i]; for (int k = 0; k <= t; k++) p[i + 1][k] -= c * p[t][k];
        }
    }
    return p[n];
}
 
template <class T> vector<T> faddeevLeVerrier(const vector<vector<T>>& A, vector<vector<T>>* inverse = nullptr) {   // c[0..n] highest first: [1, c_{n-1}, ..., c_0]
    int n = sz(A); vector<T> c(n + 1, T(0)); c[0] = T(1);
    vector<vector<T>> M(n, vector<T>(n, T(0))), prev; for (int i = 0; i < n; i++) M[i][i] = T(1);
    for (int k = 1; k <= n; k++) {
        vector<vector<T>> AM(n, vector<T>(n, T(0)));
        for (int i = 0; i < n; i++) for (int l = 0; l < n; l++) for (int j = 0; j < n; j++) AM[i][j] += A[i][l] * M[l][j];
        T tr = T(0); for (int i = 0; i < n; i++) tr += AM[i][i];
        c[k] = -tr / T(k);
        prev = M; M = AM; for (int i = 0; i < n; i++) M[i][i] += c[k];   // M_k = A M_{k-1} + c I; M_n = 0 (Cayley-Hamilton)
    }
    if (inverse && n && c[n] != T(0)) { *inverse = prev; for (auto& r : *inverse) for (auto& x : r) x = -x / c[n]; }
    return c;
}
 
vector<mint> powerTimesVector(const vector<vector<mint>>& A, vector<mint> v, ll k) {   // A^k v in O(n^2 log k + n^3)
    int n = sz(A); if (n == 0) return v;
    vector<mint> chi = characteristicPolynomial(A);         // monic, degree n
    auto mulMod = [&](const vector<mint>& a, const vector<mint>& b) {   // a * b mod chi
        vector<mint> r(2 * n, 0);
        for (int i = 0; i < n; i++) if (a[i] != mint(0)) for (int j = 0; j < n; j++) r[i + j] += a[i] * b[j];
        for (int d = 2 * n - 1; d >= n; d--) if (r[d] != mint(0)) { mint f = r[d]; for (int t = 0; t <= n; t++) r[d - n + t] -= f * chi[t]; }
        r.resize(n); return r;
    };
    vector<mint> res(n, 0), base(n, 0); res[0] = 1;
    if (n == 1) base[0] = -chi[0]; else base[1] = 1;        // x mod chi
    for (; k; k >>= 1, base = mulMod(base, base)) if (k & 1) res = mulMod(res, base);
    vector<mint> out(n, 0), cur = v;                        // sum_i res_i A^i v
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) out[j] += res[i] * cur[j];
        vector<mint> nxt(n, 0); for (int r = 0; r < n; r++) for (int c = 0; c < n; c++) nxt[r] += A[r][c] * cur[c];
        cur = nxt;
    }
    return out;
}
 
vector<mint> detAffine(const vector<vector<mint>>& A, const vector<vector<mint>>& B) {   // coefficients of det(A + xB), O(n^4), p > n
    int n = sz(A); vector<mint> ys(n + 1), res(n + 1, 0);
    auto det = [&](vector<vector<mint>> a) {
        mint d = 1;
        for (int i = 0; i < n; i++) {
            int p = -1; for (int r = i; r < n; r++) if (a[r][i] != mint(0)) { p = r; break; }
            if (p < 0) return mint(0);
            if (p != i) swap(a[p], a[i]), d = -d;
            d *= a[i][i]; mint inv = a[i][i].inv();
            for (int r = i + 1; r < n; r++) { mint f = a[r][i] * inv; for (int c = i; c < n; c++) a[r][c] -= f * a[i][c]; }
        }
        return d;
    };
    for (int x = 0; x <= n; x++) { auto M = A; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) M[i][j] += mint(x) * B[i][j]; ys[x] = det(M); }
    for (int i = 0; i <= n; i++) {                          // Lagrange: prod_{j != i} (x - j) / (i - j)
        vector<mint> poly = {1}; mint denom = 1;
        for (int j = 0; j <= n; j++) if (j != i) {
            vector<mint> q(sz(poly) + 1, 0);
            for (int t = 0; t < sz(poly); t++) q[t + 1] += poly[t], q[t] -= mint(j) * poly[t];
            poly = q; denom *= mint(i - j);
        }
        mint f = ys[i] / denom; for (int t = 0; t <= n; t++) res[t] += f * poly[t];
    }
    return res;
}

Proof: Hessenberg reduction and its recurrence

Similarity. Row $i-=c\cdot$row $j+1$ is left multiplication by $E=I-c{e}_i{e}_{j+1}^T$, and column $j+1+=c\cdot$column $i$ is right multiplication by $E^{-1}=I+c{e}_i{e}_{j+1}^T$; swapping rows and the same columns is $PA{P}^{-1}$. Similar matrices have the same characteristic polynomial: $\det (xI-SA{S}^{-1})=\det (S(xI-A)S^{-1})=\det (xI-A)$. Choosing $c=a_{ij}/a_{j+1,j}$ clears $a_{ij}$ for $i\ge j+2$; the column operation only touches column $j+1$, so earlier cleared columns stay zero, and the result is upper Hessenberg (zero below the first subdiagonal). Recurrence. Let $H_i$ be the leading $i\times i$ block and $p_i=\det (xI-H_i)$. Expand $\det (xI-H_{i+1})$ along its last column, whose entries are $x-h_{ii}$ and $-h_{ti}$ ($t<i$). Deleting row $t$ and the last column leaves a block-triangular matrix with blocks $xI-H_t$ and an upper-triangular block whose diagonal is $-h_{t+1,t},\dots ,-h_{i,i-1}$. With the cofactor sign $(-1{)}^{t+i}$ the term is $(-h_{ti})(-1{)}^{t+i}{p}_t{\prod }_{s=t}^{i-1}(-h_{s+1,s})=-h_{ti}({\prod }_{s=t}^{i-1}{h}_{s+1,s})p_t$, which is the code’s update. Total cost $O(n^3)$ for the reduction and $O(n^3)$ for the recurrence.

Proof: Cayley–Hamilton, Faddeev–LeVerrier, powers and det(A + xB)

Cayley–Hamilton. $\mathrm{adj}(xI-A)(xI-A)={\chi }_A(x)I$, and $\mathrm{adj}(xI-A)={\sum }_{k<n}{B}_k{x}^k$ has matrix coefficients. Comparing coefficients of $x^k$ in ${\sum }_k{B}_k{x}^{k+1}-{\sum }_k{B}_{k}A{x}^k={\sum }_k{c}_k{x}^{k}I$ and multiplying the $x^k$ equation by $A^k$ on the right, the left sides telescope to 0, so ${\chi }_A(A)={\sum }_k{c}_k{A}^k=0$. Faddeev–LeVerrier. Define $M_0=I$, $M_k=A{M}_{k-1}+c_{n-k}I$, so $M_k={\sum }_{j=0}^k{c}_{n-j}{A}^{k-j}$ and $A{M}_{k-1}={\sum }_{j<k}{c}_{n-j}{A}^{k-j}$. Newton’s identities for the eigenvalues give $k{c}_{n-k}=-{\sum }_{j<k}{c}_{n-j}\mathrm{tr}{A}^{k-j}=-\mathrm{tr}(A{M}_{k-1})$, which is the update. Then $M_n={\chi }_A(A)=0$ (a free correctness check), so $A{M}_{n-1}=-c_{0}I$ and $A^{-1}=-M_{n-1}/c_0$. Dividing by $k\le n$ requires characteristic 0 or $p>n$; modulo a small prime use the Hessenberg method. ($\det A=(-1{)}^n{c}_0$ and the traces of all powers appear along the way.)

Powers. Write $x^k=q(x){\chi }_A(x)+r(x)$ with $\mathrm{deg}r<n$; substituting $A$ and using ${\chi }_A(A)=0$ gives $A^k=r(A)$, and $r$ is computed by binary exponentiation in the ring ${\mathbb{F}}_p[x]/({\chi }_A)$. $\det (A+xB)$ is a polynomial of degree $\le n$ in $x$ (Leibniz), so $n+1$ values determine it (Lagrange interpolation needs the points $\mathrm{0..}n$ distinct mod $p$); if $B$ is invertible, $\det (A+xB)=\det B\det (xI+B^{-1}A)=\det B{\chi }_{-B^{-1}A}(x)$. All functions were checked against brute-force determinants of $xI-A$ evaluated at many points, direct matrix powers and $A{A}^{-1}=I$.

Related notes

• Linear Recurrences (Kitamasa, Bostan–Mori, Berlekamp–Massey)
• Sparse and Structured Matrices (Wiedemann: minimal polynomial of sparse matrices)
• Determinants (Bareiss for exact values)