矩阵快速幂及其应用

矩阵快速幂广泛应用于将线性递推关系（如斐波那契数列）、图论中的路径计数、带固定状态的动态规划问题（如铺砖或字符串构造）等场景，转化为矩阵的高次幂运算，从而在 $ O (\log n)$ 时间内高效求解第 $ n $ 项或长度为 $ n $ 的方案数。

矩阵快速幂

由于矩阵乘法满足结合律，我们仍可类比标量快速幂的思想：只需将快速幂中的乘法运算替换为矩阵乘法，并将初始结果设为单位矩阵 $ I $，即可高效实现矩阵的快速幂运算。

矩阵操作封装模板

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

矩阵加速递推

如果有一道题目让你求斐波那契数列第 $ n $ 项的值，最简单的方法莫过于直接递推了。但是如果 $ n $ 的范围达到了 $ 10^{18}$ 级别，递推就不行了，此时我们可以考虑矩阵加速递推。

1. 斐波那契数列

链接：https://www.luogu.com.cn/problem/P1962

思路与推导：

在斐波那契数列当中，$ F_1 = F_2 = 1 $，$ F_n = F_{n − 1} + F_{n − 2}(i ≥ 3)$。

设 $\text {Fib}(n)$ 表示一个 $ 1 \times 2 $ 的矩阵 $\begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$。我们希望根据 $\text {Fib}(n-1) = \begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix}$ 推出 $\text {Fib}(n)$。

试推导一个矩阵 $\text {base}$，使得 $\text {Fib}(n-1) \times \text {base} = \text {Fib}(n)$，即 $\begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix} \times \text {base} = \begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$。

怎么推呢？因为 $ F_n = F_{n-1} + F_{n-2}$，所以 $\text {base}$ 矩阵的第一列应该是 $\begin {bmatrix} 1 \\ 1 \end {bmatrix}$，这样在进行矩阵乘法运算的时候才能令 $ F_{n-1}$ 与 $ F_{n-2}$ 相加，从而得出 $ F_n $。

此时得到 $\text {base} = \begin {bmatrix} 1 & x \\ 1 & y \end {bmatrix}$，代入原式得 $\begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix} \times \begin {bmatrix} 1 & x \\ 1 & y \end {bmatrix} = \begin {bmatrix} F_{n-1} + F_{n-2} & x \times F_{n-1} + y \times F_{n-2} \end {bmatrix} = \begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$，

解得 $ x = 1, y = 0 $，因此 $base$ 矩阵就为 $\begin {bmatrix} 1 & 1 \\ 1 & 0 \end {bmatrix}$。

综上所述：原式化为 $\begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix} \times \begin {bmatrix} 1 & 1 \\ 1 & 0 \end {bmatrix} = \begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$。转化为代码，应该怎么求呢？

定义初始矩阵 $\text {start} = \begin {bmatrix} F_2 & F_1 \end {bmatrix} = \begin {bmatrix} 1 & 1 \end {bmatrix}$， $\text {base} = \begin {bmatrix} 1 & 1 \\ 1 & 0 \end {bmatrix}$。那么， $ F_n $ 就等于 $\text {start} \times \text {base}^{n-2}$ 这个矩阵的第一行第一列元素，也就是 $\begin {bmatrix} 1 & 1 \end {bmatrix} \times \begin {bmatrix} 1 & 1 \\ 1 & 0 \end {bmatrix}^{n-2}$ 的第一行第一列元素。

注意：矩阵乘法不满足交换律，所以一定不能写成 $\begin {bmatrix} 1 & 1 \\ 1 & 0 \end {bmatrix}^{n-2} \times \begin {bmatrix} 1 & 1 \end {bmatrix}$ 的第一行第一列元素。

另外，对于 $ n \leq 2 $ 的情况，直接输出 1 即可，不需要执行矩阵快速幂。

为什么要乘上 $\text {base}$ 矩阵的 $ n-2 $ 次方而不是 $ n $ 次方呢？因为 $ F_1, F_2 $ 是不需要进行矩阵乘法就能求出的。也就是说，如果只进行一次乘法，就已经求出 $ F_3 $ 了。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    ll n;
    cin >> n;

    if (n <= 2)
    {
        cout << 1 << endl;
        return;
    }

    Matrix base(2, 2);
    base.clear();
    base.M[1][1] = base.M[2][1] = base.M[1][2] = 1;

    base = base ^ (n - 2);

    Matrix start(1, 2);
    start.M[1][1] = start.M[1][2] = 1;

    auto ans = start * base;

    cout << ans.M[1][1] << endl;
}

