【学习小记】KD-Tree

Preface

听说KD树实在是个大神器
可以解决多维空间多维偏序点权和，可以求某个点的空间最近最远点

就二维平面上的来说，复杂度在\(O(n\log n)\)到\(O(n\sqrt n)\)不等

嫌KD树不平衡了还可以来一个替罪羊树式的暴力重构
再也不用担心写不出树套树了！（狗头）

Text

这是个什么东西呢？

在一维数轴上，可以简单的比较两个值的大小，我们就很容易将它们分治，建出一个二叉搜索树
但是拓展到高维往往就很困难

KD树解决了这一麻烦，KD树也是一个二叉分治结构，不同的是KD树的分治维是会变的。

比如说对于1号节点，它是按照X坐标分治的，即它的左子树的X坐标比他小，Y坐标比他大
对于另一个2号节点，他可能又按照Y坐标分治了。

可以根据以下构树的过程来加深理解。

构造

以二维平面为例，我们现在有一堆点。

按什么坐标分治其实是可以随机的，而另一种比较好的方法是每一维轮流，父亲按X分治，当前就按Y分治这样子下去

实现上来说，我们当前的节点可以表示成节点序列上的一段区间[l,r]。
现在我们想要知道节点序列上区间[l,r]中第D维的中位数，然后把比它小的放到它左边，大的放到右边。

这个东西排序似乎要\(n\log\)？
然而c++STL库提供了一个std::nth_element()的函数,找出序列第k大，并且把小的放左边，大的放右边，内部的相对顺序不保证，并且它的复杂度是线性的。（原理好像是基于分治的）
它一般可以传四个参，依次是序列头，序列第k个位置，序列维，比较函数
调用完以后序列第k个位置就已经放好了。

这样就搞定了，递归建树即可，注意如果是左边或右边没有元素了就不递归那边。

显然，平衡树能记录的值它都能记录。子树和，子树最大值等等

为了方便，我们一般还记录子树中每一维坐标的范围。

可以看出，每个子树就是一个矩形内的所有点。

动态插入

类似splay一样向下搜索到叶子，直接加入即可。

然而多次以后可能会出现不平衡的情况。（即某一侧子树过大）
那么我们可以类似替罪羊树的做法，插入以后找最高的一个不平衡点，对它的子树进行重构。

当然，由于KD树的本质是一个暴力，你完全可以每插入7000次或者10000次之类的把整棵树重构一遍（好写简洁，慢不了多少）

动态删除

类似替罪羊树，删除一个节点时，找它左子树中这一维最大的那个节点替代它。
不平衡就重构。

基本操作/询问

矩形加/矩形赋值/矩形查询，类似线段树那样做就好了。
对于K维操作，n个点，总的复杂度大概是\(O(kn^{1-{1\over k}})\)（当然1维不是）

有一个它能做而树套树完成不了的问题——K维最近/最远点
这个距离可以是曼哈顿距离，欧几里得距离，甚至切比雪夫距离。

它的本质就是一个搜索。
由于我们知道了每个子树的每维坐标的范围，我们就可以算出询问节点到这个子树的矩形的距离至少（至多）是多少

以欧几里得最近点为例，
令矩形为\((x1,y1,x2,y2)\)，询问点为\((p,q)\)，那么这个点到这个矩形内的点距离至少是\(\sqrt{max(0,max(x1-p,p-x2))^2+max(0,max(y1-q,q-y2))^2}\)

简单来说，这就是一个估价函数。
如果估价函数值比当前已经搜出来的答案还大，那么直接退出。

随机数据下复杂度是\(O(n\log n)\)的，构造下大概是\(O(n\sqrt n)\)
k维复杂度大概是\(O(kn^{1-{1\over k}})\)

这可以用来优化某些DP

先说这么多，如果有别的应用以后再补。

Code

一个相当丑的模板

namespace KDT
{
	int D,rt,t[N][2],sz[N],dw[N],tr[N],le;
	int a[N][3],rg[N][3][2],ans;
	bool cmp(int x,int y)
	{
		return a[x][D]<a[y][D];
	}
	void upd(const int &k,const int &y)
	{
		if(y) 
		{
			fo(j,0,2) 
			{
				rg[k][j][0]=min(rg[k][j][0],rg[y][j][0]);
				rg[k][j][1]=max(rg[k][j][1],rg[y][j][1]);
			}
			sz[k]=sz[k]+sz[y];
		}
	}
	void up(int k)
	{
		sz[k]=1;fo(j,0,2) rg[k][j][0]=rg[k][j][1]=a[k][j];
		if(t[k][0]) upd(k,t[k][0]);
		if(t[k][1]) upd(k,t[k][1]);
	}
	int make(int l,int r)
	{
		int mid=(l+r)>>1;
		D=(D+1)%3;
		std::nth_element(tr+l,tr+mid,tr+r+1,cmp);
		int k=tr[mid];dw[k]=D;
		t[k][0]=(l<mid)?make(l,mid-1):0;
		t[k][1]=(r>mid)?make(mid+1,r):0;
		up(k);
		return k;
	}
	void ins(int &k,int p)
	{	
		if(!k) 
		{
			k=p;sz[p]=1;dw[p]=(D+1)%3;
			fo(j,0,2) rg[p][j][0]=rg[p][j][1]=a[p][j];
			return;
		}
		D=dw[k];
		if(a[p][D]<=a[k][D]) ins(t[k][0],p);
		else ins(t[k][1],p);
		up(k);
	}
	int pri(int k,int x,int y,int z)
	{
		return rg[k][0][0]+sqr(max(0,max(rg[k][1][0]-y,y-rg[k][1][1])))+sqr(max(0,max(rg[k][2][0]-z,z-rg[k][2][1])));
	}
	int calc(int k,int x,int y,int z)
	{
		return a[k][0]+sqr(a[k][1]-y)+sqr(a[k][2]-z);
	}
	void query(int k,int x,int y,int z)
	{
		if(pri(k,x,y,z)>=ans) return;
		ans=min(ans,calc(k,x,y,z));
		if(t[k][0]) query(t[k][0],x,y,z);
		if(t[k][1]) query(t[k][1],x,y,z);
	}
}