目录见右上
站内文章欢迎各位转载,但是请注明出处(毕竟码字很辛苦哈)。

主题

寻找最近点对

这里的最近点对指的是两点之间的欧氏距离(哈哈哈事实上就是我们现在采用的平面几何中的距离,也就是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.


评论
热度(1)
©主题 —— Powered by LOFTER