分类 技术世界 下的文章

一:背景

SPFA(Shortest Path Faster Algorithm)算法,是西南交通大学段凡丁于 1994 年发表的,其在 Bellman-ford 算法的基础上加上一个队列优化,减少了冗余的松弛操作,是一种高效的最短路算法。

二:算法过程

设立一个队列用来保存待优化的顶点,优化时每次取出队首顶点 u,并且用 u 点当前的最短路径估计值dist[u]对与 u 点邻接的顶点 v 进行松弛操作,如果 v 点的最短路径估计值dist[v]可以更小,且 v 点不在当前的队列中,就将 v 点放入队尾。这样不断从队列中取出顶点来进行松弛操作,直至队列空为止。(所谓的松弛操作,简单来说,对于顶点 i,把dist[i]调整更小。更多解释请参考百科:松弛操作)

而其检测负权回路的方法也很简单,如果某个点进入队列的次数大于等于 n,则存在负权回路,其中 n 为图的顶点数。

三:代码

#include <iostream>    
#include <queue>
#include <stack>

using namespace std;

int  matrix[100][100]; // 邻接矩阵
bool visited[100];     // 标记数组
int  dist[100];        // 源点到顶点 i 的最短距离
int  path[100];        // 记录最短路的路径
int  enqueue_num[100]; // 记录入队次数
int  vertex_num;       // 顶点数
int  edge_num;         // 边数
int  source;           // 源点

bool SPFA()
{
    memset(visited, 0, sizeof(visited));
    memset(enqueue_num, 0, sizeof(enqueue_num));
    for (int i = 0; i < vertex_num; i++)
    {
        dist[i] = INT_MAX;
        path[i] = source;
    }

    queue<int> Q;
    Q.push(source);
    dist[source] = 0;
    visited[source] = 1;
    enqueue_num[source]++;

    while (!Q.empty())
    {
        int u = Q.front();
        Q.pop();
        visited[u] = 0;

        for (int v = 0; v < vertex_num; v++)
        {
            if (matrix[u][v] != INT_MAX)  // u 与 v 直接邻接
            {
                if (dist[u] + matrix[u][v] < dist[v])
                {
                    dist[v] = dist[u] + matrix[u][v];
                    path[v] = u;

                    if (!visited[v])
                    {
                        Q.push(v);
                        enqueue_num[v]++;
                        if (enqueue_num[v] >= vertex_num)
                            return false;
                        visited[v] = 1;
                    }
                }
            }
        }
    }

    return true;
}

void Print()
{
    for (int i = 0; i < vertex_num; i++)
    {
        if (i != source)
        {
            cout << source << " 到 " << i << " 的最短距离是:" << dist[i] << ",路径是:" << i;
            int t = path[i];
            while (t != source)
            {
                cout << "--" << t;
                t = path[t];
            }
            cout << "--" << source << endl;
        }
    }
}

int main()
{

    cout << "请输入图的顶点数,边数,源点:";
    cin >> vertex_num >> edge_num >> source;

    for (int i = 0; i < vertex_num; i++)
        for (int j = 0; j < vertex_num; j++)
            matrix[i][j] = (i != j) ? INT_MAX : 0;  // 初始化 matrix 数组

    cout << "请输入 " << edge_num << " 条边的信息:\n";
    int u, v, w;
    for (int i = 0; i < edge_num; i++)
    {
        cin >> u >> v >> w;
        matrix[u][v] = w;
    }

    if (SPFA())
        Print();
    else
        cout << "存在负权回路!\n";

    return 0;
}

运行如下:

/* Test 1 */
请输入图的顶点数,边数,源点:5 7 0
请输入 7 条边的信息:
0 1 100
0 2 30
0 4 10
2 1 60
2 3 60
3 1 10
4 3 50
0 到 1 的最短距离是:70,路径是:1--3--4--0
0 到 2 的最短距离是:30,路径是:2--0
0 到 3 的最短距离是:60,路径是:3--4--0
0 到 4 的最短距离是:10,路径是:4--0

/* Test 2 */
请输入图的顶点数,边数,源点:4 6 0
请输入 6 条边的信息:
0 1 20
0 2 5
3 0 -200
1 3 4
3 1 4
2 3 2
存在负权回路!

四:判断负权回路的证明

如果某个点进入队列的次数大于等于 n,则存在负权回路。为什么偏偏是 n?

对于一个不存在负权回路的图,设其顶点数为 n,我们把图稍微“转换”下,如下图 A:

  • 与源点 0 邻接的点{ 1, 2, 3 }作为第一批次;
  • 与第一批次邻接的点{ 4, 5, 6, 7, 8, 9 }作为第二批次;
  • ......
  • 与第 k-1 批次邻接的点{ ...... }作为第k批次。

其中 k≤n-1,当 k=n-1 时,即为上图 B。