2. 广义斐波那契数列

链接：https://www.luogu.com.cn/problem/P1349

思路：

本题与斐波那契数列的区别为该题中 $ F_1 = a_1 $, $ F_2 = a_2 $，$ F_n = p \times F_{n − 1} + q \times F_{n − 2} (i ≥ 3)$。

定义 $base$ 为二阶方阵，由 $ F_n = p \times F_{n − 1} + q \times F_{n − 2}$，可知 $base$ 第一列为 $\begin {bmatrix} p \\ q \end {bmatrix}$。

此时得到 $\text {base} = \begin {bmatrix} p & x \\ q & y \end {bmatrix}$，代入原式得 $\begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix} \times \begin {bmatrix} p & x \\ q & y \end {bmatrix} = \begin {bmatrix} p \times F_{n-1} + q \times F_{n-2} & x \times F_{n-1} + y \times F_{n-2} \end {bmatrix} = \begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$，解得 $ x = 1, y = 0 $，因此 $base$ 矩阵就为 $\begin {bmatrix} p & 1 \\ q & 0 \end {bmatrix}$。

综上所述：原式化为 $\begin {bmatrix} F_{n-1} & F_{n-2} \end {bmatrix} \times \begin {bmatrix} p & 1 \\ q & 0 \end {bmatrix} = \begin {bmatrix} F_n & F_{n-1} \end {bmatrix}$。剩余步骤与普通斐波那契数列相同。

代码：

int MOD;

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    ll p, q, a1, a2, n, m;
    cin >> p >> q >> a1 >> a2 >> n >> m;

    MOD = m;

    if (n == 1 || n == 2)
    {
        if (n == 1)
            cout << a1 % MOD << endl;
        else
            cout << a2 % MOD << endl;

        return;
    }

    Matrix base(2, 2);
    base.M[1][1] = p;
    base.M[2][1] = q;
    base.M[1][2] = 1;
    base.M[2][2] = 0;

    base = base ^ (n - 2);

    Matrix start(1, 2);
    start.M[1][1] = a2;
    start.M[1][2] = a1;

    auto ans = start * base;

    cout << ans.M[1][1] << endl;
}

3. 矩阵加速（数列）

链接：https://www.luogu.com.cn/problem/P1939

思路：

由递推数列

$a_x = \begin{cases} 1 & x \in {1,2,3} \\ a_{x-1} + a_{x-3} & x \geq 4 \end{cases}$

可得 $base$ 矩阵第一列为 $\begin {bmatrix} 1 \\ 0 \\ 1 \end {bmatrix}$，此时得到 $\text {base} = \begin {bmatrix} 1 & a & d \\ 0 & b & e \\ 1 & c & f \end {bmatrix}$，

由 $\begin {bmatrix} a_{x-1} & a_{x-2} & a_{x-3} \end {bmatrix} \times \begin {bmatrix} 1 & a & d \\ 0 & b & e \\ 1 & c & f \end {bmatrix} = \begin {bmatrix} a_{x} & a_{x-1} & a_{x-2} \end {bmatrix}$，

可得 $\begin {bmatrix} a_{x-1} + a_{x-3} & a_{x-1} \times a + a_{x-2} \times b + a_{x-3} \times c & a_{x-1} \times d + a_{x-2} \times e + a_{x-3} \times f \end {bmatrix} = \begin {bmatrix} a_{x} & a_{x-1} & a_{x-2} \end {bmatrix}$

