本来这个算法在笔者电脑里无人问津过一段时间了,但今天正好做HDU 1007见到了这个问题,今天就来把代码分享出来吧!
我们首先将所有点按照坐标x排序一下,再做一条直线l当作“分割线”,方便我们递归。
然后,我们就可以把这些点按照x轴的坐标分为左半部分和右半部分。那么最短距离一定在左半部分、右半部分、跨越左右的点对中的一个。
那么你可能会有疑问了:本来最近点对也一定在这三个区域内,这不还是相当于什么都没干吗?
还真不是。我们可以假设通过递归得到了左边最小距离为d1,右边最小距离为d2,令δ = min(d1,d2)
如图所示,如果跨越左右的点对可能是最短距离,那么它也必然比δ小。而在以l为中心、最大距离为2δ的区域中,最多有O(n)个如图所示的矩形。另外,可以证明对于每个矩形区域,最多尝试8个点对一定能找到最短距离(算法导论第33.4节有详细的证明,这里不再赘述)。
因此,我们可以写出递归式:T(n)=2T(n/2)+O(n),可以用主项定理(master method)解得时间复杂度T(n)=O(nlogn)。加上排序一次的时间O(nlogn),因此整个算法的运行时间T(n)' = T(n)+O(nlogn) = O(nlogn)。
下面,通过这个算法,我们就可以写出一份代码来:
/**
* Find closest distance in N points.
* Time cost: O(nlogn)
* Author: Zheng Chen / Arclabs001
* Copyright 2015 Xi'an University of Posts & Telecommunications
*/
/**
* Algorithm:
* First of all, sort the points in ascending order by x.
* Divide the array into 2 parts, and find the smallest distance within two parts.
* Let min as the smaller one, and find the smallest split distance.
*/
#include <iostream>
#include <algorithm>
#include <cmath>
#define INF 0x6FFFFFFF
using namespace std;
struct Node{
double x, y;
friend bool operator < (const Node &a, const Node &b){
if(a.x == b.x)
return a.y < b.y;
return a.x < b.x;
}
};
Node* Point = NULL;
/**
* Calculate the distance between two points.
* @return [double, the distance between a and b]
*/
double _distance(const Node a, const Node b)
{
return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}
double smaller(double p, double q)
{
return (p > q) ? q : p;
}
/**
* [Find the closest distance, divide & conquer]
* @param left [search from where]
* @param right [to where]
* @return [double, the smallest distance in the whole set of points]
*/
double Closest_distance(int left, int right)
{
double d = INF;
double distance_tmp;
if(left == right)
return 0;
if(right == left+1)
return _distance( Point[left], Point[right] );
int mid = (left + right) / 2;
d = smaller( Closest_distance(left,mid) , Closest_distance(mid,right) );
for(int i=mid-1; i>=left && Point[mid].x - Point[i].x < d; i--){
for(int j = mid+1; j<=right && Point[j].x - Point[mid].x < d && fabs( Point[i].y - Point[j].y) < d; j++){
distance_tmp = _distance( Point[i], Point[j] );
if(distance_tmp < d)
d = distance_tmp;
}
}
return d;
}
int main()
{
int n;
cin>>n;
Point = new Node[n];
for(int i=0; i<n ; i++){
cin>>Point[i].x>>Point[i].y;
}
sort(Point,Point+n);
cout<<Closest_distance(0,n-1)<<endl;
return 0;
}
当然,直接套用这个代码提交到HDOJ会导致Exceed Time Limit,你懂的^_^