矩阵加速图上问题学习笔记

whiteqwq

2020-08-05 22:07:19

Algo. & Theory

矩阵加速图上问题学习笔记

更好的阅读体验

前置知识

引入

我们看一道题：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比较大的数据范围，你想到了什么？

我们直接重构代码，上矩阵快速幂： ``` #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] 迷路（拆点，普通矩阵加速）

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

题意：给定一个带边权的无向图，求从点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的状态，这样就好办了。

这样，我们就可以直接上板子了：

#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] 迷路和P6190 [NOI Online #1 入门组]魔法（民间数据）。

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

总结

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

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

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

其他

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

参考文献：