每操作完一个批次的点,至少有一个点的最短路径被确定。这里读者只需从 Dijkstra 算法方面来考虑即可。Dijkstra 每次循环都找出dist[]里的最小值,可以对应到这里的每个批次。

一个不存在负权回路的图,最多有 n-1 个批次,每做完一个批次至少有一个点的最短路径被确定,即一个点的入队次数不超过 n-1。因为若一个顶点要入队列,则必存在一条权值之和更小的路径,而在最多做完 n-1 个批次后,所有顶点的最短路径都被确定。(这里需要注意的是,如果一个批次中,有多条路径对某顶点进行更新,则该顶点只会被入队一次,这从代码就可以看出)

五:时间复杂度

对于一个不存在负权回路的图,我们假设其顶点数为 n,边数为 m。

引自 SPFA 论文:考虑一个随机图,运用均摊分析的思想,每个点的平均出度为 $O(\frac m n)$,而每个点的平均入队次数为 2,因此时间复杂度为 $O(n⋅\frac m n⋅2)=O(2m)=O(m)$。

关于上述的“平均入队次数为 2”,2 这个数字从何得来,我也找不到证明,从网上各位朋友对此的一致态度:尚待商榷。但是可以确定的是,SPFA 算法在随机图中的平均性能是优于 Bellman_Ford 算法的。

SPFA 的最佳时间复杂度为 $O(n)$。比如上图 B,每个点只入队一次。

接着再看下 SPFA 的最差时间复杂度,它发生在一个完全图中,如下图(为突出重点,其余边未画出),

我们约定,0 点为源点,每次更新完 k 点出队后,k+1​ 点都可以再次对 k 点进行更新并入队,其中。​1≤ k≤ n-2​ 那么我们得出:

0 点,入队 1 次;
1 点,入队 n-1 次;
2 点,入队 n-2 次;
3 点,入队 n-3 次;
.
n-2 点,入队 2 次;
n-1 点,入队 1 次;

因为是完全图,所以每个点的出度为 n-1,因此总的时间复杂度为:

$$ (n-1)⋅[1+1+2+3+...+(n-2)+(n-1)]=O(n^3) $$

由于是完全图,也可以表达成 $O(nm)$。很容易看出,SPFA 算法的时间复杂度很不稳定。

一:背景

Dijkstra 算法是处理单源最短路径的有效算法,但它对存在负权回路的图就会失效。这时候,就需要使用其他的算法来应对这个问题,Bellman-Ford(中文名:贝尔曼-福特)算法就是其中一个。

Bellman-Ford 算法不仅可以求出最短路径,也可以检测负权回路的问题。该算法由美国数学家理查德•贝尔曼(Richard Bellman, 动态规划的提出者)和小莱斯特•福特(Lester Ford)发明。

二:算法过程分析

对于一个不存在负权回路的图,Bellman-Ford 算法求解最短路径的方法如下:

设其顶点数为 n,边数为 m。设其源点为 source,数组dist[i]记录从源点 source 到顶点i的最短路径,除了dist[source]初始化为 0 外,其它dist[]皆初始化为 INT_MAX。以下操作循环执行 n-1 次:

  • 对于每一条边 arc(u, v),如果 dist[u] + w(u, v) < dist[v],则使 dist[v] = dist[u] + w(u, v),其中 w(u, v) 为边 arc(u, v) 的权值。

n-1 次循环,Bellman-Ford 算法就是利用已经找到的最短路径去更新其它点的dist[]

接下来再看看 Bellman-Ford 算法是如何检测负权回路的?

检测的方法很简单,只需在求解最短路径的 n-1 次循环基础上,再进行第 n 次循环:

  • 对于所有边,只要存在一条边 arc(u, v) 使得 dist[u] + w(u, v) < dist[v],则该图存在负权回路,其中 w(u, v) 为边 arc(u, v) 的权值。
循环次数dist[0]dist[1]dist[2]
第1次0-5-3
第2次-2-5-3
第3次-2-7-5

三:完整代码

#include <iostream>
#include <stack>

using namespace std;

struct Edge
{
    int u;
    int v;
    int w;
};

Edge edge[10000]; // 记录所有边
int  dist[100];   // 源点到顶点 i 的最短距离
int  path[100];   // 记录最短路的路径
int  vertex_num;  // 顶点数
int  edge_num;    // 边数
int  source;      // 源点

bool BellmanFord()
{
    // 初始化
    for (int i = 0; i < vertex_num; i++)
        dist[i] = (i == source) ? 0 : INT_MAX;

    // n-1 次循环求最短路径
    for (int i = 1; i <= vertex_num - 1; i++)
    {
        for (int j = 0; j < edge_num; j++)
        {
            if (dist[edge[j].v] > dist[edge[j].u] + edge[j].w)
            {
                dist[edge[j].v] = dist[edge[j].u] + edge[j].w;
                path[edge[j].v] = edge[j].u;
            }
        }
    }

    bool flag = true;  // 标记是否有负权回路

    // 第 n 次循环判断负权回路
    for (int i = 0; i < edge_num; i++)
    {
        if (dist[edge[i].v] > dist[edge[i].u] + edge[i].w)
        {
            flag = false;
            break;
        }
    }

    return flag;
}