手动代入前几项得 $\text {base} = \begin {bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end {bmatrix}$，那么就有 $\begin {bmatrix} a_{x-1} & a_{x-2} & a_{x-3} \end {bmatrix} \times \begin {bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end {bmatrix} = \begin {bmatrix} a_{x} & a_{x-1} & a_{x-2} \end {bmatrix}$

定义初始矩阵 $\text {start} = \begin {bmatrix} x_3 & x_2 & x_1 \end {bmatrix} = \begin {bmatrix} 1 & 1 & 1\end {bmatrix}$。那么 $ x_n $ 就等于 $\text {start} \times \text {base}^{n-3}$ 这个矩阵的第一行第一列元素，也就是 $\begin {bmatrix} 1 & 1 & 1\end {bmatrix} \times \begin {bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end {bmatrix}^{n-3}$ 的第一行第一列元素。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    ll n;
    cin >> n;

    if (n <= 3)
    {
        cout << 1 << endl;
        return;
    }

    Matrix base(3, 3);
    base.clear();
    base.M[1][1] = base.M[1][2] = base.M[3][1] = base.M[2][3] = 1;

    base = base ^ (n - 3);

    Matrix start(1, 3);
    start.M[1][1] = start.M[1][2] = start.M[1][3] = 1;

    auto ans = start * base;

    cout << ans.M[1][1] << endl;
}

定长路径统计

模型

问题描述

给一个包含 $ n $ 个节点的有向图，每条边的边权均为 $1$，然后给一个整数 $ k $，你的任务是对于所有点对 ($ u, v $) 求出从 $ u $ 到 $ v $ 长度为 $ k $ 的路径的数量（不一定是简单路径，即路径上的点或者边可能走多次)。

解法

我们将这个图用邻接矩阵 $ G $ 表示（对于图中的边 $(u \to v)$，令 $ G [u,v] = 1 $，其余为 $0$ 的矩阵；如果有重边，则设 $ G [u,v]$ 为重边的数量），表示这个有向图。下述算法同样适用于图有自环的情况。

显然，该邻接矩阵对应 $ k = 1 $ 时的答案。

假设我们知道长度为 $ k $ 的路径条数构成的矩阵，记为矩阵 $ C_k $，我们想要求 $ C_{k+1}$。显然有 $DP$ 转移方程：$ C_{k+1}[i,j] = \sum_{p=1}^{n} C_k [i,p] \cdot G [p,j]$

我们可以把它看作矩阵乘法的运算，于是上述转移可以描述为：

$C_{k+1} = C_k \cdot G$

那么把这个递推式展开可以得到：$ C_k = \underbrace {G \cdot G \cdots G}_{k\text { 次}} = G^k $

直接使用矩阵快速幂，在 $ O (n^3 \log k)$ 的复杂度内计算结果。

1. 公交车路线

链接：https://www.luogu.com.cn/problem/P2233

问题描述：

在一个由 $8$ 个站点（$A–H$）组成的环上，从 $A$ 出发，每次只能移动到相邻站点，求恰好经过 $ n $ 次移动（即换乘 $ n $ 次车）后首次到达 $E$ 的路径数量，结果对 $1000$ 取模。

思路：

容易想到 $ G = \begin {bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end {bmatrix} = C_{1}$，那么答案就为 $ C_{n}[1, 5]$。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    int n;
    cin >> n;

    Matrix base(8, 8);
    base.clear();
    base.M[1][2] = base.M[1][8] = 1;
    base.M[2][1] = base.M[2][3] = 1;
    base.M[3][2] = base.M[3][4] = 1;
    base.M[4][3] = base.M[4][5] = 1;
    base.M[6][5] = base.M[6][7] = 1;
    base.M[7][6] = base.M[7][8] = 1;
    base.M[8][1] = base.M[8][7] = 1;

    auto ans = base ^ n;

    cout << ans.M[1][5] << endl;
}

2. 可乐

链接：https://www.luogu.com.cn/problem/P5789

问题描述：

给定一个 $ N $ 个点的无向图，可乐机器人从 $1$ 号点出发，每秒可选择：留在原地、移动到相邻点、或自爆（终止）；求 $ t $ 秒内所有可能的行为方案总数（包括中途自爆），结果对 $2017$ 取模。注意：一旦自爆就不再行动，且每秒必须执行一种行为。

思路：

考虑建图：

对于移动到相邻点，只需建一个无向邻接矩阵 $ G $；
对于留着原地，就是每个点都有一个自环；
对于自爆，定义一个超级终点 $ n + 1 $，可乐机器人在 $ u $ 点自爆，就是移动到了点 $ n + 1 $；
注意，超级终点 $ n + 1 $ 也要有条自环。

关于 $ n + 1 $ 点为什么也要有条自环：

在矩阵相乘中，如果不建自环，机器人不但不会停在那里，反而会直接凭空消失；

当 $ t = 2 $ 并选择在起点直接自爆时，在题目中路径为 $ 1 \to \text {爆炸 } (n + 1)$，但是真实的路径应该为 $ 1 \to \text {爆炸 } (n + 1) \to \text {爆炸 } (n + 1)$，虽然在时刻 $1$ 时已经自爆了，但是在真实的路径上仍然要考虑时刻 $2$ 的情况，因此 $ n + 1 $ 点也要有自环。

综上答案为 $\sum_{i = 1}^{n + 1} C_t [1, i] \text { mod 2017}$。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 110;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    int n, m, t;
    cin >> n >> m;
    vector a(n + 2, vector<int>(n + 2, 0));
    for (int i = 0; i < m; ++i)
    {
        int u, v;
        cin >> u >> v;
        a[u][v] = a[v][u] = 1;
    }

    for (int i = 1; i <= n + 1; ++i)
    {
        a[i][i] = 1;
        a[i][n + 1] = 1;
    }

    cin >> t;

    Matrix base(n + 1, n + 1);
    base.clear();
    for (int i = 1; i <= n + 1; ++i)
        for (int j = 1; j <= n + 1; ++j)
            base.M[i][j] = a[i][j];

    base = base ^ t;

    int ans = 0;
    for (int i = 1; i <= n + 1; ++i)
        ans = (ans + base.M[1][i]) % MOD;

    cout << ans << endl;
}

