P7883

发表于 2022-04-03 更新于 2022-05-04 分类于模板

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

把原来分治的代码修了一下。

再讲一遍分治的做法。

首先把点按照横坐标为第一关键字，纵坐标为第二关键字排序。

找到中间位置的点，以这个点的横坐标为界限，分出了两个区域的点，分别进入两个区域解决子问题。

于是我们得到了左区域的答案 $d_1$ 和右区域的答案 $d_2$，令 $d=\min\{d_1,d_2\}$

然后我们直接在 $[l,r]$ 范围内找距离中位数所在直线的距离 $\le d$ 的点，共计有 $cnt$ 个。

然后我们将这 $cnt$ 个点按纵坐标排序，再在这 $cnt$ 个点里打暴力，枚举一个点 $i$，再枚举另一个点 $j$，如果说 $i$，$j$ 之间的纵坐标距离 $\ge d$ 就及时退出 $j$ 的枚举。

然后我们把取到的最小值返回即可。

正确性应该是显然的，复杂度大概会有些迷惑。

实际上对于每个 $i$，我们比较的 $j$ 只有 $O(1)$ 个。

这样想，对于一个 $i$，我们枚举到的 $j$ 被严格限制在了一个 $d\times 2d$ 的小矩形里，我们将其拆分为左右两个 $d\times d$ 的小正方形，不难得出每个小正方形内部的点的距离相互之间必须 $\ge d$。

感性的画个图就会发现这个点画不了太多（我还是不会严谨证明。。。哪位巨神教教我啊 QAQ）。

时间复杂度 $O(n\log n)$。

代码（分治）：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
namespace Ehnaev{
  inline ll read() {
    ll ret=0,f=1;char ch=getchar();
    while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
    while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
    return ret*f;
  }
  inline void write(ll x) {
    static char buf[22];static ll len=-1;
    if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
    else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
    while(len>=0) putchar(buf[len--]);
  }
}using Ehnaev::read;using Ehnaev::write;

const ll N=1e6,inf=1e15;

ll n;
ll tmp[N+5];

struct node{ll x,y;}a[N+5];

inline bool _Cmp(ll x,ll y) {return a[x].y<a[y].y;}

inline ll Dis(ll x,ll y) {
  return (a[x].x-a[y].x)*(a[x].x-a[y].x)+(a[x].y-a[y].y)*(a[x].y-a[y].y);
}

inline ll Msort(ll l,ll r) {
  ll d=inf;
  if(l==r) return d;if(l+1==r) return Dis(l,r);
  ll mid=(l+r)>>1;ll d1=Msort(l,mid),d2=Msort(mid+1,r);
  d=min(d1,d2);ll cnt=0;
  for(ll i=l;i<=r;i++) {
    if((a[mid].x-a[i].x)*(a[mid].x-a[i].x)<d) {tmp[++cnt]=i;}
  }
  sort(tmp+1,tmp+cnt+1,_Cmp);
  for(ll i=1;i<=cnt;i++) {
    for(ll j=i+1;j<=cnt
    &&(a[tmp[j]].y-a[tmp[i]].y)*(a[tmp[j]].y-a[tmp[i]].y)<d;j++) {
      ll d3=Dis(tmp[i],tmp[j]);if(d>d3) {d=d3;}
    }
  }
  return d;
}

inline bool Cmp(node x,node y) {return x.x==y.x?x.y<y.y:x.x<y.x;}

int main() {
  
  n=read();
  
  for(ll i=1;i<=n;i++) {a[i].x=read();a[i].y=read();}
  
  sort(a+1,a+n+1,Cmp);
  
  write(Msort(1,n));
  
  return 0;
}

以此作为 K-D 树基础的记录。。。

我们将平面使用 K-D 树的方法划分。

OI wiki 上给出的 K-D 树的划分是这样的：

现在有一个超长方体内有 $n$ 个点。
找出方差最大的维度，并选择该维度的中位数，将所有的点根据该维度的位置分布分为两种。
然后在两边的点中继续按这种方法划分。

然后我们维护 K-D 树上每个节点（被划分出来的点集）的 K 维界限（形成的超矩形），就产生了各种各样的乱搞做法。

比如说，我们这题就可以用个 2-D 树，然后枚举 $n$ 个点，在 K-D 树上查询最近点。

查询方法很 sb，就是直接遍历树上节点（各种被划分的点集），然后把显然不可能成为答案的点集（上下左右界产生的矩形与该点的最短距离都比现有答案劣，且该点在矩形外）直接剪枝掉。

然后你就能跑过 P1429（笔者乱搞能力尚欠，过不了 P7883）。

这个复杂度听上去就很错误。

