KD - Tree 学习笔记

（这本是 P4169 的一篇无法提交的题解）

而 KD - Tree 本身在竞赛中的应用就比较少，又几乎不作为正解出现，所以这里只拿一道例题来举例讲解

题意：

给定一个二维点集 S ，每次有两种操作：

操作 1 ： 往点集 S 中加入一个点

操作 2 ： 给定一个点 $(x,y)$ ，询问该点到点集 S 中每个点曼哈顿距离的最小值

两点 $(x_1,y_1),(x_2,y_2)$ 之间的曼哈顿距离定义为：

$$|x_1-x_2|-|y_1-y_2|$$

题解：

这题正解是一个 $O((n+m){\log }_2^2(n+m))$ 复杂度的做法，需要用到 CDQ 分治

不过既然是 KD - Tree 的学笔，我们就用 KD - Tree 来做这道题

首先来介绍一下 KD - Tree ：

---

1. 定义 & 前置知识

KD - Tree ，原名 K - Dimension Tree ，是一种可以高效处理 k 维空间信息的数据结构

（然而 OI 中很少考到高维的 DS 题目）

（微信步数显然不是 DS 题目吧）

所以 KD - Tree 在 OI 中的应用大多数是二维（例如这道题），少数会有三维

而且 KD - Tree 据说可以配合其它的一些算法 比如乱搞 来获得较高分数

在学习 KD - Tree 之前，你需要了解如下东西：

• 平衡树（这里特指替罪羊树）

• STL 中的 nth_element 函数

没了，就这些，所以 KD - Tree 的入门是比较简单的，学会了上手就会很快

2. 建树过程

类似平衡树，KD - Tree 是一种 BST

但是和 BST 不同，在 k 维空间中，KD - Tree 每个结点参与比较的键值是有 k 个的（即每一维的坐标）

不过 KD - Tree 在检索的时候仍然能保证一定的次序

我们只需要在每一个深度都定义一个参与比较的维度，那么这样在检索时，每次究竟是往左子树走还是往右子树走，只要看这个点在当前检索到的深度的“比较维度”上的键值是多少即可

建树的算法已经给你了，关键是如何选择恰当的方式优化，避免像 BST 一样退化成链

那么有两个问题：

1. 每个深度怎样选维度？

2. 每个深度选好维度后怎样选点？

对于第一个问题，我们可以选择 1 ~ k 维轮流来

但是有一种更优的选择方式：选择方差最大的一维，具体证明不明，而且实现据说极为复杂

（识相点就选第一种吧，反正很少有出题人会想到要卡你）

对于第二个问题，显然选这个维度的中位数即可，这样左右两边的结点数是相等的

以上两种方式基本可以保证建树时间复杂度为 $O(n{\log }_{2}n)$，深度为 $O({\log }_{2}n)$

具体实现时，中位数可以使用 STL 中的 nth_element 函数，不会影响复杂度，用法请上网查找

这样讲会很抽象，给一张例图：

首先我们以横坐标为关键字，找到这一维的中位数

排序后即为：

$$A,H,C,B,E,D,G,F$$

我们选择中位数 B 作为根节点，此时可以把所有点分成两个部分：

那么我们递归处理左右子树，此时关键字变成纵坐标

左子树排好序后为： $A,C,H$，取中位数 C 作为 B 点的左儿子，继续递归 $\cdots \cdots$

右子树排好序后为： $D,F,G,E$，取中位数 F 做为 B 点的右儿子，继续递归 $\cdots \cdots$

那么建出来的 KD - Tree 如下图：

3. 插入 & 删除

插入操作很简单，类似 BST，我们直接通过检索找到新结点应该插入的位置即可

但是只要出题人想卡你，把新结点都放在一个位置，这样 KD - Tree 的深度越来越高，整颗树越来越不平衡，怎么办？

这个时候，替罪羊树的重构思想就派上用场了

