前置知识：平面最近点对

考虑以下分治算法：

设平面上的点都在点集S中，为了将S线性分割为大小大致相等的2个子集S1和S2，我们选取一垂直线l（方程：x=m）来作为分割直线。其中m为S中各点x坐标的中位数。由此将S分割为S1={p∈S|px≤m}和S2={p∈S|px>m}。从而使S1和S2分别位于直线l的左侧和右侧，且S=S1∪S2 。由于m是S中各点x坐标值的中位数，因此S1和S2中的点数大致相等。 递归地在S1和S2上解最接近点对问题，我们分别得到S1和S2中的最小距离δ1和δ2。现设δ=min (δ1,δ2)。

若S的最接近点对(p,q)之间的距离d(p,q)<δ则p和q必分属于S1和S2。不妨设p∈S1，q∈S2。那么p和q距直线l的距离均小于δ。因此，我们若用P1和P2分别表示直线l的左边和右边的宽为δ的2个垂直长条，则p∈S1，q∈S2

此时，P1中所有点与P2中所有点构成的点对均为最接近点对的候选者。在最坏情况下有n^2/4对这样的候选者。但是P1和P2中的点具有以下的稀疏性质，它使我们不必检查所有这n2/4对候选者。考虑P1中任意一点p,它若与P2中的点q构成最接近点对的候选者，则必有d(p,q)<δ。满足这个条件的P2中的点有多少个呢？容易看出这样的点一定落在一个δ×2δ的矩形R中，

由δ的意义可知P2中任何2个S中的点的距离都不小于δ。由此可以推出矩形R中最多只有6个S中的点。事实上，我们可以将矩形R的长为2δ的边3等分，将它的长为δ的边2等分，由此导出6个（δ/2）×（2δ/3）的矩形。

因此d(u,v)≤5δ/6<δ 。这与δ的意义相矛盾。也就是说矩形R中最多只有6个S中的点。由于这种稀疏性质，对于P1中任一点p，P2中最多只有6个点与它构成最接近点对的候选者。因此，在分治法的合并步骤中，我们最多只需要检查6×n/2=3n对候选者，而不是n^2/4对候选者。

我们只知道对于P1中每个S1中的点p最多只需要检查P2中的6个点，但是我们并不确切地知道要检查哪6个点。为了解决这个问题，我们可以将p和P2中所有S2的点投影到垂直线l上。由于能与p点一起构成最接近点对候选者的S2中点一定在矩形R中，所以它们在直线l上的投影点距p在l上投影点的距离小于δ。由上面的分析可知，这种投影点最多只有6个。因此，若将P1和P2中所有S的点按其y坐标排好序，则对P1中所有点p，对排好序的点列作一次扫描，就可以找出所有最接近点对的候选者，对P1中每一点最多只要检查P2中排好序的相继6个点。

至此，我们用分治法求出平面最接近点对。

清爽代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define N 200010
#define INF 0x3f3f3f3f

using namespace std;

struct node
{
    double x,y;
    int f;
};
node a[N],b[N];

bool cmp(node x,node y)
{
    if (x.x==y.x) return x.y<y.y;
    return x.x<y.x;
}

double calc(int x,int y)
{
    if (a[x].f==a[y].f) return INF;//与平面最近点对的唯一区别
    double x1=(a[x].x-a[y].x)*(a[x].x-a[y].x);
    double y1=(a[x].y-a[y].y)*(a[x].y-a[y].y);
    return sqrt(x1+y1);
}

int temp[N];

double devide(int l,int r)
{
    double d=INF;
    if (l==r) return d;
    if (l+1==r) return calc(l,r);
    int mid=(l+r)>>1;
    double d1=devide(l,mid);
    double d2=devide(mid+1,r);
    if (d1>d2) d=d2;
    else d=d1;
    int k=0;
    for (int i=l;i<=r;i++)
        if (fabs(a[i].x-a[mid].x)<=d)
            temp[++k]=i;
            //这边存到一个临时数组，直接n^2查找最近值，因为排序和没排差不多
            //(可以理性感受最多区域内只有6个点)
    for (int i=1;i<=k;i++)
        for (int j=i+1;j<=k;j++)
        {
            double d1=calc(temp[i],temp[j]);
            if (d>d1) d=d1;
        }
    return d;
}

int main()
{
    int T;
    cin>>T;
    while (T--)
    {
        int n;
        cin>>n;
        for (int i=1;i<=n;i++)
        {
            cin>>a[i].x>>a[i].y;
            a[i].f=0;
        }
        for (int i=n+1;i<=n*2;i++)
        {
            cin>>a[i].x>>a[i].y;
            a[i].f=1;
        }
        n*=2;
        sort(a+1,a+n+1,cmp);
        printf("%.3lf\n",devide(1,n));
    }
    return 0;
}

Reference: https://www.luogu.com.cn/blog/user277/solution-p1429