下面先给出比较简单的KD树的做法——
根据圆心建一棵KD树,然后模拟题目的过程,考虑搜索一个圆
剪枝:如果当前圆[与包含该子树内所有圆的最小矩形]都不相交就退出
然而这样的理论复杂度是$o(n^2)$,所以会被出题人卡了
但是如果将坐标系旋转45度,即对于$(x,y)$,变为$((x-y)/\sqrt{2},(x+y)/\sqrt{2})$,就不会被卡了
然后因为可能相切,精度要求高,eps大约要取$10^{-3}$
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 3000005 4 #define k (l+r>>1) 5 #define s (sqrt(2)) 6 #define eps (1e-3) 7 #define sqr(x) (1LL*(x)*(x)) 8 int n,t,id[N],ans[N],ls[N],rs[N]; 9 struct P{ 10 double x,y; 11 }mn[N],mx[N]; 12 struct ji{ 13 int r,id; 14 P a; 15 bool operator < (const ji &p){ 16 if (t)return a.x<p.a.x; 17 return a.y<p.a.y; 18 } 19 }a[N]; 20 bool cmp(int x,int y){ 21 return (a[x].r>a[y].r)||(a[x].r==a[y].r)&&(a[x].id<a[y].id); 22 } 23 P rotate(P x){ 24 return P{(x.x-x.y)/s,(x.x+x.y)/s}; 25 } 26 double dis(P x,P y){ 27 return sqr(x.x-y.x)+sqr(x.y-y.y); 28 } 29 double calc(int t,P x){ 30 P y; 31 if (x.x<mn[t].x)y.x=mn[t].x; 32 else 33 if (x.x<=mx[t].x)y.x=x.x; 34 else y.x=mx[t].x; 35 if (x.y<mn[t].y)y.y=mn[t].y; 36 else 37 if (x.y<=mx[t].y)y.y=x.y; 38 else y.y=mx[t].y; 39 return dis(x,y); 40 } 41 void up(int l,int r){ 42 mn[k].x=min(mn[ls[k]].x,mn[rs[k]].x); 43 mn[k].y=min(mn[ls[k]].y,mn[rs[k]].y); 44 mx[k].x=max(mx[ls[k]].x,mx[rs[k]].x); 45 mx[k].y=max(mx[ls[k]].y,mx[rs[k]].y); 46 if (ans[a[k].id])return; 47 mn[k].x=min(mn[k].x,a[k].a.x-a[k].r); 48 mn[k].y=min(mn[k].y,a[k].a.y-a[k].r); 49 mx[k].x=max(mx[k].x,a[k].a.x+a[k].r); 50 mx[k].y=max(mx[k].y,a[k].a.y+a[k].r); 51 } 52 int build(int l,int r,int p){ 53 if (l>r)return 0; 54 t=p; 55 nth_element(a+l,a+k,a+r+1); 56 ls[k]=build(l,k-1,p^1); 57 rs[k]=build(k+1,r,p^1); 58 up(l,r); 59 return k; 60 } 61 void query(int l,int r,int x){ 62 if (l>r)return; 63 if (calc(k,a[x].a)>sqr(a[x].r)+eps)return; 64 if ((!ans[a[k].id])&&(dis(a[x].a,a[k].a)<=sqr(a[x].r+a[k].r)+eps))ans[a[k].id]=a[x].id; 65 query(l,k-1,x); 66 query(k+1,r,x); 67 up(l,r); 68 } 69 int main(){ 70 scanf("%d",&n); 71 for(int i=1;i<=n;i++){ 72 scanf("%lf%lf%d",&a[i].a.x,&a[i].a.y,&a[i].r); 73 a[i].a=rotate(a[i].a); 74 a[i].id=id[i]=i; 75 } 76 mn[0].x=mn[0].y=2e9+1; 77 mx[0].x=mx[0].y=-2e9-1; 78 build(1,n,0); 79 sort(id+1,id+n+1,cmp); 80 for(int i=1;i<=n;i++) 81 if (!ans[a[id[i]].id])query(1,n,id[i]); 82 for(int i=1;i<=n;i++)printf("%d ",ans[i]); 83 }
当然,上面的做法并不是正解,下面给出更巧妙的做法——
对于圆$c_{i}$,设圆心为$(x_{i},y_{i})$,半径为$r_{i}$,显然$c_{i}$和$c_{j}$相交当且仅当$(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}\le (r_{i}+r_{j})^{2}$
考虑Subtask4,即所有圆半径都相同的情况(假设半径都为$r$)
建立新坐标系,原坐标系中$(x,y)$对应新坐标系中$(\lfloor\frac{x}{r}\rfloor,\lfloor\frac{y}{r}\rfloor)$
记$(x'_{i},y'_{i})$为$(x_{i},y_{i})$在新坐标系中对应的点,即$(\lfloor\frac{x_{i}}{r}\rfloor,\lfloor\frac{y_{i}}{r}\rfloor)$
此时,考虑两个圆$c_{i}$和$c_{j}$的相交与$(x'_{i},y'_{i})$的关系:
1.若$|x'_{i}-x'_{j}|>2$或$|y'_{i}-y'_{j}|>2$,则该维在原坐标系中相差大于$2r$,显然这两个圆不交
2.若$(x'_{i},y'_{i})=(x'_{j},y'_{j})$,则两维坐标在原坐标系中相差都小于等于$r$,显然这两个圆相交
当搜索$c_{i}$时,根据第1个性质,其仅需要遍历所有$|x'_{i}-x'_{j}|,|y'_{i}-y'_{j}|\le 2$的圆$c_{j}$,不妨先枚举$(x'_{j},y'_{j})$,再枚举所有对应该位置的$j$,这个将所有点按照$(x'_{i},y'_{i})$这个二元组排序即可做到
(删除可以用类似指针的方式维护,但实际上重复访问已经被删除的圆不影响复杂度)
分析此时的复杂度,不妨考虑每一个$j$会被搜索到的次数,注意到搜索$c_{i}$后,根据第2个性质,$(x'_{i},y'_{i})$上所有点都会被选,那么之后就不会从$(x'_{i},y'_{i})$搜到$j$了,因此总复杂度为$o(n)$
另外由于前面的排序,总复杂度为$o(n\log n)$,可以通过
考虑原题,实际上若$r\ge \max_{i=1}^{n}r_{i}$,上面的第1个性质仍然成立,但第2个性质并不一定成立
但是,有一个类似的性质:若$(x'_{i},y'_{i})=(x'_{j},y'_{j})$且$r_{i},r_{j}>\frac{r}{2}$,则这两个圆相交
初始令$r=\max_{i=1}^{n}r_{i}$,若当前$\max_{i=1}^{n}r_{i}\le \frac{r}{2}$(不包括删去的圆),则将$r$变为$\frac{r}{2}$并重新计算$(x'_{i},y'_{i})$
来分析此时的复杂度:对于每一个$r$,同样可以证明每一个$j$不会重复从同一个$(x'_{i},y'_{i})$访问$r$(若两次访问,由于是同一个$r$,这两个圆半径都大于$\frac{r}{2}$,根据上面的性质应该被删除了),复杂度仍为$o(n\log n)$
现在,复杂度的瓶颈已经是排序了,暴力sort复杂度为$o(n\log^{2}n)$,我们将$r$的初值改为$2^{30}$,那么每一次都可以在上一次基础上拆分,即单次排序变为$o(n)$,总排序复杂度为$o(n\log n)$
最终,总复杂度为$o(n\log n)$,但常数很大,反正被卡了QAQ
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 300005 4 #define ll long long 5 int n,R,ans[N]; 6 struct Circle{ 7 int x,y,r,id; 8 int xx()const{ 9 if (x>=0)return x/R; 10 return (x-R+1)/R; 11 } 12 int yy()const{ 13 if (y>=0)return y/R; 14 return (y-R+1)/R; 15 } 16 bool operator < (const Circle &k)const{ 17 return (xx()<k.xx())||(xx()==k.xx())&&(yy()<k.yy()); 18 } 19 }a[N],b[N]; 20 vector<Circle>v[2]; 21 ll sqr(int k){ 22 return 1LL*k*k; 23 } 24 bool cmp1(Circle x,Circle y){ 25 return (x.xx()<y.xx())||(x.xx()==y.xx())&&(x.y<y.y); 26 } 27 bool cmp2(Circle x,Circle y){ 28 return (x.r>y.r)||(x.r==y.r)&&(x.id<y.id); 29 } 30 bool check(Circle x,Circle y){ 31 return sqr(x.x-y.x)+sqr(x.y-y.y)<=sqr(x.r+y.r); 32 } 33 void find(int x,int y,Circle k){ 34 Circle o; 35 o.x=min(abs(x),(1<<30)/R)*R; 36 if (x<0)o.x=-o.x; 37 o.y=min(abs(y),(1<<30)/R)*R; 38 if (y<0)o.y=-o.y; 39 int l=lower_bound(a+1,a+n+1,o)-a; 40 int r=upper_bound(a+1,a+n+1,o)-a-1; 41 for(int i=l;i<=r;i++) 42 if ((!ans[a[i].id])&&(check(k,a[i])))ans[a[i].id]=k.id; 43 } 44 int main(){ 45 scanf("%d",&n); 46 for(int i=1;i<=n;i++){ 47 scanf("%d%d%d",&a[i].x,&a[i].y,&a[i].r); 48 a[i].id=i; 49 } 50 R=(1<<30); 51 sort(a+1,a+n+1,cmp1); 52 memcpy(b,a,sizeof(b)); 53 sort(b+1,b+n+1,cmp2); 54 for(int i=1;i<=n;i++) 55 if (!ans[b[i].id]){ 56 while (b[i].r<=R/2){ 57 for(int j=1,lst=1;j<=n;j++) 58 if ((j==n)||(a[j].xx()<a[j+1].xx())){ 59 int x=a[lst].xx(); 60 v[0].clear(),v[1].clear(); 61 R/=2; 62 for(int k=lst;k<=j;k++)v[a[k].xx()-x*2].push_back(a[k]); 63 for(int k=0;k<v[0].size();k++)a[lst++]=v[0][k]; 64 for(int k=0;k<v[1].size();k++)a[lst++]=v[1][k]; 65 R*=2; 66 } 67 R/=2; 68 } 69 for(int j=-2;j<=2;j++) 70 for(int k=-2;k<=2;k++)find(b[i].xx()+j,b[i].yy()+k,b[i]); 71 } 72 for(int i=1;i<=n;i++)printf("%d ",ans[i]); 73 }