离线静态四维数点题解
前言
本文同步自洛谷专栏，题目传送门。
怎么打了 KDT 标签却没有 KDT 题解？这是一篇复杂度 \(O(n\sqrt{n})\) 的 KDT 题解！
前置知识：【模板】KDT，【模板】三维偏序 / 陌上花开。
参考资料：KDT，处理高维数据的利器 by xde。
注：写此题时，本人仅写过 \(2\) 次 KDT，可能有些优化并不了解或实现不优秀，望各路大手子指正。
KDT 实现方式均为二进制分组优化。\(n,m\) 同阶，下文将不作区分。
分析
记原题面中给出的四维分别为 \(x,y,a,b\)，其中 \(a,b\) 维通过取负，也调整为 \(\le\) 型。（个人习惯）
题意很明确，离线静态四维数点，采用 BIT 套 KDT 的方式实现。
对 \(x\) 维排序，对 \(y\) 维离散化，这里离散化注意实现方式，把操作维的要区分开，避免重复。（和一篇莫队题解有些相似）
$\red{\text{离散化实现与一些讨论}}$
这是离散化部分的代码。
#define inf 1000000001
struct ope{//(x,y,a,b)->(x1,y1,x2,y2)
int x,y,a,b,id;
}Q[N],A[N];//查询，修改
n=read(),m=read();
for(int i=1;i<=n;i++){
int x=read(),y=read(),a=read();
A[i]={x,y,inf-a,inf-read(),0},E[i]=y;
}
sort(E+1,E+1+n);//注意离散化手法使之可以严格区分
for(int i=1;i<=n;i++){//避免往单点中大量插入
int p=lower_bound(E+1,E+1+n,A[i].y)-E;
cnt[p]++,A[i].y=p+cnt[p]-1;//cnt 是一个辅助数组
}
for(int i=1;i<=m;i++){
int x=read(),y=read(),a=read();
y=upper_bound(E+1,E+1+n,y)-E-1;
Q[i]={x,y,inf-a,inf-read(),i};
}
事实上，排序并去重的离散化也能够通过，
然而理论上在 BIT 上准确地找出 \(\log{n}\) 个位置各插入 \(\frac{n}{\log{n}}\) 个数，查询复杂度可以达到 \(O(\sqrt{n\log{n}})\)，而非 \(O(\sqrt{n})\)。
如果所有询问的查询点也相同，可以将复杂度卡到 \(O(n\sqrt{n\log{n}})\)。
可能是 lxl 没有卡？个人认为卡此方法难度较大，并且可以通过随机扰动（上述 hack 依赖于 BIT 中的位置）的方式解决，不作进一步展开。
现在要做的就是动态三维偏序，则用 BIT 解决 \(y\) 维小于等于的限制，BIT 中的每一位开一棵 KDT 解决动态的 \(a,b\) 维，做完了。
$\red{\text{外层 BIT 的实现}}$
sort(A+1,A+1+n),sort(Q+1,Q+1+m);//cnt 数组用来存答案
for(int i=1,pl=1,p,tmp;i<=m;cnt[Q[i].id]=tmp,i++){
for(;pl<=n&&A[pl].x<=Q[i].x;pl++){//归并
for(p=A[pl].y;p<=n;p+=p&-p) KDT[p].insert(A[pl].a,A[pl].b);
}
for(tmp=0,p=Q[i].y;p;p-=p&-p) tmp+=KDT[p].query(Q[i].a,Q[i].b);
}
for(int i=1;i<=m;i++) write(cnt[i]);
复杂度分析
空间复杂度是显然的 \(O(n\log{n})\)。
此方法的第一难点在于分析时间复杂度，首先给出结论：复杂度为 \(O(n\sqrt{n})\)。
第一眼看上去会很疑惑：BIT 套 KDT，内层 \(O(\sqrt{n})\)，外层 \(O(\log{n})\)，总计不应该是 \(O(n\sqrt{n}\log{n})\) 吗？
那么我们就要从二进制分组优化出发来寻找答案了。
我们根据 \(T(n)=2T(\frac{n}{4})+O(1)\) 分析得到 KDT 的复杂度是 \(O(\sqrt{n})\)，其中 \(n\) 为点数，下面写为 \(siz\)。
二进制分组优化下，记 \(k=\log{siz}\)，则其总时间复杂度可以写为：
\(\sum_{i=0}^k\sqrt{2^i}=\frac{2^{\frac{k+1}{2}}-1}{\sqrt{2}-1}\le(2+\sqrt{2})2^{\frac{k}{2}}=O(\sqrt{siz})\)。
再应用在 BIT 上，注意 BIT 中，每次询问的区间长度至少减半，又因为离散化方式，每个位置只对应一个点，也就是每次询问的 KDT 内的总点数至少减半，复杂度也可以分析到 \(O(\sqrt{n})\)。
总时间复杂度自然是 \(O(n\sqrt{n})\)。
实现方式和常数优化
此方法的第二个难点就是实现方式与卡常了。
如果同一个 KDT 中的点对应的下标不连续，不仅代码难度增大，而且常数也会更劣，用 vector 自然也有常数的问题。那么可以考虑给这 \(n\) 个 KDT 预分配空间。
根据离散化方式，对于位置 \(i\) 上的 KDT，其最多有 \(\operatorname{lowbit}(i)\) 个结点，基于此给每个 KDT 分配一个下标区间，全局开一个类似 pool 的数组即可。
$\red{\text{一个优化}}$
每次在 KDT 中查询时，查询的是一个紧贴 \(a,b\) 轴的矩形，可以省掉两个参数的判断，降低常数。
当你发现 \(\sqrt{n}=632<6444=\log^3{n}\) 时，不免大喜：碾压 cdq 套 cdq，最优解我来了！
兴致冲冲地写了一份代码交上去后，你会发现：TLE，被卡常了！
一看出题人 lxl 就不感觉奇怪了。
来不严谨地讨论一下此方法的常数，便于直观理解。
最劣情况下，二进制分组优化和 BIT 各自产生一个 \(\sqrt{2}+2\) 的常数，两者相乘大概是 \(11.6\)，\(11.6n\sqrt{n}=3\times 10^9\)，然而 KDT 内部实现为递归查询，
$\red{\text{注}}$
其实上面的分析中，由于一些下取整的原因，每次多出了一个 \(\sqrt{2}\)，但是本人不太能证明到底什么情况下最劣，所以分析的比较粗糙。
或许可以把常数证到 \(\sqrt{2}+1\)，两个一乘就是 \(6\) 左右。但常数分析毕竟玄学，仅供感性理解。
\(3\times 10^9\) 次递归，二十秒时限都不一定能跑过去，怎么办？
众所周知，循环比递归快，\(1\times 10^9\) 次简单的循环可以跑进一秒，那 KDT 的递归能不能写成循环呢？
显然是可以的，只需要把类似深搜的递归改成类似广搜的循环即可，手写队列就能通过。
代码
KDT 是交替排序，尝试过方差优化，然而好像没有什么用？
$\red{\text{code}}$
#include<bits/stdc++.h>
using namespace std;//判了 EOF
#define gc() (rp1==rp2&&(rp2=(rp1=buf)+fread(buf,1,IO,stdin))==rp1?EOF:*rp1++)
#define pc(a) ((wrp==obuf+IO&&(fwrite(obuf,1,IO,stdout),wrp=obuf)),(*wrp++)=a)
#define N 400005
#define inf 1000000001
struct node{
int x,y,ls,rs,sum;
int lx,ly,rx,ry;
}tr[N*20];
struct ope{
int x,y,a,b,id;
inline bool operator <(const ope &z)const{
return x<z.x;
}
}Q[N],A[N];
const int IO=1<<22;
char buf[IO+1],obuf[IO+1],*wrp=obuf,*rp1,*rp2;
int cnt[N],n,m,E[N],qp[N];
inline void cmax(int &a,int b){(a<b)&&(a=b);}
inline void cmin(int &a,int b){(a>b)&&(a=b);}
inline int read(){
int a=0,c=gc();
while(!isdigit(c)) c=gc();
while(isdigit(c)) a=10*a+c-'0',c=gc();
return a;
}
inline void write(int x){
int sk[20],top=0;
do{
sk[++top]=x%10,x/=10;
}while(x);
while(top) pc(sk[top--]+'0');
pc('\n');
}
class KDTree{
private:
int st,now,rt[20],ret;
static inline bool cmpx(const struct node &a,const struct node &b){return a.x<b.x;}
static inline bool cmpy(const struct node &a,const struct node &b){return a.y<b.y;}
inline void update(node &a,node &b){
cmin(a.lx,b.lx),cmax(a.rx,b.rx);
cmin(a.ly,b.ly),cmax(a.ry,b.ry),a.sum+=b.sum;
}
int build(int l,int r,int o){
int p=(l+r+1)>>1,a,b;
if(o) nth_element(tr+l,tr+p,tr+r+1,cmpx);
else nth_element(tr+l,tr+p,tr+r+1,cmpy);
a=tr[p].ls=l<=p-1?build(l,p-1,o^1):0;
b=tr[p].rs=p+1<=r?build(p+1,r,o^1):0,tr[p].sum=1;
tr[p].lx=tr[p].rx=tr[p].x,tr[p].ly=tr[p].ry=tr[p].y;
a&&(update(tr[p],tr[a]),0),b&&(update(tr[p],tr[b]),0);
return p;
}
void qry(int p,int qxr,int qyr){
int head=1,tail=1;qp[tail++]=p;
for(int a,b;head<tail;){
p=qp[head++],a=tr[p].ls,b=tr[p].rs;
if(tr[p].x<=qxr&&tr[p].y<=qyr) ret++;
if(a&&!(tr[a].lx>qxr||tr[a].ly>qyr)){
if(tr[a].rx<=qxr&&tr[a].ry<=qyr) ret+=tr[a].sum;
else qp[tail++]=a;
}
if(b&&!(tr[b].lx>qxr||tr[b].ly>qyr)){
if(tr[b].rx<=qxr&&tr[b].ry<=qyr) ret+=tr[b].sum;
else qp[tail++]=b;
}
}
}
public:
inline void init(int x){st=x,now=0;}
void insert(int x,int y){
now++,tr[st+now].x=x,tr[st+now].y=y;
int p=31^__builtin_clz(now&-now);
rt[p]=build(st+1+now-(1<<p),st+now,0);
}
int query(int x,int y){
ret=0;
for(int tmp=now,p;tmp;tmp-=(1<<p)){
p=31^__builtin_clz(tmp),qry(rt[p],x,y);
}
return ret;
}
}KDT[N];
int main(){
n=read(),m=read(),tr[0].sum=0;
for(int i=1;i<=n;i++){
int x=read(),y=read(),a=read();
A[i]={x,y,inf-a,inf-read(),0},E[i]=y;
}
sort(E+1,E+1+n);//注意离散化手法使之可以严格区分
for(int i=1;i<=n;i++){//避免往单点中大量插入
int p=lower_bound(E+1,E+1+n,A[i].y)-E;
cnt[p]++,A[i].y=p+cnt[p]-1;
}
for(int i=1;i<=m;i++){
int x=read(),y=read(),a=read();
y=upper_bound(E+1,E+1+n,y)-E-1;
Q[i]={x,y,inf-a,inf-read(),i};
}
for(int i=1,s=0;i<=n;i++) KDT[i].init(s),s+=i&-i;
sort(A+1,A+1+n),sort(Q+1,Q+1+m);
for(int i=1,pl=1,p,tmp;i<=m;cnt[Q[i].id]=tmp,i++){
for(;pl<=n&&A[pl].x<=Q[i].x;pl++){
for(p=A[pl].y;p<=n;p+=p&-p) KDT[p].insert(A[pl].a,A[pl].b);
}
for(tmp=0,p=Q[i].y;p;p-=p&-p) tmp+=KDT[p].query(Q[i].a,Q[i].b);
}
for(int i=1;i<=m;i++) write(cnt[i]);
return fwrite(obuf,1,wrp-obuf,stdout),0;
}
KDT：二维数点我是 \(O(n\sqrt{n})\)，三维数点我是 \(O(n\sqrt{n})\)，四维数点我还是 \(O(n\sqrt{n})\)！毒瘤围困万千重，我自岿然不动！
浙公网安备 33010602011771号