矩阵加速图上问题题型总结

题型总结 学习笔记

---

分类：$\bigstar$题型总结$\rightarrow$$\mathcal{A}$矩阵加速图上问题，简称$\bigstar.\mathcal{A}$

更好的阅读体验

前置知识

• 矩阵乘法与矩阵快速幂
• 简单的图上dp
• 一颗愿意学习的大脑
• 基本的码力

引入

我们看一道题：P2233 [HNOI2002]公交车路线

题意：有$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{gather*} \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} \end{gather*}$$

经过平方后，$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] 迷路（拆点，普通矩阵加速）

现在，我们就可以做一道题目来实践刚刚的结论了！

题意：给定一个带边权的无向图，求从点$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去散步（点边互换，普通矩阵加速）

同样，我们再做一道比较新颖的题目来巩固一下：P2151 [SDOI2009]HH去散步

这道题运用了点边互换的小trick，分析详见P2151 [SDOI2009]HH去散步解题报告。

例题3：CF576D Flights for Regular Customers（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浅谈进行学习）

代码：

#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]沼泽鳄鱼（分组矩阵加速）

有时候，单个的矩阵已经无法满足我们表示状态的需求了，此时我们需要处理多个矩阵，并同样地利用这些矩阵进行加速，我们以P2579 [ZJOI2005]沼泽鳄鱼为例。

题意：给定一个$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（矩阵乘法重定义）

上文说过，会对矩阵加速Floyd进行讨论，我们便以此题为例来探讨矩阵乘法重定义后对Floyd的加速。

我们以这道题为例子，来了解如何对矩阵乘法重定义，并证明矩阵加速的正确性。

题意：给定一张$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 提高组]魔法值

在这里，我们来根据P6569 [NOI Online #3 提高组]魔法值谈谈矩阵乘法的倍增预处理与二进制拆分优化。

题意：给定一个$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]美食家

学了矩阵加速图上问题，结果同步赛写挂了，只有$40$分，沦为暴力分/kk。

这道题是矩阵加速图上问题的一些技巧综合题，前置知识：拆点，矩阵乘法转向量乘法，即P4159 [SCOI2009] 迷路和P6569 [NOI Online #3 提高组]魔法值。

加油，在前面的学习中，你已经锻炼了基本的矩阵乘法思维，现在是实践的时候了！

总结

矩阵加速是一个比较神奇的科技，复杂度是$O(|S|^3\log n)$（其中$|S|$为矩阵的边长，$n$为题目中询问的值，常为时间，距离等）。

因此，矩阵加速图上问题的题目大多会有这样的规律：图的点数很小，或者边数很小，同时$n$非常大，这个时候你基本上就可以确定使用矩阵快速幂了。

再加上矩阵乘法可以有重定义，bitset优化等黑科技，所以矩阵加速图上问题的应用其实非常广（但是我太菜了，并不知道很多应用），希望大家可以探索出更多的毒瘤操作，危害社会（狗头）。

其他

这里给出所有本人已知的此类型题目，有兴趣的读者可以随意切：

• P2886 [USACO07NOV]Cow Relays G（矩阵乘法重定义）
• P3758 [TJOI2017]可乐（加强版：P5789 [TJOI2017]可乐（数据加强版））（普通矩阵加速）
• P2233 [HNOI2002]公交车路线（普通矩阵加速）
• P2151 [SDOI2009]HH去散步（点边互换，普通矩阵加速）
• P4159 [SCOI2009] 迷路（拆点，普通矩阵加速）
• P2579 [ZJOI2005]沼泽鳄鱼（预处理多个矩阵，分组快速幂）
• P6569 [NOI Online #3 提高组]魔法值
  （矩阵乘法重定义，二进制拆分）
• P6190 [NOI Online #1 入门组]魔法（矩阵乘法重定义）
• P3502 [POI2010]CHO-Hamsters（把字符串转换为图，用矩阵加速Floyd）
• P3821 Isaac（与沼泽鳄鱼很像，预处理多个矩阵，分组快速幂）
• CF576D Flights for Regular Customers（bitset优化矩阵乘法）
• P6772 [NOI2020]美食家（拆点，矩阵乘法重定义，二进制拆分）

参考文献：

• 俞华程《矩阵乘法在信息学中的应用》
• 扶苏的bitset浅谈