KD樹小結 【bzoj3489】 A simple rmq problem

好久以前我就想過怎麼快速在二維平面上查找一個區域的信息,思考許久無果,只能想到幾種優秀一點的暴力。php

KD樹就是幹上面那件事的。node

別的很少說,趕忙把本身的理解寫下來,省得涼了。ios

 

KD樹的組成

以維護k維空間(x,y,……)內的KD樹爲例,主要由一下三部分組成:數組

  1. p[k],表明樹上這個結點所儲存的點(在題目中給出的/你本身加上的點集中的一個點)。
  2. ch[2],表示它的子結點(沒錯,KD樹是一棵二叉樹)
  3. mi[k]與mx[k],mi/mx[i]表明KD樹這個結點統轄的全部點的第i-1範圍。好比說mi[1]=2,mx[1]=4,就表明這棵樹統轄的點的y座標都在[2,4]內。

 

不看mi和mx,長得就和splay/trie樹同樣,一個p維護當前節點,一個ch[2]記錄左右兒子。數據結構

不看p[k],長得就和線段樹同樣,有左右兒子和區間信息。ide

沒錯,KD樹功能如線段樹,結點維護區域信息;形態如splay/trie樹,每一個結點有實際的值和意義。函數

 

KD樹的構建

通常題目都是二維平面。下面就以二維平面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樹的精髓以後你就會發現……嘿嘿嘿……

(就是爲了把這個暴力數據結構剪枝剪更多跑更快)

 

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')最近的點(曼哈頓距離,|x-x'|+|y-y'|):
    • 首先咱們看暴力的剪枝:按某一維排序,若是該維的差過大就無論了。
    • 而令咱們期待的KD樹呢?呃很差意思,它也是這麼作的……
    • 咱們維護過兩個叫作mi[]和mx[]的東西吧……這個時候就是它派上用場了。
    • 具體還請看代碼吧:
      //查詢的點(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]);
        }
      }
    • dis():若是當前點在這個區間內就是0,不然就是最極的點到它的距離。
    • 聰明絕頂的你已經發現了……這TM就是個暴力。
    • 其實這個數據結構就是一個暴力……
    • 當暴力有了時間複雜度證實……還叫暴力麼?讀書人的事,能叫偷麼?
    • 這麼暴力有幾個好處:不用枚舉全部點;剪枝有效及時。
    • 複雜度有保障,大概在O(√n)級別下,主要看數據。
  • 若是是區間查詢,以區間查詢點權和爲例(以前就有維護好):
    • 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;
}
簡單題
相關文章
相關標籤/搜索