题目链接
P7883 平面最近点对（加强加强版）
题目背景
P1429 平面最近点对（加强版）里最高赞题解写道：
我们充分发扬人类智慧：
将所有点全部绕原点旋转同一个角度，然后按 \(x\) 坐标排序
根据数学直觉，在随机旋转后，答案中的两个点在数组中肯定不会离得太远
所以我们只取每个点向后的 \(5\) 个点来计算答案
这样速度快得飞起，在 \(n=1000000\) 时都可以在 1 s 内卡过
当然，这是错的。
题目描述
给定 \(n\) 个二维欧几里得平面上的点 \(p_1, p_2, \dots, p_n\)，请输出距离最近的两个点的距离。
输入格式
输入第一行为一个正整数 \(n\)，表示点数。
接下来 \(n\) 行，第 \(i\) 行为用空格隔开的整数 \(x_i, y_i\)，表示 \(p_i = (x_i, y_i)\)。
输入保证：没有两个坐标完全相同的点。
输出格式
输出一行，包含一个整数 \(D^2\)，表示距离最近的两个点的距离的平方。
由于输入的点为整点，因此这个值一定是整数。
样例 #1
样例输入 #1
2
-10000000 -10000000
10000000 10000000
样例输出 #1
800000000000000
样例 #2
样例输入 #2
5
1 1
1 9
9 1
9 9
0 10
样例输出 #2
2
提示
对于第二组样例，\((1, 9)\)、\((0, 10)\) 两个点最近，距离为 \(\sqrt 2\)，因此你需要输出 \(2\)。
数据范围
对于 \(100 \%\) 的数据，\(2 \leq n \leq 4 \times 10^5\)，\(-10^7 \leq x_i, y_i \leq 10^7\)。
本题目标复杂度是 \(O(n \log ^2 n)\)。设置 350ms 时限的原因是：
- \(O(n \log ^2 n)\) 参考代码使用
cin不会 TLE。最快的 std 能 \(<\) 100ms。 - @wlzhouzhuan 的程序能恰好在 350ms 内跑 \(1000n\) 次检查。
- 150 组测试数据，为了防止卡评测。
解题思路
kdtree
kdtree 是一种能够高效处理 \(k\) 维空间的数据结构，当节点数远大于 \(2^k\) 时，kdtree 往往能取得很好的效果
以 \(2\) 维空间为例，按照某种划分原则以及选择某个维度（\(x\) 或 \(y\)），将某个点作为 kdtree 的根节点，假设选择的维度是 \(x\)，则所有小于 \(x\) 的点作为左子树，其余点作为右子树，同理递归划分，一般选择中点进行划分，这样可以使左子树和右子树的节点数量尽量相等，使得树高为 \(O(logn)\) 级别，同时为了使 kdtree 上的节点的某一个维度不那么集中（主要是为了加快查询等操作），通常会选择方差较大的维度进行划分
本题可以枚举每个点，找到与其距离最近的点，由于 kdtree 上每个以节点为根的子树都可以表示为一个矩形，所以求出查询点和矩形的最近距离，如果该最近距离仍大于需要更新的答案的话，这棵子树中肯定没有更优解，即可以不用遍历该子树，另外可以进行启发式搜索：两棵子树中选择最近距离较近的子树搜索
对于矩阵查询，2dtree 单次时间复杂度最优 \(O(logn)\)，最坏 \(O(\sqrt{n})\)，扩展到 \(k\) 维，最坏时间复杂度：\(O(n^{1-\frac{1}{k}})\)
- 最坏时间复杂度：\(O(n^2)\)
代码
// Problem: P7883 平面最近点对（加强加强版）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P7883
// Memory Limit: 512 MB
// Time Limit: 350 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=4e5+5;
int n;
LL res=2e18;
struct Point
{
int x,y;
}point[N];
struct Tr
{
int s[2];
int U,D,L,R;
}tr[N];
inline bool cmp1(Point &x,Point &y)
{
return x.x<y.x;
}
inline bool cmp2(Point &x,Point &y)
{
return x.y<y.y;
}
void pushup(int x)
{
tr[x].L=tr[x].R=point[x].x;
tr[x].U=tr[x].D=point[x].y;
if(tr[x].s[0])
{
tr[x].L=min(tr[x].L,tr[tr[x].s[0]].L);
tr[x].D=min(tr[x].D,tr[tr[x].s[0]].D);
tr[x].R=max(tr[x].R,tr[tr[x].s[0]].R);
tr[x].U=max(tr[x].U,tr[tr[x].s[0]].U);
}
if(tr[x].s[1])
{
tr[x].L=min(tr[x].L,tr[tr[x].s[1]].L);
tr[x].D=min(tr[x].D,tr[tr[x].s[1]].D);
tr[x].R=max(tr[x].R,tr[tr[x].s[1]].R);
tr[x].U=max(tr[x].U,tr[tr[x].s[1]].U);
}
}
LL dist(int x,int y)
{
return (LL)(point[x].x-point[y].x)*(point[x].x-point[y].x)+(LL)(point[x].y-point[y].y)*(point[x].y-point[y].y);
}
LL get(int x,int y)
{
LL res=0;
if(point[y].x<tr[x].L)res+=(LL)(tr[x].L-point[y].x)*(LL)(tr[x].L-point[y].x);
if(point[y].x>tr[x].R)res+=(LL)(tr[x].R-point[y].x)*(LL)(tr[x].R-point[y].x);
if(point[y].y<tr[x].D)res+=(LL)(tr[x].D-point[y].y)*(LL)(tr[x].D-point[y].y);
if(point[y].y>tr[x].U)res+=(LL)(tr[x].U-point[y].y)*(LL)(tr[x].U-point[y].y);
return res;
}
void ask(int l,int r,int x)
{
if(l>r)return ;
int mid=l+r>>1;
if(mid!=x)res=min(res,dist(x,mid));
if(l==r)return ;
LL distL=get(tr[mid].s[0],x),distR=get(tr[mid].s[1],x);
if(distL<res&&distR<res)
{
if(distL<distR)
{
ask(l,mid-1,x);
if(distR<res)ask(mid+1,r,x);
}
else
{
ask(mid+1,r,x);
if(distL<res)ask(l,mid-1,x);
}
}
else
{
if(distL<res)ask(l,mid-1,x);
if(distR<res)ask(mid+1,r,x);
}
}
int build(int l,int r)
{
if(l>r)return 0;
if(l==r)
{
pushup(l);
return l;
}
int mid=l+r>>1;
double avg[2]={0},va[2]={0};
for(int i=l;i<=r;i++)avg[0]+=point[i].x,avg[1]+=point[i].y;
avg[0]/=(double)(r-l+1),avg[1]/=(double)(r-l+1);
for(int i=l;i<=r;i++)
va[0]+=(point[i].x-avg[0])*(point[i].x-avg[0]),
va[1]+=(point[i].y-avg[1])*(point[i].y-avg[1]);
if(va[0]>=va[1])
nth_element(point+l,point+mid,point+r+1,cmp1);
else
nth_element(point+l,point+mid,point+r+1,cmp2);
tr[mid].s[0]=build(l,mid-1),tr[mid].s[1]=build(mid+1,r);
pushup(mid);
return mid;
}
int main()
{
tr[0].L=tr[0].D=res,tr[0].R=tr[0].U=-res;
read(n);
for(int i=1;i<=n;i++)
read(point[i].x),read(point[i].y);
sort(point+1,point+1+n,[](const Point &a,const Point &b){
if(a.x!=b.x)return a.x<b.x;
return a.y<b.y;
});
build(1,n);
for(int i=1;i<=n;i++)ask(1,n,i);
printf("%lld",LL(res));
return 0;
}