好久以前我就想過怎麼快速在二維平面上查找一個區域的信息,思考許久無果,只能想到幾種優秀一點的暴力。php
KD樹就是幹上面那件事的。node
別的很少說,趕忙把本身的理解寫下來,省得涼了。ios
以維護k維空間(x,y,……)內的KD樹爲例,主要由一下三部分組成:數組
不看mi和mx,長得就和splay/trie樹同樣,一個p維護當前節點,一個ch[2]記錄左右兒子。數據結構
不看p[k],長得就和線段樹同樣,有左右兒子和區間信息。ide
沒錯,KD樹功能如線段樹,結點維護區域信息;形態如splay/trie樹,每一個結點有實際的值和意義。函數
通常題目都是二維平面。下面就以二維平面KD樹的構建爲例。post
讀入把點存進結構體數組a中,座標分別爲a[x].p[i]。大數據
inline void build(int &x,int l,int r,int type){ x=(l+r)>>1;now=type; nth_element(a+l,a+x,a+r+1,cmp); nd=a[x];newnode(x); if(l<x)build(ch[x][0],l,x-1,type^1);else ch[x][0]=0; if(x<r)build(ch[x][1],x+1,r,type^1);else ch[x][1]=0; pushup(x); } build(kd.root,1,n,0);
很是優美……對type、now做用不明的同窗請繼續閱讀……你要如今就明白就奇怪了ui
系統函數nth_element(a+l,a+x,a+r+1),頭文件algorithm,需定義<或cmp函數。
做用:把排序後第x大的放到第x位,比它小的放進左邊,比它大的放進右邊(兩邊無序)。
注意區間開閉:左閉右開,中間也是閉合的。
複雜度:平均,指望是O(n)?能夠接受。
下面給出cmp、newnode、pushup代碼。
struct Node{int p[2],mi[2],mx[2];}a[N]; inline bool cmp(const Node &a,const Node &b){return a.p[now]<b.p[now];} inline void Min(int &x,int y){x=x<y?x:y;} inline void Max(int &x,int y){x=x>y?x:y;} inline void pushup(int x){ int ls=ch[x][0],rs=ch[x][1]; if(ls){ Min(T[x].mi[0],T[ls].mi[0]);Max(T[x].mx[0],T[ls].mx[0]); Min(T[x].mi[1],T[ls].mi[1]);Max(T[x].mx[1],T[ls].mx[1]); } if(rs){ Min(T[x].mi[0],T[rs].mi[0]);Max(T[x].mx[0],T[rs].mx[0]); Min(T[x].mi[1],T[rs].mi[1]);Max(T[x].mx[1],T[rs].mx[1]); } } inline void newnode(int x){ T[x].p[0]=T[x].mi[0]=T[x].mx[0]=nd.p[0]; T[x].p[1]=T[x].mi[1]=T[x].mx[1]=nd.p[1]; }
不要問我爲何辣麼長,爲了減常衝榜,把循環展開了……
聰明的讀者已經發現KD樹的構建巧妙之處。它不是純粹按照x維,或者某一維排序,而是先按x維,再按y維,再按z維,再……最後又回到x維……
這樣分割的區域更加整齊劃一更加均勻緊縮,不像上面的按照某一維劃分,到最後變成一條條長條,KD樹劃分到底的圖案仍是很好看的。
這樣分割有什麼好處呢?等你真正領悟了KD樹的精髓以後你就會發現……嘿嘿嘿……
(就是爲了把這個暴力數據結構剪枝剪更多跑更快)
1.往KD樹上插點
插點能夠分爲插新點和插老點。若是有老點,特判一句,把信息覆蓋便可。
inline void insert(int &x,int type){ if(!x){x=++cnt,newnode(cnt);return;} if(nd.p[0]==T[x].p[0] && nd.p[1]==T[x].p[1]){ ……(自行維護);return; } if(nd.p[type]<T[x].p[type])insert(ch[x][0],type^1); else insert(ch[x][1],type^1); pushup(x); }
依然很是的美妙……等等有什麼不對?
咱們能估計出一棵剛建好的KD樹深度是O(log)的。
但你這麼隨便亂插……有道題叫HNOI2017 spaly 插入不旋轉的單旋spaly見過?T成苟。
這都不是問題!知不知道有一種數據結構叫作替罪羊樹哇?
知道替罪羊樹怎麼保證複雜度的嗎?
重構!大力重構!自信重構!不爽就重構!
爲了省事大概每插入10000次就重構一次好了……
if(kd.cnt==sz){ for(int i=1;i<=sz;++i)a[i]=kd.T[i]; kd.rebuild(kd.root,1,sz,0);sz+=10000; }
2.在KD樹上查詢
//查詢的點(x',y')儲存在nd中。 //這裏的l,r就是mi,mx的意思。 inline int dis(Node p,int x,int ans=0){ for(int i=0;i<2;++i) ans+=max(0,t[x].l[i]-p.p[i])+max(0,p.p[i]-t[x].r[i]); return ans; } inline void query(int x){ Ans=min(Ans,abs(t[x].p[0]-nd.p[0])+abs(t[x].p[1]-nd.p[1])); int dl=ch[x][0]?dis(nd,ch[x][0]):Inf; int dr=ch[x][1]?dis(nd,ch[x][1]):Inf; if(dl<dr){ if(dl<Ans)query(ch[x][0]); if(dr<Ans)query(ch[x][1]); } else{ if(dr<Ans)query(ch[x][1]); if(dl<Ans)query(ch[x][0]); } }
inline bool in(int l,int r,int xl,int xr){return l<=xl && xr<=r;} inline bool out(int l,int r,int xl,int xr){return xr<l || r<xl;} inline int query(int x,int x1,int y1,int x2,int y2){ int ans=0;if(!x)return ans; if(in(x1,x2,T[x].mi[0],T[x].mx[0])) if(in(y1,y2,T[x].mi[1],T[x].mx[1])) return T[x].sum; if(out(x1,x2,T[x].mi[0],T[x].mx[0]))return 0; if(out(y1,y2,T[x].mi[1],T[x].mx[1]))return 0; if(in(x1,x2,T[x].p[0],T[x].p[0])) if(in(y1,y2,T[x].p[1],T[x].p[1])) ans+=T[x].val; return ans+query(ch[x][0],x1,y1,x2,y2)+query(ch[x][1],x1,y1,x2,y2); }
KD樹的基本姿式大概就是這個樣子……好寫很差寫錯,基本上都是個板子。
附上學長的一(三)句話,從各方面進行了深度總結:
「能不寫最好仍是不要寫吧,輕鬆被卡→_→,也許能夠出奇制勝?若是要寫,從新構樹是個不錯的選擇。發現大數據跑不過,多半是剪枝掛了。」
仍是給個連接……MashiroSky大爺。
upd:以當前座標差最大的來作type應該比輪換type更優秀……
例題有"SJY擺棋子"、"簡單題"等,在此就不作贅述了。
比較有意思的應用就是【bzoj3489】 A simple rmq problem,正如上面所言,KD樹解決傳統數據結構題出奇制勝。
附上"BZOJ4066簡單題"代碼一份,操做是單點修改+矩形求和在線。
#include <iostream> #include <cstdio> #include <cstdlib> #include <algorithm> #include <vector> #include <cstring> #include <queue> #include <complex> #include <stack> #define LL long long int #define dob double #define FILE "bzoj_4066" //#define FILE "簡單題" using namespace std; const int N = 200010; int n,lst,now,sz=10000; inline int gi(){ int x=0,res=1;char ch=getchar(); while(ch>'9'||ch<'0'){if(ch=='-')res*=-1;ch=getchar();} while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar(); return x*res; } inline void Min(int &x,int y){x=x<y?x:y;} inline void Max(int &x,int y){x=x>y?x:y;} struct Node{int p[2],mi[2],mx[2],val,sum;}a[N]; inline bool cmp(const Node &a,const Node &b){return a.p[now]<b.p[now];} struct KD_Tree{ int ch[N][2],root,cnt; Node T[N],nd; inline void pushup(int x){ int ls=ch[x][0],rs=ch[x][1]; if(ls){ Min(T[x].mi[0],T[ls].mi[0]);Max(T[x].mx[0],T[ls].mx[0]); Min(T[x].mi[1],T[ls].mi[1]);Max(T[x].mx[1],T[ls].mx[1]); } if(rs){ Min(T[x].mi[0],T[rs].mi[0]);Max(T[x].mx[0],T[rs].mx[0]); Min(T[x].mi[1],T[rs].mi[1]);Max(T[x].mx[1],T[rs].mx[1]); } T[x].sum=T[x].val; if(ls)T[x].sum+=T[ls].sum; if(rs)T[x].sum+=T[rs].sum; } inline void newnode(int x){ T[x].p[0]=T[x].mi[0]=T[x].mx[0]=nd.p[0]; T[x].p[1]=T[x].mi[1]=T[x].mx[1]=nd.p[1]; T[x].sum=T[x].val=nd.val; } inline void insert(int &x,int type){ if(!x){x=++cnt,newnode(cnt);return;} if(nd.p[0]==T[x].p[0] && nd.p[1]==T[x].p[1]){ T[x].val+=nd.val;T[x].sum+=nd.val; return; } if(nd.p[type]<T[x].p[type])insert(ch[x][0],type^1); else insert(ch[x][1],type^1); pushup(x); } inline void rebuild(int &x,int l,int r,int type){ x=(l+r)>>1;now=type; nth_element(a+l,a+x,a+r+1,cmp); nd=a[x];newnode(x); if(l<x)rebuild(ch[x][0],l,x-1,type^1);else ch[x][0]=0; if(x<r)rebuild(ch[x][1],x+1,r,type^1);else ch[x][1]=0; pushup(x); } inline bool in(int l,int r,int xl,int xr){return l<=xl && xr<=r;} inline bool out(int l,int r,int xl,int xr){return xr<l || r<xl;} inline int query(int x,int x1,int y1,int x2,int y2){ int ans=0;if(!x)return ans; if(in(x1,x2,T[x].mi[0],T[x].mx[0])) if(in(y1,y2,T[x].mi[1],T[x].mx[1])) return T[x].sum; if(out(x1,x2,T[x].mi[0],T[x].mx[0]))return 0; if(out(y1,y2,T[x].mi[1],T[x].mx[1]))return 0; if(in(x1,x2,T[x].p[0],T[x].p[0])) if(in(y1,y2,T[x].p[1],T[x].p[1])) ans+=T[x].val; return ans+query(ch[x][0],x1,y1,x2,y2)+query(ch[x][1],x1,y1,x2,y2); } }kd; int main() { freopen(FILE".in","r",stdin); freopen(FILE".out","w",stdout); n=gi(); while(1){ int type=gi();if(type==3)break; int x1=gi()^lst,y1=gi()^lst; if(type==1){ int A=gi()^lst; kd.nd.p[0]=x1;kd.nd.p[1]=y1; kd.nd.sum=kd.nd.val=A; kd.insert(kd.root,0); if(kd.cnt==sz){ for(int i=1;i<=sz;++i)a[i]=kd.T[i]; kd.rebuild(kd.root,1,sz,0);sz+=10000; } } if(type==2){ int x2=gi()^lst,y2=gi()^lst; lst=kd.query(kd.root,x1,y1,x2,y2); printf("%d\n",lst); } } fclose(stdin);fclose(stdout); return 0; }