void Print()
{
    for (int i = 0; i < vertex_num; i++)
    {
        if (i != source)
        {
            cout << source << " 到 " << i << " 的最短距离是:" << dist[i] << ",路径是:" << i;
            int t = path[i];
            while (t != source)
            {
                cout << "--" << t;
                t = path[t];
            }
            cout << "--" << source << endl;
        }
    }
}

int main()
{

    cout << "请输入图的顶点数,边数,源点:";
    cin >> vertex_num >> edge_num >> source;

    cout << "请输入 " << edge_num << " 条边的信息:\n";
    for (int i = 0; i < edge_num; i++)
        cin >> edge[i].u >> edge[i].v >> edge[i].w;

    if (BellmanFord())
        Print();
    else
        cout << "存在负权回路!\n";

    return 0;
}

测试如下:

/* Test 1 */
请输入图的顶点数,边数,源点:5 7 0
请输入 7 条边的信息:
0 1 100
0 2 30
0 4 10
2 1 60
2 3 60
3 1 10
4 3 50
0 到 1 的最短距离是:70,路径是:1--3--4--0
0 到 2 的最短距离是:30,路径是:2--0
0 到 3 的最短距离是:60,路径是:3--4--0
0 到 4 的最短距离是:10,路径是:4--0

/* Test 2 */
请输入图的顶点数,边数,源点:4 6 0
请输入 6 条边的信息:
0 1 20
0 2 5
3 0 -200
1 3 4
3 1 4
2 3 2
存在负权回路!

四:算法优化

以下除非特殊说明,否则都默认是不存在负权回路的。

先来看看 Bellman-Ford 算法为何需要循环 n-1 次来求解最短路径?

读者可以从 Dijkstra 算法来考虑,想一下,Dijkstra 从源点开始,更新dist[],找到最小值,再更新dist[] ,,,每次循环都可以确定至少一个点的最短路。Bellman-Ford 算法同样也是这样,它的每次循环也可以确定至少一个点的最短路,故需要 n-1 次循环。

Bellman-Ford 算法的时间复杂度为 $O(nm)$,其中 n 为顶点数,m 为边数。每一次循环都需要对 m 条边进行操作,$O(nm)$ 的时间,其实大多数都浪费了。

大家可以考虑一个随机图(点和边随机生成),除了已确定最短路的顶点与尚未确定最短路的顶点之间的边,其它的边所做的都是无用的,大致描述为下图(分割线以左为已确定最短路的顶点),

其中红色部分为所做无用的边,蓝色部分为实际有用的边。

既然只需用到中间蓝色部分的边,那算法优化的方向就找到了,请接着看本系列第三篇文章:spfa 算法。

一:背景

Dijkstra 算法(中文名:迪杰斯特拉算法)是由荷兰计算机科学家 Edsger Wybe Dijkstra 提出。该算法常用于路由算法或者作为其他图算法的一个子模块。举例来说,如果图中的顶点表示城市,而边上的权重表示城市间开车行经的距离,该算法可以用来找到两个城市之间的最短路径。

二:算法过程

我们用一个例子来具体说明迪杰斯特拉算法的流程。

定义源点为 0,dist[i]为源点 0 到顶点i的最短路径。其过程描述如下:

步骤dist[1]dist[2]dist[3]dist[4]已找到的集合
第1步812+∞{ 2 }
第2步8×24{ 2, 3 }
第3步5××4{ 2, 3, 4 }
第4步5×××{ 2, 3, 4, 1 }
第5步××××{ 2, 3, 4, 1 }

第 1 步:从源点 0 开始,找到与其邻接的点:1,2,3,更新dist[]数组,因 0 不与 4 邻接,故dist[4]为正无穷。在dist[]中找到最小值,其顶点为 2,即此时已找到 0 到 2 的最短路。

第 2 步:从 2 开始,继续更新dist[]数组:2 与 1 不邻接,不更新;2 与 3 邻接,因0→2→3dist[3]大,故不更新dist[3] ;2 与 4 邻接,因0→2→4dist[4]小,故更新dist[4]为 4。在dist[]中找到最小值,其顶点为 3,即此时又找到 0 到 3 的最短路。

第 3 步:从 3 开始,继续更新dist[]数组:3 与 1 邻接,因0→3→1dist[1]小,更新dist[1]为 5;3 与 4 邻接,因0→3→4dist[4]大,故不更新。在dist[]中找到最小值,其顶点为 4,即此时又找到 0 到 4 的最短路。

第 4 步:从 4 开始,继续更新dist[]数组:4 与 1 不邻接,不更新。在dist[]中找到最小值,其顶点为 1,即此时又找到 0 到 1 的最短路。