我们是先钦定一个“平衡因子” $\alpha$（一般取 0.75 最优），当插入某个结点后，某个子树的大小超过了以它父亲为根的子树大小乘上 $\alpha$ 时，我们直接把整颗子树拿出来重构

表达式：

$$size(father)\ast \alpha <size(leftchildorrightchild)$$

但是每次可能会有多个满足条件的结点，取出深度最小的结点重构即可

（实际上本题对于每个满足条件的点都可以选择暴力重构，洛谷上不会超时）

删除操作稍微有点复杂，因为 KD - Tree 是比较偏静态的数据结构

我们可以采用惰性删除的方法，删除时先给这个点打标记，并更新它的祖先，到重构时再真正从树的结构里删去

然而如果一个子树中删除掉的结点过多，使得未被删去的结点数小于子树总结点数乘上 $\alpha$，我们同样需要重构（这样就可以爽快地删掉那些本该被删的结点了）

带重构的 KD - Tree，采用势能分析法，均摊时间复杂度仍然是 $O(n{\log }_{2}n)$ 的

---

回到这题，我们有一个显然的做法，就是先建好 KD - Tree，然后对于每次询问遍历整颗 KD - Tree，求出最短哈密顿距离

显然 TLE

显然要剪枝

我们需要考虑一种最优性剪枝

比如现在已经求得的答案是 $ans$，现在所在的结点是 $now$

我们一定要进入 $now$ 的左右子树吗？并不一定

类似 $A\ast$，我们可以用一个估价函数，求出一个理论上一定小于等于实际距离的最短距离

什么意思？就是说当前询问的点对于某颗子树估计的最短距离，一定要小于等于真实的最短距离

如果这个估计出来的最短距离大于等于 $ans$，那么显然实际的最短距离一定会大于等于这个估计距离，也就大于等于 $ans$，因此，我们进入这颗子树，对答案不会有任何贡献

如何估计呢？我们对于一颗子树，维护一个最小的，可以包括子树内所有点的矩形

那么这个点到这个矩形的最小曼哈顿距离，一定小于实际上的最小曼哈顿距离

这样的算法已经很优秀了，但是即使吸氧也会 T 掉一个点

我们还可以进行一个优化：

假设左子树的估价为 $leftvalue$，右子树的估价为 $rightvalue$

如果 $leftvalue<rightvalue$，我们实际上可以先递归左子树，再递归右子树

这样，当进入左子树时，我们有望可以获得更小的 $ans$ ，然后就有一定概率不会进入右子树

以上所有操作的理论复杂度为 $O((n+m)\sqrt{n+m})$，实际则不是很稳定，要注意常数不要过大

代码：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cctype>
#include<algorithm>

#define MAXN 600005
#define INF 0x3F3F3F3F

#define up(i, l, r) for(int i = (l); i <= (r); i = -~i)
#define down(i, l, r) for(int i = (l); i >= (r); i--)

using i64 = long long;
using namespace std;

namespace IO
{
	int f; char c;

	template<typename T> inline void _Read(T &v)
	{
		v = 0; f = 1; c = getchar();
		while(!isdigit(c)) { if(c == '-') f = -1; c = getchar(); }
		while(isdigit(c)) { v = (v << 3) + (v << 1) + (int)(c - '0'); c = getchar(); }
		v *= f;
		return;
	}

	template<typename T> inline void _Write(T k)
	{
		if(k < 0) { putchar('-'); k = -k; }
		if(k > 9) _Write(k / 10);
		putchar((char)(k % 10 + '0'));
		return;
	}

	inline int Read() { int v; _Read(v); return v; }
	inline i64 Read64() { i64 v; _Read(v); return v; }
	inline void Write(int v, char ed = '\n') { _Write(v); putchar(ed); return; }
	inline void Write64(i64 v, char ed = '\n') { _Write(v); putchar(ed); return; }
}

using IO :: Read;
using IO :: Read64;
using IO :: Write;
using IO :: Write64;

