矩阵加速图上问题学习笔记
Verdandi
2020-08-05 22:07:19
# 矩阵加速图上问题学习笔记
[更好的阅读体验](https://zybuluo.com/xiaoziyao/note/1730580)
## 前置知识
- 矩阵乘法与矩阵快速幂
- 简单的图上dp
- 一颗愿意学习的大脑
- 基本的码力
## 引入
我们看一道题：[P2233 [HNOI2002]公交车路线](https://www.luogu.com.cn/problem/P2233)
题意：有$8$个车站$A$~$H$围成一圈，相邻的车站可以互相到达（比如$A$和$H$，$C$和$D$），从车站$A$出发换$n$次车到达车站$E$有多少种方案，注意，到达车站$E$后将不会继续行动。
数据范围：$1\leqslant n\leqslant 10^7$。
分析：
这道题是一道图上dp裸题，我们可以定义$tot_i$为到达车站$i$的方案数，并枚举换车次数$t$来更新它：每一次更新的时候枚举出发车站与到达车站，只要这两个车站可以相互到达，且出发车站不是$E$，那么我们就对到达车站的$tot$进行更新。
核心代码：
```
inline int abs(int x){
return x<0? -x:x;
}
inline int check(int a,int b){//判断能否到达
return (a==1&&b==8)||(a==8&&b==1)||abs(a-b)==1;
}
tot[1]=1;//从车站A出发
for(t=1;t<=n;t++){
for(i=1;i<=8;i++)
tmp[i]=0;
for(i=1;i<=8;i++)//枚举出发车站
for(j=1;j<=8;j++)//枚举到达车站
if(i!=5&&check(i,j))
tmp[j]=(tmp[j]+tot[i])%mod;
for(i=1;i<=8;i++)
tot[i]=tmp[i];
}
```
但是这份代码交上去后会T，获得$80pts$的好成绩。
啊这……
怎么办呢？
我们可以尝试给这份代码换一个形式，再看看有什么特殊的地方：
```
a[1][2]=a[1][8]=1;//a为可以相互到达点组成的邻接矩阵
a[2][3]=a[2][1]=1;
a[3][2]=a[3][4]=1;
a[4][3]=a[4][5]=1;
//这里没有a[5][4]=a[5][6]=1的原因是到了车站E就不能继续行动了
a[6][5]=a[6][7]=1;
a[7][6]=a[7][8]=1;
a[8][7]=a[8][1]=1;
tot[1]=1;
for(t=1;t<=n;t++){
for(i=1;i<=8;i++)
tmp[i]=0;
for(i=1;i<=8;i++)
for(j=1;j<=8;j++)
if(a[j][i])
tmp[i]=(tmp[i]+tot[j])%mod;
for(i=1;i<=8;i++)
tot[i]=tmp[i];
}
```
诶诶诶诶诶，好像发现了什么不得了的东西。
让我们再换个形式，让它更加明了：
```
a[1][2]=a[1][8]=1;
a[2][3]=a[2][1]=1;
a[3][2]=a[3][4]=1;
a[4][3]=a[4][5]=1;
a[6][5]=a[6][7]=1;
a[7][6]=a[7][8]=1;
a[8][7]=a[8][1]=1;
//注意，这里tmp和tot多加了一维
tot[1][1]=1;
for(t=1;t<=n;t++){
for(i=1;i<=1;i++)//这个i的值始终为1，不用管
for(j=1;j<=8;j++)
tmp[i][j]=0;
for(i=1;i<=1;i++)//这个i的值始终为1，不用管
for(j=1;j<=8;j++)
for(k=1;k<=8;k++)
tmp[i][j]=(tmp[i][j]+tot[i][k]*a[k][j]%mod)%mod;//其实就相当于a[k][j]等于1就更新tot
for(i=1;i<=8;i++)
tot[1][i]=tmp[1][i];
}
```
我们对比一下矩阵乘法的代码：
```
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=(res.mt[i][j]+a.mt[i][k]*b.mt[k][j]%mod)%mod;
return res;
}
```
！！！
加上$n$比较大的数据范围，你想到了什么？
$O(\log n)$！矩阵快速幂！
我们直接重构代码，上矩阵快速幂：
```
#include<stdio.h>
const int maxk=10,mod=1000;
int i,j,k,m,n;
struct matrix{
int len,wid;
int mt[maxk][maxk];
};
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=(res.mt[i][j]+a.mt[i][k]*b.mt[k][j]%mod)%mod;
return res;
}
int calc(int b){
matrix a,res;
a.len=a.wid=8;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=a.wid;j++)
a.mt[i][j]=0;
a.mt[1][2]=a.mt[1][8]=1;
a.mt[2][3]=a.mt[2][1]=1;
a.mt[3][2]=a.mt[3][4]=1;
a.mt[4][3]=a.mt[4][5]=1;
a.mt[6][5]=a.mt[6][7]=1;
a.mt[7][6]=a.mt[7][8]=1;
a.mt[8][7]=a.mt[8][1]=1;
res.len=1,res.wid=8;
res.mt[1][1]=1,res.mt[1][2]=res.mt[1][3]=res.mt[1][4]=res.mt[1][5]=res.mt[1][6]=res.mt[1][7]=res.mt[1][8]=0;
while(b){
if(b&1)
res=mul(res,a);
a=mul(a,a),b>>=1;
}
return res.mt[1][5];
}
int main(){
scanf("%d",&n);
printf("%d\n",calc(n));
return 0;
}
```
然后就AC了！
复杂度：$O(\log n)$，带一个$8^2$的大常数。
## 正确性证明
那，为什么这是对的呢？
我们考虑上面构造的矩阵$a$，看它的实际含义：
$$
\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}
$$
我们发现$a_{i,j}$的含义是除了车站$E$，车站$i$是否可以到达车站$j$，等价于除了车站$E$，车站$i$经过一条边（即换乘一次）是否能到达车站$j$。因为里面的数由$0$和$1$组成，所以等价于即除了车站$E$，车站$i$经过一条边到达车站$j$的方案数。
那把它自乘一次（即平方一次）呢？
$$
\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}^2=
\begin{bmatrix}
2&0&1&0&0&0&1&0\\
0&2&0&1&0&0&0&1\\
1&0&2&0&1&0&0&0\\
0&1&0&1&0&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&1&0&1\\
1&0&0&0&1&0&2&0\\
0&1&0&0&0&1&0&2
\end{bmatrix}
$$
经过平方后，$a_{i,j}^2$其意义其实就是除了车站$E$，车站$i$经过二条边到达车站$j$的方案数。
原因是矩阵平方时的代码是这样的：
```
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=(res.mt[i][j]+a.mt[i][k]*b.mt[k][j]%mod)%mod;
return res;
}
```
即有$a^2_{i,j}=\sum_{k=1}^n a_{i,k}\cdot a_{k,j}$
这可以看为一个变相的Floyd（矩阵加速Floyd会在后面进行讨论）：枚举起点$i$与终点$j$，枚举决策点$k$，使用乘法原理，利用$i$到$k$和$k$到$j$的路径数量更新$i$到$j$的路径数量。
再看看$a^3$，它有$a^3_{i,j}=\sum_{x=1}^n a_{i,x}\cdot(\sum_{y=1}^n a_{x,y}\cdot a_{y,j})=\sum_{x=1}^n\sum_{y=1}^n a_{i,x}\cdot a_{x,y}\cdot a_{y,j}$，刚好是每一条从$i$到$j$的长度为$3$路径数量相乘之和。
以此类推，我们就知道这个矩阵的$k$次方代表的是除了车站$E$，每个车站经过$k$条边到达其他车站的方案数。
这刚好与我们的题目不谋而合！
**推广到一般情况，给定邻接矩阵$a$，$a$的$k$次幂代表每个点经过$k$条边到达其他车站的方案数**，这样，我们就可以用矩阵快速幂实现许多黑科技！（这个结论可以用归纳法证明，~~留为作业~~）
## 例题1：[P4159 [SCOI2009] 迷路](https://www.luogu.com.cn/problem/P4159)（拆点，普通矩阵加速）
现在，我们就可以做一道题目来实践刚刚的结论了！
题意：给定一个带边权的无向图，求从点$1$到点$n$，长度为$t$的路径数量。
数据范围：$1\leqslant n\leqslant 10$，$1\leqslant t\leqslant10^9$，所有边权为$1$到$9$内的整数。
分析：
这道题如果去掉边权就是我们上面讨论的结论了，现在考虑怎么加上边权。
可以拆点，将每个点拆做$9$个状态（$1$到$9$），第$k$个状态的$i$点用$(i,k)$表示，代表还要经过$i-1$个长度的边才能到达点$i$的状态，这样就好办了。
- 对于原图中的边$(i,j,k)$，$(i,0)$连向$(j,k-1)$（从点$i$出发时必须到达了点$i$，所以是$(i,0)$，在转移的时候默认边长为$1$，所以$k$要减一）
- 对于每个点$i$，它所有状态需要相互连边，$(i,k)$连向$(i,k-1)$（转移时走过一个长度的边）
这样，我们就可以直接上板子了：
```
#include<stdio.h>
#include<iostream>
using namespace std;
const int maxn=105,mod=2009;
int i,j,k,m,n;
string s;
struct matrix{
int len,wid;
int mt[maxn][maxn];
}st;
inline int getnum(int x,int k){
return (k-1)*n+x;
}
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=(res.mt[i][j]+a.mt[i][k]*b.mt[k][j]%mod)%mod;
return res;
}
int calc(int b){
matrix a=st,res;
res.len=1,res.wid=9*n;
res.mt[1][getnum(1,1)]=1;
for(int i=2;i<=9*n;i++)
res.mt[1][i]=0;
while(b){
if(b&1)
res=mul(res,a);
a=mul(a,a),b>>=1;
}
return res.mt[1][getnum(n,1)];
}
int main(){
scanf("%d%d",&n,&m);
st.len=st.wid=9*n;
for(i=1;i<=n;i++){
cin>>s;
for(j=0;j<n;j++)
if(s[j]!='0')
st.mt[getnum(i,1)][getnum(j+1,s[j]-48)]=1;
}
for(i=1;i<=n;i++)
for(j=1;j<=8;j++)
st.mt[getnum(i,j+1)][getnum(i,j)]=1;
printf("%d\n",calc(m));
return 0;
}
```
## 例题2：[P2151 [SDOI2009]HH去散步](https://www.luogu.com.cn/problem/P2151)（点边互换，普通矩阵加速）
同样，我们再做一道比较新颖的题目来巩固一下：[P2151 [SDOI2009]HH去散步](https://www.luogu.com.cn/problem/P2151)
这道题运用了点边互换的小trick，分析详见[P2151 [SDOI2009]HH去散步解题报告](https://zybuluo.com/xiaoziyao/note/1730724)。
## 例题3：[CF576D Flights for Regular Customers](https://www.luogu.com.cn/problem/CF576D)（bitset优化矩阵乘法）
有时候，$O(n^3\log k)$的复杂度不足以我们通过题目，我们以这一题为例，来了解bitset对矩阵乘法的优化。
题意：给定一张$n$个点$m$条边的有向图，起点为$1$，终点为$n$，第$i$条边可以通行当且仅当已经走过$d_i$条边，求起点到终点最少要走多少条边，若无法到达输出“Impossible”。
数据范围：$1\leqslant n,m\leqslant 150$，$1\leqslant d_i\leqslant 10^9$。
分析：
这道题比较经典，值得一做。
考虑边的解锁顺序，先套路地对边进行排序（按$d$值从小到达）。然后每次添加进一条边，判断是否能到达，如果能到达则求答案最小值，否则输出“Impossible”。
考虑怎么添加边：假设已经枚举到了第$i$条边，我们维护了一个矩阵$g$表示当前的邻接矩阵（$g$现在还没有更新），那么一定有矩阵$now=g^{d_i-d_{i-1}}$，表示在解锁第$i-1$条边后，走了$d_i-d_{i-1}$条边后，能到达的点的邻接矩阵。然后我们加入第$i$条边，更新$g$，并用Floyd求$dis$数组，其中$dis_{i,j}$表示加入第$i$条边后的最短路。
那么更新答案也比较简单了：因为我们需要经过$d_i$条边才能解锁第$i$条边，所以$d_i$的边数是必须要的；在解锁了第$i$条边后，如果我们在$x$点，那么我们到达终点还需要的边数就是$dis_{x,n}$，已经用Floyd求出来了，因此我们答案就可以求出来了（注意，$x$点可以更新答案的原因是可以在$d_i-d_{i-1}$的步数内由上一个$now$转移来，所以需要判断一下是否满足$now_{1,x}=1$）：$ans=\min_{i=1}^m\{\min_{j=1,now_{1,j}=1}^n\{d_i+dis_{j,n}\}\}$。
但即使我们使用了矩阵快速幂，我们的复杂度也是$O(n^3\cdot m\cdot \log dmax)$的，不足以通过本题（其中$dmax$指$d$的最大值）。
因为我们复杂度的瓶颈在于求$now$矩阵上，因此我们考虑如何优化矩阵快速幂，也就是矩阵乘法的复杂度。因为$now$只需要判断是否可以到达，所以只需要存储$01$数据就可以了。什么东西处理$01$数据最快呢？bitset！因此我们可以用bitset维护矩阵，复杂度$O(\frac{n^3\cdot m\cdot\log dmax}{\omega})$（注，如果不知道bitset可以看[扶苏的bitset浅谈](https://www.luogu.com.cn/blog/fusu2333/fu-su-di-bitset-qian-tan)进行学习）
代码：
```
#include<stdio.h>
#include<bitset>
#include<algorithm>
#define inf 0x3f3f3f3f
using namespace std;
const int maxn=155,maxm=155;
int i,j,k,m,n,t,ans;
int dis[maxn][maxn];
struct matrix{
int len,wid;
bitset<maxn>mt[maxn];
}now,g;
struct edge{
int x,y,d;
}e[maxm];
inline int cmp(edge a,edge b){
return a.d<b.d;
}
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int k=1;k<=a.wid;k++)
for(int i=1;i<=a.len;i++)
if(a.mt[i][k])
res.mt[i]|=b.mt[k];
return res;
}
matrix calc(matrix a,int b,matrix res){
while(b){
if(b&1)
res=mul(res,a);
a=mul(a,a),b>>=1;
}
return res;
}
int main(){
scanf("%d%d",&n,&m);
for(i=1;i<=m;i++)
scanf("%d%d%d",&e[i].x,&e[i].y,&e[i].d);
sort(e+1,e+1+m,cmp);
ans=inf;
now.len=n,now.wid=n;
for(i=1;i<=n;i++)
now.mt[i][i]=1;
g.len=n,g.wid=n;
for(t=1;t<=m;t++){
now=calc(g,e[t].d-e[t-1].d,now);
g.mt[e[t].x][e[t].y]=1;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++){
if(g.mt[i][j]==0)
dis[i][j]=i==j? 0:inf;
else dis[i][j]=g.mt[i][j];
}
for(k=1;k<=n;k++)
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
dis[i][j]=min(dis[i][j],dis[i][k]+dis[k][j]);
for(i=1;i<=n;i++)
if(now.mt[1][i])
ans=min(ans,e[t].d+dis[i][n]);
}
if(ans==inf)
puts("Impossible");
else printf("%d\n",ans);
return 0;
}
```
## 例题4：[P2579 [ZJOI2005]沼泽鳄鱼](https://www.luogu.com.cn/problem/P2579)（分组矩阵加速）
有时候，单个的矩阵已经无法满足我们表示状态的需求了，此时我们需要处理多个矩阵，并同样地利用这些矩阵进行加速，我们以[P2579 [ZJOI2005]沼泽鳄鱼](https://www.luogu.com.cn/problem/P2579)为例。
题意：给定一个$n$个点$m$条边的无向图，其中会有一些食人鱼按照$p1\rightarrow p2\rightarrow p1\rightarrow\cdots$，或$p1\rightarrow p2\rightarrow p3\rightarrow p1\rightarrow\cdots$，或$p1\rightarrow p2\rightarrow p3\rightarrow p4\rightarrow p1\rightarrow\cdots$的顺序进行行动（每次行动时间需要$1$单位时间），每条边需要花费$1$单位时间，求从$s$到$t$，长度为$k$，且到的每个点都没有食人鱼的路径数量。
数据范围：$1\leqslant n\leqslant 50$，$1\leqslant k\leqslant 2\cdot10^9$，食人鱼的数量不超过$20$。
分析：
如果不考虑食人鱼，这就是一个显然的矩阵快速幂题，但是有了食人鱼的限制，我们就需要把一些点禁掉，不可以到某些点。
但是，食人鱼会动啊？这怎么处理呢？
其实比较简单，我们发现食人鱼的运动周期比较小，只有$2,3,4$，这样，我们就可以枚举每一个状态（最多$12=lcm(2,3,4)$种状态），处理这个状态时的邻接矩阵，然后判断一下$k$的大小：如果小于$12$就暴力乘，如果大于等于$12$就计算有几个循环，用$12$个状态相乘得到的矩阵来矩阵快速幂，余数暴力乘。
有一点需要注意的是，由于**矩阵乘法不满足交换律**，所以我们在把$12$个状态相乘的时候，需要按$2,3\cdots12,1$的顺序**顺次**相乘（因为时间$1$还没有出发）。
代码：
```
#include<stdio.h>
const int maxn=55,mod=10000;
int i,j,k,l,m,n,s,t,p,fish;
struct matrix{
int len,wid;
int mt[maxn][maxn];
}r[13];
matrix mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=0;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=(res.mt[i][j]+a.mt[i][k]*b.mt[k][j]%mod)%mod;
return res;
}
int calc(int x,int y,int b){
int rst=b%12;
matrix a,res;
res.len=1,res.wid=n;
for(int i=1;i<=n;i++)
res.mt[1][i]=0;
res.mt[1][x]=1;
if(b<12){
for(int i=2;i<=b+1;i++)
res=mul(res,r[i]);
return res.mt[1][y];
}
a=r[2];
for(int i=3;i<=12;i++)
a=mul(a,r[i]);
a=mul(a,r[1]);
b/=12;
while(b){
if(b&1)
res=mul(res,a);
a=mul(a,a),b>>=1;
}
for(int i=2;i<=rst+1;i++)
res=mul(res,r[i]);
return res.mt[1][y];
}
int main(){
scanf("%d%d%d%d%d",&n,&m,&s,&t,&p);
s++,t++;
for(i=1;i<=12;i++)
r[i].len=n,r[i].wid=n;
for(i=1;i<=m;i++){
int x,y;
scanf("%d%d",&x,&y);
x++,y++;
for(j=1;j<=12;j++)
r[j].mt[x][y]=r[j].mt[y][x]=1;
}
scanf("%d",&fish);
for(i=1;i<=fish;i++){
int x,y;
scanf("%d",&x);
for(j=1;j<=x;j++){
scanf("%d",&y);
y++;
for(k=j;k<=12;k+=x)
for(l=1;l<=n;l++)
r[k].mt[l][y]=0;
}
}
printf("%d\n",calc(s,t,p));
return 0;
}
```
## 例题5：[P2886 [USACO07NOV]Cow Relays G](https://www.luogu.com.cn/problem/P2886)（矩阵乘法重定义）
上文说过，会对矩阵加速Floyd进行讨论，我们便以此题为例来探讨**矩阵乘法重定义后对Floyd的加速**。
我们以[这道题](https://www.luogu.com.cn/problem/P2886)为例子，来了解如何对矩阵乘法重定义，并证明矩阵加速的正确性。
题意：给定一张$m$条边的无向图，求$s$到$e$经过$k$条边的最短路径长度。
数据范围：$2\leqslant m\leqslant 100$，$1\leqslant k\leqslant 10^6$。
分析：
我们可以先想一想这道题的暴力解法，虽然复杂度不能通过本题，但可以对正解产生一定的启发：
我们可以直接建图，设$dis_{i,j}$然后做$k$遍Floyd，每次Floyd用更新最短路长度
（注，文中的Floyd可能枚举顺序与大部分的Floyd顺序不一样，但是仍然是正确的）
```
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
dis[i][j]=i==j? 0:inf;
for(p=1;p<=q;p++){
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
tmp[i][j]=inf;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
for(k=1;k<=n;k++)
tmp[i][j]=min(tmp[i][j],dis[i][k]+g[k][j]);
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
dis[i][j]=tmp[i][j];
}
```
（暴力比正解还难写）
其实，Floyd中核心内容就一句话：
$$dis_{i,j}'=\min_{k=1}^n (dis_{i,k}+g_{k,j})$$
而矩阵乘法的核心是这样：
$$c_{i,j}=\sum_{k=1}^n(a_{i,k}\cdot b_{k,j})$$
这两个式子惊人的相似，我们考虑是否可以把矩阵乘法重定义，实现Floyd的转移，同时也要满足结合律（为了让矩阵快速幂可以对其加速）。
我们定义矩阵乘法为：
$$c_{i,j}=\min_{k=1}^n(a_{i,k}+b_{k,j})$$
我们证明一下这个矩阵乘法满足结合律：
设$D=(A\times B)\times C$，那么有$D_{i,j}=\min_{b=1}^n(\min_{a=1}^n(A_{i,a}+B_{a,b})+C_{b,j})$，因为外层的$\min$要让$A_{i,a}+B_{a,b}+C_{b,j}$最小，里层的$\min$要让$A_{i,a}+B_{a,b}$最小，所以可以直接拆开括号，得到$D_{i,j}=\min_{a=1}^n\min_{b=1}^n(A_{i,a}+B_{a,b}+C_{b,j})$。
设$E=A\times(B\times C)$，同理有$E=\min_{a=1}^n\min_{b=1}^n(A_{i,a}+B_{a,b}+C_{b,j})$。
因此，我们证明了这个矩阵乘法满足结合律，这样就可以利用矩阵快速幂对邻接矩阵的Floyd转移进行加速了！
时间复杂度是$O(n^3\log k)$，可以通过本题。
代码：
```
#include<stdio.h>
#include<map>
#define inf 1000000000
using namespace std;
const int maxn=105;
int i,j,k,m,n,s,t,tot;
map<int,int>mp;
struct matrix{
int len,wid;
int mt[maxn][maxn];
}st;
matrix _mul(matrix a,matrix b){
matrix res;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=res.len;i++)
for(int j=1;j<=res.wid;j++)
res.mt[i][j]=inf;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]=min(res.mt[i][j],a.mt[i][k]+b.mt[k][j]);
return res;
}
int calc(int x,int y,int b){
matrix a=st,res;
res.len=1,res.wid=n;
for(int i=1;i<=n;i++)
res.mt[1][i]=inf;
res.mt[1][x]=0;
while(b){
if(b&1)
res=_mul(res,a);
a=_mul(a,a),b>>=1;
}
return res.mt[1][y];
}
int main(){
scanf("%d%d%d%d",&k,&m,&s,&t);
for(i=1;i<=m;i++){
int x,y,z;
scanf("%d%d%d",&z,&x,&y);
if(mp.count(x)==0)
mp[x]=++tot;
if(mp.count(y)==0)
mp[y]=++tot;
x=mp[x],y=mp[y];
st.mt[x][y]=st.mt[y][x]=z;
}
s=mp[s],t=mp[t],n=tot;
st.len=n,st.wid=n;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
if(st.mt[i][j]==0)
st.mt[i][j]=inf;
printf("%d\n",calc(s,t,k));
return 0;
}
```
## 例题6：[P6569 [NOI Online #3 提高组]魔法值](https://www.luogu.com.cn/problem/P6569)
在这里，我们来根据[P6569 [NOI Online #3 提高组]魔法值](https://www.luogu.com.cn/problem/P6569)谈谈矩阵乘法的倍增预处理与二进制拆分优化。
题意：给定一个$n$个点，$m$条边的无向图（带点权），每一天，每个节点的点权会异或上与它相邻的所有点昨天的点权，$q$次询问，每次询问节点$1$在第$k$天的点权。
看到这道题目，很显然我们可以重定义一下矩阵乘法：
$$c_{i,j}=\oplus_{k=1}^n(a_{i,k}\times b_{k,j})$$
然后，我们可以让$a$矩阵为原点权矩阵（大小$1\times n$），$b$矩阵为邻接矩阵（大小$n\times n$），此时$c$矩阵就是下一天的点权。
然后我们每一次询问做一次矩阵快速幂可以做到$O(qn^3\log k)$的复杂度。
但是这个复杂度不能通过本题，需要考虑优化。
我们发现做矩阵快速幂的时候我们的矩阵大小都是$n\times n$的，十分浪费，而我们最后的答案只需要求节点$1$的点权，应该可以省去很多不需要求的信息。
我们预处理所有$2^p$次幂的邻接矩阵，然后在询问的时候对于$k$进行二进制拆分，乘上所有对应的邻接矩阵。
但这样的复杂度为什么是对的呢？因为询问的时候进行的矩阵乘法是$1\times n$与$n\times n$矩阵的乘法，而根据矩阵乘法的过程可以知道它的复杂度为$O(n^2)$。
这样，时间复杂度就优化到了：$O(n^3\log k+qn^2\log k)$。
```
#include<stdio.h>
#include<string.h>
#define int long long
const int maxn=105,maxk=40;
int n,m,q;
int s[maxn];
struct matrix{
int len,wid;
int mt[maxn][maxn];
}zero,ans[maxk];
matrix mul(matrix a,matrix b){
matrix res=zero;
res.len=a.len,res.wid=b.wid;
for(int i=1;i<=a.len;i++)
for(int j=1;j<=b.wid;j++)
for(int k=1;k<=a.wid;k++)
res.mt[i][j]^=(a.mt[i][k]*b.mt[k][j]);
return res;
}
int query(int x){
matrix res;
res.len=1,res.wid=n;
for(int i=1;i<=n;i++)
res.mt[1][i]=s[i];
for(int i=0;i<maxk;i++)
if((x>>i)&1)
res=mul(res,ans[i]);
return res.mt[1][1];
}
signed main(){
zero.len=zero.wid=0;
memset(zero.mt,0,sizeof(zero.mt));
scanf("%lld%lld%lld",&n,&m,&q);
for(int i=1;i<=n;i++)
scanf("%lld",&s[i]);
ans[0].len=ans[0].wid=n;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
ans[0].mt[i][j]=0;
for(int i=1;i<=m;i++){
int x,y;
scanf("%lld%lld",&x,&y);
ans[0].mt[x][y]=ans[0].mt[y][x]=1;
}
for(int i=1;i<maxk;i++)
ans[i]=mul(ans[i-1],ans[i-1]);
for(int i=1;i<=q;i++){
int x;
scanf("%lld",&x);
printf("%lld\n",query(x));
}
return 0;
}
```
## 终极挑战：[P6772 [NOI2020]美食家](https://www.luogu.com.cn/problem/P6772)
学了矩阵加速图上问题，结果同步赛写挂了，只有$40$分，沦为暴力分/kk。
这道题是矩阵加速图上问题的一些技巧综合题，前置知识：拆点，预处理多个矩阵进行二进制拆分，即[P4159 [SCOI2009] 迷路](https://www.luogu.com.cn/problem/P4159)和[P6190 [NOI Online #1 入门组]魔法（民间数据）](https://www.luogu.com.cn/problem/P6190)。
加油，在前面的学习中，你已经锻炼了基本的矩阵乘法思维，现在是实践的时候了！
## 总结
矩阵加速是一个比较神奇的科技，复杂度是$O(|S|^3\log n)$（其中$|S|$为矩阵的边长，$n$为题目中询问的值，常为时间，距离等）。
因此，矩阵加速图上问题的题目大多会有这样的规律：**图的点数很小，或者边数很小，同时$n$非常大**，这个时候你基本上就可以确定使用矩阵快速幂了。
再加上矩阵乘法可以有重定义，bitset优化等黑科技，所以矩阵加速图上问题的应用其实非常广（但是我太菜了，并不知道很多应用），希望大家可以探索出更多的毒瘤操作，~~危害社会~~（狗头）。
## 其他
这里给出所有本人已知的此类型题目，有兴趣的读者可以随意切：
- [P2886 [USACO07NOV]Cow Relays G](https://www.luogu.com.cn/problem/P2886)（矩阵乘法重定义）
- [P3758 [TJOI2017]可乐](https://www.luogu.com.cn/problem/P3758)（加强版：[P5789 [TJOI2017]可乐（数据加强版）](https://www.luogu.com.cn/problem/P5789)）（普通矩阵加速）
- [P2233 [HNOI2002]公交车路线](https://www.luogu.com.cn/problem/P2233)（普通矩阵加速）
- [P2151 [SDOI2009]HH去散步](https://www.luogu.com.cn/problem/P2151)（点边互换，普通矩阵加速）
- [P4159 [SCOI2009] 迷路](https://www.luogu.com.cn/problem/P4159)（拆点，普通矩阵加速）
- [P2579 [ZJOI2005]沼泽鳄鱼](https://www.luogu.com.cn/problem/P2579)（预处理多个矩阵，分组快速幂）
- [P6569 [NOI Online #3 提高组]魔法值
](https://www.luogu.com.cn/problem/P6569)（矩阵乘法重定义，二进制拆分）
- [P6190 [NOI Online #1 入门组]魔法](https://www.luogu.com.cn/problem/P6190)（矩阵乘法重定义）
- [P3502 [POI2010]CHO-Hamsters](https://www.luogu.com.cn/problem/P3502)（把字符串转换为图，用矩阵加速Floyd）
- [P3821 Isaac](https://www.luogu.com.cn/problem/P3821)（与沼泽鳄鱼很像，预处理多个矩阵，分组快速幂）
- [CF576D Flights for Regular Customers](https://www.luogu.com.cn/problem/CF576D)（bitset优化矩阵乘法）
- [P6772 [NOI2020]美食家](https://www.luogu.com.cn/problem/P6772)（拆点，矩阵乘法重定义，二进制拆分）
参考文献：
- 俞华程《矩阵乘法在信息学中的应用》
- [扶苏的bitset浅谈](https://www.luogu.com.cn/blog/fusu2333/fu-su-di-bitset-qian-tan)