第5步:所有点都已找到,停止。

对于上述步骤,你可能存在以下的疑问:

若 A 作为源点,与其邻接的只有 B,C,D 三点,其dist[]最小时顶点为 C,即就可以确定A→C为 A 到 C 的最短路。但是我们存在疑问的是:是否还存在另一条路径使 A 到 C 的距离更小? 用反证法证明。

假设存在如上图的红色虚线路径,使A→D→C的距离更小,那么A→D作为A→D→C的子路径,其距离也比A→C小,这与前面所述“dist[]最小时顶点为 C”矛盾,故假设不成立。因此这个疑问不存在。

根据上面的证明,我们可以推断出,Dijkstra 每次循环都可以确定一个顶点的最短路径,故程序需要循环 n-1 次。

三:完整代码

#include <iostream>

using namespace std;

int  matrix[100][100]; // 邻接矩阵
bool visited[100];     // 标记数组
int  dist[100];        // 源点到顶点 i 的最短距离
int  path[100];        // 记录最短路的路径
int  source;           // 源点
int  vertex_num;       // 顶点数
int  edge_num;         // 边数

void Dijkstra(int source)
{
    memset(visited, 0, sizeof(visited));  // 初始化标记数组
    visited[source] = true;
    for (int i = 0; i < vertex_num; i++)
    {
        dist[i] = matrix[source][i];
        path[i] = source;
    }

    int min_cost;        // 权值最小
    int min_cost_index;  // 权值最小的下标

    for (int i = 1; i < vertex_num; i++)  // 找到源点到另外 vertex_num-1 个点的最短路径
    {
        min_cost = INT_MAX;

        for (int j = 0; j < vertex_num; j++)
        {
            if (visited[j] == false && dist[j] < min_cost)  // 找到权值最小
            {
                min_cost = dist[j];
                min_cost_index = j;
            }
        }

        visited[min_cost_index] = true;  // 该点已找到,进行标记

        for (int j = 0; j < vertex_num; j++)  // 更新 dist 数组
        {
            if (visited[j] == false &&
                matrix[min_cost_index][j] != INT_MAX &&  // 确保两点之间有边
                matrix[min_cost_index][j] + min_cost < dist[j])
            {
                dist[j] = matrix[min_cost_index][j] + min_cost;
                path[j] = min_cost_index;
            }
        }
    }
}

int main()
{
    cout << "请输入图的顶点数(<100):";
    cin >> vertex_num;
    cout << "请输入图的边数:";
    cin >> edge_num;

    for (int i = 0; i < vertex_num; i++)
        for (int j = 0; j < vertex_num; j++)
            matrix[i][j] = (i != j) ? INT_MAX : 0;  // 初始化 matrix 数组

    cout << "请输入边的信息:\n";
    int u, v, w;
    for (int i = 0; i < edge_num; i++)
    {
        cin >> u >> v >> w;
        matrix[u][v] = matrix[v][u] = w;
    }

    cout << "请输入源点(<" << vertex_num << "):";
    cin >> source;
    Dijkstra(source);

    for (int i = 0; i < vertex_num; i++)
    {
        if (i != source)
        {
            cout << source << " 到 " << i << " 的最短距离是:" << dist[i] << ",路径是:" << i;
            int t = path[i];
            while (t != source)
            {
                cout << "--" << t;
                t = path[t];
            }
            cout << "--" << source << endl;
        }
    }

    return 0;
}

输入数据,结果为:

请输入图的顶点数(<100):5
请输入图的边数:7
请输入边的信息:
0 1 3
0 2 1
0 3 2
1 3 3
2 3 2
3 4 3
2 4 3
请输入源点(<5):0
0 到 1 的最短距离是:3,路径是:1--0
0 到 2 的最短距离是:1,路径是:2--0
0 到 3 的最短距离是:2,路径是:3--0
0 到 4 的最短距离是:4,路径是:4--2--0

四:时间复杂度

设图的边数为 m,顶点数为 n。

Dijkstra 算法最简单的实现方法是用一个数组来存储所有顶点的dist[](即本程序采用的策略),所以搜索dist[]中最小元素的运算需要线性搜索 $O(n)$。这样的话算法的运行时间是 $O(n^2)$。

对于边数远少于 $n^2$ 的稀疏图来说,我们可以用邻接表来更有效的实现该算法,同时需要将一个二叉堆或者斐波纳契堆用作优先队列来查找最小的顶点。当用到二叉堆的时候,算法所需的时间为 $O((m+n)logn)$,斐波纳契堆能稍微提高一些性能,让算法运行时间达到 $O(m+nlogn)$。然而,使用斐波纳契堆进行编程,常常会由于算法常数过大而导致速度没有显著提高。

