P7883 平面最近点对（加强加强版）

题目链接

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 时限的原因是：

1. \(O(n \log ^2 n)\) 参考代码使用 cin 不会 TLE。最快的 std 能 \(<\) 100ms。
2. @wlzhouzhuan 的程序能恰好在 350ms 内跑 \(1000n\) 次检查。
3. 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)\)

• 均摊时间复杂度：\(O(n\times \sqrt{n})\)

代码

// 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()
{
	read(n);
	for(int i=1;i<=n;i++)
    	read(point[i].x),read(point[i].y);
    build(1,n);
    for(int i=1;i<=n;i++)ask(1,n,i);
    printf("%lld",LL(res));
    return 0;
}