3. 迷路

链接：https://www.luogu.com.cn/problem/P4159

问题描述：

给定一个 $ n $ 个节点的有向图，每条边的权重是 $1$ 到 $9$ 的整数（字符形式给出），求从节点 $1$ 出发恰好用时 $ t $ 到达节点 $ n $ 的路径总数（不能停留，必须走完总时间 $ t $），结果对 $2009$ 取模。

思路：

由于权重不只为 $ 0,1 $，不能直接套模板。

考虑使用类似分层图的方法，既将点 $ u $ 拆分成 $<u, i>(i \in [1, 9])$，此时边权全变为 1。

建图：

由于每个点都拆分成了 $9$ 份，因此邻接矩阵大小也变成的原有的 $9$ 倍；
当存在边 $(u, v, w)$ 时，将 $<u, 1>$ 到 $<u, w>$ 以及中间节点全部连成链，$<u, w>$ 到 $<v, 1>$ 建边；
例：设存在边 $(u, v, 3)$，建 $<u, 1> \to <u, 2> \to <u, 3> \to <v, 1>$。

综上答案为 $ C_t [<1, 1>, <n, 1>]$。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 1;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                if (A.M[i][k] == 0)
                    continue;

                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = (ans.M[i][j] + A.M[i][k] * B.M[k][j]) % MOD;
            }
        }
        return ans;
    }

    // 重载 + 运算符 (两个矩阵维度完全一致)
    Matrix friend operator+(const Matrix &A, const Matrix &B)
    {
        Matrix ans(A.ROW, A.COL);

        for (int i = 1; i <= A.ROW; ++i)
            for (int j = 1; j <= A.COL; ++j)
                ans.M[i][j] = (A.M[i][j] + B.M[i][j]) % MOD;

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

void miyan()
{
    int n, t;
    cin >> n >> t;

    vector a(n + 1, vector<int>(n + 1, 0));
    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= n; ++j)
        {
            char c;
            cin >> c;
            int w = c - '0';

            a[i][j] = w;
        }
    }

    int idx = 1;
    map<pii, int> mp;
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= 9; ++j)
            mp[{i, j}] = idx++;

    vector b(n * 9 + 1, vector<int>(n * 9 + 1, 0));
    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= n; ++j)
        {
            if (!a[i][j])
                continue;

            int x = a[i][j];
            for (int k = 1; k < x; ++k)
                b[mp[{i, k}]][mp[{i, k + 1}]] = 1;

            b[mp[{i, x}]][mp[{j, 1}]] = 1;
        }
    }

    Matrix base(n * 9, n * 9);
    base.clear();
    for (int i = 1; i <= n * 9; ++i)
        for (int j = 1; j <= n * 9; ++j)
            base.M[i][j] = b[i][j];

    base = base ^ t;
    cout << base.M[mp[{1, 1}]][mp[{n, 1}]] << endl;
}

定长最短路

模型

问题描述

给一个包含 $ n $ 个节点的加权有向图和一个整数 $ k $。对于每个点对 ($ u, v $) 找到从 $ u $ 到 $ v $ 的恰好包含 $ k $ 条边的最短路的长度。（不一定是简单路径，即路径上的点或者边可能走多次)

解法

