寻找最近点对
这里的最近点对指的是两点之间的欧氏距离(哈哈哈事实上就是我们现在采用的平面几何中的距离,也就是sqrt((x1-x2)^2+(y1-y2)^2))
讨论一个很多人写了但是还是要写的n(logn)^2的分治算法。
在《算法导论》中学会的,我尝试以尽量通俗的话来说明这个算法。
考虑二分。
我们视xi为横坐标,yi为纵坐标。
(放张图避免看得很累)
假设当前我们要处理区间[l..r]的问题。
如果[l..r]只包含了3个或以下的点——就直接暴力求解就好了。
如果没有,那么——首先画一条垂直x轴的线,把[l..r]的点分为两份,左边的点集为[l..m],记为XL。右边的点集为[m+1..r],记为XR。
设p1与p2是[l..r]内构成最短距离的两个点。那么只会有三种情况:
p1与p2均在XL中
p1与p2均在XR中
p1与p2一个在XL中,一个在XR中。
前两种情况比较好解决,递归求解。现在关键是第三种情况。
假设已经解出XL这一部分点集的最短距离minXL与XR这一部分的最短距离minXR。令minLR=min{minXL,minXR}
设点p1(与之前p1的定义不同了)在XL内,p2在XR内。p1与p2的距离是所有 “一个在XL内,另一个在XR内” 的点对中 距离最小的点对的距离。
设p1与p2的距离为dis(p1,p2)
那么现在答案只有两种情况:
dis(p1,p2)
minLR
假设答案是dis(p1,p2)
那么dis(p1,p2)<minLR,所以p1与p2横坐标之差不会超过minLR,纵坐标亦然。所以p1与p2会在一个如图所示的矩形内(这是最大边界)。中间那条线是我们之前画出的分离XL与XR的线。
思考便可知,左边这个正方形里面每一个之间的距离都应该不小于minLR。所以左边这个正方形里面最多只会有4个点。右边这个正方形亦然,也最多只有4个点。所以整个矩形内不会多于8个点。
我们是知道了这个矩形的(长)——况且画出的线是这个矩形的中线。所以我们是确定了矩形内的点的横坐标的范围。把在这个范围内的点划为另一个点集YM。
如果YM内的点是按照纵坐标(升序)排序的,由于一个矩形最多只有8个点,所以我们认为存在dis(p1,p2)<minLR,那么只需要考虑p1与p1之后的7个点的距离是否小于minLR。只有这7个点可能成为p2。
接下来分析时间复杂度?
首先按照横坐标x轴排一趟序O(nlogn)
对于画线,我们可以用一个二分查找(O(logn))。
接下来划分元素,这个也是O(n)可以做到的(具体说可能不用真正去“分”,所以时间复杂度可以视为O(1)不管)。
把横坐标在矩形内的点加入YM也是可以O(n)做到的。
最后一步按照纵坐标Y找7个点——意味着要重新排序。许多博客认为是采用了重新排序的方法,这样使得算法复杂度达到了O(n*(logn)^2)。
显然这是一个O(n*(logn)^2)的算法。这是很优的了。
我相信很多人都是博客教的(maybe学长教的)不是原著教的:
根据原著来看,算法是一个O(nlogn)可以达到的。所以最后的问题就在于如何把最后一步——按纵坐标Y找7个点。
原著中只是提及了预排序——但是并没有提及如何操作。个人想到了两种方法——
第一种方法:
首先按照横坐标x轴排一趟序O(nlogn),记录在A中。
接着按照纵坐标Y轴排一趟序O(nlogn),记录在B中。
我们按照横坐标X轴来进行分治操作。每一次在划分元素的时候讲元素按照在B中位置放置在YM中。这时YM是一个中间有空格的数组。之后扫一遍把空格填好。每一层操作复杂度为O(n)。
之后就在YM中每个点找7个点去计算答案就好,每一层操作复杂度O(n)。
后来个人发现——事实上这不是O(nlogn)的算法。因为这不是每一层O(n),而是每一个递归程序O(n),一共有2^n个递归程序,所以总的时间复杂度为O(n*2^n)。
所以还是抛弃它。
第二种方法:
开始不按照X轴排序,只按照Y轴排序O(nlogn)。
我们按照纵坐标Y轴来进行分治操作。每一次(层)划分元素先O(n)找到要画的线,之后把元素分出去,记录在YM中。
之后(每一层)再归并一次O(n)变为按照Y轴排序,将其储存在MYY中。
之后就在MYY中每个点找7个点去计算答案,每一层操作复杂度O(n)。
注意当这一层求出答案之后,要把数组重新变为为按照X轴排序的(可以用一个辅助数组记录),否则第三步不能O(n)归并。这个操作也是O(n)的。
所以总操作复杂度为O(nlogn)。
这个方法关键在于如何确定画的线。如果O(n)计算平均数显然可以会被卡得很死,而寻找中位数也不太方便(离散数组)。
可以参考随机快排的思想,随机取一个点作为线,之后划分。不过对于边界情况要思考一下,可能会有所有的点横坐标都是一个数,为了算法性能平均性,不能单纯的用“<”或者“>”。
接着还有一个很无奈的就是空间复杂度。每一个递归程序我们需要一个还原数组(最后一步要还原数组)。所以为了防止递归栈溢出,只能手动模拟栈。
代码复杂度极高。建议还是不要写了。
所以最后我也是向时间复杂度是O(n(logn)^2)低头了。不过这个算法也很优了。10^4的数据,比起暴力大概快10倍。
代码如下:
type
point=record
x,y:longint;//记录点的信息
end;
var
A,Temp:array[1..100000] of point;//A数原数组,Temp是辅助数组(文中的YM数组)。
n:longint;
ans:double;//答案
function Dis(x1,x2,y1,y2:longint):Double;//计算距离(不用sqrt,最后输出答案的时候用sqrt可以少一点常数)
begin
exit(sqrt(sqr(x1-x2)+sqr(y1-y2)));
end;
function Doublemin(a,b:double):double;//两个double取较小值
begin
if a<b then exit(a);
exit(b);
end;
procedure swap(var a,b:longint);//交换两个整数
var
c:longint;
begin
c:=a;
a:=b;
b:=c;
end;
procedure Linit;//初始化
var
vi:longint;
begin
read(n);
for vi:=1 to n do
read(A[vi].x,A[vi].y);
randomize;
end;
procedure Xsort(l,r:longint);//将A按照x轴排序
var
tl,tr,x:longint;
begin
tl:=l;
tr:=r;
x:=A[l+random(r-l+1)].x;
repeat
while A[tl].x<x do inc(tl);
while A[tr].x>x do dec(tr);
if tl<=tr then
begin
swap(A[tl].x,A[tr].x);
swap(A[tl].y,A[tr].y);
inc(tl);
dec(tr);
end;
until tl>tr;
if l<tr then Xsort(l,tr);
if tl<r then Xsort(tl,r);
end;
procedure Ysort(l,r:longint);//将Temp按照Y轴排序
var
tl,tr,y:longint;
begin
tl:=l;
tr:=r;
y:=Temp[l+random(r-l+1)].y;
repeat
while Temp[tl].y<y do inc(tl);
while Temp[tr].y>y do dec(tr);
if tl<=tr then
begin
swap(Temp[tl].x,Temp[tr].x);
swap(Temp[tl].y,Temp[tr].y);
inc(tl);
dec(tr);
end;
until tl>tr;
if l<tr then Ysort(l,tr);
if tl<r then Ysort(tl,r);
end;
function Division(l,r:longint):double;//分治
var
vi,vj,midline,Arrline,tl,tr:longint;//画的线是取的是点Arrline的横坐标,midline是点Arrline的横坐标的值。tl与tr记录递归之后需要考虑的跨边界的点对的范围
now:double;//当前这一个递归程序的答案
begin
now:=1 shl 30;//先把now变为一个很大的数
if r-l+1<=3 then//如果元素个数不多于3,就直接暴力求解
begin
for vi:=l to r do
for vj:=vi+1 to r do
now:=doublemin(now,dis(A[vi].x,A[vj].x,A[vi].y,A[vj].y));
exit(now);
end;
Arrline:=l+random(r-l+1);//随便找一个画的线(二分也可以)
Midline:=A[Arrline].x;
now:=doublemin(now,Division(l,Arrline));
now:=doublemin(now,Division(Arrline+1,r));//分治,答案放在now中
tl:=l;
tr:=r;
for vi:=Arrline downto l do
if A[vi].x>=Midline-now then//找到tl与tr,确定考虑 哪些点成为跨边界的点对 去更新答案
begin
Temp[vi].x:=A[vi].x;
Temp[vi].y:=A[vi].y;
tl:=vi;
end
else
break;
for vi:=Arrline+1 to r do
if A[vi].x<=Midline+now then
begin
Temp[vi].x:=A[vi].x;
Temp[vi].y:=A[vi].y;
tr:=vi;
end
else break;
Ysort(tl,tr);//将Temp按Y轴排序
for vi:=tl to tr do
for vj:=1 to 7 do
begin
if vi+vj>tr then break;
now:=doublemin(now,dis(A[vi].x,A[vi+vj].x,A[vi].y,A[vi+vj].y));//更新now
end;
exit(now);//返回答案
end;
begin
Linit;//初始化
Xsort(1,n);//将A按照X轴排序
ans:=Division(1,n);//计算答案
writeln(ans:0:6);//输出
end.