KDT小记
command_block
2020-10-23 15:17:23
$\rm KD-tree$ 是一种能对 $k$ 维空间中的**点**进行快速检索的树形数据结构。
相对另一个传统方法《树套树》而言，有以下优势 ：
- 部分情况下代码实现更简单
- 空间复杂度 $O(n)$
- 可以搭配传统的整体标记系统。
$\rm KDT$ 搭配剪枝乱搞后会在某些问题中取得奇效，特别是在考场上收益率特别高。
# $\rm KD-tree$ 的构造
$\rm KDT$ 是一棵二叉排序树，每个节点有一个管辖(超)矩形。保证 $u$ 子树中的所有点都在 $u$ 的管辖矩形内。
左右儿子的管辖矩形是由父亲的管辖矩形水平划分而得。$\rm KDT$ 的核心思想就是对矩形的水平分割。
在满足上述要求的前提下，$\rm KDT$ 还要使得高度严格为 $\log_2n+O(1)$ 。
这是一个 $\rm KDT$ 以及平面划分的例子。
- 具体构造流程 :
挑选一个维度，使用与该维度正交的超平面来分割当前区域。
对于二维情况，即选择 $x,y$ 其中一维，使用水平直线 $y=c$ 或竖直直线 $x=c$ 来划分当前矩形。
不妨设我们在 $x$ 坐标上划分。
为了分割得尽量均匀，我们令 $c$ 为当前区域中，点的 $x$ 坐标的中位数，这样恰能将点集分成两半。
将当前区域以及点集划分成两部分后，继续分治。
如上图，第一次在 $x$ 坐标划分选中了 $①$ 号点，划分线为图中红线。
关于分割维度的选择，有如下两种常见的方案 ：
- 轮换划分 ： 每个维度轮着来。
- 方差划分 ： 优先划分方差较大的维度。
后者实现略复杂，但是传说有奇妙加速，具体原理不明。
寻找中位数时，使用 `nth_element` 指令可以做到线性。
这样，建树的总复杂度就是 $T(n)=2T(n/2)+O(n)=O(n\log n)$。
# $\rm KD-tree$ 的 $\rm range\ query$
在不利用差分的前提下，支持将一个矩形区域内的点，划分成 $O(\sqrt{n})$ 个点负责的区域。
做法非常暴力，若查询范围包含当前范围，则取当前范围。
- **复杂度分析**
假设我们查询矩形 $R$ ,将树上的节点按照管辖矩形分类 ：
1. 与 $R$ 无交集。
2. 完全包含在 $R$ 内。
3. 部分包含在 $R$ 内。
观察查询算法的行为，其一路经过 $3$ 类点，且一旦遇到 $1,2$ 类点就停下。
而且，$1,2$ 类点的子树中显然无三类点，这样，查询的复杂度和 $3$ 类点的个数直接相关。
对于矩形的每条边，分别计算其“穿过”的节点个数。
我们以轮换划分来分析。在相邻的**两轮**中，分别对 $x,y$ 进行了划分。
能够发现，在划分出的四个区域中，每个区域的大小都是原来的 $1/4$ ，任意的一条直线只能穿过其中两个。
这样，设 $T(n)$ 为在大小为 $n$ 的 $\rm KDT$ 中一条直线所经过的区域数，那么有
$T(n)=2T(n/4)+O(1)$ 即 $T(n)=O(\sqrt{n})$。
所以，一条直线经过的区域个数为 $O(\sqrt{n})$。
**附** : 有一点细节，我们的划分出的矩形区域并非严格不交，而在边界上是相交的。
如果一条直线恰好于边界重合，其可能经过多于 $2$ 个区域。
但是，如果出现这种情况，递归向儿子时，直线总是在矩形的一侧，则不会再次经过多于 $2$ 个区域，复杂度仍是正确的。
如果想避免这种情况，将坐标乘以 $2$ ，并将矩形略微放大，使得没有点恰在矩形边上即可。
**扩展①** : 若在 $k$ 维空间中，连续的 $k$ 轮划分会剖分出 $2^k$ 个部分。
一个用于划分的超平面最多穿过其中 $2^{k-1}$ 个部分（可以归纳证明）。
这样，复杂度递推式就为 $T(h)=2^{k-1}T(h-k)+O(1)$ ,即 $T(h)=2^{h-(h/k)}$
所以，超平面经过的区域个数为 $O(n^{1-(1/k)})$。
**扩展②** : 有时候，两个维度的操作个数不同，我们可以改变 $\rm KDT$ 的结构来改进复杂度。
若在相邻的两轮划分中， $x$ 上分成 $a$ 段，$y$ 上分成 $b$ 段，则有 $ab$ 份，每一份大小为 $n/ab$。
对于垂直于 $x$ 轴的一条直线，其经过区域数为 $T_1(n)=bT_1(n/ab)+O(1)$。
对于垂直于 $y$ 轴的一条直线，其经过区域数为 $T_2(n)=aT_2(n/ab)+O(1)$。
将两者相乘，可得 $T_1(n)T_2(n)=abT_1(n/ab)T_2(n/ab)$ （和 $O(1)$ 相关的部分可不计）
这样，不难看出 $T_1(n)T_2(n)=O(n)$ ，即两个维度操作复杂度乘积为 $O(n)$。
（注意，若 $a,b$ 较大，复杂度也会退化）
对于某个问题，若 $x$ 维度上有 $O(n)$ 个修改， $y$ 维度上有 $O(n\log^2n)$ 个查询。
则可以让 $x$ 上复杂度为 $O(\sqrt{n}\log n)$ ，$y$ 上复杂度 $O(\sqrt{n}/\log n)$ ，可得总复杂度为 $O(n\sqrt{n}\log n)$。
具体操作举例 ： 相邻的三次划分中， $x$ 占两次， $y$ 占一次，则有 ：
$x:T_1(n)=2T(n/8)+O(1)=O(n^{1/3})$
$y:T_1(n)=4T(n/8)+O(1)=O(n^{2/3})$
- **例题** ： [CF44G Shooting Gallery](https://www.luogu.com.cn/problem/CF44G)
**题意** ： 给出一些矩形，每个矩形都有一个高度，保证高度互不相同。
有若干次从下到上的射击，会击中并摧毁射击点上最低的那个矩形。
求出每次射击击中的矩形，或指出未射中。
考虑使用矩形来匹配射击点。
可以把矩形按高度排序后，依次寻找范围内的编号最小的射击点，并将其删除，正确性显然。
先对射击点建立 $\rm KDT$ ，维护子树内未被删除的最小射击点。`range query` 之后单点修改即可。
剪枝 : 若子树内编号最小点仍然大于当前答案，则不进入。优先进入最小点编号更小的子树。(大幅优化)
由于我个人比较偏好 $\rm Leafy$ 数据结构，下面是“线段树式”的 $\rm KDT$ 实现。
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 100500
using namespace std;
struct Node
{int xl,yl,xr,yr,p;}a[MaxN<<2];
struct Point{int x,y,p;}p[MaxN];
bool cmpX(const Point &A,const Point &B)
{return A.x<B.x;}
bool cmpY(const Point &A,const Point &B)
{return A.y<B.y;}
void build(int l,int r,int u,bool kd)
{
if (l==r){
a[u].xl=a[u].xr=p[l].x;
a[u].yl=a[u].yr=p[l].y;
a[u].p=p[l].p;
return ;
}int mid=(l+r)>>1;
nth_element(p+l,p+mid,p+r+1,kd ? cmpX : cmpY);
build(l,mid,u<<1,!kd);
build(mid+1,r,u<<1|1,!kd);
int ls=u<<1,rs=u<<1|1;
a[u].p=min(a[ls].p,a[rs].p);
a[u].xl=min(a[ls].xl,a[rs].xl);
a[u].yl=min(a[ls].yl,a[rs].yl);
a[u].xr=max(a[ls].xr,a[rs].xr);
a[u].yr=max(a[ls].yr,a[rs].yr);
}
int to,n;
void upd(int l=1,int r=n,int u=1)
{
if (l==r){a[u].p=n+1;return ;}
int mid=(l+r)>>1;
if (to<=mid)upd(l,mid,u<<1);
else upd(mid+1,r,u<<1|1);
a[u].p=min(a[u<<1].p,a[u<<1|1].p);
}
Node wf;
void qry(int u=1)
{
if (a[u].p>=wf.p||wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl|a[u].yr<wf.yl)return ;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
{wf.p=min(wf.p,a[u].p);return ;}
if (a[u<<1].p<a[u<<1|1].p)
{qry(u<<1);qry(u<<1|1);}
else {qry(u<<1|1);qry(u<<1);}
}
struct Squ{int xl,yl,xr,yr,p,h;}s[MaxN];
bool cmpS(const Squ &A,const Squ &B)
{return A.h<B.h;}
int m,tp[MaxN],ans[MaxN];
int main()
{
scanf("%d",&m);
for (int i=1;i<=m;i++){
scanf("%d%d%d%d%d",&s[i].xl,&s[i].xr,&s[i].yl,&s[i].yr,&s[i].h);
s[i].p=i;
}
sort(s+1,s+m+1,cmpS);
scanf("%d",&n);
for (int i=1;i<=n;i++){
scanf("%d%d",&p[i].x,&p[i].y);
p[i].p=i;
}
build(1,n,1,0);
for (int i=1;i<=n;i++)tp[p[i].p]=i;
for (int i=1;i<=m;i++){
wf=(Node){s[i].xl,s[i].yl,s[i].xr,s[i].yr,n+1};
qry();
if (wf.p==n+1)continue;
ans[wf.p]=s[i].p;
to=tp[wf.p];upd();
}
for (int i=1;i<=n;i++)
printf("%d\n",ans[i]);
return 0;
}
```
# $\rm KD-tree$ 配合标记
$\rm KDT$ 和树套树有本质的不同，其能把 $k$ 维空间内的一个区域映射成 $O(n^{1-(1/k)})$ 个区间的并。（类似于树剖？）
所以，传统的（线段树）标记体系可以应用于 $\rm KDT$ 上。
为了减少细节，可以把 $\rm KDT$ 写成 $\rm Leafy$ 的。
- **例题①** : BZOJ4154 Generating Synergy
**题意** ： 给定一棵有根树,支持下列操作 ：
- 将距离节点 $u$ 不超过 $l$ 的子节点染成颜色 $c$。
- 或询问点 $u$ 的颜色。
不难发现，子树的限制是 $dfs$ 序上的一个区间，而深度的限制是 $dep$ 上的一个前缀。
所以，问题相当于矩形覆盖，单点查询，使用 $\rm KDT$ 即可。
------------
- **例题②** ： [P6349 [PA2011]Kangaroos](https://www.luogu.com.cn/problem/P6349)
区间 $[L,R]$ 与区间 $[l,r]$ 有交集的充要条件是 $[l\leq R][L\leq r]$ ，是二维偏序。
考虑离线，将每个询问看作二位平面上的点 $(l,r)$ ，建立 $\rm KDT$。
设 $f[i][j]$ 为 ： 第 $j$ 个询问区间，从第 $i$ 个序列区间开始向前，能连续有交的区间个数。
每加入一个序列区间 $L,R$ ，则满足 $[l\leq R][L\leq r]$ 的询问区间 $f$ 值加一，其余清零。
$f[?][j]$ 的（历史）最大值即为询问 $j$ 的答案。
- **思考题** ：[[DS记录]P6783 [Ynoi2008]rrusq](https://www.luogu.com.cn/blog/command-block/ds-ji-lu-p6783-ynoi2008rrusq)
# $\rm KD-tree$ 的插入
- **例题①** ：[P4148 简单题](https://www.luogu.com.cn/problem/P4148)
- **局部重构**
一种广为流传的做法是 ： 首先暴力插入，然后模仿替罪羊树，当左右子树不平衡时重构部分子树。
大部分人把 $\alpha$ 设为 $0.75$ 左右，并分析出此时重构的复杂度是 $O(n\log^2n)$。
这时不能想当然地认为查询复杂度仍然为 $O(\sqrt{n})$。
树的高度并非 $\log_2 n+O(1)$ ,而是 $\log_{4/3}n+O(1)$ ，为原来的两倍还多。
$h>2\log_2 n\Rightarrow$ 查询复杂度 $O(2^{h/2})>O(2^{\log_2n})=O(n)$ ，似乎基本上没保证了。
实际上还是有更紧的界。不妨假设整颗树都是最不平衡的情况，即每个点的大儿子大小都是总大小的 $\alpha$ 倍。
$\rm KDT$ 在向下的两层产生的四个分支中，只会选择两个分支。
四个分支和原树的大小比例如图 :
这样，最坏情况下只能选出大小为 $\alpha^2,\alpha(1-\alpha)$ 的两个分支。
设 $T(n)$ 为大小为 $n$ 的树的查询复杂度，有递推式 $T(n)=T(\alpha^2n)+T(\alpha(1-\alpha)n)+O(1)$
当 $\alpha=0.75$ 时, $T(n)=T(9n/16)+T(3n/16)+O(1)$
设 $T(n)=O(n^c)$ ，则有 $n^c=(9n/16)^c+(3n/16)^c$ ,即 $(9/16)^c+(3/16)^c=1$ 解得 $c≈0.68$。
也就是说，实际上是大约 $O(n^{0.68})$ 的复杂度，而且很难卡满。
或者 ,我们取更小的 $\alpha=3/5$ ,此时有 $T(n)=T(9n/25)+T(6n/25)+O(1)$ ,解得 $T(n)≈O(n^{0.57})$。
这个复杂度仍然很难卡满，已经能够令人满意了。
重构的时间消耗则约为 $O(\frac{1-3/5}{6/5-1}n\log n\log_{5/3}n)$ 大概是 $2\sim 4$ 倍常数的样子。
附 : [替罪羊树复杂度简略分析](https://www.luogu.com.cn/paste/v5fxeowc)
当然，有条件的话，根据实际情况调 $\alpha$ 当然是最好的。
下面是 $\alpha=3/5$ 的代码。[评测记录](https://www.luogu.com.cn/record/40478991)
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 200500
using namespace std;
int read(){}
struct Point{int x,y,v;}p[MaxN];
int tot;
bool cmpX(const Point &A,const Point &B)
{return A.x<B.x;}
bool cmpY(const Point &A,const Point &B)
{return A.y<B.y;}
struct Node{
int xl,yl,xr,yr,l,r,s,c;
void leaf(Point &P){
xl=xr=P.x;yl=yr=P.y;
s=P.v;c=1;
}
}a[MaxN<<1];
int tn,rt,rub[MaxN<<1],tb;
int cre(){
if (!tb)return ++tn;
int u=rub[tb--];
a[u]=(Node){0,0,0,0,0,0};
return u;
}
void pia(int u)
{
if (a[u].l==a[u].r){
p[++tot]=(Point){a[u].xl,a[u].yl,a[u].s};
rub[++tb]=u;return ;
}pia(a[u].l);pia(a[u].r);
rub[++tb]=u;
}
bool chk(int u)
{return a[u].c>3&&max(a[a[u].l].c,a[a[u].r].c)*5>a[u].c*3;}
inline void up(int u)
{
int l=a[u].l,r=a[u].r;
a[u].s=a[l].s+a[r].s;
a[u].c=a[l].c+a[r].c;
a[u].xl=min(a[l].xl,a[r].xl);
a[u].yl=min(a[l].yl,a[r].yl);
a[u].xr=max(a[l].xr,a[r].xr);
a[u].yr=max(a[l].yr,a[r].yr);
}
void build(int l,int r,int &u,bool kd)
{
u=cre();
if (l==r){a[u].leaf(p[l]);return ;}
int mid=(l+r)>>1;
nth_element(p+l,p+mid,p+r+1,kd ? cmpX : cmpY);
build(l,mid,a[u].l,!kd);
build(mid+1,r,a[u].r,!kd);
up(u);
}
Point wfp,reb;
void ins(int u,bool kd)
{
if (a[u].l==a[u].r){
int v0=cre(),v1=cre();
a[v0]=a[u];a[v1].leaf(wfp);
a[u].l=v0;a[u].r=v1;up(u);
return ;
}int ls=a[u].l;
if (a[ls].xl<=wfp.x&&wfp.x<=a[ls].xr
&&a[ls].yl<=wfp.y&&wfp.y<=a[ls].yr)
ins(ls,!kd);
else ins(a[u].r,!kd);
up(u);
if (chk(u)){reb.x=u;reb.y=kd;}
}
void insert()
{
if (!tn){a[tn=1].leaf(wfp);return ;}
reb.x=0;ins(1,0);
if (reb.x){
tot=0;pia(reb.x);
int tmp;build(1,tot,tmp,reb.y);
}
}
Node wf;
void qry(int u=1)
{
if (wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl|a[u].yr<wf.yl)return ;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
{wf.s+=a[u].s;return ;}
qry(a[u].l);qry(a[u].r);
}
int n,las;
int main()
{
while(1){
int op=read();
if (op==1){
wfp.x=read()^las;
wfp.y=read()^las;
wfp.v=read()^las;
insert();
}if (op==2){
wf.xl=read()^las;wf.yl=read()^las;
wf.xr=read()^las;wf.yr=read()^las;
wf.s=0;qry();
printf("%d\n",las=wf.s);
}if (op==3)break;
}return 0;
}
```
- **根号分治**
设一个阈值 $B$ ，当积攒了 $B$ 个插入后，重构整棵树。
重构复杂度为 $O(n^2\log n/B)$ 查询复杂度为 $O(B+\sqrt{n})$ ,取 $B=\sqrt{n\log n}$ 可得复杂度为 $O(n\sqrt{n\log n})$。
如果不想带 $\log$，可以开两棵 $\rm KDT$。
插入的点先暂存，积攒了 $\sqrt{n}$ 个插入后，和第一棵树合并。
积攒了 $B$ 个插入后，和第二棵树合并。
查询需要在大小 $\rm KDT$ 中查询，以及查看散点，复杂度为 $O(\sqrt{n})$。
合并的复杂度为 $O(B\sqrt{n}\log n+n^2\log n/B)$ ，取 $b=n^{3/4}$ 可得复杂度为 $O(n^{5/4}\log n)$。
这样一来就得到了一个均摊 $O(n\sqrt{n})$ 的算法。
常数较小，跑得比局部重构要快得多（似乎也更好写）。
$O(n\sqrt{n\log n})$ [评测记录](https://www.luogu.com.cn/record/40488192)
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 100500
using namespace std;
int read(){}
struct Node
{int xl,yl,xr,yr,s;}a[MaxN<<2];
struct Point{int x,y,v;}p[MaxN];
bool cmpX(const Point &A,const Point &B)
{return A.x<B.x;}
bool cmpY(const Point &A,const Point &B)
{return A.y<B.y;}
void build(int l,int r,int u,bool kd)
{
if (l==r){
a[u].xl=a[u].xr=p[l].x;
a[u].yl=a[u].yr=p[l].y;
a[u].s=p[l].v;
return ;
}int mid=(l+r)>>1;
nth_element(p+l,p+mid,p+r+1,kd ? cmpX : cmpY);
build(l,mid,u<<1,!kd);
build(mid+1,r,u<<1|1,!kd);
int ls=u<<1,rs=u<<1|1;
a[u].s=a[ls].s+a[rs].s;
a[u].xl=min(a[ls].xl,a[rs].xl);
a[u].yl=min(a[ls].yl,a[rs].yl);
a[u].xr=max(a[ls].xr,a[rs].xr);
a[u].yr=max(a[ls].yr,a[rs].yr);
}
Node wf;
void qry(int u=1)
{
if (wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl||a[u].yr<wf.yl)return ;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
{wf.s+=a[u].s;return ;}
qry(u<<1|1);qry(u<<1);
}
int main()
{
int pl=0,pr=0,n=read(),las=0;
while(1){
int op=read();
if (op==1){
pr++;
p[pr].x=read()^las;
p[pr].y=read()^las;
p[pr].v=read()^las;
if ((pr-pl)*(pr-pl)>pr*40)build(1,pl=pr,1,0);
}if (op==2){
wf.xl=read()^las;wf.yl=read()^las;
wf.xr=read()^las;wf.yr=read()^las;
wf.s=0;qry();
for (int i=pl+1;i<=pr;i++)
if (wf.xl<=p[i].x&&p[i].x<=wf.xr&&wf.yl<=p[i].y&&p[i].y<=wf.yr)
wf.s+=p[i].v;
printf("%d\n",las=wf.s);
}if (op==3)break;
}return 0;
}
```
$O(n\sqrt{n}+n^{5/4}\log n)$ [评测记录](https://www.luogu.com.cn/record/40510598)
```cpp
#include<algorithm>
#include<cstdio>
#include<cmath>
#define MaxN 100500
using namespace std;
int read(){
int X=0;char ch=0,w=0;
while(ch<48||ch>57)ch=getchar(),w|=(ch=='-');
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return w?-X:X;
}
struct Point{int x,y,v;}p[MaxN];
bool cmpX(const Point &A,const Point &B)
{return A.x<B.x;}
bool cmpY(const Point &A,const Point &B)
{return A.y<B.y;}
struct Node{int xl,yl,xr,yr,s;};
Node wf;
struct KDT
{
Node a[MaxN<<2];
void build(int l,int r,int u,bool kd)
{
if (l==r){
a[u].xl=a[u].xr=p[l].x;
a[u].yl=a[u].yr=p[l].y;
a[u].s=p[l].v;
return ;
}int mid=(l+r)>>1;
nth_element(p+l,p+mid,p+r+1,kd ? cmpX : cmpY);
build(l,mid,u<<1,!kd);
build(mid+1,r,u<<1|1,!kd);
int ls=u<<1,rs=u<<1|1;
a[u].s=a[ls].s+a[rs].s;
a[u].xl=min(a[ls].xl,a[rs].xl);
a[u].yl=min(a[ls].yl,a[rs].yl);
a[u].xr=max(a[ls].xr,a[rs].xr);
a[u].yr=max(a[ls].yr,a[rs].yr);
}
void qry(int u=1)
{
if (wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl||a[u].yr<wf.yl)return ;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
{wf.s+=a[u].s;return ;}
qry(u<<1|1);qry(u<<1);
}
}T0,T1;
int main()
{
int p0=0,p1=0,p2=0,n=read(),las=0;
while(1){
int op=read();
if (op==1){
p2++;
p[p2].x=read()^las;
p[p2].y=read()^las;
p[p2].v=read()^las;
if ((p2-p1)*(p2-p1)>p2*1.5)
T1.build(p0+1,(p1=p2),1,0);
else if (p1-p0>1.2*pow(p2,0.75)){
T0.build(1,p0=p1,1,0);
T1.build(p0+1,(p1=p2),1,0);
}
}if (op==2){
wf.xl=read()^las;wf.yl=read()^las;
wf.xr=read()^las;wf.yr=read()^las;
wf.s=0;T0.qry();T1.qry();
for (int i=p1+1;i<=p2;i++)
if (wf.xl<=p[i].x&&p[i].x<=wf.xr&&wf.yl<=p[i].y&&p[i].y<=wf.yr)
wf.s+=p[i].v;
printf("%d\n",las=wf.s);
}if (op==3)break;
}return 0;
}
```
------------
- **例题②** ： [P4848 崂山白花蛇草水](https://www.luogu.com.cn/problem/P4848)
- **做法1** :
使用 $\rm KDT$ 套权值线段树来维护，查询时多树二分。
建树的复杂度可以用线段树合并做到 $O(q\log w)$，查询的复杂度是 $O(\sqrt{q}\log w)$ ，且常数较大。
使用上文的根号分治方法，总时间复杂度为 $O(q\sqrt{q}\log w)$，空间复杂度为 $O(q\log w)$。
实现较为繁琐，且不能通过本题。
- **做法2** :
可以改用权值线段树套 $\rm KDT$ ，这样时间复杂度仍然是 $O(q\sqrt{q}\log w)$。
注意到，在线段树上二分时，我们可以选择查询左儿子的 $\rm KDT$ 还是右儿子的 $\rm KDT$。(具体做法见代码)
我们选择较小的那一侧来查询，这样时间开销较低。
最坏情况下，我们每次都会进入点较多的子树。
这相当于把点集分成 $O(\log w)$ 个部分，然后分别查询，查询复杂度降为 $O(\log w*\sqrt{q/\log w})=O(\sqrt{q\log w})$。
重构的复杂度是 $O(q^{5/4}\log q\log w)$。
由于本题数据随机，很多代码都将复杂度向查询倾斜。实际上，这在极限数据下表现并不好。
同时，由于本题时限较紧（ $\log^3$ 跑 $10^5$ ）,需利用数据特性，将重构次数减少才能通过。
```cpp
#include<algorithm>
#include<cstdio>
#include<vector>
#include<cmath>
#define pb push_back
#define MaxN 100500
using namespace std;
int read(){}
struct Point{int x,y;}sp[MaxN];
bool cmpX(const Point &A,const Point &B)
{return A.x<B.x;}
bool cmpY(const Point &A,const Point &B)
{return A.y<B.y;}
struct Node{int xl,yl,xr,yr,c,l,r;}a[MaxN*80];
int tn,rub[MaxN*80],tb;
int cre(){
if (!tb)return ++tn;
int u=rub[tb--];
a[u]=(Node){0,0,0,0,0,0};
return u;
}
Node wf;
struct KDT
{
int rt;
void del(int u){
if (!u)return ;
rub[++tb]=u;
del(a[u].l);del(a[u].r);
}
void build(int l,int r,int &u,bool kd)
{
a[u=cre()].c=r-l+1;
if (l==r){
a[u].xl=a[u].xr=sp[l].x;
a[u].yl=a[u].yr=sp[l].y;
return ;
}int mid=(l+r)>>1;
nth_element(sp+l,sp+mid,sp+r+1,kd ? cmpX : cmpY);
build(l,mid,a[u].l,!kd);
build(mid+1,r,a[u].r,!kd);
int ls=a[u].l,rs=a[u].r;
a[u].xl=min(a[ls].xl,a[rs].xl);
a[u].yl=min(a[ls].yl,a[rs].yl);
a[u].xr=max(a[ls].xr,a[rs].xr);
a[u].yr=max(a[ls].yr,a[rs].yr);
}
void qry(int u)
{
if (wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl||a[u].yr<wf.yl)return ;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
{wf.c+=a[u].c;return ;}
qry(a[u].l);qry(a[u].r);
}
};
struct SquDS
{
vector<Point> p;
KDT T0,T1;int p0,p1;
void insert(Point tp)
{
p.pb(tp);
if ((p.size()-p1)*(p.size()-p1)>p.size()*48){
p1=p.size();
for (int i=p0;i<p1;i++)sp[i]=p[i];
T1.del(T1.rt);T1.build(p0,p1-1,T1.rt,0);
}else if (p1-p0>pow(p.size()*7,0.75)){
for (int i=0;i<p.size();i++)sp[i]=p[i];
T0.del(T0.rt);T0.build(0,(p0=p1)-1,T0.rt,0);
T1.del(T1.rt);T1.build(p0,(p1=p.size())-1,T1.rt,0);
}
}
int qry()
{
wf.c=0;T0.qry(T0.rt);T1.qry(T1.rt);
for (int i=p1;i<p.size();i++)
wf.c+=(wf.xl<=p[i].x&&p[i].x<=wf.xr&&wf.yl<=p[i].y&&p[i].y<=wf.yr);
return wf.c;
}
};
struct SGT_Node
{int l,r;SquDS s;}t[MaxN<<5];
int tot,to;Point wfp;
#define Lim 1000000000
void add(int l,int r,int &u)
{
if (!u)u=++tot;
t[u].s.insert(wfp);
if (l==r)return ;
int mid=(l+r)>>1;
if (to<=mid)add(l,mid,t[u].l);
else add(mid+1,r,t[u].r);
}
int qry(int l,int r,int u,int c,int k)
{
if (l==r)return l;
int mid=(l+r)>>1,ls=t[u].l,rs=t[u].r,lc,rc;
if (t[ls].s.p.size()<t[rs].s.p.size())
{lc=t[ls].s.qry();rc=c-lc;}
else {rc=t[rs].s.qry();lc=c-rc;}
if (k>lc)
return qry(mid+1,r,rs,rc,k-lc);
return qry(l,mid,ls,lc,k);
}
int n,q,las,SGTrt;
int main()
{
n=read();q=read();
for (int i=1,op;i<=q;i++){
op=read();
if (op==1){
wfp.x=read()^las;wfp.y=read()^las;
to=read()^las;add(1,Lim,SGTrt);
}else {
wf.xl=read()^las;wf.yl=read()^las;
wf.xr=read()^las;wf.yr=read()^las;
int c=t[SGTrt].s.qry(),k=read()^las;
if (k>c){
puts("NAIVE!ORZzyz.");
las=0;continue;
}printf("%d\n",las=qry(1,Lim,SGTrt,c,c-k+1));
}
}return 0;
}
```
- **做法3** :
若我们的线段树恰好在每个叶子有一个点(分布均匀)，那么我们一次询问所遇到的 $\rm KDT$ 点数分别是 $q,q/2,q/4...$
此时查询复杂度 $\sum\limits_{i=0}O(\sqrt{q/2^i})=O(\sqrt{q})$，插入亦然。
当然，实际数据可能往邻近值域塞入较多点，这样我们每次都会遇到很大棵的 $\rm KDT$ 。
外层可以模仿替罪羊树的重构，这样保证了深度加深时划分的较为平均，但是需要支付 $O(q\log^3q)$ 的重构代价，以及带来巨大常数。(所以本做法并不实用)
.
- **做法4** :
考虑不定叉的权值线段树，给定常数 $k$ ,第零层分 $2^1$ 叉，第一层分 $2^{k-1}$ 叉，第二层分 $2^{k^2-k}$ 叉……
以此类推,第 $t$ 层分 $2^{k^t-k^{t-1}}$ 叉。
若共有 $n$ 个叶子，对于第 $m$ 层的点 $u$ :
- 兄弟数量为 $2^{\sum\limits_{i=0}^{m-1}k^i-k^{i-1}}=O(2^{k^{m-1}})$
- 有 $2^{k^{m}-k^{m-1}}$ 个分支
- 一个分支的大小为 $n/2^{k^m}$ 级别。
这样，在第 $m$ 层消耗的查询的复杂度为 $O\Big(2^{k^{m}-k^{m-1}}*\sqrt{n/2^{k^m}}\Big)=O(\sqrt{n}*2^{k^{m-1}(k/2-1)})$。
总复杂度也就是 $O\Big(\sqrt{n}*\sum\limits_{i=0}^{2^{k^m}\leq n}2^{k^{i-1}(k/2-1)}\Big)$
若取 $k=2$ ,则有 $k/2-1=0$ ,后半部分是 $\sum\limits_{i=0}^{2^{2^m}\leq n}2^0=O(\log\log n)$
这颗树一共有 $O(\log\log n)$ 层，所以查询复杂度为 $O(\sqrt{n}\log\log n)$ ，重构复杂度为 $O(n\log n\log\log n)$。
使用根号分治，以 $B$ 为阈值，查询的复杂度为 $O(\sqrt{q}\log\log n+B)$ ,重构的复杂度为 $O(q^2\log q\log\log q/B)$
取 $B=\sqrt{q\log q\log\log q}$ 可以做到 $O(q\sqrt{q\log q\log\log q})$ 的复杂度。空间低至 $O(q\log\log q)$。
本做法只需要写静态 $\rm KDT$ ，而且瓶颈在散点处理，常数较小。
若取 $1<k<2$ ，可能可以获得更好的理论复杂度。
```cpp
#include<algorithm>
#include<cstdio>
#include<vector>
#include<cmath>
#define pb push_back
#define MaxN 100500
using namespace std;
int read(){}
struct Point{int x,y,v;}p[MaxN],sp[MaxN];
bool cmpX(const Point &A,const Point &B){return A.x<B.x;}
bool cmpY(const Point &A,const Point &B){return A.y<B.y;}
bool cmpV(const Point &A,const Point &B){return A.v<B.v;}
struct Node{int xl,yl,xr,yr,c;}a[MaxN*4*4];
Node *tn,wf;
struct KDT
{
Node *a;
void build(int l,int r,int u,bool kd)
{
a[u].c=r-l+1;
if (l==r){
a[u].xl=a[u].xr=sp[l].x;
a[u].yl=a[u].yr=sp[l].y;
return ;
}int mid=(l+r)>>1;
nth_element(sp+l,sp+mid,sp+r+1,kd ? cmpX : cmpY);
build(l,mid,u<<1,!kd);
build(mid+1,r,u<<1|1,!kd);
int ls=u<<1,rs=u<<1|1;
a[u].xl=min(a[ls].xl,a[rs].xl);
a[u].yl=min(a[ls].yl,a[rs].yl);
a[u].xr=max(a[ls].xr,a[rs].xr);
a[u].yr=max(a[ls].yr,a[rs].yr);
}
int qry(int u=1)
{
if (wf.xr<a[u].xl||a[u].xr<wf.xl||wf.yr<a[u].yl||a[u].yr<wf.yl)return 0;
if (wf.xl<=a[u].xl&&a[u].xr<=wf.xr&&wf.yl<=a[u].yl&&a[u].yr<=wf.yr)
return a[u].c;
return qry(u<<1)+qry(u<<1|1);
}
};
struct SGT_Node{int mid[12];KDT s;}t[205];
const int d[6]={0,2,4,10,1251};
void build(int l,int r,int u,int dep)
{
if (l>r)return ;
int k=d[dep];
t[u].s.a=tn;tn+=(r-l+1)<<2;
for (int i=l;i<=r;i++)sp[i]=p[i];
t[u].s.build(l,r,1,0);
if (k>r-l+1)return ;
for (int i=0;i<=k;i++)
t[u].mid[i]=l+(r-l+1)*i/k;
for (int i=0;i<k;i++)
build(t[u].mid[i],t[u].mid[i+1]-1,u*k+i,dep+1);
}
inline bool in(const Point &p,const Node &r)
{return r.xl<=p.x&&p.x<=r.xr&&r.yl<=p.y&&p.y<=r.yr;}
int tp,p0,p1,wfk,ret;
int ltg(int lim){
int ret=0;
while(p[tp].v<=lim&&tp<=p1)
ret+=in(p[tp++],wf);
return ret;
}
int ltk()
{
while(tp<=p1){
int c=in(p[tp],wf);
if (wfk>c)wfk-=c;
else return p[tp].v;
tp++;
}return -1;
}
int qry(int l,int r,int u,int dep)
{
if (l>r)return -1;
int k=d[dep];
if (k>r-l+1){
ret=l;
for (int i=l;i<=r;i++){
int stp=tp,pc=ltg(p[i].v),tc=in(p[i],wf)+pc;
if (wfk>tc)wfk-=tc;
else {
if (wfk>pc)return p[i].v;
else {tp=stp;break;}
}
}return -1;
}
for (int i=0;i<k;i++){
int stp=tp,
tc=t[u*k+i].s.qry()+ltg(p[t[u].mid[i+1]-1].v);
if (wfk>tc)wfk-=tc;
else {
tp=stp;
return qry(t[u].mid[i],t[u].mid[i+1]-1,u*k+i,dep+1);
}
}return -1;
}
int n,q,las;
int main()
{
n=read();q=read();
for (int i=1;i<=200;i++)t[i].s.a=a;
for (int i=1,op;i<=q;i++){
op=read();
if (op==1){
Point sav;
sav.x=read()^las;
sav.y=read()^las;
sav.v=1000000000-(read()^las);
bool fl=0;
for (int i=p1;i>p0;i--)
if (sav.v<=p[i].v)p[i+1]=p[i];
else {p[i+1]=sav;fl=1;break;}
if (!fl)p[p0+1]=sav;
p1++;
if ((p1-p0)*(p1-p0)>=p1*120){
sort(p+1,p+(p0=p1)+1,cmpV);
tn=a+1;build(1,p1,1,1);
}
}else {
wf.xl=read()^las;wf.yl=read()^las;
wf.xr=read()^las;wf.yr=read()^las;
wfk=read()^las;
tp=p0+1;
las=qry(1,p0,1,1);
if (las==-1)las=ltk();
if (las==-1){
puts("NAIVE!ORZzyz.");
las=0;continue;
}printf("%d\n",las=(1000000000-las));
}
}return 0;
}
```
# $\rm KDT$ 乱搞
人，要有梦想（？）
这部分的代码实现先咕了。
- **邻近点查询** : [HDU2966 In case of failure](http://acm.hdu.edu.cn/showproblem.php?pid=2966)
**题意** ：给平面图上若干点（无重点），对每个点，找到欧几里德距离最近点。
一个非常直觉的做法 ： 在树上暴力搜索，若已经得到答案 $d$ ，以询问点为圆心， $d$ 为半径画一个圆，与该圆无交的区域不去搜索。
如果把点安排在圆周上并微扰，询问圆心的复杂度可以卡到 $O(n)$ 级别。
但是，该算法在随机数据下表现良好。
在本题中，由于询问点同时也是贡献点，难以构造出 $\rm Hack$ 数据。
- [P4475 巧克力王国](https://www.luogu.com.cn/problem/P4475)
**题意** ： 给出若干个点 $(x,y)$ ，带有点权。
每次给出 $a,b,c$ ,询问满足 $ax+by\leq c$ 的点的点权和。
不难发现等价于半平面数点，正经的做法可见 [$\rm EI$ 的讨论](https://www.luogu.com.cn/discuss/show/186834) 。
数据随机的情况下，当然是 $\rm KDT$ 大法好。
直觉做法 ： 若半平面包含当前矩形，直接返回，否则继续递归。
数据随机的情况下复杂度可以估计为 $O(n^{1/2+eps})$ 或 $O(n^{\log_2\sqrt{3}})≈O(n^{3/4})$。
- [P4631 [APIO2018] Circle selection 选圆圈](https://www.luogu.com.cn/problem/P4631)
记下一个区域内最大的半径，如果区域内不可能有圆和当前圆有交，则不搜索。
这样可以直接过掉 $\rm luogu$ 的数据，$\rm Loj$ 上随机旋转后也可通过。
# $\rm KDT$ 分治
类似线段树分治。
- **移动迷宫**
**题意** : 给出一张无向图，每条边有两个权值 $a,b$，支持如下操作。
- 修改某条边的权值。
- 给出常数 $k$ ,查询从某个点出发，经过权值满足 $k\in[a,b]$ 的边，能到达的点数。
允许离线。
我们把询问按照 $($时间$,k)$ 看作平面上的点。
一条边有权值 $a,b$ ,以及存在区间 $l,r$ ,能够贡献的询问是一个矩形。
对询问建立 $\rm KDT$ ，将边挂上去，然后模仿线段树分治即可。
时间复杂度为 $O(n\sqrt{n}\log n)$ 空间为 $O(n\sqrt{n})$。
空间复杂度可以做到更低，我们不使用提前 $O(\sqrt{n})$ 挂一个操作的方式，而把所有操作一起移动（这才是真正的分治）
每次到达一个点时，将当前节点上滞留的操作 :
- 若完全包含当前区域，则执行该操作。
- 否则，分到下面相应的孩子中。
最后，将该节点上滞留标记占用的空间回收。(需要手写链表)
考虑我们在 $\rm KDT$ 上 $dfs$ 的过程，若处理到点 $u$ ，只有 $u$ 到根路径的侧链共 $O(\log n)$ 个节点可能挂着操作。
这样，一个操作最多同时存在于 $O(\log n)$ 个节点中，空间复杂度也就是 $O(n\log n)$ 了。
- $\rm Sub\ Problem1$ ：[P5443 [APIO2019]桥梁](https://www.luogu.com.cn/problem/P5443)
由于正解是 $O(n\sqrt{n\log n})$ 的，且常数较小，所以 $\rm KDT$ 分治并不能通过，只能卡到 $73'$。
[Loj评测记录](https://loj.ac/submission/968381)
- $\rm Sub\ Problem2$ ：[P3247 [HNOI2016]最小公倍数](https://www.luogu.com.cn/problem/P3247)
由于时限较松，可以通过。[评测记录](https://www.luogu.com.cn/record/40637608)