KDT小记

command_block

2020-10-23 15:17:23

算法·理论

相对另一个传统方法《树套树》而言，有以下优势： - 部分情况下代码实现更简单 - 空间复杂度 $O(n)

可以搭配传统的整体标记系统。

# $\rm KD-tree$ 的构造 $\rm KDT$ 是一棵二叉排序树，每个节点有一个管辖(超)矩形。保证 $u$ 子树中的所有点都在 $u$ 的管辖矩形内。左右儿子的管辖矩形是由父亲的管辖矩形水平划分而得。$\rm KDT$ 的核心思想就是对矩形的水平分割。在满足上述要求的前提下，$\rm KDT$ 还要使得高度严格为 $\log_2n+O(1)$ 。 ![](https://cdn.luogu.com.cn/upload/image_hosting/n204tnzp.png) 这是一个 $\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 ,将树上的节点按照管辖矩形分类：

与 R 无交集。
完全包含在 R 内。
部分包含在 R 内。

观察查询算法的行为，其一路经过 3 类点，且一旦遇到 1,2 类点就停下。

而且，1,2 类点的子树中显然无三类点，这样，查询的复杂度和 3 类点的个数直接相关。

对于矩形的每条边，分别计算其“穿过”的节点个数。

我们以轮换划分来分析。在相邻的两轮中，分别对 x,y 进行了划分。

能够发现，在划分出的四个区域中，每个区域的大小都是原来的 1/4 ，任意的一条直线只能穿过其中两个。

这样，设 T(n) 为在大小为 n 的 \rm KDT 中一条直线所经过的区域数，那么有

所以，一条直线经过的区域个数为 $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

题意：给出一些矩形，每个矩形都有一个高度，保证高度互不相同。

有若干次从下到上的射击，会击中并摧毁射击点上最低的那个矩形。

求出每次射击击中的矩形，或指出未射中。

考虑使用矩形来匹配射击点。

可以把矩形按高度排序后，依次寻找范围内的编号最小的射击点，并将其删除，正确性显然。

先对射击点建立 \rm KDT ，维护子树内未被删除的最小射击点。range query 之后单点修改即可。

剪枝 : 若子树内编号最小点仍然大于当前答案，则不进入。优先进入最小点编号更小的子树。(大幅优化)

由于我个人比较偏好 \rm Leafy 数据结构，下面是“线段树式”的 \rm KDT 实现。

#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$ 上。为了减少细节，可以把 $\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$ 在向下的两层产生的四个分支中，只会选择两个分支。四个分支和原树的大小比例如图 : ![](https://cdn.luogu.com.cn/upload/image_hosting/zwzl3kjj.png) 这样，最坏情况下只能选出大小为 $\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 倍常数的样子。

附 : 替罪羊树复杂度简略分析

当然，有条件的话，根据实际情况调 \alpha 当然是最好的。

下面是 \alpha=3/5 的代码。评测记录

#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}) 的算法。

常数较小，跑得比局部重构要快得多（似乎也更好写）。

```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 ，可能可以获得更好的理论复杂度。

#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

题意：给平面图上若干点（无重点），对每个点，找到欧几里德距离最近点。

一个非常直觉的做法：在树上暴力搜索，若已经得到答案 d ，以询问点为圆心， d 为半径画一个圆，与该圆无交的区域不去搜索。

如果把点安排在圆周上并微扰，询问圆心的复杂度可以卡到 O(n) 级别。

但是，该算法在随机数据下表现良好。

在本题中，由于询问点同时也是贡献点，难以构造出 \rm Hack 数据。

P4475 巧克力王国

题意：给出若干个点 (x,y) ，带有点权。

每次给出 a,b,c ,询问满足 ax+by\leq c 的点的点权和。

不难发现等价于半平面数点，正经的做法可见 \rm EI 的讨论。

数据随机的情况下，当然是 \rm KDT 大法好。

直觉做法：若半平面包含当前矩形，直接返回，否则继续递归。

数据随机的情况下复杂度可以估计为 O(n^{1/2+eps}) 或 O(n^{\log_2\sqrt{3}})≈O(n^{3/4})。

P4631 [APIO2018] Circle selection 选圆圈

记下一个区域内最大的半径，如果区域内不可能有圆和当前圆有交，则不搜索。

这样可以直接过掉 \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) 了。

由于正解是 $O(n\sqrt{n\log n})$ 的，且常数较小，所以 $\rm KDT$ 分治并不能通过，只能卡到 $73'$。 [Loj评测记录](https://loj.ac/submission/968381)
由于时限较松，可以通过。[评测记录](https://www.luogu.com.cn/record/40637608)