「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]);
}