Sorry, your browser cannot access this site
This page requires browser support (enable) JavaScript
Learn more >

P4169-CDQ分治 K-D tree（三维偏序）-天使玩偶

题解

发布于：Apr 30, 2022

更新于：Aug 15, 2024

次浏览

P4169-CDQ分治/K-D tree（三维偏序）-天使玩偶

这是一篇两种做法都有的题解

题外话

我写吐了……

本着不看题解的原则，没写（不会）K-D tree，就写了个cdq分治的做法。下面是我的写题步骤：

  1. 想着树状数组维护不了区间最值，于是写了线段树，因为一个**的错误调了几个小时；

  2. cdq只写了两个方向。显然是错的，因为没考虑修改。所以挂了；

  3. 加上另外两个方向，正确性终于ok，兴高采烈地交上去然后TLE；

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #include<cctype>
    #include<algorithm>
    #include<cmath>
    using namespace std;
    char buf[1<<23],*p1=buf,*p2=buf,obuf[1<<23],*O=obuf;
    #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
    inline int read(){
    	int w=0,x=0;char c=getchar();
    	while(!isdigit(c))w|=c=='-',c=getchar();
    	while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
    	return w?-x:x;
    }
    void print(long long x) {
        if(x>9) print(x/10);
        *O++=x%10+'0';
    }
    namespace star
    {
    	const int maxn=6e5+10,INF=0x3f3f3f3f;
    	int n,m,ans[maxn],LN=0,RN=1000006;
    	struct query{
    		int x,y,id,op;
    	}q[maxn<<1];
    	inline bool cmp1(query a,query b){return a.x<b.x;}
    	inline bool cmp2(query a,query b){return a.x>b.x;}
    	struct SegmentTree{
    		#define ls (ro<<1)
    		#define rs (ro<<1|1)
    		struct tree{
    			int l,r,mx;
    		}e[16000005];
    		void build(int ro,int l,int r){
    			e[ro].l=l,e[ro].r=r;
    			if(l==r){
    				e[ro].mx=-INF;return;
    			}
    			int mid=l+r>>1;
    			build(ls,l,mid);
    			build(rs,mid+1,r);
    			e[ro].mx=max(e[ls].mx,e[rs].mx);
    		}
    		void update(int ro,int k,int v){
    			int l=e[ro].l,r=e[ro].r;
    			if(l==r){
    				e[ro].mx=max(e[ro].mx,v);return;
    			}
    			int mid=l+r>>1;
    			if(k<=mid)update(ls,k,v);
    			else update(rs,k,v);
    			e[ro].mx=max(e[ls].mx,e[rs].mx);
    		}
    		void update2(int ro,int k,int v){
    			int l=e[ro].l,r=e[ro].r;
    			if(l==r){
    				e[ro].mx=v;return;
    			}
    			int mid=l+r>>1;
    			if(k<=mid)update2(ls,k,v);
    			else update2(rs,k,v);
    			e[ro].mx=max(e[ls].mx,e[rs].mx);
    		}
    		int query(int ro,int x,int y){
    			int l=e[ro].l,r=e[ro].r;
    			if(l==x and r==y)return e[ro].mx;
    			int mid=l+r>>1;
    			if(y<=mid)return query(ls,x,y);
    			else if(x>mid)return query(rs,x,y);
    			else return max(query(ls,x,mid),query(rs,mid+1,y));
    		}
    		#undef ls
    		#undef rs
    	}T;
    	void cdq(int l,int r){
    		if(l==r)return;
    		int mid=l+r>>1;
    		cdq(l,mid),cdq(mid+1,r);
    		int i,j;
    		sort(q+l,q+mid+1,cmp1),sort(q+mid+1,q+r+1,cmp1);
    		for(i=mid+1,j=l;i<=r;i++){
    			while(q[j].x<=q[i].x and j<=mid){
    				if(!q[j].op) T.update(1,q[j].y,q[j].x+q[j].y);j++;
    			}
    			if(q[i].op)ans[q[i].id]=min(ans[q[i].id],q[i].x+q[i].y-T.query(1,LN,q[i].y));
    		}
    		for(i=l;i<j;i++) T.update2(1,q[i].y,-INF);
    		
    		for(i=mid+1,j=l;i<=r;i++){
    			while(q[j].x<=q[i].x and j<=mid){
    				if(!q[j].op) T.update(1,q[j].y,q[j].x-q[j].y);j++;
    			}
    			if(q[i].op)ans[q[i].id]=min(ans[q[i].id],q[i].x-q[i].y-T.query(1,q[i].y,RN));
    		}
    		for(i=l;i<j;i++) T.update2(1,q[i].y,-INF);
    		
    		sort(q+l,q+mid+1,cmp2),sort(q+mid+1,q+r+1,cmp2);
    		for(i=mid+1,j=l;i<=r;i++){
    			while(q[j].x>=q[i].x and j<=mid){
    				if(!q[j].op) T.update(1,q[j].y,q[j].y-q[j].x);j++;
    			}
    			if(q[i].op)ans[q[i].id]=min(ans[q[i].id],q[i].y-q[i].x-T.query(1,LN,q[i].y));
    		}
    		for(i=l;i<j;i++) T.update2(1,q[i].y,-INF);
    		
    		for(i=mid+1,j=l;i<=r;i++){
    			while(q[j].x>=q[i].x and j<=mid){
    				if(!q[j].op) T.update(1,q[j].y,-q[j].y-q[j].x);j++;
    			}
    			if(q[i].op)ans[q[i].id]=min(ans[q[i].id],-q[i].y-q[i].x-T.query(1,q[i].y,RN));
    		}
    		for(i=l;i<j;i++) T.update2(1,q[i].y,-INF);
    	}
    	inline void work(){
    		n=read(),m=read();
    		memset(ans,INF,sizeof ans);
    		for(int i=1;i<=n;i++) q[i].x=read(),q[i].y=read(),LN=min(LN,q[i].y),RN=max(RN,q[i].y),q[i].op=0;
    		for(int i=n+1;i<=n+m;i++){
    			q[i].op=(read()-1);
    			q[i].id=q[i-1].id+q[i].op;
    			q[i].x=read(),q[i].y=read();
    		}
    		int mx=q[n+m].id;
    		T.build(1,LN,RN);
    		cdq(1,n+m);
    		for(int i=1;i<=mx;i++)print(ans[i]),*O++='\n';
    		fwrite(obuf,O-obuf,1,stdout);
    	}
    }
    signed main(){
    	star::work();
    	return 0;
    }

  4. 原以为是线段树常数过大，学习了树状数组解法然后WA了；

  5. 改回了线段树，把sort换成了merge，又WA；

  6. 只有最后一个办法了：改成树状数组+merge，要是再挂我就当场去学KD-tree.

  7. 写挂了，我滚去学K-D tree了。

