KDT小记

command_block · 2020-10-23 15:17:23 · 算法·理论

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

在不利用差分的前提下，支持将一个矩形区域内的点，划分成 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 中一条直线所经过的区域数，那么有

所以，超平面经过的区域个数为 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 占一次，则有 ：

• 例题 ： 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 配合标记

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

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

此时查询复杂度 \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)