划分 2-D 树时，求方差是 $O(n)$ 的，找中位数并把点分到两边可以使用 nth_element 函数（属于 algorithm），时间复杂度 $O(n)$，递归的话用一下主定理得到总时间复杂度 $O(n\log n)$。

但是查询的时候，枚举点就是 $O(n)$ 的，在 2-D 树上查询的时候，最坏情况是把所有节点都遍历了，时间复杂度仍是 $O(n)$，算上去的话总复杂度应该是 $O(n^2)$！

但是它就是跑得飞快，我也不知道为什么。

可能又有什么高深的期望复杂度理论吧，我不会我不会（

代码（K-D 树，P1429）：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define ll long long
#define ld long double
using namespace std;
namespace Ehnaev{
  inline ll read() {
    ll ret=0,f=1;char ch=getchar();
    while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
    while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
    return ret*f;
  }
  inline void write(ll x) {
    static char buf[22];static ll len=-1;
    if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
    else {putchar(45);do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
    while(len>=0) putchar(buf[len--]);
  }
}using Ehnaev::read;using Ehnaev::write;

const ll N=2e6;
const ld inf=2e18;

ll n;
ld ans=inf;
ll d[N+5],lc[N+5],rc[N+5];
ld L[N+5],R[N+5],D[N+5],U[N+5];

struct node{ld x,y;}a[N+5];

inline ld Dist(ll u,ll v) {
  return (a[u].x-a[v].x)*(a[u].x-a[v].x)+(a[u].y-a[v].y)*(a[u].y-a[v].y);
}

inline bool Cmp1(node a,node b) {return a.x<b.x;}
inline bool Cmp2(node a,node b) {return a.y<b.y;}

inline void Maintain(ll x) {
  L[x]=R[x]=a[x].x;D[x]=U[x]=a[x].y;
  if(lc[x]) {
    L[x]=min(L[x],L[lc[x]]);R[x]=max(R[x],R[lc[x]]);
    D[x]=min(D[x],D[lc[x]]);U[x]=max(U[x],U[lc[x]]);
  }
  if(rc[x]) {
    L[x]=min(L[x],L[rc[x]]);R[x]=max(R[x],R[rc[x]]);
    D[x]=min(D[x],D[rc[x]]);U[x]=max(U[x],U[rc[x]]);
  }
}

inline ll Build(ll l,ll r) {
  if(l>=r) return 0;ll mid=(l+r)>>1;
  ld avx=0,avy=0,vax=0,vay=0; // Average & Variance.
  for(ll i=l;i<=r;i++) {avx+=a[i].x;avy+=a[i].y;}
  avx/=(ld)(r-l+1);avy/=(ld)(r-l+1);
  for(ll i=l;i<=r;i++) {
    vax+=(a[i].x-avx)*(a[i].x-avx);vay+=(a[i].y-avy)*(a[i].y-avy);
  }
  if(vax>=vay) {d[mid]=1;nth_element(a+l,a+mid,a+r+1,Cmp1);}
  else {d[mid]=2;nth_element(a+l,a+mid,a+r+1,Cmp2);}
  lc[mid]=Build(l,mid-1);rc[mid]=Build(mid+1,r);
  Maintain(mid);return mid;
}

inline ld F(ll u,ll v) {
  ld ret=0;
  if(L[v]>a[u].x) ret+=(L[v]-a[u].x)*(L[v]-a[u].x);
  if(R[v]<a[u].x) ret+=(a[u].x-R[v])*(a[u].x-R[v]);
  if(D[v]>a[u].y) ret+=(D[v]-a[u].y)*(D[v]-a[u].y);
  if(U[v]<a[u].y) ret+=(a[u].y-U[v])*(a[u].y-U[v]);
  return ret;
}

inline void Query(ll l,ll r,ll x) {
  if(l>r) return;ll mid=(l+r)>>1;
  if(mid!=x) ans=min(ans,Dist(x,mid));
  if(l==r) return;ld distl=F(x,lc[mid]),distr=F(x,rc[mid]);
  if(distl<ans&&distr<ans) {
    if(distl<distr) {
      Query(l,mid-1,x);if(distr<ans) Query(mid+1,r,x);
    }
    else {
      Query(mid+1,r,x);if(distl<ans) Query(l,mid-1,x);
    }
  }
  else {
    if(distl<ans) Query(l,mid-1,x);if(distr<ans) Query(mid+1,r,x);
  }
}

int main() {

  n=read();

  for(ll i=1;i<=n;i++) scanf("%Lf %Lf",&a[i].x,&a[i].y);

  Build(1,n);

  for(ll i=1;i<=n;i++) Query(1,n,i);

  printf("%.4Lf",sqrt(ans));

  return 0;
}