「APIO2018」选圆圈(K-D Tree/CDQ+Set)

「APIO2018」选圆圈(K-D Tree/CDQ+Set)

Part1 K-D Tree做法

K-D Tree经常用来优化大暴力。。

把圆((x,y,r))视为矩形((x-r,y-r,x+r,y+r))，依据((x,y))构建K-D Tree

维护K-D Tree每个节点所有矩形最小和最大的(x,y)，通过判断当前圆与其是否有交来剪枝

删去的节点(x,y)不算进矩形范围即可

很显然这是一个最坏(O(n^2))的算法，直接(x,y)轮换建树这样写APIO的数据已经卡过了。。

比较好的办法是按照(x,y)的方差大的一维建树，当然旋转角度也是可以的

实际运行速度堪比(O(nlog n))

Code1:旋转+(x,y)轮换建树

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double db;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) f=IO=='-';
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=3e5+10,INF=1e9+10;
const db eps=1e-7;
const db co=cos(1123),si=sin(1123);

int n,typ;
struct Node{
	db x,y; int r,id;
}A[N],B[N];
db Dis(db x,db y,Node z){ return (x-z.x)*(x-z.x)+(y-z.y)*(y-z.y); }
db Dis(Node x,Node y){ return Dis(x.x,x.y,y); }
db In(db x,db y,Node z){ return Dis(x,y,z)-eps<=(db)z.r*z.r; }
int cmp(Node x,Node y){ return typ?x.x<y.x:x.y<y.y; }
int c[N],ch[N][2],rt;
db lx[N],rx[N],ly[N],ry[N];
void Up(int u) {
	if(c[B[u].id]) lx[u]=ly[u]=1e18,rx[u]=ry[u]=-1e18;
	else lx[u]=B[u].x-B[u].r,rx[u]=B[u].x+B[u].r,ly[u]=B[u].y-B[u].r,ry[u]=B[u].y+B[u].r;
	for(int v:ch[u]) if(v) cmin(lx[u],lx[v]),cmax(rx[u],rx[v]),cmin(ly[u],ly[v]),cmax(ry[u],ry[v]);
}
int Build(int l,int r) {
	if(l>r) return 0;
	int u=(l+r)>>1;
	nth_element(B+l,B+u,B+r+1,cmp);
	typ^=1; ch[u][0]=Build(l,u-1),ch[u][1]=Build(u+1,r); typ^=1;
	return Up(u),u;
}
int Cross(int x,Node y){
	if(lx[x]>rx[x]) return 0;
	if(In(lx[x],ly[x],y) || In(lx[x],ry[x],y) || In(rx[x],ly[x],y) || In(rx[x],ry[x],y)) return 1;
	if(lx[x]-eps<=y.x && y.x<=rx[x]+eps && ly[x]-eps<=y.y && y.y<=ry[x]+eps) return 1;
	if(lx[x]-eps<=y.x && y.x<=rx[x]+eps) if(In(y.x,ly[x],y) || In(y.x,ry[x],y)) return 1;
	if(ly[x]-eps<=y.y && y.y<=ry[x]+eps) if(In(lx[x],y.y,y) || In(rx[x],y.y,y)) return 1;
	return 0;
}
void Del(int u,Node x){
	if(!u || !Cross(u,x)) return;
	if(!c[B[u].id] && Dis(x,B[u])-eps<=(db)(x.r+B[u].r)*(x.r+B[u].r)) c[B[u].id]=x.id;
	Del(ch[u][0],x),Del(ch[u][1],x);
	Up(u);
}

int main(){
	n=rd();
	rep(i,1,n) {
		db x=rd(),y=rd();
		A[i]=B[i]=(Node){x*co-y*si,x*si+y*co,rd(),i};
	}
	sort(A+1,A+n+1,[&](Node x,Node y){ return x.r!=y.r?x.r>y.r:x.id<y.id;}),rt=Build(1,n);
	rep(i,1,n) if(!c[A[i].id]) Del(rt,A[i]);
	rep(i,1,n) printf("%d ",c[i]);
}

[ ]

Code2:方差建树

#include<bits/stdc++.h>
using namespace std;
using ll=long long;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
char buf[200000],*p1,*p2;
#define getchar() (((p1==p2)&&(p2=(p1=buf)+fread(buf,1,200000,stdin))),*p1++)
char IO;
int rd(){
	int s=0,f=0;
	while(!isdigit(IO=getchar())) f=IO=='-';
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}
void wt(int x){
	static int buf[10],l=0;
	while(x) buf[++l]=x%10+'0',x/=10;
	drep(i,l,1) putchar(buf[i]);
	l=0;
}

const int N=3e5+10,INF=1e9+10;

int n,typ;
struct Node{
	ll x,y,r;
	int id;
}A[N],B[N];
ll Dis(ll x,ll y,Node z){ return (x-z.x)*(x-z.x)+(y-z.y)*(y-z.y); }
ll Dis(Node x,Node y){ return Dis(x.x,x.y,y); }
int In(ll x,ll y,Node z){ return Dis(x,y,z)<=z.r*z.r; }
int cmp(Node x,Node y){ return typ?x.x<y.x:x.y<y.y; }
int c[N],ch[N][2],rt;
int lx[N],rx[N],ly[N],ry[N];
void Up(int u) {
	if(c[B[u].id]) lx[u]=ly[u]=INF,rx[u]=ry[u]=-INF;
	else lx[u]=B[u].x-B[u].r,rx[u]=B[u].x+B[u].r,ly[u]=B[u].y-B[u].r,ry[u]=B[u].y+B[u].r;
	for(int v:ch[u]) if(v) cmin(lx[u],lx[v]),cmax(rx[u],rx[v]),cmin(ly[u],ly[v]),cmax(ry[u],ry[v]);
}