关于 $O((m+n)logn)$ 的由来,我简单的证明了下:

  • dist[]数组调整成最小堆,需要 $O(n)$ 的时间;
  • 因为是最小堆,所以每次取出最小值只需 $O(1)$ 的时间,接着把数组尾的数放置堆顶,并花费 $O(logn)$ 的时间重新调整成最小堆;
  • 我们需要 n-1 次操作才可以找出剩下的 n-1 个点,在这期间,大约需要访问 m 次边,每次访问都可能造成dist[]的改变,因此还需要 $O(logn)$ 的时间来进行最小堆的重新调整(从当前dist[]改变的位置往上调整)。

综上所述:总的时间复杂度为:$O(n)+O(nlogn)+O(mlogn)=O((m+n)logn)$

最后简单说下 Dijkstra 优化时二叉堆的两种实现方式:

  • 优先队列,把每个顶点的序号和其dist[]压在一个结构体再放进队列里;
  • 自己建一个小顶堆heap[],存储顶点序号,再用一个数组pos[]记录第 i 个顶点在堆中的位置。

相比之下,前者的编码难度较低,因此在平时编程甚至算法竞赛中,都是首选。

五:该算法的缺陷

Dijkstra 算法有个巨大的缺陷,请考虑下面这幅图:

u→v间存在一条负权回路(所谓的负权回路,维基和百科并未收录其名词,但从网上的一致态度来看,其含义为:如果存在一个环(从某个点出发又回到自己的路径),而且这个环上所有权值之和是负数,那这就是一个负权环,也叫负权回路),那么只要无限次地走这条负权回路,便可以无限制地减少它的最短路径权值,这就变相地说明最短路径不存在。一个不存在最短路径的图,Dijkstra 算法无法检测出这个问题,其最后求解的dist[]也是错的。

那么对于上述的“一个不存在最短路径的图”,我们该用什么方法来解决呢?请接着看本系列第二篇文章。

首先介绍下什么是凸包?如下图:

在一个二维坐标系中,有若干点杂乱排列着,将最外层的点连接起来构成的凸多边型,它能包含给定的所有的点,这个多边形就是凸包。

寻找凸包的算法有很多种,Graham Scan 算法是一种十分简单高效的二维凸包算法,能够在 $O(nlogn)$ 的时间内找到凸包。

Graham Scan 算法的做法是先确定一个起点(一般是最左边的点和最右边的点),然后一个个点扫过去,如果新加入的点和之前已经找到的点所构成的"壳"凸性没有变化,就继续扫,否则就把已经找到的最后一个点删去,再比较凸性,直到凸性不发生变化。分别扫描上下两个"壳",合并在一起,凸包就找到了。这么说很抽象,我们看图来解释:

先找"下壳",上下其实是一样的。首先加入两个点 A 和 B。

然后插入第三个点 C,并计算 $\overrightarrow{AB}×\overrightarrow{BC}$ 的向量积,却发现向量积系数小于(等于)0,也就是说 $\overrightarrow{BC}$ 在 $\overrightarrow{AB}$ 的顺时针方向上。

于是删去 B 点。

按照这样的方法依次扫描,找完"下壳"后,再找"上壳"。

关于扫描的顺序,有坐标序和极角序两种,本文采用前者。坐标序是比较两个点的 x 坐标,小的先被扫描(扫描上凸壳的时候反过来),如果两个点 x 坐标相同,那么就比较 y 坐标,同样的也是小的先被扫描(扫描上凸壳的时候也是反过来)。极角序使用atan2函数的返回值进行比较,读者可以自己尝试写下。

下面贴下代码:

struct Point
{
    double x, y;

    Point operator-(Point & p)
    {
        Point t;
        t.x = x - p.x;
        t.y = y - p.y;
        return t;
    }

    double cross(Point p) // 向量叉积
    {
        return x * p.y - p.x * y;
    }
};

bool cmp(Point & p1, Point & p2)
{
    if (p1.x != p2.x)
        return p1.x < p2.x;

    return p1.y < p2.y;
}

Point point[1005];  // 无序点
int   convex[1005]; // 保存组成凸包的点的下标
int   n;            // 坐标系的无序点的个数

int GetConvexHull()
{
    sort(point, point + n, cmp);
    int temp;
    int total = 0;

    for (int i = 0; i < n; i++) // 下凸包
    {
        while (total > 1 && 
              (point[convex[total - 1]] - point[convex[total - 2]]).cross(point[i] - point[convex[total - 1]]) <= 0)
            total--;

        convex[total++] = i;
    }

    temp = total;

    for (int i = n - 2; i >= 0; i--) // 上凸包
    {
        while (total > temp && 
              (point[convex[total - 1]] - point[convex[total - 2]]).cross(point[i] - point[convex[total - 1]]) <= 0)
            total--;

        convex[total++] = i;
    }

    return total - 1; // 返回组成凸包的点的个数,实际上多了一个,就是起点,所以组成凸包的点个数是 total - 1
}

参考文献:

