「学习笔记」矩阵快速幂
一.快速幂
P1226 【模板】快速幂||取余运算
快速幂用来解决这样的式子 \(a^b\)，大家应该都会，我就放个代码。
ll ksm (ll a, ll b) {//a^b
ll res = 1;
while (b) {
if (b & 1) {
res = res * a;
}
b >>= 1;
a = a * a;
}
return res;
}
Ybtoj【例题1】序列的第k个数
这道题就可以直接根据前三项判断是等差数列还是等比数列。
等差数列 \(a_n=a_1+(n-1)×q\)，等比数列 \(a_n=a_1×q_{n-1}\)，\(q\) 为公差。
二.矩阵
概念1.单位矩阵
主对角线上的元素都为 \(1\)，其余元素都为 \(0\) 的 \(n\) 阶矩阵称为 \(n\) 阶单位矩阵。
单位矩阵的重要点在于：单位矩阵的任意次方都等于它本身。
如图，这就是一个三阶单位矩阵：
在矩阵快速幂转移时，可以将其设为初始矩阵。
概念2.矩阵加/减法
两个矩阵相加/减，就将其对应位置上的数相加/减。
设 \(A,B\) 为长宽相等的矩阵，矩阵 \(C=A±B\)，那么 \(A_{i,j}±B_{i,j}=C_{i,j}\)。
下图就是一个加法的例子。
概念3.矩阵乘法
设矩阵 \(C=A×B\)，有几个要求。
\(1.\) \(A\) 的列数等于 \(B\) 的行数。
\(2.\) \(A\) 是 \(n×r\) 的矩阵， \(B\) 是 \(r×m\) 的矩阵，那么 \(C\) 就是 \(n×m\) 的矩阵。
\(3.\) \(C_{i,j}=\sum ^{r}_{k=1}A_{i,k}B_{k,j}\)。
如图就是矩阵乘法的例子：
我们显然可以得到一个结论：
矩阵乘法满足结合律，不满足交换律。
这里给出 \(O(n^3)\) 的矩阵乘法代码：
inline void Mul (int f[][], int x[][], int y[][]) {//计算矩阵x乘y=f。
int t[N][N];//临时记录答案矩阵。
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j ++) {
t[i][j] = 0;//清空答案矩阵。
}
}
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j ++) {
for (int k = 0; k < n; k ++) {
t[i][j] = t[i][j] + x[i][k] * y[k][j];//矩阵乘法法则。
}
}
}
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j ++) {
f[i][j] = t[i][j];//复制答案。
}
}
}
概念4.方阵乘幂
方阵是一个长和宽相等的矩阵。
将方阵 \(A\) 自乘 \(n\) 次，记作 \(A^n\)。
上文说到矩阵乘法满足结合律，那么我们就可以将方阵乘幂利用快速幂的思想优化。
三.矩阵快速幂的引入
P3390 【模板】矩阵快速幂
先解决这道模板题。
给定矩阵 \(A\)，求 \(A^k\)。
只需要将快速幂中的数字变为矩阵即可，乘法变为矩阵乘法。
给出一个模板，这是以后做题的基础。
inline void qpow (int T) {
initans ();//这里是将ans初始化为单位矩阵。
while (T) {
if (T & 1) {
Mul (ans, ans, a);
//ans=ans*a。
}
Mul (a, a, a);
//a=a*a。
T >>= 1;
}
}
观察一下，是不是与原来的快速幂模板一模一样？
Ybtoj【例题2】斐波那契数列
P1962 斐波那契数列
求斐波那契数列的第 \(n\) 项 \(\bmod10^9+7，n<2^{63}\)。
\(O(n)\) 递推显然不可做，我们就需要采用矩阵快速幂来进行优化。
我们的思路在于：
我们将最后求的答案变为一个矩阵 \(B\)，构建一个初始矩阵 \(A\) 和一个转移矩阵 \(base\)，使 \(A×base=B\)。
那么，我们设答案矩阵 \(fib(n)=\begin{bmatrix} f_{n} & f_{n-1} \end{bmatrix}\)，我们希望通过 \(fib(n-1)=\begin{bmatrix} f_{n-1} & f_{n-2} \end{bmatrix}\) 得到 \(fib(n)\)，也就是说让 \(fib(n)=fib(n-1)×base\)。
由 \(f_n=f_{n-1}+f_{n-2}\) 得：
转移矩阵的第一列为 \(\begin{bmatrix} 1 \\ 1 \end{bmatrix}\)，因为 \(f_{n-1}×1+f_{n-2}×1=f_n\)。
那么可以推出转移矩阵的第二列为 \(\begin{bmatrix} 1 \\ 0 \end{bmatrix}\)，因为 \(f_{n-1}×1+f_{n-2}×0=f_{n-1}\)。
那么转移矩阵就是 \(\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\)。
因为 \(f_1=f_2=1\)，那么设初始矩阵 \(fib(2)=\begin{bmatrix} 1 & 1 \end{bmatrix}\)，乘一次 \(base\) 都能计算一项，也就是说我们要从 \(fib(2)\) 推到 \(fib(n)\) 需要乘 \(n-2\) 次 \(base\)。
也就说 \(A×base^{n-2}=B\)。
核心代码如下：
inline void qpow (int b) {
while (b) {
if (b & 1) {
Mul (f, f, a);///计算f=f×a。
}
Mul (a, a, a);//a=a^2。
b >>= 1;
}
}
最后输出 \(f_{0,0}\)，答案矩阵的第一位就是 \(f_n\)。
这里的矩阵快速幂仅仅是把快速幂中的数字改为了矩阵而已。
P1939 【模板】矩阵加速（数列）
设答案矩阵为 \(F(n)=\begin{bmatrix} a_{n} & a_{n-1} & a_{n-2} \end{bmatrix}\)，有一个矩阵 \(F(n-1)=\begin{bmatrix} a_{n-1}&a_{n-1} & a_{n-3} \end{bmatrix}\)，我们要构建矩阵 \(base\) 使 \(F(n-1)×base=F(n)\)。
然后可以推出矩阵 \(base=\begin{bmatrix} 1 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}\)。
套模板即可。
四.矩阵快速幂的运用
P1349 广义斐波那契数列
易得转移矩阵 \(base=\begin{bmatrix} p & 1 \\ q & 0 \end{bmatrix}\)，初始矩阵 \(A=\begin{bmatrix} a_{2} & a_{1} \end{bmatrix}\)，套矩阵快速幂模板转移 \(n-2\) 次即可。
P2044 [NOI2012] 随机数生成器
题目给出了递推式，发现要加上 \(c\)，我们的初始矩阵就设为 \([x_0,c]\)。
然后就根据递推式得出转移矩阵 \(\begin{bmatrix} a & 0 \\ 1 & 1 \end{bmatrix}\)，矩阵快速幂 \(n\) 次转移即可。
但是我们发现数据非常大，那么就要使用龟速乘。
Ybtoj【例题3】行为方案
P3758 [TJOI2017]可乐
P5789 [TJOI2017]可乐（数据加强版）
发现这个自爆操作不太好搞，可以创建一个虚拟点 \(0\)，自爆操作就相当于走到 \(0\) 点。
可以设出一个状态 \(f_{i,j,k}\) 为 \(i\) 走到 \(j\) 走 \(k\) 步的方案数，那么就有一个暴力方程 \(f_{i,j,k}=\sum\) \(f_{i,p,k-1} × f_{p,j,1}\)。
然后我们发现这个第三维可以压掉，那么 \(f_{i,j}=\sum f_{i,p}×f_{p,j}\)。
这样的复杂度是 \(O(tn^2)\) 的，无法通过怎么办？
发现这个转移方程非常符合矩阵快速幂的形式，我们可以用矩阵快速幂优化到 \(O(\log tn^2)\)。
最后的答案就是 \(\sum^{n}_{i=0} ans_{1,i}\)。
代码如下：
#include <bits/stdc++.h>
using namespace std;
const int N = 105;
const int mod = 2017;
int n, m, t, base[N][N], ans[N][N];
inline void Mul (int f[N][N], int x[N][N], int y[N][N]) {
int t[N][N];
for (int i = 0; i <= n; i ++) {
for (int j = 0; j <= n; j ++) {
t[i][j] = 0;
}
}
for (int i = 0; i <= n; i ++) {
for (int j = 0; j <= n; j ++) {
for (int k = 0; k <= n; k ++) {
t[i][j] = (t[i][j] + x[i][k] * y[k][j] % mod) % mod;
}
}
}
for (int i = 0; i <= n; i ++) {
for (int j = 0; j <= n; j ++) {
f[i][j] = t[i][j];
}
}
}
inline void qpow (int T) {
while (T) {
if (T & 1) {
Mul (ans, ans, base);
}
Mul (base, base, base);
T >>= 1;
}
}
int main () {
scanf ("%d%d", &n, &m);
for (int i = 1; i <= m; i ++) {
int u, v;
scanf ("%d%d", &u, &v);
base[u][v] = 1;
base[v][u] = 1;//有连边的方案为1。
}
for (int i = 0; i <= n; i ++) {
base[i][i] = 1;//自己到自己方案为1。
base[i][0] = 1;//到爆炸点方案为1。
ans[i][i] = 1;//初始矩阵设为单位矩阵。
}
scanf ("%d", &t);
qpow (t);//因为初始的时刻为0，所以总共要转移t次。
int sum = 0;
for (int i = 0; i <= n; i ++) {
sum = (sum + ans[1][i]) % mod;
}
printf ("%d", sum);
return 0;
}
代码中初始矩阵和答案矩阵是同一个矩阵，因为我在快速幂时不断覆盖初始矩阵。
初始矩阵因为是没操作的情况，那么就设为单位矩阵。
Ybtoj【例题4】矩阵求和
这题我就稍微讲一下思路。（因为我也不太会做）
我们构造矩阵 \(B=\begin{bmatrix} A & I \\ 0 & I \end{bmatrix}\)，其中 \(A\) 为给出的初始矩阵， \(I\) 为 \(n×n\) 的单位矩阵， \(0\) 为 \(n×n\) 的零矩阵。
那么 \(B^{k+1}=\begin{bmatrix} A^{k+1} & I+\Sigma _{i=1}^{k}A^{i} \\ 0 & 1 \end{bmatrix}\)，就发现 \(B^{k+1}\) 右上角的矩阵减去一个单位矩阵就是我们所要求的的答案。
Ybtoj【例题5】最短路径
P2886 [USACO07NOV]Cow Relays G
我们发现，有 \(100\) 条边，但有 \(1000\) 个点，那我们先离散化一下。
抛开边数限制，我们可以采用 \(Floyd\) 算法。
图的邻接矩阵存储本质上就是一个矩阵，那我们可以想到用矩阵快速幂优化。
由于矩阵乘法满足结合律，我们直接转化，将乘法变为取 \(\min\)。
每次乘幂都是相当于走了一条边，那么乘 \(n\) 次转移矩阵就相当于走了 \(n\) 条边。
模板转化即可。
inline void Mul (int f[N][N], int x[N][N], int y[N][N]) {
int t[N][N];
for (int i = 0; i < tot; i ++) {
for (int j = 0; j < tot; j ++) {
t[i][j] = inf;//floyd算法初始为inf。
}
}
for (int i = 0; i < tot; i ++) {
for (int j = 0; j < tot; j ++) {
for (int k = 0; k < tot; k ++) {
t[i][j] = min (t[i][j], x[i][k] + y[k][j]);//将乘法改为min。
}
}
}
for (int i = 0; i < tot; i ++) {
for (int j = 0; j < tot; j ++) {
f[i][j] = t[i][j];
}
}
}
多说一句，由于 \(Ybtoj\) 上的数据有问题，所以有一个点需要数据特判才能 \(AC\)。
P2233 [HNOI2002]公交车路线
这种题一看就是线性递推，那么就容易想到用矩阵快速幂进行优化。
阅读题面先写出几个信息：
\(1.\) \(n\) 为偶数时方案数为 \(0\)。
\(2.\) \(n<4\) 时方案数为 \(0\)。
那么我们先特判 \(n\)，然后把 \(n\) 除以二，这样每一个数字都是相邻的，便于处理。
我们可以马上秒一个 \(dfs\) 爆搜，得到前几位：
\(f_1=0,f_2=2,f_3=8,f_4=28\)。
猜测 \(f_i\) 与 \(f_{i-1}\) 和 \(f_{i-2}\) 有关系，那我们列方程组，解得递推式如下：
\(f_i=4×f_{i-1}-2×f_{i-2}\)。
设答案矩阵为 \(F(n)\)，那么这样构建初始矩阵 \(F(2)=[0, 2]\)，转移矩阵 \(base=\begin{bmatrix} 4 & 1 \\ -2 & 0 \end{bmatrix}\)。
然后转移 \(n-2\) 次就可以了。
答案要乘 \(2\) ，因为我们只计算走向一边的方案。
又因为转移矩阵中有负数，所以我们要在答案输出时防止出现负数，像这样操作：
cout << (2 * ans[0][0] % mod + mod) % mod;
完整代码不放了就套模板即可。
