矩阵快速幂
矩阵
不严谨地讲，可以把一个存储数的二维数组当成一个矩阵。
矩阵加法：计算矩阵 \(C = A + B\)，每个 $c_{i,j} = a_{i,j} + b_{i,j} $，需要 \(A\) 矩阵和 \(B\) 矩阵行列数对应相等才能加。
矩阵数乘：计算矩阵 \(C=kA\) 或 \(C=Ak\)，每个 \(c_{i,j} = k a_{i,j}\)，其中 \(k\) 是数。
矩阵乘法：计算矩阵 \(C=AB\)，每个 \(c_{i,j} = \sum \limits_{k=1}^n a_{i,k} \cdot b_{kj}\)，三层循环，时间复杂度 \(O(n^3)\)。
两个矩阵能乘当且仅当 \(A\) 矩阵是 \(a\) 行 \(b\) 列的，\(B\) 矩阵是 \(b\) 行 \(c\) 列的，这样乘完会得到一个 \(a\) 行 \(c\) 列的矩阵。
矩阵乘法具有结合律 \((AB)C = A(BC)\)，但没有交换律（\(AB \ne BA\)），甚至都不一定能保证两个乘法都能做，行数列数可能都不满足能乘的条件，不过就算能乘，结果也不一定是相等的。
矩阵的幂：\(A^k = \overbrace{A \cdot A \cdots A}^\text{k 个}\)，需要是方阵才能做幂运算。
图的邻接矩阵的幂：设邻接矩阵为 \(A\)，则 \(A^k\) 包含了图上任意两点恰好经过 \(k\) 条边的路径数量。因为邻接矩阵 \(A^1\) 显然包含图上任意两点是否能直接相连。而 \({A^2}_{i,j} = \sum \limits_{k=1}^n A_{i,k} \times A_{k,j}\)，实际上就等于枚举中间点为 \(k\)，从 \(i\) 到 \(j\) 恰好经过 \(2\) 条边的路径数。同理，可以推广到 \(k\) 次方。
矩阵快速幂
因为矩阵乘法有结合律，所以计算一个 \(n\) 行 \(n\) 列的矩阵 \(A\) 的 \(b\) 次方时可以用快速幂加速。
比如可以把 \(A^{13} = A^1 A^4 A^8\)，然后用快速幂的写法做这个过程。
while (b > 0) {
if (b & 1) ans = mul(ans, A); // 这里是矩阵乘法
A = mul(A, A); // 这里是矩阵乘法，将A自己乘自己
b >>= 1;
}
时间复杂度 \(O(n^3 \log b)\)，其中 ans 矩阵初始应为单位矩阵（\(ans_{i,i}=1\)，其他地方为 \(0\)），这个矩阵乘任意矩阵 \(A\) 都会得到 \(A\) 矩阵本身。
例题：P3390 【模板】矩阵快速幂
解题思路
矩阵快速幂模板，中间计算乘加的时候注意取模，另外注意输入有负数。
矩阵乘法就通过三层循环计算。时间复杂度为 \(O(n^3 \log k)\)。
参考代码
#include <cstdio>
using ll = long long;
const int N = 105;
const int MOD = 1000000007;
int n, a[N][N], ans[N][N], c[N][N];
void mul(int a[N][N], int b[N][N]) {
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
c[i][j] = 0;
for (int k = 1; k <= n; k++) {
c[i][j] += 1ll * a[i][k] * b[k][j] % MOD;
c[i][j] %= MOD;
}
}
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
a[i][j] = c[i][j];
}
int main()
{
ll k; scanf("%d%lld", &n, &k);
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
scanf("%d", &a[i][j]);
if (a[i][j] < 0) a[i][j] = (a[i][j] + MOD) % MOD;
}
for (int i = 1; i <= n; i++) ans[i][i] = 1;
while (k > 0) {
if (k & 1) mul(ans, a);
mul(a, a);
k >>= 1;
}
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) printf("%d ", ans[i][j]);
printf("\n");
}
return 0;
}
也可以将矩阵封装成一种类型。可以重载 * 运算符，并使用三重循环进行矩阵乘法的操作，这样进行矩阵乘法就像计算整数乘法一样方便，增强可读性。
参考代码
#include <cstdio>
#include <vector>
#include <algorithm>
using ll = long long;
using std::vector;
using std::min;
const int MOD = 1000000007;
struct Matrix {
int n, m; // shape
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() { // I
int l = min(n, m);
for (int i = 0; i < l; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += 1ll * mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
int main()
{
int n; ll k; scanf("%d%lld", &n, &k);
Matrix a; a.set(n, n);
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++) {
scanf("%d", &a.mat[i][j]);
if (a.mat[i][j] < 0) a.mat[i][j] += MOD;
}
Matrix ans; ans.set(n, n); ans.eye();
while (k > 0) {
if (k & 1) ans = ans * a;
a = a * a;
k >>= 1;
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) printf("%d ", ans.mat[i][j]);
printf("\n");
}
return 0;
}
例题：P1962 斐波那契数列
由于 \(n\) 的数据范围极大，直接一项一项递推会超时。这个时候可以考虑使用矩阵加速递推过程。
考虑当 \(n \ge 3\) 时，每一项都和前面两项有关，所以考虑 \(\begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix}\) 和 \(\begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix}\) 的关系。
可以使用待定系数法，假设 \(\begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix} = \begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix} \begin{bmatrix} a & b \\ c & d \end{bmatrix}\)，则容易解出这个矩阵是 \(A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\)。
那么 \(\begin{bmatrix} F_n & F_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \end{bmatrix} A^{n-2}\)，矩阵乘法满足结合律，可以使用矩阵快速幂加速。
参考代码
#include <cstdio>
#include <vector>
#include <algorithm>
using ll = long long;
using std::vector;
using std::min;
const int MOD = 1000000007;
struct Matrix {
int n, m; // shape
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() { // I
int l = min(n, m);
for (int i = 0; i < l; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += 1ll * mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
Matrix qpow(Matrix a, ll k) {
Matrix res; res.set(a.n, a.n); res.eye();
while (k > 0) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
int main()
{
ll n; scanf("%lld", &n);
if (n < 3) {
printf("1\n");
} else {
n -= 2;
Matrix a; a.set(2, 2); a.mat[0][0] = a.mat[0][1] = a.mat[1][0] = 1;
a = qpow(a, n);
Matrix ans; ans.set(1, 2); ans.mat[0][0] = ans.mat[0][1] = 1;
ans = ans * a;
printf("%d\n", ans.mat[0][0]);
}
return 0;
}
习题：P1349 广义斐波那契数列
解题思路
当 \(n \ge 3\) 时，\(\begin{bmatrix} a_n & a_{n-1} \end{bmatrix} = \begin{bmatrix} a_2 & a_1 \end{bmatrix} A^{n-2}\)，其中 \(A = \begin{bmatrix} p & 1 \\ q & 0 \end{bmatrix}\)。
注意本题的模数是 int 范围内任意一个数，所以在矩阵乘法计算过程中可能溢出。
参考代码
#include <cstdio>
#include <vector>
#include <algorithm>
using ll = long long;
using std::vector;
using std::min;
int mod;
struct Matrix {
int n, m;
vector<vector<ll>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() {
for (int i = 0; i < n; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += mat[i][k] * other.mat[k][j] % mod;
res.mat[i][j] %= mod;
}
}
return res;
}
};
Matrix qpow(Matrix a, int k) {
Matrix res; res.set(a.n, a.m); res.eye();
while (k > 0) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
int main()
{
int p, q, a1, a2, n, m;
scanf("%d%d%d%d%d%d", &p, &q, &a1, &a2, &n, &m);
mod = m;
if (n == 1) {
printf("%d\n", a1 % m);
} else if (n == 2) {
printf("%d\n", a2 % m);
} else {
n -= 2;
Matrix a; a.set(2, 2);
a.mat[0][0] = p; a.mat[0][1] = 1; a.mat[1][0] = q;
a = qpow(a, n);
Matrix ans; ans.set(1, 2); ans.mat[0][0] = a2; ans.mat[0][1] = a1;
ans = ans * a;
printf("%lld\n", ans.mat[0][0]);
}
return 0;
}
习题：P5789 [TJOI2017] 可乐（数据加强版）
解题思路
可以将自爆看成一个节点，则在原始图的基础上 \(1 \sim n\) 每个点再向 \(0\) 号点连一条单向边，而停在原地相当于每个点加一个自环，设这个图的邻接矩阵为 \(A\)。
设 \(dp_{i,j}\) 表示经过 \(i\) 秒，落在城市 \(j\) 的方案数。那么 \(dp_{i,j} = \sum \limits_{A_{k,j} = 1} dp_{i-1,k}\)，而这个递推式可以用矩阵加速，即 \(\begin{bmatrix} dp_{i,0} & dp_{i,1} & \cdots & dp_{i,n} \end{bmatrix} = \begin{bmatrix} dp_{i-1,0} & dp_{i-1,1} & \cdots & dp_{i-1,n} \end{bmatrix} A\)，也就是 \(\begin{bmatrix} dp_{t,0} & dp_{t,1} & \cdots & dp_{t,n} \end{bmatrix} = \begin{bmatrix} 0 & 1 & \cdots & 0 \end{bmatrix} A^{t}\)，最终 \(\sum \limits_{i=0}^n dp_{t,i}\) 为总方案数。
注意，\(0\) 号节点需要有一个自环，因为如果没有这个自环，最终 \(dp_{t,0}\) 表示的是 \(t\) 秒后恰好自爆的方案数，但是根据题意，在 \(t\) 秒之前自爆的情况也要算进总的方案数中。
参考代码
#include <cstdio>
#include <vector>
using std::vector;
const int MOD = 2017;
struct Matrix {
int n, m;
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() {
for (int i = 0; i < n; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
Matrix qpow(Matrix a, int k) {
Matrix res; res.set(a.n, a.m); res.eye();
while (k > 0) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
int main()
{
int n, m; scanf("%d%d", &n, &m);
Matrix a; a.set(n + 1, n + 1); a.eye();
for (int i = 0; i <= n; i++) a.mat[i][0] = 1;
for (int i = 1; i <= m; i++) {
int u, v; scanf("%d%d", &u, &v);
a.mat[u][v] = a.mat[v][u] = 1;
}
int t; scanf("%d", &t);
a = qpow(a, t);
Matrix ans; ans.set(1, n + 1); ans.mat[0][1] = 1;
ans = ans * a;
int sum = 0;
for (int i = 0; i <= n; i++) sum = (sum + ans.mat[0][i]) % MOD;
printf("%d\n", sum);
return 0;
}
习题：P2579 [ZJOI2005] 沼泽鳄鱼
解题思路
如果没有食人鱼，那么这里就可以用邻接矩阵快速幂去加速方案数递推问题。而考虑食人鱼之后相当于原始邻接矩阵上有些边会用不了。
注意到题目中食人鱼的周期只有三种，分别为 \(2,3,4\)，因此第 \(t\) 个时间单位时可以转移的邻接矩阵的状态实际上和 \(t \bmod 12\) 个时间单位时是一样的。所以可以先预处理 \(12\) 个矩阵 \(A_0, A_1, A_2, \dots, A_{11}\) 分别表示对应单位时间下的邻接矩阵。
令矩阵 \(B = A_1 A_2 \cdots A_{11} A_0\)，则这个矩阵代表一个 \(12\) 的周期下要乘的矩阵，而 \(A_1 A_2 \cdots A_{k \bmod 12}\) 则是最后一个不满 \(12\) 的周期要乘的矩阵。
参考代码
#include <cstdio>
#include <vector>
using std::vector;
const int MOD = 10000;
struct Matrix {
int n, m;
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() {
for (int i = 0; i < n; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
Matrix qpow(Matrix a, int k) {
Matrix res; res.set(a.n, a.m); res.eye();
while (k > 0) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
int main()
{
int n, m, start, end, k;
scanf("%d%d%d%d%d", &n, &m, &start, &end, &k);
vector<Matrix> a(12);
for (int i = 0; i < 12; i++) a[i].set(n, n);
for (int i = 1; i <= m; i++) {
int x, y; scanf("%d%d", &x, &y);
for (int j = 0; j < 12; j++) a[j].mat[x][y] = a[j].mat[y][x] = 1;
}
int nfish; scanf("%d", &nfish);
for (int i = 1; i <= nfish; i++) {
int t; scanf("%d", &t);
for (int j = 0; j < t; j++) {
int p; scanf("%d", &p);
for (int w = j; w < 12; w += t) {
for (int u = 0; u < n; u++) {
a[w].mat[u][p] = 0;
}
}
}
}
Matrix b; b.set(n, n); b.eye();
for (int i = 1; i < 12; i++) b = b * a[i];
b = b * a[0];
b = qpow(b, k / 12);
Matrix ans; ans.set(1, n); ans.mat[0][start] = 1;
ans = ans * b;
for (int i = 1; i <= k % 12; i++) {
ans = ans * a[i];
}
printf("%d\n", ans.mat[0][end]);
return 0;
}
例题：P1357 花园
这是一个在环上进行满足特定约束的计数问题。
- 约束特性：问题的核心约束是“任意连续 \(m\) 个花圃中……”，这表明，在决定第 \(i\) 个位置种什么花时，只需要知道前面 \(m-1\) 个位置的种花情况，就可以判断新形成的 \(m\) 长窗口是否合法，这是一种典型的动态规划模式。
- 数据规模：\(n\) 的范围非常大（高达 \(10^{15}\)），常规的线性 DP（复杂度为 \(O(n \cdot \dots)\)）会超时。而 \(m\) 非常小（\(m \le 5\)），这意味着 DP 的状态数会很少。这种 \(n\) 巨大但状态数少的线性递推问题，是矩阵快速幂优化的经典应用场景。
- 环形处理：对于环形问题，一个常用的技巧是将其转化为线性问题求解。一个长度为 \(n\) 的环形排列可以看作是一个长度为 \(n\) 的线性排列，且其首尾部分能够“无缝”拼接而不违反约束。在基于状态转移的 DP 中，这对应于从某个状态出发，经过 \(n\) 步转移后，恰好又回到了自身状态，所有这类路径的总和就是环形方案的总数。
由于约束窗口大小为 \(m\)，决定下一个位置的关键信息是前面 \(m-1\) 个位置的模式。因此，可以用一个长度为 \(m-1\) 的二进制数 \(\text{mask}\) 作为 DP 状态，其中 1 代表 C，0 代表 P，总的状态数 \(S = 2^{m-1}\)。
可以构建一个 \(S \times S\) 大小的转移矩阵 \(T\)，其中 \(T_{j,i}\) 表示从状态 \(i\) 经过一步（即增加一个花圃）转移到状态 \(j\) 的方案数。
构建好转移矩阵 \(T\) 后，它描述了状态向量经过一步的变化，经过 \(n\) 步的变化就是 \(T^n\)。可以使用矩阵快速幂算法，在 \(O(S^3 \log n)\) 的时间内计算出 \(T^n\)。
如前所述，环形方案对应于从一个状态出发，经过 \(n\) 步转移后恰好回到自身状态的路径。在矩阵 \(T^n\) 中，对角线上的元素 \({T^n}_{i,i}\) 正好表示从状态 \(i\) 出发，经过 \(n\) 步后回到状态 \(i\) 的总路径数。
因此，最终答案就是 \(T^n\) 的迹，即主对角线上所有元素的和。
参考代码
#include <cstdio>
#include <vector>
using namespace std;
using ll = long long;
const ll MOD = 1e9 + 7;
ll n;
int m, k;
// 矩阵结构体，用于矩阵快速幂
struct Matrix {
int size;
vector<vector<ll>> mat;
Matrix(int s) : size(s), mat(s, vector<ll>(s, 0)) {}
// 静态成员函数，用于计算两个矩阵的乘积
static Matrix multiply(const Matrix& a, const Matrix& b) {
Matrix result(a.size);
for (int i = 0; i < a.size; ++i) {
for (int j = 0; j < a.size; ++j) {
for (int l = 0; l < a.size; ++l) {
result.mat[i][j] = (result.mat[i][j] + a.mat[i][l] * b.mat[l][j]) % MOD;
}
}
}
return result;
}
};
// 矩阵快速幂，用于计算 base^exp
Matrix matrix_pow(Matrix base, ll exp) {
Matrix result(base.size);
// 先创建一个单位矩阵作为初始结果
for (int i = 0; i < base.size; ++i) {
result.mat[i][i] = 1;
}
while (exp > 0) {
if (exp % 2 == 1) {
result = Matrix::multiply(result, base);
}
base = Matrix::multiply(base, base);
exp /= 2;
}
return result;
}
// 计算一个整数的二进制表示中 '1' 的个数
int popcount(int x) {
int count = 0;
while (x > 0) {
x &= (x - 1); // Brian Kernighan 算法
count++;
}
return count;
}
int main() {
scanf("%lld%d%d", &n, &m, &k);
// DP的状态是 m-1 位的二进制数，表示最后 m-1 个花圃的模式
// 状态总数 state_size = 2^(m-1)
int state_size = 1 << (m - 1);
Matrix transition(state_size);
// 构建转移矩阵 T
// T[j][i] 表示从状态 i 转移到状态 j 的方案数
for (int current_state = 0; current_state < state_size; ++current_state) {
// 尝试在当前 m-1 位状态后添加一个 'P' (0)
// 形成一个 m 位的窗口，其模式为 current_state 左移一位
int window_p = current_state << 1;
if (popcount(window_p) <= k) {
// 如果这个 m 位窗口合法，则计算转移到的下一个 m-1 位状态
int next_state_p = window_p & (state_size - 1); // 取低 m-1 位
transition.mat[next_state_p][current_state]++;
}
// 尝试在当前 m-1 位状态后添加一个 'C' (1)
int window_c = (current_state << 1) | 1;
if (popcount(window_c) <= k) {
// 如果这个 m 位窗口合法，计算转移到的下一个 m-1 位状态
int next_state_c = window_c & (state_size - 1); // 取低 m-1 位
transition.mat[next_state_c][current_state]++;
}
}
// 使用矩阵快速幂计算 T^n
Matrix result_matrix = matrix_pow(transition, n);
// 对于环形问题，答案是 T^n 的迹 (Trace)
// T^n[i][i] 表示从状态 i 出发，经过 n 步后恰好回到状态 i 的路径数
// 所有这些路径的和就是所有合法的环形方案数
ll ans = 0;
for (int i = 0; i < state_size; ++i) {
ans = (ans + result_matrix.mat[i][i]) % MOD;
}
printf("%lld\n", ans);
return 0;
}
例题：P1939 矩阵加速（数列）
这个递推过程是常规的，\(\begin{bmatrix} a_{n} & a_{n-1} & a_{n-2} \end{bmatrix} = \begin{bmatrix} a_3 & a_2 & a_1 \end{bmatrix} \begin{bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}^{n-3}\)。
参考代码
#include <cstdio>
#include <vector>
using std::vector;
const int MOD = 1000000007;
struct Matrix {
int n, m;
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() {
for (int i = 0; i < n; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += 1ll * mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
Matrix qpow(Matrix a, int k) {
Matrix res; res.set(a.n, a.m); res.eye();
while (k > 0) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
void solve() {
int n; scanf("%d", &n);
if (n <= 3) {
printf("1\n"); return;
}
n -= 3;
Matrix a; a.set(3, 3);
a.mat[0][0] = a.mat[0][1] = a.mat[1][2] = a.mat[2][0] = 1;
a = qpow(a, n);
Matrix ans; ans.set(1, 3); ans.mat[0][0] = ans.mat[0][1] = ans.mat[0][2] = 1;
ans = ans * a;
printf("%d\n", ans.mat[0][0]);
}
int main()
{
int t; scanf("%d", &t);
for (int i = 1; i <= t; i++) solve();
return 0;
}
分析一下这里的时间复杂度，假设递推式中用到的矩阵为 \(A\)，规模为 \(m \times m\)，对于本题来说 \(m=3\)，时间复杂度为 \(O(T m^3 \log n)\)。
但是实际上因为本题有 \(T\) 次询问，可以发现 \(A^2, A^4, A^8, \dots\) 可能在多次询问中被反复计算。
考虑预处理 \(A^1, A^2, A^4, A^8, \dots\) 的结果，则这样的预计算矩阵结果需要 \(O(\log n)\) 个。预处理的时间复杂度为 \(O(m^3 \log n)\)。
例如当询问 \(a_{16}\) 的结果时，需要计算 \(\begin{bmatrix} a_3 & a_2 & a_1 \end{bmatrix} A^{13} = \begin{bmatrix} a_3 & a_2 & a_1 \end{bmatrix} A^1 A^4 A^8\)。
从左往右计算，由于最左边的是一个 \(1 \times m\) 的矩阵（行向量），则它乘一个 \(m \times m\) 的矩阵的时间复杂度为 \(O(m^2)\)，因此如果持续这个过程，计算的时间复杂度为 \(O(m^2 \log n)\)。
这样处理整体的时间复杂度为 \(O(m^3 \log n + T m^2 \log n)\)，在多次询问的场景下优于一开始的计算方式。
参考代码
#include <cstdio>
#include <vector>
using std::vector;
const int MOD = 1000000007;
struct Matrix {
int n, m;
vector<vector<int>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
void eye() {
for (int i = 0; i < n; i++) mat[i][i] = 1;
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] += 1ll * mat[i][k] * other.mat[k][j] % MOD;
res.mat[i][j] %= MOD;
}
}
return res;
}
};
Matrix a[31];
void solve() {
int n; scanf("%d", &n);
if (n <= 3) {
printf("1\n"); return;
}
n -= 3;
Matrix ans; ans.set(1, 3);
ans.mat[0][0] = ans.mat[0][1] = ans.mat[0][2] = 1;
int i = 0;
while (n > 0) {
if (n & 1) ans = ans * a[i];
i++; n >>= 1;
}
printf("%d\n", ans.mat[0][0]);
}
void init() {
a[0].set(3, 3);
a[0].mat[0][0] = a[0].mat[0][1] = a[0].mat[1][2] = a[0].mat[2][0] = 1;
for (int i = 1; i < 31; i++) a[i] = a[i - 1] * a[i - 1];
}
int main()
{
init();
int t; scanf("%d", &t);
for (int i = 1; i <= t; i++) solve();
return 0;
}
习题：P6569 [NOI Online #3 提高组] 魔法值
解题思路
根据题意，定义两个矩阵的“异或”运算 \(C = A \oplus B\) 为 \(C_{i,j} = \bigoplus \limits_{k=1}^n A_{i,k} \times B_{k,j}\)。
设图的邻接矩阵为 \(A\)，则有 \(\begin{bmatrix} f_{1,j} & \cdots & f_{n,j} \end{bmatrix} = \begin{bmatrix} f_{1,j-1} & \cdots & f_{n,j-1} \end{bmatrix} \oplus A\)，所以 \(\begin{bmatrix} f_{1,a_i} & \cdots & f_{n,a_i} \end{bmatrix} = \begin{bmatrix} f_{1,0} & \cdots & f_{n,0} \end{bmatrix} \overbrace{\oplus A \oplus A \cdots \oplus A}^{a_i 个}\)。
如果要加速这个计算过程，需要证明这种“异或”运算有结合律。
对于本题中的运算，\(A \oplus B\) 的 \(B\) 因为是邻接矩阵，里面的取值只有 \(0,1\) 两种。
\(((A \oplus B) \oplus C)_{i,j} = \bigoplus \limits_{x=1}^n (\bigoplus \limits_{y=1}^n A_{i,y} \times B_{y,x}) \times C_{x,j}\)，而不管 \(C_{x,j}\) 等于 \(1\) 还是 \(0\)，都可以把 \(C_{x,j}\) 乘进去，即 \(((A \oplus B) \oplus C)_{i,j} = \bigoplus \limits_{x=1}^n \bigoplus \limits_{y=1}^n A_{i,y} \times B_{y,x} \times C_{x,j}\)。可以理解为对于 \(A_{i,y}\) 这一项，看 \(B\) 的第 \(y\) 行和 \(C\) 的第 \(j\) 列有多少个对应位置两者都是 \(1\)，如果有奇数个位置则相当于奇数个 \(A_{i,y}\) 相异或，也就是 \(A_{i,y}\) 要参与到最终的异或和式中，如果偶数个位置那么抵消掉为 \(0\)。
\((A \oplus (B \oplus C))_{i,j} = \bigoplus \limits_{y=1}^n A_{i,y} \times (\bigoplus \limits_{x=1}^n B_{y,x} \times C_{x,j})\)。如果 \(B\) 的第 \(y\) 行和 \(C\) 的第 \(j\) 列有奇数个对应位置两者都是 \(1\)，则 \(\bigoplus \limits_{x=1}^n B_{y,x} \times C_{x,j} = 1\)，偶数个则 \(\bigoplus \limits_{x=1}^n B_{y,x} \times C_{x,j} = 0\)。因此 \((A \oplus B) \oplus C = A \oplus (B \oplus C)\)。
因为有结合律，则可以对 \(a_i\) 次异或做二进制拆分，类似于 P1939 矩阵加速（数列），先预处理邻接矩阵的 \(1,2,4,8,\dots\) 次异或的结果，时间复杂度为 \(O(n^3 \log a_i)\)。针对每一次询问的 \(a_i\)，做二进制拆分，然后依次进行矩阵“异或”运算。由于每一次“异或”左侧的矩阵是 \(1 \times n\) 的，而“异或”右侧的矩阵是 \(n \times n\) 的，因此计算一次询问结果的时间复杂度为 \(O(n^2 \log a_i)\)。总的时间复杂度为 \(O(n^3 \log a_i + q n^2 \log a_i)\)。
参考代码
#include <cstdio>
#include <vector>
using ll = long long;
using std::vector;
struct Matrix {
int n, m;
vector<vector<ll>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m);
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = 0;
for (int k = 0; k < m; k++) {
res.mat[i][j] ^= mat[i][k] * other.mat[k][j];
}
}
return res;
}
};
Matrix f0, a[32];
void init() {
for (int i = 1; i < 32; i++) a[i] = a[i - 1] * a[i - 1];
}
void solve() {
ll k; scanf("%lld", &k);
Matrix ans = f0;
int i = 0;
while (k > 0) {
if (k & 1) ans = ans * a[i];
i++; k >>= 1;
}
printf("%lld\n", ans.mat[0][0]);
}
int main()
{
int n, m, q; scanf("%d%d%d", &n, &m, &q);
f0.set(1, n);
for (int i = 0; i < n; i++) scanf("%lld", &f0.mat[0][i]);
a[0].set(n, n);
for (int i = 1; i <= m; i++) {
int u, v; scanf("%d%d", &u, &v);
u--; v--;
a[0].mat[u][v] = a[0].mat[v][u] = 1;
}
init();
for (int i = 1; i <= q; i++) solve();
return 0;
}
习题：P6772 [NOI2020] 美食家
解题思路
测试点 \(1 \sim 8\)
设 \(dp_{i,j}\) 表示第 \(i\) 天时，到达城市 \(j\) 的最大愉悦值。则有 \(dp_{i,j} = \max \{ dp_{i-w, k] + c_j \}\)，其中 \(k\) 满足存在 \(k \rightarrow j\) 的边，其消耗的时间为 \(w\)。
而美食节则相当于对于计算过程中可达的状态 \(dp_{t, x}\)，可以额外获得 \(y\) 的愉悦值。
这个朴素的 \(DP\) 过程的时间复杂度为 \(O(Tm)\)。
测试点 \(9 \sim 10\)
输入的图是一个环，因此可以求出走一圈的时间，只有当 \(T\) 是一圈时间的倍数时才能回到城市 \(1\)，收获的愉悦值包括圈数乘上一圈的愉悦值总和，以及美食节产生的额外愉悦值（需要判断每个美食节举办城市及时间是否和行程吻合），注意还要加上刚开始在城市 \(1\) 时就获得的 \(c_1\) 愉悦值。
参考代码
#include <cstdio>
#include <algorithm>
#include <vector>
using ll = long long;
using std::vector;
using std::sort;
using std::max;
const int N = 55;
const int M = 505;
const int T = 52505;
const int K = 205;
int c[N], u[M], v[M], w[M];
ll dp[T][N];
struct Festival {
int t, x, y;
};
Festival f[K];
int main()
{
int n, m, t, k;
scanf("%d%d%d%d", &n, &m, &t, &k);
for (int i = 1; i <= n; i++) scanf("%d", &c[i]);
bool cycle = true;
for (int i = 1; i <= m; i++) {
scanf("%d%d%d", &u[i], &v[i], &w[i]);
if (v[i] != u[i] % n + 1) {
cycle = false;
}
}
for (int i = 1; i <= k; i++) scanf("%d%d%d", &f[i].t, &f[i].x, &f[i].y);
sort(f + 1, f + k + 1, [](Festival f1, Festival f2) {
return f1.t < f2.t;
});
if (cycle) {
vector<int> cost(n + 1), arr(n + 1);
for (int i = 1; i <= m; i++) cost[u[i]] = w[i];
ll sum = 0;
for (int i = 1; i <= n; i++) {
arr[i] = arr[i - 1] + cost[i - 1]; // 计算第一次到达城市i的时间
sum += c[i];
}
int cycle_time = arr[n] + cost[n];
ll bonus = 0;
for (int i = 1; i <= k; i++) {
if (f[i].t % cycle_time == arr[f[i].x])
bonus += f[i].y;
}
if (t % cycle_time == 0)
printf("%lld\n", 1ll * t / cycle_time * sum + bonus + c[1]);
else
printf("-1\n");
return 0;
}
dp[0][1] = c[1];
int idx = 1;
for (int i = 1; i <= t; i++) {
for (int j = 1; j <= m; j++) {
if (i - w[j] >= 0 && dp[i - w[j]][u[j]] > 0) {
dp[i][v[j]] = max(dp[i][v[j]], dp[i - w[j]][u[j]] + c[v[j]]);
}
}
while (idx <= k && f[idx].t == i) {
if (dp[i][f[idx].x] > 0) {
dp[i][f[idx].x] += f[idx].y;
}
idx++;
}
}
printf("%lld\n", dp[t][1] == 0 ? -1 : dp[t][1]);
return 0;
}
测试点 \(11 \sim 13\)
上面的 \(DP\) 递推过程瓶颈在于第一维时间 \(T\) 非常大，容易想到可能与矩阵加速递推有关。
分析状态转移方程，定义一种新的矩阵“乘法”运算：\(C = A \otimes B\)。\(C_{i,j} = \max \limits_{k=1}^n \{ A_{i,k} + B_{k,j} \}\)。利用邻接矩阵构造转移矩阵，可以设每个元素默认值为 \(-1\)，而有边的地方设为 \(c_j\)，\(A_{i,k} + B_{k,j}\) 参与计算需要满足 \(A_{i,k}\) 和 \(B_{k,j}\) 都不等于 \(-1\)。
可以证明这个新定义的“乘法”依然具有结合律。
先考虑 \(k=0\) 也就是没有美食节的情况。如果边权只有 \(w = 1\)，设上述方式构造的转移矩阵为 \(A\)，则 \(dp_{i, 1 \dots n}\) 可以由 \(dp_{i-1, 1 \dots n} \otimes A\) 得到。
但是边权的范围是 \(1 \le w \le 5\)，注意到 \(w\) 的范围不大，可以把每个原本的节点 \(u\) 拆成五个节点 \(u_1, u_2, u_3, u_4, u_5\)，此时原本输入的一条边权为 \(w\) 的边 \(u \rightarrow v\) 可以做等效处理，比如 \(w = 3\) 相当于 \(u_1 \rightarrow u_2 \rightarrow u_3 \rightarrow v_1\)，而这个新图上相当于边权全部看作 \(1\)，对应的转移矩阵的取值包括 \(-1\)（不连边），\(0\)（\(u_1 \rightarrow u_2 \rightarrow u_3 \rightarrow \cdots\)），\(c_v\)（\(u_w \rightarrow v\)）。这个转移矩阵的大小为 \(wn \times wn\)，因此矩阵加速递推的时间复杂度为 \(O(w^3n^3 \log T)\)。
测试点 \(14 \sim 20\)
处理多个美食节的情况。
美食节的存在相当于需要矩阵加速过程在中间暂停下来。例如有在中间两个美食节，分别先后发生在 \(t_1\) 和 \(t_2\) 时间，那么递推过程需要 \(dp_{1, 1 \dots n} \rightarrow dp_{t_1, 1 \dots n}\)，然后根据此时的转移结果判断是否需要加对应的额外愉悦值，然后继续递推 \(dp_{t_1, 1 \dots n} \rightarrow dp_{t_2, 1 \dots n}\)，再判断是否需要加 \(t_2\) 美食节产生的额外愉悦值，再完成最后到 \(T\) 天的递推过程。
也就是将 \(k\) 个美食节按时间先后依次处理，每次矩阵递推过程乘的是 \(A\) 的 \(\Delta t\) 次方，这里的实现类似于 P6569 [NOI Online #3 提高组] 魔法值，先预处理 \(A^1, A^2, A^4, A^8, \cdots\)，这样之后每次乘 \(A^{\Delta t}\) 时时间复杂度为 \(O(w^2n^2 \log \delta t)\)。整个预处理和矩阵递推过程总的时间复杂度为 \(O(w^3n^3\log T + (k + \log T) w^2 n^2)\)。
参考代码
#include <cstdio>
#include <algorithm>
#include <vector>
using ll = long long;
using std::sort;
using std::max;
using std::vector;
const int N = 55;
const int K = 205;
int c[N];
struct Festival {
int t, x, y;
};
Festival f[K];
struct Matrix {
int n, m;
vector<vector<ll>> mat;
void set(int x, int y) {
n = x; m = y;
mat.resize(n);
for (int i = 0; i < n; i++) mat[i].resize(m, -1);
}
Matrix operator*(const Matrix &other) const {
Matrix res; res.set(n, other.m);
for (int i = 0; i < n; i++)
for (int j = 0; j < other.m; j++) {
res.mat[i][j] = -1;
for (int k = 0; k < m; k++) {
if (mat[i][k] == -1 || other.mat[k][j] == -1) continue;
res.mat[i][j] = max(res.mat[i][j], mat[i][k] + other.mat[k][j]);
}
}
return res;
}
};
Matrix a[31];
void init() {
for (int i = 1; i < 31; i++) a[i] = a[i - 1] * a[i - 1];
}
int main()
{
int n, m, t, k;
scanf("%d%d%d%d", &n, &m, &t, &k);
for (int i = 1; i <= n; i++) scanf("%d", &c[i]);
a[0].set(n * 5, n * 5);
for (int i = 0; i < n; i++) {
for (int j = 0; j < 4; j++){
int u = j * n + i, v = u + n;
a[0].mat[u][v] = 0;
}
}
for (int i = 1; i <= m; i++) {
int u, v, w;
scanf("%d%d%d", &u, &v, &w);
u--; v--;
a[0].mat[u + (w - 1) * n][v] = c[v + 1];
}
init();
for (int i = 1; i <= k; i++) scanf("%d%d%d", &f[i].t, &f[i].x, &f[i].y);
sort(f + 1, f + k + 1, [](Festival f1, Festival f2) {
return f1.t < f2.t;
});
int ct = 0;
Matrix ans; ans.set(1, n * 5); ans.mat[0][0] = c[1];
f[k + 1].t = t;
for (int i = 1; i <= k + 1; i++) {
int j = 0, p = f[i].t - ct;
while (p > 0) {
if (p & 1) ans = ans * a[j];
j++; p >>= 1;
}
if (i <= k) {
if (ans.mat[0][f[i].x - 1] != -1)
ans.mat[0][f[i].x - 1] += f[i].y;
}
ct = f[i].t;
}
printf("%lld\n", ans.mat[0][0]);
return 0;
}
浙公网安备 33010602011771号