int Get(int l,int r){
	long double _x=0,_y=0,x=0,y=0;
	rep(i,l,r) _x+=B[i].x,_y+=B[i].y;
	_x/=r-l+1,_y/=r-l+1;
	rep(i,l,r) x+=(B[i].x-_x)*(B[i].x-_x),y+=(B[i].y-_y)*(B[i].y-_y);
	return x>y;
}

int Build(int l,int r) {
	if(l>r) return 0;
	int u=(l+r)>>1;
	typ=Get(l,r),nth_element(B+l,B+u,B+r+1,cmp);
	ch[u][0]=Build(l,u-1),ch[u][1]=Build(u+1,r);
	return Up(u),u;
}
int Cross(int x,Node y){
	if(lx[x]>rx[x]) return 0;
	if(In(lx[x],ly[x],y) || In(lx[x],ry[x],y) || In(rx[x],ly[x],y) || In(rx[x],ry[x],y)) return 1;
	if(lx[x]<=y.x && y.x<=rx[x] && ly[x]<=y.y && y.y<=ry[x]) return 1;
	if(lx[x]<=y.x && y.x<=rx[x]) if(In(y.x,ly[x],y) || In(y.x,ry[x],y)) return 1;
	if(ly[x]<=y.y && y.y<=ry[x]) if(In(lx[x],y.y,y) || In(rx[x],y.y,y)) return 1;
	return 0;
}
void Del(int u,Node x){
	if(!u || !Cross(u,x)) return;
	if(!c[B[u].id] && Dis(x,B[u])<=(x.r+B[u].r)*(x.r+B[u].r)) c[B[u].id]=x.id;
	Del(ch[u][0],x),Del(ch[u][1],x);
	Up(u);
}

int main(){
	n=rd();
	rep(i,1,n) {
		int x=rd(),y=rd();
		A[i]=B[i]=(Node){x,y,rd(),i};
	}
	sort(A+1,A+n+1,[&](Node x,Node y){ return x.r!=y.r?x.r>y.r:x.id<y.id;}),rt=Build(1,n);
	rep(i,1,n) if(!c[A[i].id]) Del(rt,A[i]);
	rep(i,1,n) printf("%d ",c[i]);
}

[ ]

Part2 CDQ+Set

这是一个稳定(O(nlog ^2n))的算法

按照(r)递减，(id)递增的顺序对于圆排序后，(CDQ)考虑([l,mid])对([mid+1,r])的贡献

先处理([l,mid])的部分，就能知道哪些圆可以对([mid+1,r])产生贡献

处理贡献时，依然把圆视为矩形，按照(x)插入、删除和查询矩形的左右边界((x-r,y),(x+r,y))

插入、删除和查询均是在(set)中维护(y)的前驱后继

同时还需要交换(x,y)重新进行一遍

正确性：

与每个圆交的圆一定在(x)或(y)上与它相邻

如果这个圆在x,y上都不与它相邻还与它相交，则必然会跨过一个相邻的圆，这个圆不会被加入set

故不存在这种情况

实际运行常数很大，被K-D Tree吊起来打

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(i=a;i<=b;++i)
using P=pair<int,int>;
#define M make_pair
#define X first
#define Y second
#define S(x) 1ll*(x)*(x)
const int N=1e6+10;
int n,c[N],i,j,D[N];
struct C{ int x,y,r,i; } A[N]; 
void Upd(int i,int j) { if(S(A[i].x-A[j].x)+S(A[i].y-A[j].y)<=S(A[i].r+A[j].r)&&D[c[A[j].i]]>i) c[A[j].i]=A[i].i; }
set<P>st;
P I[N],E[N],Q[N];
void Work(int l,int r){
	int mid=(l+r)>>1,n=0,m=0,x=0,y=0,t;
	rep(i,l,mid) if(c[A[i].i]==A[i].i) I[m]=M(A[i].x-A[i].r,i),E[m++]=M(A[i].x+A[i].r+1,i);
	rep(i,mid+1,r) Q[n++]=M(A[i].x-A[i].r,i),Q[n++]=M(A[i].x+A[i].r,i);
	sort(I,I+m),sort(E,E+m),sort(Q,Q+n),st.clear();
	rep(i,0,n-1) {
		while(x<m&&I[x].X<=Q[i].X) st.insert(M(A[t=I[x++].Y].y,t));
		while(y<m&&E[y].X<=Q[i].X) st.erase(M(A[t=E[y++].Y].y,t));
		auto j=st.lower_bound(M(A[t=Q[i].Y].y,t));
		if(j!=st.end()) Upd(j->Y,t);
		if(j!=st.begin()) Upd((--j)->Y,t);
	}
}
void Solve(int l,int r) {
	if(r-l+1<=80) {
		rep(i,l,r)if(c[A[i].i]==A[i].i)rep(j,i+1,r) Upd(i,j);
		return;
	}
	int mid=(l+r)>>1;
	Solve(l,mid);
	Work(l,r);rep(i,l,r) swap(A[i].x,A[i].y);
	Work(l,r);rep(i,l,r) swap(A[i].x,A[i].y);
	Solve(mid+1,r);
}
int main(){
	scanf("%d",&n);
	rep(i,1,n) scanf("%d%d%d",&A[i].x,&A[i].y,&A[i].r),A[i].i=i;
	sort(A+1,A+n+1,[&](C x,C y){return M(-x.r,x.i)<M(-y.r,y.i);});
	rep(i,1,n) D[A[i].i]=c[i]=i;
	Solve(1,n);
	rep(i,1,n) printf("%d ",c[i]);
}