推薦：

k-d tree演算法

對於D維的點若干，多次查詢距離某個點第K大的點是什麼。

處理這一類問題的一個數據結構，叫K-D Tree

基本思想是對點進行區域分塊處理。

圖示：

K-D Tree是一個二叉樹。

每個點維護的資訊是，

split ：分裂座標軸

ls、rs：左右兒子

node：該節點儲存的真實點

建樹：

遞迴建樹。類似線段樹（但是每個點有實際的點）

選擇當前區域的點的各個維度的方差最大的維度（傳說如果方差大，資料分散，複雜度或者精度有所保證？？），把這個維度當做split

這個節點的真實點就是c[mid]

然後，把這個維度[s]小於c[mid][s]的放在左邊，大於的放在右邊。

（實現時，用一個nth_element，再過載小於號，可以O(n)實現把中間的放在mid位置上，並且這個維度[s]小於c[mid][s]的放在左邊，大於的放在右邊。）

然後遞迴建樹即可。

x，y是split

這樣，整個K-D Tree就把一些點分成了若干個塊。

我們一塊一塊處理會比較容易剪枝。

查詢：最近的點（即K=1）

本質是爆搜+剪枝。。。

設查詢距離點st的最近的點to

設距離為now

法一：

不斷通過當前split維度和st這個維度的大小比較，

我們先走st所屬的塊，

回溯回來之後，

由於可能在另一半有更近的點。

如果分界線到st的該維度距離小於now

那麼再走另外一個塊搜尋。

法二：

上面那個剪枝比較粗糙。

我們發現，一個塊的所有點，其實可以用一個矩形框住。

那麼，如果st到這個矩形可能的最近點距離小於now的話，再搜下去。

具體來說，我們每個節點維護這個節點代表的塊內，最大最小的x，y座標。（其實就是矩形四個頂點）

這個最短距離：

x差為：dx=max(st.x-x.mxx,0)+max(x.min-st.x,0)

畫個圖理解下，如果st的x在mix,mxx之間的話，那麼x差就認為是0

如果在mix左邊，那麼就是x-mix,

如果在mxx右邊，那麼就是mxx-x

y同理。

dis=sqrt(dx*dx+dy*dy)

發現，當st在所屬的塊內時，dis一定是0

然後就可以剪枝了。

對於兩個兒子，選擇估價距離較小的那個先搜，回溯時，如果另一個距離還比now小的話，再搜另一個。

理論上，應該比法一多減一些枝。

總之，複雜度不明。

傳說最差O(k*n^(1-1/k))每次（k是維度）

例題：

平面最近點對（加強版）

這個題用分治是最好的。

我們用KD樹來試試。

列舉所有的點，找到與它最近的點距離，然後所有距離取min即可

直接做就好了。

注意：

1.如果寫法二，那麼對於0號節點的哨兵必須mi=inf,mx=-inf。否則剪枝就掛了。

2.建樹的時候，build返回節點編號，不能返回計數器tot。。。。

這個題法一更快？？

可能資料水，然後法二常數大吧。。。

