二维情形下的最接近点对问题

在上一篇中我们讨论了一维情形下的最近点对问题，本篇我们类比一维的思路讨论二维情形下的最接近点对问题。
将以上过程推广到二维最接近点对问题，设S中的点为平面上的点，它们都有2个坐标值x和y。同样我们将平面上的点集S线性分割为大小大致相等的2个子集S1和S2，我们选取一垂直线l:x=m来作为分割直线。递归地在S1和S2上解最接近点对问题，我们分别得到S1和S2中的最小距离d1和d2。令d=min(d1,d2)若S的最接近点对(p,q)之间的距离小于d,那么这两个点必定位于直线l的两侧，此时，S1中的所有点与S2中的所有点构成的点对均为最接近点对的候选者。在最坏情况下有(n*n)/4对这样的候选者。但是S1和S2中的点具有以下的稀疏性质，它使我们不必检查所有这n^2/4对候选者。考虑S1中任意一点p,它若与S2中的点q构成最接近点对的候选者，则必有d(p,q)小于d,满足这个条件的P2中的点有多少个呢?
容易知道这样的点一定落在一个d×2d的矩形R中。由d的意义可知S2中任何2个S中的点的距离都不小于d。由此可以推出矩形R中最多只有6个S中的点。
因此，在分治法的合并步骤中，我们最多只需要检查6×n/2=3n对候选者，而不是n^2/4对候选者。这是否就意味着我们可以在O(n)时间内完成分治法的合并步骤呢？现在还不能作出这个结论，因为我们只知道对于S1中的每个点p最多只需要检查S2中的6个点，但是我们并不确切地知道要检查哪6个点。为了解决这个问题，我们可以将p和P2中所有S2的点投影到垂直线l上。由于能与p点一起构成最接近点对候选者的S2中点一定在矩形R中，所以它们在直线l上的投影点距p在l上投影点的距离小于d。由上面的分析可知，这种投影点最多只有6个。因此，若将S1和S2中所有S的点按其y坐标排好序，则对S1中所有点p，对排好序的点列作一次扫描，就可以找出所有最接近点对的候选者，对S1中每一点最多只要检查S2中排好序的相继6个点。
下面贴代码：

#include <algorithm>#include <cmath>#include <cstdio>#include <cstdlib>#include <cstring>#include <ctime>#include <ctype.h>#include <iostream>#include <map>#include <queue>#include <set>#include <stack>#include <string>#include <vector>#define eps 1e-8#define INF 0x7fffffff#define PI acos(-1.0)#define seed 31//131,1313typedef long long LL;typedef unsigned long long ULL;using namespace std;const int maxn=100005;//分治算法求二维最近点对struct Point{  double x,y;}p[maxn];int a[maxn];int cmpx(Point a,Point b){  return a.x<b.x;}int cmpy(int a,int b){  return p[a].y<p[b].y;}inline double dis(Point a,Point b){   return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}double closest(int low,int high){   int i,j,k;   if(low+1==high){  //只有两个点    return dis(p[low],p[high]);   }   if(low+2==high){ //只有三个点     return min(dis(p[low],p[high]),min(dis(p[low],p[low+1]),dis(p[low+1],p[high])));   }   int mid=(low+high)/2; //求中点即左右子集的分界线   double d=min(closest(low,mid),closest(mid+1,high));   for(i=low,k=0;i<=high;i++){ //把x坐标在p[mid].x-d  ~  p[mid].x+d范围内的点筛选出来     if(p[i].x>=p[mid].x-d&&p[i].x<=p[mid].x+d){        a[k++]=i; //保存这些点的下标索引     }   }   sort(a,a+k,cmpy); //按y坐标进行升序排序   for(i=0;i<k;i++){    for(j=i+1;j<k;j++){       if(p[a[j]].y-p[a[i]].y>=d) //注意下标索引            break;       d=min(d,dis(p[a[i]],p[a[j]]));    }   }   return d;}int main(){    int i,n;    while(scanf("%d",&n)!=EOF){        if(n==0) break;        for(i = 0 ; i < n ; ++i)            scanf("%lf %lf",&p[i].x,&p[i].y);        sort(p , p + n , cmpx);//按x坐标进行升序排序        printf("%.2lf\n",closest(0 , n - 1)/2);//最近点对间的距离    }    return 0;}