我们仍构造这个图的邻接矩阵 $ G $，其中 $ G [i,j]$ 表示从 $ i $ 到 $ j $ 的边权。如果 $ i, j $ 两点之间没有边，那么 $ G [i,j] = \infty $。（有重边的情况取边权的最小值）

显然上述矩阵对应 $ k = 1 $ 时问题的答案。我们仍假设我们知道 $ k $ 的答案，记为矩阵 $ L_k $。现在我们想求 $ k+1 $ 的答案。显然有转移方程：

$L_{k+1}[i,j] = \min_{1 \leq p \leq n} { L_k[i,p] + G[p,j] }$

事实上我们可以类比矩阵乘法，你发现上述转移只是把矩阵乘法的乘积求和变成相加取最小值，于是我们定义这个运算为 $\odot $，即

$A \odot B = C \quad \Longleftrightarrow \quad C[i,j] = \min_{1 \leq p \leq n} { A[i,p] + B[p,j] }$

于是得到

$L_{k+1} = L_k \odot G$

展开递推式得到

$ L_k = \underbrace {G \odot G \cdots G}_{k\text { 次}} = G^{\odot k}$

我们仍然可以用矩阵快速幂的方法计算上式，因为它显然是具有结合律的。时间复杂度 $ O (n^3 \log k)$。

模板

广义矩阵乘法（也叫 $(\min, +)$ 乘法）。

运算规则的对比

运算类型	算术乘法 ×	算术加法 +	单位元 (Identity)
普通矩阵乘法	乘法 (\times)	加法 (+)	1 (乘 1 不变，加 0 不变)
最短路矩阵乘法	加法 (+)	取最小值 (\min)	0 (加 0 不变，\min 无穷大不变)

单位元性质

单位矩阵 $ I $ 必须满足：对于任何矩阵 $ A $，都有 $ A \times I = A $。

让我们看广义乘法的计算公式：

$$(A \times I)_{i,j} = \min_{k=1}^{n} \{ A_{i,k} + I_{k,j} \}$$

为了让结果等于 $ A_{i,j}$，我们需要这个 $\min $ 序列中：

当 $ k = j $ 时：$ A_{i,k} + I_{k,j}$ 必须等于 $ A_{i,j}$。

这要求 $ I_{j,j} = 0 $（对角线为 $0$）。

当 $ k \neq j $ 时：$ A_{i,k} + I_{k,j}$ 必须大于或等于 $ A_{i,j}$，确保它不会影响 $\min $ 的结果。

为了在任何情况下都不影响最小值，最保险的值就是 $ I_{k,j} = \infty $（非对角线为无穷大）。

代码

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0x3f, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 0;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);
        ans.clear();

        for (int i = 1; i <= A.ROW; ++i)
            for (int k = 1; k <= A.COL; ++k)
                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = min(ans.M[i][j], A.M[i][k] + B.M[k][j]);

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

const ll inf = 0x3f3f3f3f3f3f3f3f;

1. Exactly K Steps 2

链接：https://atcoder.jp/contests/abc445/tasks/abc445_f

问题描述：

给定一个 $ N $ 个点的完全有向图（任意两点间都有单向边，边权为 $ C_{i,j}$），对每个起点 $ s $，求一条从 $ s $ 出发、恰好走 $ K $ 步后回到 $ s $ 的路径，使得总边权最小。

思路：

题目已经给出邻接矩阵 $ G $，对该矩阵跑快速幂得到 $ L_k $，那么最终答案就为对于每个 $ s = 1, 2, \dots, N $，都输出 $ L_k [s, s]$。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 100;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0x3f, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 0;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);
        ans.clear();

        for (int i = 1; i <= A.ROW; ++i)
            for (int k = 1; k <= A.COL; ++k)
                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = min(ans.M[i][j], A.M[i][k] + B.M[k][j]);

        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

const ll inf = 0x3f3f3f3f3f3f3f3f;

void miyan()
{
    int n, k;
    cin >> n >> k;
    vector a(n + 1, vector<ll>(n + 1, inf));
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            cin >> a[i][j];

    Matrix base(n, n);
    base.clear();
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= n; ++j)
            base.M[i][j] = a[i][j];

    base = base ^ k;

    for (int i = 1; i <= n; ++i)
        cout << base.M[i][i] << endl;
}

2. Cow Relays G

链接：https://www.luogu.com.cn/problem/P2886

问题描述：

在无向图中，求从起点 $ S $ 到终点 $ E $ 恰好经过 $ N $ 条边的最短路径长度（允许重复经过点和边）。

