分治法求解平面最近点对问题

问题描述

给出一个平面内的若干个点，问这些点两两之间的最近距离是多少
弱数据题目链接

• 暴力自不必说，两个$for$挨个看，时间复杂度是$O(n^2)$，可以通过数据较小的测试

分治求解

数据加强题目链接

• 这里我们可以回忆一下归并排序是怎么做的，如何达到$O(nlogn)$的时间复杂度，采用的即是分治思想，也就是把一组数据对半分成两部分，继续划分，直至最小，排序拼回原数组，再排序，再拼回，直至全部排序完成

• 那么借助这个想法，我们能不能把平面上的点分成两部分分而治之呢？
  

• 我们想求上图任意两个点之间的最近点对，首先，要对这些点进行递归划分，根据分治的思想，我们先将其分成两部分，然后递归下去，按横坐标从小到大进行划分，直到最后的一个点或者两个点，可以得到如下的划分图，线条长度代表递归的深度，长度越长递归深度越小
  

• 我们可以先考虑一下第一次划分之后，如何求最近点对，因为我们知道，这两个点要么都在左面，要么都在右面，要么一个在左面，另一个在右面，在同一个方向的时候只需要递归处理即可，关键是位居两侧时候怎么办，因为这时候我们已经求出左面和右面的最小值$dis$，那么位居两侧的可能点有哪些呢？我们应该把这些点筛选出来，如果存在这样的一个点对，则点对中的两个点和划分线之间的距离都应该在$dis$以内，也就是说如果$fabs(x_{mid}-x_1)<dis$，这里的$mid$是靠近中间划分线的点（黑书中写的是$\le$，似乎应该是小于，因为如果$x$的差值等于$dis$，那么这两个点之间距离不可能小于$dis$）,那么就把这个点加入到待选数组中

• 接下来待选数组中都是可能的选项，刚才已经按照横坐标筛选一遍了，如果暴力枚举两个点，有$T(n)=T(\frac{n}{2})+O(n^2)$，时间复杂度退化为$O(n^2)$，仍然会超时，黑书中介绍的方法是现在按照纵坐标筛选，先按纵坐标升序排列，然后如果纵坐标差值超过$dis$，那么剪枝，否则更新$dis$

• 算法导论中讲在枚举$tmp\_p$数组中元素时，只需要枚举每个点之后的七个点即可，这是从两个边长均为$dis$的正方形拼成长方形中点的最多数量这个角度考虑的（具体参考算法导论），那么这一次枚举时间复杂度就是线性的
  

• 例题程序

#include <iostream>
#include <cstring>
#include <algorithm>
#include <queue>
#include <stack>
#include <vector>
#include <cmath>
#include <cstdio>
#include <map>
#include <iomanip>
using namespace std;
typedef long long ll;
const int MAXN = 2e5 + 100;
const double INF = 1e20;
const double eps = 1e-8;
int Data[MAXN];
struct Point{
    double x, y;
};
Point p[MAXN], tmp_p[MAXN];
int sgn(double x){
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
double distance(Point A, Point B){
    return hypot(A.x - B.x, A.y - B.y);
}
double Closest_Pair(int l, int r, int n){
    double dis = INF;
    if(l == r) return dis;
    if(l + 1 == r) return distance(p[l], p[r]);
    int mid = ((r - l) >> 1) + l;
    double d1 = Closest_Pair(l, mid, n);
    double d2 = Closest_Pair(mid + 1, r, n);
    dis = min(d1, d2);
    int k = 0;
    for(int i=l;i<=r;i++){
        if(fabs(p[mid].x - p[i].x) < dis) tmp_p[k++] = p[i];
    }
    sort(tmp_p, tmp_p + k, [](Point A, Point B){
        return A.y < B.y;
    });
    for(int i=0;i<k;i++){
        for(int j=i+1;j<k;j++){
            if(distance(tmp_p[j], tmp_p[i]) >= dis) break;
            dis = min(dis, distance(tmp_p[j], tmp_p[i]));
        }
    }
    return dis;
}
int main(){
    ios::sync_with_stdio(false);
    int n;
    cin >> n;
    for(int i=0;i<n;i++){
        cin >> p[i].x >> p[i].y;
    }
    sort(p, p + n, [](Point A, Point B){
        if(A.x == B.x){
            return A.y < B.y;
        }return A.x < B.x;
    });
    cout << fixed << setprecision(4) << Closest_Pair(0, n - 1, n);
    return 0;
}

时间复杂度

• 参考算法导论，运行时间可以使用递归式$T(n)=T(\frac{n}{2})+O(n)$描述，所以$T(n)=O(nlogn)$，其中$T(n)$表示每一步递归所用的运行时间，算法每次递归都需要进行一次$nlogn$的排序，设$T^{\prime }(n)$表示整个算法的运行时间，则$T^{\prime }(n)=T(n)+O(nlogn)$，总的时间复杂度都是$O(nlogn)$