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

在上一篇中我们讨论了一维情形下的最近点对问题，本篇我们类比一维的思路讨论二维情形下的最接近点对问题。
将以上过程推广到二维最接近点对问题，设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,1313
typedef 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;
}