一:什么是向量积?

向量积,也称(向量)叉积,(向量)叉乘,外积,是一种在向量空间中对向量进行的二元运算。常见于物理学力学、电磁学、光学和计算机图形学等理工学科中,是一种很重要的概念。

设向量 $\overrightarrow{c}$ 由两个向量 $\overrightarrow{a}$ 和 $\overrightarrow{b}$ 按如下公式定出:$\overrightarrow{c}$ 的模 $|\overrightarrow{c}|=|\overrightarrow{a}||\overrightarrow{b}|sinθ$,其中 $θ$ 为 $\overrightarrow{a}$ 和 $\overrightarrow{b}$ 间的夹角;$\overrightarrow{c}$ 的方向垂直于 $\overrightarrow{a}$ 和 $\overrightarrow{b}$ 所决定的平面,指向按右手规则从 $\overrightarrow{a}$ 转向 $\overrightarrow{b}$ 来确定,如下图:

那么,向量 $\overrightarrow{c}$ 叫做向量 $\overrightarrow{a}$ 与 $\overrightarrow{b}$ 的向量积,记作 $\overrightarrow{a}×\overrightarrow{b}$。

由上述的定义,我们很容易总结出两条性质:

$$ \begin{align} \overrightarrow{a}×\overrightarrow{b}&=\overrightarrow{0} \tag{其中 $\overrightarrow{a}$ 平行$\overrightarrow{b}$}\\ \overrightarrow{a}×\overrightarrow{b}&=- \overrightarrow{b}×\overrightarrow{a} \tag{不满足交换律} \end{align} $$

下面来推导向量积的坐标表达式,以二维向量为例。设 $\overrightarrow{a}=(a_x, a_y),\overrightarrow{b}=(b_x, b_y)$,得:

$$ \begin{align} \overrightarrow{a}×\overrightarrow{b}&=(a_x\overrightarrow{i}+a_y\overrightarrow{j})×(b_x\overrightarrow{i}+b_y\overrightarrow{j}) \tag{其中 $\overrightarrow{i}$ 和 $\overrightarrow{j}$ 分别是 $xy$ 轴上的单位向量}\\ &=a_x\overrightarrow{i}×(b_x\overrightarrow{i}+b_y\overrightarrow{j})+a_y\overrightarrow{j}×(b_x\overrightarrow{i}+b_y\overrightarrow{j}) \tag{分解}\\ &=a_xb_x(\overrightarrow{i}×\overrightarrow{i})+a_yb_y(\overrightarrow{j}×\overrightarrow{j})+(a_xb_y-a_yb_x)(\overrightarrow{i}×\overrightarrow{j})\tag{合并}\\ &=(a_xb_x+a_yb_y)\overrightarrow{0}+(a_xb_y-a_yb_x)(\overrightarrow{i}×\overrightarrow{j}) \tag{消去平行向量}\\ &=(a_xb_y-a_yb_x)(\overrightarrow{i}×\overrightarrow{j}) \tag{消去 0 向量} \end{align} $$

仔细观察上式,得出:

  1. $a_xb_y-a_yb_x>0$,则 $\overrightarrow{b}$ 在 $\overrightarrow{a}$ 的逆时针方向上(参照 $\overrightarrow{i}$ 和 $\overrightarrow{j}$ 的位置);
  2. $a_xb_y-a_yb_x<0$,则 $\overrightarrow{b}$ 在 $\overrightarrow{a}$ 的顺时针方向上;
  3. $a_xb_y-a_yb_x=0$,则 $\overrightarrow{a}$ 和 $\overrightarrow{b}$ 共线,但是否同向不确定。

二:向量积与计算几何(线段相交)

在计算几何中,向量积也有很大的应用,最经典的就是"线段相交"问题:给定两条线段 $AB$ 和 $CD$,其中 $A(a_x, a_y),B(b_x, b_y),C(c_x, c_y),D(d_x, d_y)$,判断它们是否相交。

我们只需利用两个实验即可完成判断,快速排斥实验和跨立实验。

先看跨立实验,若两条线段相交,则两线段必然会相互跨立对方。

代码如下:

struct Point
{
    double x, y;
};

/* 求出向量 AB 与向量 AC 的向量积,返回 0 代表共线 */
double cross(Point A, Point B, Point C)
{
    return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x);
}

if (cross(A, B, C) * cross(A, B, D) <= 0 &&
    cross(C, D, A) * cross(C, D, B) <= 0)
    pass test;

但是请考虑一个特殊情况,若两条线段共线呢?此时仅仅通过上面的跨立实验是无法检测的,因此还需要借助快速排斥实验。

所谓的快速排斥实验,就是以线段 $AB$,$CD$ 为对角线作两个矩形,若两个矩形外离,则两条线段必定不相交。

代码如下:

