AcWing 345 牛站

题目传送门

一、前言

本题是求恰好经过\(k\)条边的最短路径，\(y\)总对这题的讲解比较粗糙 ，很多先修知识点()并没有讲解(菜是原罪)，这题是可以将很多重要的知识点串起来的好题，仔细分析会大有收获的。

二、BellmanFord算法

有边数限制的最短路径问题，第一反应自然是\(BellmanFord\)算法，只需对\(BellmanFord\)算法的逻辑稍加修改，再吸吸氧应该就可以\(AC\)了
为防止正解不突出，另开一篇文章写一下如何使用\(BellmanFord\)算法进行解决此问题：

三、矩阵快速幂

下面详细介绍下用矩阵快速幂的方法来求解本题。

1、图论中的矩阵乘法

设g[i][j] = 1表示i到j之间存在关联边，g[i][j] = 0表示不存在关联边，t[i][j]表示从i恰好经过两条边到达j的路径条数，很容易得出状态转移方程t[i][j] += g[i][k] * g[k][j]，k从1到n，如何得到这个式子的呢？只有i到k和k到j都存在边的时候，从i到k再到j才是刚好经过两条边的路径，所有这种路径条数的和就是t[i][j]了，转化为代码的形式就是:

for(int k = 1;k <= n;k++)
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n;j++)
            t[i][j] += g[i][k] * g[k][j];

观察上面的代码可以发现，上面的代码正是求矩阵t = g * g的代码，也就是说，对没有边权的图的邻接矩阵g而言，g * g中的第i行第j列就表示了从i到j存在几条路径长度为2的路径，同理g * g * g = t * g就等于从i先经过两条边再经过一条边到达j的路径条数，也就是经过三条边的路径条数，以此类推可以推出，g的k次方中存储的就是任意两点间经过k条边路径的条数。这里我们可以看出两点:

• 其一是在图论中,矩阵乘法具有求两点间经过若干条边路径条数的功能；

• 其二是对邻接矩阵做一次自乘就可以得到经过两条边的路径数，如果需要求经过更多条边的路径数，就需要求g的k次方了。

2、广义的矩阵乘法

上面的问题的实现就是矩阵乘法的代码，解决的是求出两点间路径的条数，但是现实中很多问题并不是要我们求路径的条数，比如本题就是在路径条数限定的情况下，求最短路径长度。现在设g就是图的邻接矩阵，t[i][j]表示从i恰好经过两条边到达j的最短路径长度，则t[i][j] = min(g[i][k] + g[k][j]);你可以想象一下，i,j两个点，如果之间有两条边，则需要引入一个点才行，设为k,则路径长度就是g[i][k],g[k][j],而这个k的范围就是枚举所有点得到,转化为代码的形式就是：

for(int k = 1;k <= n;k++)
    for(int i = 1;i <= n;i++)
        for(int j = 1;j <= n;j++)
            t[i][j] = min(t[i][j],g[i][k] + g[k][j]);

可以发现，这次的代码与矩阵乘法的代码很相似，可以看作广义的矩阵乘法。

3、矩阵快速幂对动态规划的优化

首先还是介绍矩阵快速幂是如何求普通矩阵的\(n\)次方的，之前写过的关于矩阵快速幂的题解见:
\(AcWing\) \(205\) 斐波那契

首先快速幂是如何求\(a^b\)的，比如求\(3^{13}\),\(13 = 1101\),也就是\(13 = 8 + 4 + 1\),\(3^{13}\) = \(3^8 * 3^4 * 3^1\),普通数的乘法具有结合律，\(a*b*c = a * (b * c)\)，所以快速幂可以将一系列的乘法拆分成几组乘法。

对于矩阵而言也具有结合律，\(A * B * C = A * (B * C)\)，矩阵t = g * g的展开式就是t[i][j] = sum(g[i][.] * g[.][j])，g * g * g = t * g = sum(t[i][.] * g[.][j]) = sum(g[i][.] * t[.][j]) = g * t，这就是矩阵乘法具有结合律的原因，有了结合律，六个矩阵相乘，就可以先算前两个矩阵的乘积和后四个矩阵的乘积，最后再乘到一起即可。

void mul(int c[][N],int a[][N],int b[][N]){
    static int t[N][N];
    memset(t,0,sizeof t);
    for(int k = 1;k <= n;k++)
        for(int i = 1;i <= n;i++)
            for(int j = 1;j <= n;j++)
                t[i][j] += a[i][k] + b[k][j];
    memcpy(c,t,sizeof t);
}
void qmi(){
    memset(f,0,sizeof f);
    for(int i = 1;i <= n;i++)   f[i][i] = 1;
    while(k){
        if(k & 1)   mul(f,f,g);
        mul(g,g,g);
        k >>= 1;
    }
}

上面的代码就是矩阵快速幂对普通矩阵乘法的优化，也就是将快速幂代码中的普通乘法修改为矩阵乘法\(mul\)即可。

看起来很简单，但是很多细节不自己去思考是不知道为什么要这么设置的，比如\(qmi\)里面的代码，为什么对存储矩阵乘法的结果矩阵\(f\)要初始化为单位矩阵呢？我们考虑边界情况，k = 1时候，f = f * g，我们需要\(g^1 = g\)，也就是需要\(f * g = g\)，所以\(f\)的初值需要设置为单位矩阵。

