多源最短路的两种做法和详细解析

一、概述

前置知识：单源最短路请参考这篇文章。

最短路算法是图论算法中的一个分支，大致可以分为单源最短路和多源最短路两种。单源最短路指从一个点去往图中的其他点所花费的最短距离，而多源最短路指图中任意两个点的最短距离。
多源最短路可以通过运行 $N$（ $N$ 为图的点数， $M$ 为图的边数）次单元最短路（dijkstra，bellman-ford，SPFA）来实现。运行 $N$ 次单元最短路的算法复杂度如下：

算法名称	单次复杂度	$N$ 次复杂度
dijkstra	$O((N+M)\log M)$	$O(N(N+M)\log M)$
bellman-ford	$O(M^2)$	$O(N{M}^2)$
SPFA	$O(kM)$	$O(kMN)$

运行 $N$ 次多源最短路的代码复杂度较高（需要将单元最短路的算法嵌套在一个 for 循环里面，还需特别注意初始化的问题），但专门用于解决多源最短路问题、基于动态规划实现的 Floyd 算法，可以使用 $5$ 行代码，用 $O(N^3)$ 且较为稳定的时间复杂度解决这个问题。

二、Floyd 算法的原理讲解

Floyd 算法是基于动态规划实现的最短路算法，主要包括初始化和进行动态规划的过程两部分。

1. 初始化：

使用 Floyd 算法计算需要使用一个邻接矩阵 $E$ 来存图。进行输入的时候，依据题意建立无向边或者有向边。例如，令图中一条边的两个端点分别为 $u,v$ ，边权为 $w$ ，当输入格式为 $uvw$ 时，我们可以这样通过输入建立一个邻接矩阵：

• 使 $E(i,j)=\infty (i,j\in [1,N]\text{ and }i\ne j)$
• 输入为 $123$ ，使 $E(1,2)=3,E(2,1)=3$ ，这样可以建立无向边
• 输入为 $134$ ，使 $E(1,3)=4,E(3,1)=4$
• 这样我们建立了一个如下样式的图：
  

使用邻接矩阵表示为：
$$E=\left[\begin{array}{ccc}0& 3& 4\\ 3& 0& \infty \\ 4& \infty & 0\end{array}\right]$$

使用 Floyd 算法时需要注意必须使用邻接表存图，否则 Floyd 算法在运算过程中不能高效率地处理边权的问题。

2. 算法的原理以及其过程

Floyd 算法是基于动态规划算法实现的一种最短路算法，因此我们需要使用分析 DP 的思路来分析 Floyd 算法的原理和过程。

Floyd 算法的大致过程如下：
该算法可以大致的表示为中转法。
依次枚举可以中转的点 $k\in [1,N]$ ，然后看点 $i,j$ 之间的最短距离能否被上一步计算的从 $i$ 到 $k$ 的最短距离和从 $k$ 到 $j$ 的最短距离更新。如果可以，就让 $f(i,j)$ 进行一次松弛操作。例如：

• 我们先假设有一个如下样式的图。
  
• 首先，我们先按照上面的方法进行初始化，让整个邻接矩阵变成如下样式：
  $$E=\left[\begin{array}{cccc}0& 3& 2& 1\\ 3& 0& \infty & 5\\ 2& \infty & 0& 2\\ 1& 5& 2& 0\end{array}\right]$$
  其中， $\infty$ 表示还未找到长度的点。
  我们直接把 $E$ 这个邻接矩阵看作为已知的临时最短路径答案。
• 然后，我们以 $1$ 作为中转点，枚举已知的点。此时，我们可以更新 $E(2,4),E(4,2),E(3,2),E(2,3)$ 的值。此时整个邻接矩阵变成如下样式：
  $$E=\left[\begin{array}{cccc}0& 3& 2& 1\\ 3& 0& 5& 4\\ 2& 5& 0& 2\\ 1& 4& 2& 0\end{array}\right]$$
• 我们以 $2$ 作为中转点，枚举已知的点。此时，我们可以根据上面的邻接矩阵更新当前的状态，不能更新任何一项的值。
• 同理，我们使用 $3,4$ 作为中转点，均不能使 $E$ 变化。
• 因此，记 $i,j$ 之间的最短路为 $E(i,j)$ ，则最终的答案是：
  $$E=\left[\begin{array}{cccc}0& 3& 2& 1\\ 3& 0& 5& 4\\ 2& 5& 0& 2\\ 1& 4& 2& 0\end{array}\right]$$