if (min(A.x, B.x) <= max(C.x, D.x) &&
    min(C.x, D.x) <= max(A.x, B.x) &&
    min(A.y, B.y) <= max(C.y, D.y) &&
    min(C.y, D.y) <= max(A.y, B.y))
    pass test;

综合上面的分析,若要断定两条线段相交,则它们必须同时通过快速排斥实验和跨立实验。

全部代码如下:

struct Point
{
    double x, y;
};

/* 求出向量 AB 与向量 AC 的向量积,返回 0 代表共线 */
double cross(Point A, Point B, Point C)
{
    return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x);
}

/* 判断线段 AB 与线段 CD 是否相交,相交返回 true */
bool is_intersect(Point A, Point B, Point C, Point D)
{
    if (min(A.x, B.x) <= max(C.x, D.x) &&
        min(C.x, D.x) <= max(A.x, B.x) &&
        min(A.y, B.y) <= max(C.y, D.y) &&
        min(C.y, D.y) <= max(A.y, B.y) &&
        cross(A, B, C) * cross(A, B, D) <= 0 &&
        cross(C, D, A) * cross(C, D, B) <= 0)
        return true;

    return false;
}

三:参考文献

一:背景介绍

在一堆数中求其前 k 大或前 k 小的问题,简称 TOP-K 问题。而目前解决 TOP-K 问题最有效的算法即是 BFPRT 算法,又称为中位数的中位数算法,该算法由 Blum、Floyd、Pratt、Rivest、Tarjan 提出,最坏时间复杂度为 $O(n)$。

在首次接触 TOP-K 问题时,我们的第一反应就是可以先对所有数据进行一次排序,然后取其前 k 即可,但是这么做有两个问题:

  1. 快速排序的平均复杂度为 $O(nlogn)$,但最坏时间复杂度为 $O(n^2)$,不能始终保证较好的复杂度;
  2. 我们只需要前 k 大的,而对其余不需要的数也进行了排序,浪费了大量排序时间。

除这种方法之外,堆排序也是一个比较好的选择,可以维护一个大小为 k 的堆,时间复杂度为 $O(nlogk)$。

那是否还存在更有效的方法呢?我们来看下 BFPRT 算法的做法。

在快速排序的基础上,首先通过判断主元位置与k的大小使递归的规模变小,其次通过修改快速排序中主元的选取方法来降低快速排序在最坏情况下的时间复杂度

下面先来简单回顾下快速排序的过程,以升序为例:

  1. 选取主元;
  2. 以选取的主元为分界点,把小于主元的放在左边,大于主元的放在右边;
  3. 分别对左边和右边进行递归,重复上述过程。

二:算法过程及代码

BFPRT 算法步骤如下:

  1. 选取主元;
    1.1. 将 n 个元素按顺序分为 $⌊\frac n5⌋$ 个组,每组 5 个元素,若有剩余,舍去;
    1.2. 对于这 $⌊\frac n5⌋$ 个组中的每一组使用插入排序找到它们各自的中位数;
    1.3. 对于 1.2 中找到的所有中位数,调用 BFPRT 算法求出它们的中位数,作为主元;
  2. 以 1.3 选取的主元为分界点,把小于主元的放在左边,大于主元的放在右边;
  3. 判断主元的位置与 k 的大小,有选择的对左边或右边递归。

上面的描述可能并不易理解,先看下面这幅图:

BFPRT() 调用 GetPivotIndex() 和 Partition() 来求解第 k 小,在这过程中,GetPivotIndex() 也调用了 BFPRT(),即 GetPivotIndex() 和 BFPRT() 为互递归的关系。

下面为代码实现,其所求为前 k 小的数

#include <iostream>
#include <algorithm>

using namespace std;

int InsertSort(int array[], int left, int right);
int GetPivotIndex(int array[], int left, int right);
int Partition(int array[], int left, int right, int pivot_index);
int BFPRT(int array[], int left, int right, int k);

int main()
{
    int k = 8; // 1 <= k <= array.size
    int array[20] = { 11,9,10,1,13,8,15,0,16,2,17,5,14,3,6,18,12,7,19,4 };

    cout << "原数组:";
    for (int i = 0; i < 20; i++)
        cout << array[i] << " ";
    cout << endl;

    // 因为是以 k 为划分,所以还可以求出第 k 小值
    cout << "第 " << k << " 小值为:" << array[BFPRT(array, 0, 19, k)] << endl;

    cout << "变换后的数组:";
    for (int i = 0; i < 20; i++)
        cout << array[i] << " ";
    cout << endl;

    return 0;
}

/**
 * 对数组 array[left, right] 进行插入排序,并返回 [left, right]
 * 的中位数。
 */
int InsertSort(int array[], int left, int right)
{
    int temp;
    int j;

    for (int i = left + 1; i <= right; i++)
    {
        temp = array[i];
        j = i - 1;

        while (j >= left && array[j] > temp)
        {
            array[j + 1] = array[j];
            j--;
        }

        array[j + 1] = temp;
    }

    return ((right - left) >> 1) + left;
}