K-D tree

K-D tree是一种维护多维空间点的数据结构，在这道题上是两维所以也叫作2-D tree。具体实现方式是将每个节点所在的区域切割递归建成二叉树，每一个树上的节点代表一个实际的结点，但又存储着选择这个结点时的区域大小。

建树

KDT在建树时，为了让它保持平衡，我们需要尽量选择所在空间维度的中位数，所以在保证时间复杂度正确的情况下我们可以使用nth_element函数，其作用是在线性时间内将一个数摆到它排序后应该在的位置，而将比它小的放在左边，比他大的放在右边（不重载运算符情况下是从小到大），但不是有序排列。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
int build(int l,int r){
	if(l>r)return 0;
	int mid=l+r>>1;
	double av[2],va[2];
	av[0]=av[1]=va[0]=va[1]=0;
	for(int i=l;i<=r;i++)av[0]+=e[g[i]].x,av[1]+=e[g[i]].y;
	av[0]/=(r-l+1),av[1]/=(r-l+1);
	for(int i=l;i<=r;i++)
		va[0]+=(av[0]-e[g[i]].x)*(av[0]-e[g[i]].x),
		va[1]+=(av[1]-e[g[i]].y)*(av[1]-e[g[i]].y);//寻找应该切割的维度
	if(va[0]>va[1])nth_element(g+l,g+mid,g+r+1,cmpx),e[g[mid]].d=1;//d是维度
	else nth_element(g+l,g+mid,g+r+1,cmpy),e[g[mid]].d=2;
	e[g[mid]].ls=build(l,mid-1);
	e[g[mid]].rs=build(mid+1,r);
	maintain(g[mid]);//pushup函数，更新节点信息。
	return g[mid];
}