下面我们开始介绍如何使用动态规划 $\text{DP}$ 策略解决这个问题。
我们先记 $d(i,j)$ 表示 $i,j$ 之间的距离。
我们先使中转点 $k\in [1,N]$ 。
通过上面的过程，我们可以发现，每一次更新当前的最短路都会依据前一次的最短路计算结果更新当前的计算结果。因此，状态 $d$ 其实应该是 $3$ 维的，记 $d(i,j,k)$ 表示以 $p\in \forall [1,k]$ 为中转点的 $i,j$ 之间的最短距离。可以看出，每个点位有两种可能的情况：
$$d(i,j,k)=\{\begin{array}{ll}d(i,j,k-1)& \text{不使用上一步的数值更新当前的数值}\\ d(i,k,k-1)+d(k,j,k-1)& \text{使用上一步的 }d(i,k)\text{ 和 }d(k,j)\text{ 更新}\end{array}$$
显然，最后的最短距离即两者之间的最小值。因此，我们得到了 Floyd 算法的 $\text{DP}$ 状态转移方程：
d ( i , j , k ) = min ⁡ { d ( i , j , k − 1 ) , d ( i , k , k − 1 ) + d ( k , j , k − 1 ) } d(i,j,k)=\min\{d(i,j,k-1),d(i,k,k-1)+d(k,j,k-1)\} d(i,j,k)=min{d(i,j,k−1),d(i,k,k−1)+d(k,j,k−1)}

显然，当前的状态 $d(i,j,k)$ 是一个三维的状态。下面考虑优化，将三维的状态划分为二维的状态。
∵ d ( i , j , k ) = min ⁡ { d ( i , j , k − 1 ) , d ( i , k , k − 1 ) + d ( k , j , k − 1 ) } \because d(i,j,k)=\min\{d(i,j,k-1),d(i,k,k-1)+d(k,j,k-1) \} ∵d(i,j,k)=min{d(i,j,k−1),d(i,k,k−1)+d(k,j,k−1)} ，
$\therefore d(i,j,k)$ 只和 $d(p,q,k-1)(p=i\text{ or }p=k,q=k\text{ or }q=j)$ 的值有关。
$\therefore$ 根据 ⁰/₁ 背包的优化方法，我们可以让 $d(i,j,k)$ 用同样的方法变成二维的状态 $d(i,j)$ 。
$\therefore$ 此时状态转移方程变为：
d ( i , j ) = min ⁡ k ∈ [ 1 , N ] { d ( i , k ) + d ( k , j ) } d(i,j)=\min_{k \in [1,N]}\{d(i,k)+d(k,j)\} d(i,j)=k∈[1,N]min​{d(i,k)+d(k,j)}
利用上面的状态转移方程，我们可以将空间复杂度优化到 $O(N^2)$ 。

3. 例题展示和代码实例

使用 Floyd 算法和两点间距离公式 $s=\sqrt{(x_0-x_1{)}^2+(y_0-y_1{)}^2}$ 即可解决。代码如下：

#include <bits/stdc++.h>

using namespace std;

const int N = 157;

struct Node
{
    double x, y;
} a[N];

int n, f[N], k;
double d[N][N], e[N][N], ans = 1e9, ma[N], f_ma[N];

void dfs(int x)
{
    for (int i = 1; i <= n; ++i)
    {
        if (x == i) continue;
        if (e[x][i] > 1e8) continue;
        if (f[i]) continue;

        f[i] = k;
        dfs(i);
    }
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++i)
    {
        cin >> a[i].x >> a[i].y;
    }
    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= n; ++j)
        {
            int t;
            scanf("%1d", &t);
            
            if (t == 0) 
            {
                if (i != j) e[i][j] = 1e9;
                else e[i][j] = 0;
            }
            else
            {
                e[i][j] = sqrt((a[j].x - a[i].x) * (a[j].x - a[i].x) + (a[j].y - a[i].y) * (a[j].y - a[i].y) * 1.0);
            }
        }
    }

    for (int i = 1; i <= n; ++i)
    {
        if (f[i]) continue;
        
        f[i] = ++k;
        dfs(i);
    }

	//Floyd
    memcpy(d, e, sizeof(e));
    for (int k = 1; k <= n; ++k)
    {
        for (int i = 1; i <= n; ++i)
        {
            for (int j = 1; j <= n; ++j)
            {
                d[i][j] = min(d[i][j], d[i][k] + d[k][j]);
            }
        }
    }
    //Floyd-end

    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= n; ++j)
        {
            if (f[i] == f[j])
            {
                ma[i] = max(ma[i], d[i][j]);
            }
        }

        f_ma[f[i]] = max(f_ma[f[i]], ma[i]);
    }

    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j <= n; ++j)
        {
            if (f[i] == f[j]) continue;

            double t = max(max(f_ma[f[i]], f_ma[f[j]]), ma[i] + sqrt((a[j].x - a[i].x) * (a[j].x - a[i].x) + (a[j].y - a[i].y) * (a[j].y - a[i].y) * 1.0) + ma[j]);
            ans = min(ans, t);
        }
    }

    cout << fixed << setprecision(6) << ans << endl;

    return 0;
}