/**
 * 数组 array[left, right] 每五个元素作为一组,并计算每组的中位数,
 * 最后返回这些中位数的中位数下标(即主元下标)。
 *
 * @attention 末尾返回语句最后一个参数多加一个 1 的作用其实就是向上取整的意思,
 * 这样可以始终保持 k 大于 0。
 */
int GetPivotIndex(int array[], int left, int right)
{
    if (right - left < 5)
        return InsertSort(array, left, right);

    int sub_right = left - 1;

    // 每五个作为一组,求出中位数,并把这些中位数全部依次移动到数组左边
    for (int i = left; i + 4 <= right; i += 5)
    {
        int index = InsertSort(array, i, i + 4);
        swap(array[++sub_right], array[index]);
    }

    // 利用 BFPRT 得到这些中位数的中位数下标(即主元下标)
    return BFPRT(array, left, sub_right, ((sub_right - left + 1) >> 1) + 1);
}

/**
 * 利用主元下标 pivot_index 进行对数组 array[left, right] 划分,并返回
 * 划分后的分界线下标。
 */
int Partition(int array[], int left, int right, int pivot_index)
{
    swap(array[pivot_index], array[right]); // 把主元放置于末尾

    int partition_index = left; // 跟踪划分的分界线
    for (int i = left; i < right; i++)
    {
        if (array[i] < array[right])
        {
            swap(array[partition_index++], array[i]); // 比主元小的都放在左侧
        }
    }

    swap(array[partition_index], array[right]); // 最后把主元换回来

    return partition_index;
}

/**
 * 返回数组 array[left, right] 的第 k 小数的下标
 */
int BFPRT(int array[], int left, int right, int k)
{
    int pivot_index = GetPivotIndex(array, left, right); // 得到中位数的中位数下标(即主元下标)
    int partition_index = Partition(array, left, right, pivot_index); // 进行划分,返回划分边界
    int num = partition_index - left + 1;

    if (num == k)
        return partition_index;
    else if (num > k)
        return BFPRT(array, left, partition_index - 1, k);
    else
        return BFPRT(array, partition_index + 1, right, k - num);
}

运行如下:

原数组:11 9 10 1 13 8 15 0 16 2 17 5 14 3 6 18 12 7 19 4
第 8 小值为:7
变换后的数组:4 0 1 3 2 5 6 7 8 9 10 12 13 14 17 15 16 11 18 19

三:时间复杂度分析

BFPRT 算法在最坏情况下的时间复杂度是 $O(n)$,下面予以证明。令 $T(n)$ 为所求的时间复杂度,则有:

$$ T(n)≤T(\frac n 5)+T(\frac {7n}{10})+c⋅n\tag{c 为一个正常数} $$

其中:

  • $T(\frac n 5)$ 来自 GetPivotIndex(),n 个元素,5 个一组,共有 $⌊\frac n5⌋$ 个中位数;
  • $T(\frac {7n}{10})$ 来自 BFPRT(),在 $⌊\frac n5⌋$ 个中位数中,主元 x 大于其中 $\frac 12⋅\frac n5=\frac n{10}$ 的中位数,而每个中位数在其本来的 5 个数的小组中又大于或等于其中的 3 个数,所以主元 x 至少大于所有数中的 $\frac n{10}⋅3=\frac {3n}{10}$ 个。即划分之后,任意一边的长度至少为 $\frac 3{10}$,在最坏情况下,每次选择都选到了 $\frac 7{10}$ 的那一部分。
  • $c⋅n$ 来自其它操作,比如 InsertSort(),以及 GetPivotIndex() 和 Partition() 里所需的一些额外操作。

设 $T(n)=t⋅n$,其中 t 为未知,它可以是一个正常数,也可以是一个关于 n 的函数,代入上式:

$$ \begin{align} t⋅n&≤\frac {t⋅n}5+\frac{7t⋅n}{10}+c⋅n \tag{两边消去 n}\\ t&≤\frac t 5+\frac {7t}{10}+c \tag{再化简}\\ t&≤10c \tag{c 为一个正常数} \end{align} $$

其中 c 为一个正常数,故t也是一个正常数,即 $T(n)≤10c⋅n$,因此 $T(n)=O(n)$,至此证明结束。

接下来我们再来探讨下 BFPRT 算法为何选 5 作为分组主元,而不是 2, 3, 7, 9 呢?

首先排除偶数,对于偶数我们很难取舍其中位数,而奇数很容易。再者对于 3 而言,会有 $T(n)≤T(\frac n 3)+T(\frac {2n}3)+c⋅n$,它本身还是操作了 n 个元素,与以 5 为主元的 $\frac {9n}{10}$ 相比,其复杂度并没有减少。对于 7,9,... 而言,上式中的 10c,其整体都会增加,所以与 5 相比,5 更适合。

四:参考文献