int n, m, cnt, rt, top, stk[MAXN];
bool wyh;
const double Alpha = 0.75;

struct Point { int a[2]; bool operator < (const Point &t) const { return a[wyh] < t.a[wyh]; } } p[MAXN]; 

struct KD_Tree
{
	int ls[MAXN], rs[MAXN], si[MAXN], mx[MAXN][2], mi[MAXN][2];
	Point val[MAXN];

	inline int New_Node() { return top ? stk[top--] : ++cnt; }

	inline int Pushup(int now)
	{
		up(i, 0, 1)
		{
			mx[now][i] = mi[now][i] = val[now].a[i];
			if(ls[now]) mx[now][i] = max(mx[now][i], mx[ls[now]][i]), mi[now][i] = min(mi[now][i], mi[ls[now]][i]);
			if(rs[now]) mx[now][i] = max(mx[now][i], mx[rs[now]][i]), mi[now][i] = min(mi[now][i], mi[rs[now]][i]);
		}
		return si[now] = si[ls[now]] + si[rs[now]] + 1, now;
	}

	inline int Build(int l, int r, int dep)
	{
		wyh = dep ^ 1; int mid = (l + r) >> 1; nth_element(p + l, p + mid, p + r + 1);
		int now = New_Node(); val[now] = p[mid];
		ls[now] = l < mid ? Build(l, mid - 1, dep + 1) : 0, rs[now] = mid < r ? Build(mid + 1, r, dep + 1) : 0;
		return Pushup(now);
	}

	inline void Pia(int now, int tot)
	{
		if(ls[now]) Pia(ls[now], tot);
		p[tot + si[ls[now]] + 1] = val[now], stk[++top] = now;
		if(rs[now]) Pia(rs[now], tot + si[ls[now]] + 1);
		return;
	}

	inline void Check(int &now, int dep)
	{
		if((double)si[ls[now]] > Alpha * si[now] || (double)si[rs[now]] > Alpha * si[now]) Pia(now, 0), now = Build(1, si[now], dep);
		return;
	}

	inline void Insert(int &now, Point ins, int dep)
	{
		if(!now) return now = New_Node(), val[now] = ins, ls[now] = rs[now] = 0, Pushup(now), void();
		if(ins.a[dep ^ 1] <= val[now].a[dep ^ 1]) Insert(ls[now], ins, dep + 1);
		else Insert(rs[now], ins, dep + 1);
		return Pushup(now), Check(now, dep);
	}

	inline int Dist(Point p, Point q) { return abs(p.a[0] - q.a[0]) + abs(p.a[1] - q.a[1]); }

	inline int Get_Dis(Point k, int c)
	{
		int res = 0;
		up(i, 0, 1) res += max(0, k.a[i] - mx[c][i]) + max(0, mi[c][i] - k.a[i]);
		return res;
	}

	inline void Query(int now, Point qry, int &ans)
	{
		ans = min(ans, Dist(qry, val[now]));
		int L = INF, R = INF;
		if(ls[now]) L = Get_Dis(qry, ls[now]); 
		if(rs[now]) R = Get_Dis(qry, rs[now]);
		if(L < R)
		{
			if(L < ans) Query(ls[now], qry, ans);
			if(R < ans) Query(rs[now], qry, ans);
		}
		else
		{
			if(R < ans) Query(rs[now], qry, ans);
			if(L < ans) Query(ls[now], qry, ans);
		}
		return;
	}
} kdt;

int main()
{
	n = Read(), m = Read();
	up(i, 1, n) p[i].a[0] = Read(), p[i].a[1] = Read();
	rt = kdt.Build(1, n, 1);
	while(m--)
	{
		int op = Read();
		Point tmp; tmp.a[0] = Read(), tmp.a[1] = Read();
		if(op & 1) kdt.Insert(rt, tmp, 1);
		else
		{
			int ans = INF;
			kdt.Query(rt, tmp, ans);
			Write(ans);
		}
	}
	return 0;
}

$END$