插入

从root递归下去查找，按照当前节点所在维度向下插入，到达空节点时新建节点存储信息。

注意：当插入节点过多时KDT有可能失衡，此时我们需要将它拍扁重建（pia~）（因为KDT的结构，好像没有别的方法了？）\(^ ①\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
void insert(int &x,int k){
	if(!x){
		x=k,maintain(k);return;
	}
	if(e[x].d==1){
		if(e[k].x<=e[x].x)insert(e[x].ls,k);
		else insert(e[x].rs,k);
	}else{
		if(e[k].y<=e[x].y)insert(e[x].ls,k);
		else insert(e[x].rs,k);
	}
	maintain(x);
	if(bad(x)) rebuild(x);//pia
}
1
2
inline bool bad(int x){return 0.9*e[x].siz<=(double)max(e[e[x].ls].siz,e[e[x].rs].siz);}
//0.9为拍扁的阈值，越大拍扁越不频繁，但有可能失衡，按照实际情况调整。注意，实测其为0.8时会爆栈。
1
2
3
4
5
6
7
8
9
10
11
void getson(int x){
	if(!x)return;
	getson(e[x].ls);
	g[++t]=x;
	getson(e[x].rs);
}
inline void rebuild(int &x){
	t=0;
	getson(x);//找到所有被拍扁的节点（其实没必要，新建节点也行）
	x=build(1,t);
}

查询

此题要求查询距离关键点最近的点的距离。

注意：我们在查询下传的时候需要比较的是区块位置距离关键点的距离，而统计答案是按照当前节点的坐标统计。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int ans=INF;//请不要学鄙人用全局变量传参，我被人嘴了
void query(int x){
	cmin(ans,dist(x));
	int dl=INF,dr=INF;
	if(e[x].ls)dl=getdis(e[x].ls);
	if(e[x].rs)dr=getdis(e[x].rs);
	if(dl<dr){
		if(dl<ans)query(e[x].ls);
		if(dr<ans)query(e[x].rs);
	}else{
		if(dr<ans)query(e[x].rs);
		if(dl<ans)query(e[x].ls);
	}
}

完整代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<cmath>
#define cmin(i,j) (i)=min((i),(j))
#define cmax(i,j) (i)=max((i),(j))
using namespace std;
inline int read(){
	int w=0,x=0;char c=getchar();
	while(!isdigit(c))w|=c=='-',c=getchar();
	while(isdigit(c))x=(x<<3)+(x<<1)+(c^48),c=getchar();
	return w?-x:x;
}
namespace star
{
	const int maxn=6e5+10,INF=0x3f3f3f3f;
	int n,m,rt,tot,qx,qy,g[maxn],t;
	struct node{
		int x,y,ls,rs,l[2],r[2],siz,d;
	}e[maxn];
	inline bool cmpx(int a,int b){return e[a].x<e[b].x;}
	inline bool cmpy(int a,int b){return e[a].y<e[b].y;}
	inline void maintain(int x){
		e[x].siz=e[e[x].ls].siz+e[e[x].rs].siz+1;
		e[x].l[0]=e[x].r[0]=e[x].x;
		e[x].l[1]=e[x].r[1]=e[x].y;
		if(e[x].ls)
			cmin(e[x].l[0],e[e[x].ls].l[0]),cmax(e[x].r[0],e[e[x].ls].r[0]),
			cmin(e[x].l[1],e[e[x].ls].l[1]),cmax(e[x].r[1],e[e[x].ls].r[1]);
		if(e[x].rs)
			cmin(e[x].l[0],e[e[x].rs].l[0]),cmax(e[x].r[0],e[e[x].rs].r[0]),
			cmin(e[x].l[1],e[e[x].rs].l[1]),cmax(e[x].r[1],e[e[x].rs].r[1]);
	}
	int build(int l,int r){
		if(l>r)return 0;
		int mid=l+r>>1;
		double av[2],va[2];
		av[0]=av[1]=va[0]=va[1]=0;
		for(int i=l;i<=r;i++)av[0]+=e[g[i]].x,av[1]+=e[g[i]].y;
		av[0]/=(r-l+1),av[1]/=(r-l+1);
		for(int i=l;i<=r;i++)
			va[0]+=(av[0]-e[g[i]].x)*(av[0]-e[g[i]].x),
			va[1]+=(av[1]-e[g[i]].y)*(av[1]-e[g[i]].y);
		if(va[0]>va[1])nth_element(g+l,g+mid,g+r+1,cmpx),e[g[mid]].d=1;
		else nth_element(g+l,g+mid,g+r+1,cmpy),e[g[mid]].d=2;
		e[g[mid]].ls=build(l,mid-1);
		e[g[mid]].rs=build(mid+1,r);
		maintain(g[mid]);
		return g[mid];
	}
	void getson(int x){
		if(!x)return;
		getson(e[x].ls);
		g[++t]=x;
		getson(e[x].rs);
	}
	inline void rebuild(int &x){
		t=0;
		getson(x);
		x=build(1,t);
	}
	inline bool bad(int x){return 0.9*e[x].siz<=(double)max(e[e[x].ls].siz,e[e[x].rs].siz);}
	void insert(int &x,int k){
		if(!x){
			x=k,maintain(k);return;
		}
		if(e[x].d==1){
			if(e[k].x<=e[x].x)insert(e[x].ls,k);
			else insert(e[x].rs,k);
		}else{
			if(e[k].y<=e[x].y)insert(e[x].ls,k);
			else insert(e[x].rs,k);
		}
		maintain(x);
		if(bad(x)) rebuild(x);
	}
	inline int getdis(int x){return max(0,qx-e[x].r[0])+max(0,e[x].l[0]-qx)+max(0,qy-e[x].r[1])+max(0,e[x].l[1]-qy);}
	inline int dist(int x){return abs(qx-e[x].x)+abs(qy-e[x].y);}
	int ans;
	void query(int x){
		cmin(ans,dist(x));
		int dl=INF,dr=INF;
		if(e[x].ls)dl=getdis(e[x].ls);
		if(e[x].rs)dr=getdis(e[x].rs);
		if(dl<dr){
			if(dl<ans)query(e[x].ls);
			if(dr<ans)query(e[x].rs);
		}else{
			if(dr<ans)query(e[x].rs);
			if(dl<ans)query(e[x].ls);
		}
	}
	inline void work(){
		n=read(),m=read();
		while(n--){
			e[++tot].x=read(),e[tot].y=read();
			g[tot]=tot;
		}
		rt=build(1,tot);
		while(m--){
			if(read()==1){
				e[++tot].x=read(),e[tot].y=read();
				insert(rt,tot);
			}else{
				qx=read(),qy=read();ans=INF;query(rt);
				printf("%d\n",ans);
			}
		}
	}
}
signed main(){
	star::work();
	return 0;
}

注释

① K-D tree的重建方式含争议。事实上，作者是从OIwiki上学习的K-D tree写法，而在当页下方评论区中有人提出KDT不应该像替罪羊一样按照阈值重构，因为这样的话会在极端情况下多次重构导致时间复杂度退化。他提出，应该每插入\(\sqrt n\)个点后将全树进行重构，这样保证时间复杂度最劣情况下全部插入为\(O(n\sqrt n)\)，查询最劣情况下单次为\(O(\sqrt n)\)。但是大家不用担心，愚以为大多数情况下本人这种写法的均摊时间复杂度是比另一种优的，实在不行可以适当调整阈值。欢迎大家激烈对线。

更新于：Aug 15, 2024

博客园

OI

分治

给小狼留言