矩阵乘法的函数里对\(t\)的初值设置为\(0\)无可厚非，为什么需要先定义一个备份矩阵\(t\)存下矩阵乘法的结果再在计算完之后拷贝给矩阵\(c\)呢？ (黄海思考了很长时间也没想明白原因，但事实上不这么整确实不行~)

是因为很可能传进去的\(a\)、\(b\)、\(c\)都是一个矩阵，对同一个二维数组做自乘必然会改变矩阵原始的值引发运算错误，所以需要\(t\)存储中间的计算结果。 [这个解释太棒了！！终于明白了！！]

再来说下矩阵快速幂对广义矩阵乘法的优化。我们将二中实现t[i][j] = min(g[i][k] + g[k][j]);的代码，也就是广义矩阵乘法用mul表示，mul可以求g mul g，但是如何求若干个g相乘呢？考察min运算是否具有结合律，t[i][j] = min(g[i][.] + g[.][j]);，则min(t[i][.] + g[.][j]) = min(g[i][.] + t[.][j])，也就是说广义的矩阵乘法mul也具有结合律，所以也可以用矩阵快速幂对其做优化。
上面说过，广义矩阵乘法mul可以求i经过两条边到达j的最短路径长度，则对g做k次自乘就可以求经过k条边的最短路径长度了。原理很好理解，t = g mul g,t[i][k] + g[k][j]就是i经过k到达j也就是经过三条边的路径长度，在其中取个min就是经过三条边的最短路径长度了，k条边以此类推。

void mul(int c[][N],int a[][N],int b[][N]){
    static int t[N][N];
    memset(t,0x3f,sizeof t);
    for(int k = 1;k <= n;k++)
        for(int i = 1;i <= n;i++)
            for(int j = 1;j <= n;j++)
                t[i][j] = min(t[i][j],a[i][k] + b[k][j]);
    memcpy(c,t,sizeof t);
}
void qmi(){
    memset(f,0x3f,sizeof f);
    for(int i = 1;i <= n;i++)   f[i][i] = 0;
    while(k){
        if(k & 1)   mul(f,f,g);
        mul(g,g,g);
        k >>= 1;
    }
}

上面是矩阵快速幂优化本题中广义矩阵乘法的代码。对t[i][j]的求法改变是理所应当的，但是仔细观察会发现，前面矩阵快速幂对普通矩阵乘法的优化强调过的地方在这里就完全不一样了:

首先，qmi函数里结果矩阵f的初值变成了主对角线元素都是0，其他元素都是INF了!
这样做能否使得k = 1时，g * f还等于g呢？mul中求a[i][k] + b[k][j]时，也就是求f[i][k] + g[k][j]，f[i][k]除了在k = i时候是0，其他时候都是INF，不会更新最小值，而f[i][i] + g[i][j]就算i到j的距离，也就没有改变g的初值，这就是f的初值这样设置的原因。

对于mul中备份数组t的初值，设置为INF就理所应该了，因为后面需要求的是最小值。

本题的完整代码如下，值得注意的是点的编号最大是1000，最最多只有100条边，也就是最多200个节点，所以可以对节点编号做下离散化，较为简单，参考代码即可。还有个要注意的地方时，本题读取边的信息的时候是先读取边长，再读取两个节点的编号，一般题目都是最后读取边权，要注意不要写错数据的读取代码。

#include 
using namespace std;
const int N = 205;
unordered_map id;
int s;       //起点
int e;       //终点
int n;       //离散化后的节点号
int k;       //恰好k条边
int g[N][N]; //图
int f[N][N]; //最短距离

//矩阵乘法
void mul(int c[][N], int a[][N], int b[][N]) {
    static int t[N][N];
    memset(t, 0x3f, sizeof t);
    for (int k = 1; k <= n; k++)
        for (int i = 1; i <= n; i++)
            for (int j = 1; j <= n; j++)
                //求指定边数情况下的最短距离
                //广义上的矩阵乘法
                t[i][j] = min(t[i][j], a[i][k] + b[k][j]);
    memcpy(c, t, sizeof t);
}
//快速幂
void qmi() {
    memset(f, 0x3f, sizeof f);
    for (int i = 1; i <= n; i++) f[i][i] = 0; //这个初始化注意一下
    while (k) {
        if (k & 1) mul(f, f, g); //矩阵快速幂
        mul(g, g, g);
        k >>= 1;
    }
}
int main() {
    int m; // m条边
    cin >> k >> m >> s >> e;
    //起点的离散化后号为1
    id[s] = ++n;
    //如果e与s不同，那么e的新号为2,否则为1
    if (!id.count(e)) id[e] = ++n;
    //重新对s和e给定新的号码
    s = id[s], e = id[e];
    //初始化邻接矩阵
    memset(g, 0x3f, sizeof g);
    for (int i = 0; i < m; i++) {
        int a, b, c;
        cin >> c >> a >> b;
        if (!id.count(a)) id[a] = ++n; //记录点的映射关系a-> id[a]
        if (!id.count(b)) id[b] = ++n; //记录点的映射关系b-> id[b]
        a = id[a], b = id[b];          //对a,b给定新的号码
        //利用新的号码将边长c记录到邻接矩阵中
        g[a][b] = g[b][a] = min(g[a][b], c);
    }
    //快速幂+动态规划思想
    qmi();
    //输出从起点到终点，恰好经过k条边的最短路径
    cout << f[s][e] << endl;
    return 0;
}