思路：

本题与 1 题的不同仅为只需输出 $ S $ 到 $ E $ 的最短路径长度并且需要自己建邻接矩阵。

代码：

// 矩阵大小，记得改初始矩阵大小
const int SIZE = 200;
struct Matrix
{
    int ROW, COL; // 行数, 列数
    ll M[SIZE + 1][SIZE + 1];

    // 构造函数
    Matrix(int _ROW = 0, int _COL = 0)
    {
        ROW = _ROW;
        COL = _COL;
        clear();
    }

    // 清零，生成零矩阵
    void clear() { memset(M, 0x3f, sizeof(M)); }

    // 初始化，生成大小为 n 的单位矩阵 (必须是方阵)
    void init()
    {
        clear();
        for (int i = 1; i <= ROW; ++i)
            M[i][i] = 0;
    }

    // 重载 * 运算符 (A 是 N * K 矩阵，B 是 K * M 矩阵)
    Matrix friend operator*(const Matrix &A, const Matrix &B)
    {
        // 前提：A 的列数必须等于 B 的行数
        // 结果：矩阵 ans 的大小为 A.ROW * B.COL
        Matrix ans(A.ROW, B.COL);
        ans.clear();

        for (int i = 1; i <= A.ROW; ++i)
        {
            for (int k = 1; k <= A.COL; ++k)
            {
                for (int j = 1; j <= B.COL; ++j)
                    ans.M[i][j] = min(ans.M[i][j], A.M[i][k] + B.M[k][j]);
            }
        }
        return ans;
    }

    // 重载 ^ 运算符，实现矩阵快速幂
    Matrix friend operator^(Matrix base, ll exp)
    {
        Matrix ans(base.ROW, base.ROW);
        ans.init(); // 生成与 base 同阶的单位矩阵

        while (exp)
        {
            if (exp & 1)
                ans = ans * base;
            base = base * base;
            exp >>= 1;
        }

        return ans;
    }
};

const ll inf = 0x3f3f3f3f3f3f3f3f;

void miyan()
{
    int n, m, s, e;
    cin >> n >> m >> s >> e;

    vector<array<int, 3>> edge(m);
    for (int i = 0; i < m; ++i)
    {
        int w, u, v;
        cin >> w >> u >> v;

        edge[i] = {u, v, w};
    }

    int idx = 1;
    map<int, int> mp;
    if (!mp.count(s))
        mp[s] = idx++;
    if (!mp.count(e))
        mp[e] = idx++;
    for (auto [u, v, w] : edge)
    {
        if (!mp.count(u))
            mp[u] = idx++;
        if (!mp.count(v))
            mp[v] = idx++;
    }

    vector a(idx, vector<ll>(idx, inf));
    for (auto [u, v, w] : edge)
        a[mp[u]][mp[v]] = a[mp[v]][mp[u]] = w;

    Matrix base(idx - 1, idx - 1);
    base.clear();
    for (int i = 1; i < idx; ++i)
        for (int j = 1; j < idx; ++j)
            base.M[i][j] = a[i][j];

    base = base ^ n;

    cout << base.M[mp[s]][mp[e]] << endl;
}

矩阵快速幂

矩阵操作封装模板

矩阵加速递推

1. 斐波那契数列

链接：https://www.luogu.com.cn/problem/P1962

思路与推导：

代码：

2. 广义斐波那契数列

链接：https://www.luogu.com.cn/problem/P1349

思路：

代码：

3. 矩阵加速（数列）

链接：https://www.luogu.com.cn/problem/P1939

思路：

代码：

定长路径统计

模型

问题描述

解法

1. 公交车路线

链接：https://www.luogu.com.cn/problem/P2233

问题描述：

思路：

代码：

2. 可乐

链接：https://www.luogu.com.cn/problem/P5789

问题描述：

思路：

代码：

3. 迷路

链接：https://www.luogu.com.cn/problem/P4159

问题描述：

思路：

代码：

定长最短路

模型

问题描述

解法

模板

运算规则的对比

单位元性质

代码

1. Exactly K Steps 2

链接：https://atcoder.jp/contests/abc445/tasks/abc445_f

问题描述：

思路：

代码：

2. Cow Relays G

链接：https://www.luogu.com.cn/problem/P2886

问题描述：

思路：

代码：

发送评论 编辑评论

推荐文章

发送评论编辑评论