「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(n\log 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(n\log ^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]);}