矩阵快速幂加速递推
P3390 【模板】矩阵快速幂
矩阵乘法：
\(c_{i,j} = \sum a_{i,k} * b_{k,j}\)
矩阵快速幂：与普通快速幂类似。
特殊地，定义 \(A^0\) 为单位矩阵 \(I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}\)。
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
using i64 = long long;
const int N = 110, mod = 1e9 + 7;
struct matrix {
matrix() { memset(a, 0, sizeof(a)); }
int a[N][N];
};
int n;
i64 k;
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int k = 1; k <= n; k++)
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
return c;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
cin >> n >> k;
matrix ans, b;
for (int i = 1; i <= n; i++) ans.a[i][i] = 1;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
cin >> b.a[i][j];
while (k) {
if (k & 1) ans = ans * b;
b = b * b;
k >>= 1;
}
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) cout << ans.a[i][j] << ' ' ;
cout << '\n';
}
return 0;
}
POJ-3070/P1962 斐波那契数列
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
using i64 = long long;
const int N = 3, mod = 1e9 + 7;
struct matrix {
matrix() { memset(a, 0, sizeof(a)); }
int a[N][N];
};
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int i = 1; i < N; i++)
for (int j = 1; j < N; j++)
for (int k = 1; k < N; k++)
c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
return c;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
matrix ans, b;
i64 n;
cin >> n;
for (int i = 0; i < N; i++) ans.a[i][i] = 1;
b.a[1][1] = 1;
b.a[1][2] = 1;
b.a[2][1] = 1;
b.a[2][2] = 0;
while (n) {
if (n & 1) ans = ans * b;
b = b * b;
n >>= 1;
}
cout << ans.a[1][2] << '\n';
return 0;
}
P1939 【模板】矩阵加速（数列）
本题是斐波那契数列的进阶版。
因为 \(a_x\) 的计算涉及三项：\(a_{x-1}\)，\(a_{x-2}\)，\(a_{x-3}\)。
所以我们有
\(\begin{bmatrix} a_x & a_{x - 1} & a_{x - 2}\end{bmatrix} = \begin{bmatrix} a_{x - 1} & a_{x - 2} & a_{x - 3}\end{bmatrix} \times A\)
\(A\) 应为：
\(A = \begin{bmatrix} a & b & c\\d & e & f\\g & h & i\end{bmatrix}\)
按照套路，可以解得：
\(A = \begin{bmatrix} 1 & 0 & 1\\1 & 0 & 0\\0 & 1 & 0\end{bmatrix}\)
所以原问题转化为 \(\begin{bmatrix} a_{x + 2} & a_{x + 1} & a_x\end{bmatrix} = \begin{bmatrix} a_{3} & a_{2} & a_{1}\end{bmatrix} \times A^x\)
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
const int N = 4, mod = 1e9 + 7;
struct matrix {
matrix() { memset(a, 0, sizeof(a)); }
int a[N][N];
};
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int i = 1; i < N; i++)
for (int j = 1; j < N; j++)
for (int k = 1; k < N; k++)
c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
return c;
}
void solve() {
int n;
cin >> n;
matrix ans, b;
for (int i = 0; i < N; i++) ans.a[i][i] = 1;
b.a[1][1] = 1;
b.a[1][2] = 1;
b.a[2][3] = 1;
b.a[3][1] = 1;
while (n) {
if (n & 1) ans = ans * b;
b = b * b;
n >>= 1;
}
memset(b.a, 0, sizeof(b.a));
b.a[1][1] = 1;
b.a[1][2] = 1;
b.a[1][3] = 0;
ans = b * ans;
cout << ans.a[1][3] << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T;
cin >> T;
while (T--) solve();
return 0;
}
POJ 3233/UVA 11149 Power of Matrix
以UVA的题目为例，POJ上的题目就是把：solve()函数执行一遍即可。
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
const int N = 90, mod = 10;
struct matrix {
matrix() { memset(a, 0, sizeof(a)); }
int a[N][N];
};
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int i = 1; i < N; i++)
for (int j = 1; j < N; j++)
for (int k = 1; k < N; k++)
c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j]) % mod;
return c;
}
int n, k;
void solve() {
matrix ans, b;
for (int i = 1; i <= (n << 1); i++) ans.a[i][i] = 1;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
int x;
cin >> x;
x %= 10;
b.a[i][j] = x;
b.a[i][j + n] = 0;
b.a[i + n][j] = x;
if (i == j) b.a[i + n][j + n] = 1;
else b.a[i + n][j + n] = 0;
}
}
while (k) {
if (k & 1) ans = ans * b;
b = b * b;
k >>= 1;
}
for (int i = 1; i <= n; i++) {
for (int j = 1; j < n; j++) {
cout << ans.a[i + n][j] << ' ';
}
cout << ans.a[i + n][n] << '\n';
}
cout << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
while (cin >> n >> k, (n && k)) solve();
return 0;
}
P4910 帕秋莉的手环
吐槽一句：题面写得太好了！！！
建议看翻译。
本题先使用普通DP进行求解：
设 \(f[i][0/1]\) 为前 \(i - 1\) 位都已经选择好，且第 \(i\) 位染上了颜色 \(0/1\)的方案数。（设 \(1\) 为黑色，\(0\) 为白色。
于是有了：
\(f[i][1] = f[i - 1][0]\),
\(f[i][0] = f[i - 1][0] + f[i - 1][1]\)
而第一个珠子分为两种情况，黑的或白的，分别dp就可以了。
这种暴力dp可以拿 48 分，dalao们有的可以 AC。
我们通过运算，不难发现，答案为：
设 \(fib[i]\) 为斐波那契数列的第 \(i\) 项。
\(ans = fib[n - 1] \times 2 + fib[n]\)。
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
using i64 = long long;
const int N = 3, mod = 1e9 + 7;
struct matrix {
matrix() { memset(a, 0, sizeof(a)); }
int a[N][N];
};
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int i = 1; i < N; i++)
for (int j = 1; j < N; j++)
for (int k = 1; k < N; k++)
c.a[i][j] = (c.a[i][j] + 1ll * a.a[i][k] * b.a[k][j]) % mod;
return c;
}
void solve() {
i64 n;
cin >> n;
n--;
matrix a, res;
for (int i = 1; i < N; i++) res.a[i][i] = 1;
a.a[1][1] = 1;
a.a[1][2] = 1;
a.a[2][1] = 1;
a.a[2][2] = 0;
while (n) {
if (n & 1) res = res * a;
a = a * a;
n >>= 1;
}
cout << (res.a[1][2] * 2 % mod + res.a[1][1]) % mod << '\n';
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T;
cin >> T;
while (T--) solve();
return 0;
}
P2886 [USACO07NOV]Cow Relays G
题意：
给定一张 \(T\) 条边的无向连通图，求从 \(S\) 到 \(E\) 经过 \(N\) 条边的最短路长度。
用邻接矩阵 \(M\) 表示一个图，\(M(i,j)\) 是其第 \(i\) 行第 \(j\) 列元素，表示点 \(i\) 和 点 \(j\) 的直接关系，若点 \(i\) 和点 \(j\)直连，令 \(M[i, j]\) 等于 \(ij\) 的权值，否则等于无穷大（\(\infty\)）；当 \(i\) 等于 \(j\) 时令 \(M[i, i] = \infty\)。
定义矩阵的广义乘法：
\(M \times M = \min_{a = 1}^{N}[M(i, a) + M(a, j)]\)，
也就是把普通的矩阵乘法从求和改成了取最小值，把内部项相乘改为了相加。
定理： 计算邻接矩阵的广义幂 \(G = M ^ n\)，\(G(i, j)\) 的值是从 \(i\) 到 \(j\) 经过 \(n\) 条边（或 \(n - 1\) 个点）的最短路径长度。
以上摘自 《算法竞赛》。
注意：在做矩阵乘法时，一定要处理好上界(代码中为 \(idx\))。
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;
const int N = 1010, M = 110;
int n, m;
int S, T;
int ver[N], idx;
int add(int a) {
if (ver[a]) return ver[a];
ver[a] = ++idx;
return ver[a];
}
struct matrix {
matrix() { memset(a, 0x3f, sizeof(a)); }
int a[M][M];
};
matrix operator*(matrix& a, matrix& b) {
matrix c;
for (int k = 1; k <= idx; k++)
for (int i = 1; i <= idx; i++)
for (int j = 1; j <= idx; j++)
c.a[i][j] = min(c.a[i][j], a.a[i][k] + b.a[k][j]);
return c;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
matrix ans, g;
cin >> n >> m >> S >> T;
for (int i = 1; i <= m; i++)
{
int w, a, b, t1, t2;
cin >> w >> a >> b;
t1 = add(a);
t2 = add(b);
g.a[t1][t2] = g.a[t2][t1] = w;
}
n--;
ans = g;
while (n) {
if (n & 1) ans = ans * g;
g = g * g;
n >>= 1;
}
cout << ans.a[add(S)][add(T)] << '\n';
return 0;
}
[参考文献]
《算法竞赛》
浙公网安备 33010602011771号