法一：

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
namespace Miracle{
const int N=200000+5;
const double inf=2333333333.00;
int n;
struct po{
    double x,y;
    int id;
    po(){}
    po(double xx,double yy){
        x=xx;y=yy;
    }
}a[N],c[N],st,to;
bool cmp1(po a,po b){
    return a.x<b.x;
}
bool cmp2(po a,po b){
    return a.y<b.y;
}
double ans;
double now;
double dis(po a,po b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
struct tr{
    int sp;
    po O;
    int ls,rs;
}t[2*N];
int tot;
int rt;
int build(int l,int r){
    if(l>r){
        return 0;
    }
    if(l==r){
        ++tot;
        t[tot].O=c[l];
        t[tot].ls=t[tot].rs=0;
        t[tot].sp=1;
        return tot;
    }
    int mid=(l+r)>>1;
    double ax=0,ay=0;
    for(reg i=l;i<=r;++i) ax+=c[i].x,ay+=c[i].y;
    ax/=(r-l+1);ay/=(r-l+1);
    double fx=0,fy=0;
    for(reg i=l;i<=r;++i) fx+=(c[i].x-ax)*(c[i].x-ax),fy+=(c[i].y-ay)*(c[i].y-ay);
    fx/=(r-l+1);fy/=(r-l+1);
    int ret=++tot;
    if(fx>fy){//choose x;
        nth_element(c+l,c+mid,c+r+1,cmp1);
        t[ret].sp=1;
    }else{
        nth_element(c+l,c+mid,c+r+1,cmp2);
        t[ret].sp=2;
    }
    t[ret].O=c[mid];
    t[ret].ls=build(l,mid-1);
    t[ret].rs=build(mid+1,r);
    return ret;
}
void dfs(int x){
    if(!x) return;
    if(st.id!=t[x].O.id&&dis(st,t[x].O)<now){
        now=dis(st,t[x].O);
    }
    if(t[x].sp==1){
        double d=fabs(t[x].O.x-st.x);
        if(st.x<=t[x].O.x){
            dfs(t[x].ls);
            if(d<now) dfs(t[x].rs);
        }
        else{
            dfs(t[x].rs);
            if(d<now) dfs(t[x].ls);
        }
    }
    else{
        double d=fabs(t[x].O.y-st.y);
        if(st.y<=t[x].O.y){
            dfs(t[x].ls);
            if(d<now) dfs(t[x].rs);
        }
        else{
            dfs(t[x].rs);
            if(d<now) dfs(t[x].ls);
        }
    }
}
int main(){
    scanf("%d",&n);
    for(reg i=1;i<=n;++i){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        a[i].id=i;
        c[i]=a[i];
    }
    rt=build(1,n);
    ans=inf;
    for(reg i=1;i<=n;++i){
        st=a[i];
        now=inf;
        to=po(inf,inf);
        dfs(1);
        ans=min(ans,now);
    }
    printf("%.4lf",ans);
    return 0;
}

}
int main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2018/11/26 8:43:17
*/

法二：

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){
    char ch;x=0;bool fl=false;
    while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
    for(x=numb;isdigit(ch=getchar());x=x*10+numb);
    (fl==true)&&(x=-x);
}
namespace Miracle{
const int N=200000+5;
const double inf=2333333333.00;
int n;
struct po{
    double x,y;
    int id;
    po(){}
    po(double xx,double yy){
        x=xx;y=yy;
    }
}a[N],c[N],st,to;
bool cmp1(po a,po b){
    return a.x<b.x;
}
bool cmp2(po a,po b){
    return a.y<b.y;
}
double ans;
double now;
double dis(po a,po b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
struct tr{
    double mxx,mix,mxy,miy;
    int sp;
    po O;
    int ls,rs;
}t[2*N];
int tot;
int rt;
int build(int l,int r){
    if(l>r){
        return 0;
    }
    if(l==r){
        ++tot;
        t[tot].mxx=t[tot].mix=c[l].x;
        t[tot].mxy=t[tot].miy=c[l].y;
        t[tot].O=c[l];
        t[tot].ls=t[tot].rs=0;
        t[tot].sp=1;
        return tot;
    }
    int mid=(l+r)>>1;
    double ax=0,ay=0;
    for(reg i=l;i<=r;++i) ax+=c[i].x,ay+=c[i].y;
    ax/=(r-l+1);ay/=(r-l+1);
    double fx=0,fy=0;
    for(reg i=l;i<=r;++i) fx+=(c[i].x-ax)*(c[i].x-ax),fy+=(c[i].y-ay)*(c[i].y-ay);
    fx/=(r-l+1);fy/=(r-l+1);
    int ret=++tot;
    if(fx>fy){//choose x;
        nth_element(c+l,c+mid,c+r+1,cmp1);
        t[ret].sp=1;
    }else{
        nth_element(c+l,c+mid,c+r+1,cmp2);
        t[ret].sp=2;
    }
    t[ret].O=c[mid];
    t[ret].ls=build(l,mid-1);
    t[ret].rs=build(mid+1,r);
    t[ret].mxx=max(t[t[ret].rs].mxx,t[t[ret].ls].mxx);
    t[ret].mix=min(t[t[ret].rs].mix,t[t[ret].ls].mix);
    t[ret].mxy=max(t[t[ret].rs].mxy,t[t[ret].ls].mxy);
    t[ret].miy=min(t[t[ret].rs].miy,t[t[ret].ls].miy);
    //cout<<" ret "<<ret<<" "<<l<<" "<<r<<endl;
    return ret;
}
void dfs(int x){
    if(st.id!=t[x].O.id&&dis(st,t[x].O)<now){
        now=dis(st,t[x].O);
        to=t[x].O;
    }
    if(t[x].ls&&t[x].rs){
        double lx=max(st.x-t[t[x].ls].mxx,0.0)+max(t[t[x].ls].mix-st.x,0.0);
        double ly=max(st.y-t[t[x].ls].mxy,0.0)+max(t[t[x].ls].miy-st.y,0.0);
        double len1=sqrt(lx*lx+ly*ly);
        double rx=max(st.x-t[t[x].rs].mxx,0.0)+max(t[t[x].rs].mix-st.x,0.0);
        double ry=max(st.y-t[t[x].rs].mxy,0.0)+max(t[t[x].rs].miy-st.y,0.0);
        double len2=sqrt(rx*rx+ry*ry);
        if(len1<=len2&&len1<now){
            dfs(t[x].ls);
            if(len2<now) 
                dfs(t[x].rs);
        }
        else if(len2<=len1&&len2<now){
            dfs(t[x].rs);
            if(len1<now) 
                dfs(t[x].ls);
        }
    }
    else if(t[x].ls){
        double lx=max(st.x-t[t[x].ls].mxx,0.0)+max(t[t[x].ls].mix-st.x,0.0);
        double ly=max(st.y-t[t[x].ls].mxy,0.0)+max(t[t[x].ls].miy-st.y,0.0);
        double len1=sqrt(lx*lx+ly*ly);
        if(len1<now) 
        dfs(t[x].ls);
    }
    else if(t[x].rs){
        double rx=max(st.x-t[t[x].rs].mxx,0.0)+max(t[t[x].rs].mix-st.x,0.0);
        double ry=max(st.y-t[t[x].rs].mxy,0.0)+max(t[t[x].rs].miy-st.y,0.0);
        double len2=sqrt(rx*rx+ry*ry);
        if(len2<now) 
        dfs(t[x].rs);
    }
    else return;
}
int main(){
    scanf("%d",&n);
    t[0].mix=inf;t[0].mxx=-inf;
    t[0].miy=inf;t[0].mxy=-inf;
    for(reg i=1;i<=n;++i){
        scanf("%lf%lf",&a[i].x,&a[i].y);
        a[i].id=i;
        c[i]=a[i];
    }
    rt=build(1,n);
    ans=inf;
    for(reg i=1;i<=n;++i){
    //    cout<<" ii "<<i<<" : "<<a[i].x<<" "<<a[i].y<<" ------------------ "<<endl;
        st=a[i];
        now=inf;
        to=po(inf,inf);
        dfs(1);
    //    cout<<" after "<<now<<endl;
        ans=min(ans,now);
    }
    printf("%.4lf",ans);
    return 0;
}

}
int main(){
    Miracle::main();
    return 0;
}

/*
   Author: *Miracle*
   Date: 2018/11/26 8:43:17
*/

 對了，KD-Tree其實也可以不記錄左右兒子，以及代表實際點

因為，每次我們選擇的是mid位置的點，之後這個點的位置也不會再動了。

而左右兒子區間也是定值。

所以，query時記錄(l,r)即可，訪問實際點的話，直接取c[mid]就好。