首先考虑找到距离
≤
r
\le r
≤r的交点数量。发现这等价于圆上两段圆弧相交,因此将圆上的点离散化后排序,用一个主席树来求就做完了。
然后是距离求和。这看起来非常棘手。事实上,只要把所有交点都找出来就做完了。首先可以放心的将圆环从一个位置断开。其次,考虑以某种顺序将所有直线依次删掉。发现当按照长度从大到小删时,假设这条线段是
[
l
i
,
r
i
]
[l_i,r_i]
[li,ri],发现只要满足
l
j
∈
[
l
i
,
r
i
]
l_j\in [l_i,r_i]
lj∈[li,ri]或者
r
j
∈
[
l
i
,
r
i
]
r_j\in [l_i,r_i]
rj∈[li,ri],那么直线
i
,
j
i,j
i,j就一定相交。那么前面一个问题也得到了解决,只需用树状树组维护就做完了。
当然,事实上我们可以
O
(
1
)
O(1)
O(1)求出一个交点。考虑倒着做,每次删除一条线段,用链表维护就做完了。事实上交点在圆上的情况并不影响答案,所以可以少一些细节。
那么问题来了,为啥我被卡常了。
#include<bits/stdc++.h>#define ll long long#define fi first#define se second#define pb push_back#define inf 0x3f3f3f3f#define db double#define cpx complex<db>usingnamespace std;
constint N=1e5+5;
structseg{
db k,b;
}seg[N];
structpoint{
db x,y;
}p;
int n,m,cntseg,cnt;
int bit[N],sa[N],pos[N];
int le[N],ri[N],L[N],R[N];//链表db res;
pair<db,int>A[N];
boolcmp(int x,int y){
return R[x]-L[x]<R[y]-L[y];
}
ll ask(int x){
ll tot=0;
for(;x;x-=x&-x)tot+=bit[x];
return tot;
}
voidadd(int x,int y){
for(;x<=cnt;x+=x&-x)bit[x]+=y;
}
db getdist(point x,point y){
returnsqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y));
}
point calc(int x,int y){
x=pos[x],y=pos[y];
db tx=(seg[y].b-seg[x].b)/(seg[x].k-seg[y].k),ty=seg[x].k*tx+seg[x].b;
return {tx,ty};
}
voiddel(int x){
if(le[x])ri[le[x]]=ri[x];
if(ri[x])le[ri[x]]=le[x];
}
ll check(db mid){
ll res=0;
cntseg=cnt=0;
for(int i=1;i<=n;i++){
db dist=abs(seg[i].k*p.x-p.y+seg[i].b)/sqrt(seg[i].k*seg[i].k+1);
if(dist<=mid){
db a=seg[i].k*seg[i].k+1,b=2*seg[i].k*(seg[i].b-p.y)-2*p.x,c=p.x*p.x+(seg[i].b-p.y)*(seg[i].b-p.y)-mid*mid;
db delta=b*b-4*a*c;
db lx=(-b-sqrt(delta))/(2*a),ly=seg[i].k*lx+seg[i].b;
db rx=(-b+sqrt(delta))/(2*a),ry=seg[i].k*rx+seg[i].b;
cntseg++;
L[cntseg]=R[cntseg]=0;
//fixed pos[cntseg]=i;
A[++cnt]={atan2(ry-p.y,rx-p.x),cntseg};
A[++cnt]={atan2(ly-p.y,lx-p.x),cntseg};
}
}
sort(A+1,A+1+cnt);
for(int i=1;i<=cnt;i++){
if(!L[A[i].se])L[A[i].se]=i;
else R[A[i].se]=i;
}
for(int i=1;i<=cntseg;i++){
sa[i]=i;
}
sort(sa+1,sa+1+cntseg,cmp);
for(int i=1;i<=cnt;i++)le[i]=i-1,ri[i]=i+1;ri[cnt]=0;
for(int i=1;i<=cnt;i++)add(i,1);
for(int i=1;i<=cntseg;i++){
int x=sa[i];
res+=ask(R[x]-1)-ask(L[x]);
add(L[x],-1),add(R[x],-1);
del(L[x]),del(R[x]);
}
return res;
}
db getans(db mid){
ll res=0;
db tot=0;
cntseg=cnt=0;
for(int i=1;i<=n;i++){
db dist=abs(seg[i].k*p.x-p.y+seg[i].b)/sqrt(seg[i].k*seg[i].k+1);
if(dist<=mid){
db a=seg[i].k*seg[i].k+1,b=2*seg[i].k*(seg[i].b-p.y)-2*p.x,c=p.x*p.x+(seg[i].b-p.y)*(seg[i].b-p.y)-mid*mid;
db delta=b*b-4*a*c;
db lx=(-b-sqrt(delta))/(2*a),ly=seg[i].k*lx+seg[i].b;
db rx=(-b+sqrt(delta))/(2*a),ry=seg[i].k*rx+seg[i].b;
cntseg++;
L[cntseg]=R[cntseg]=0;
//fixed pos[cntseg]=i;
A[++cnt]={atan2(ry-p.y,rx-p.x),cntseg};
A[++cnt]={atan2(ly-p.y,lx-p.x),cntseg};
}
}
sort(A+1,A+1+cnt);
for(int i=1;i<=cnt;i++){
if(!L[A[i].se])L[A[i].se]=i;
else R[A[i].se]=i;
}
for(int i=1;i<=cntseg;i++)sa[i]=i;
sort(sa+1,sa+1+cntseg,cmp);
for(int i=1;i<=cnt;i++)le[i]=i-1,ri[i]=i+1;ri[cnt]=0;
for(int i=1;i<=cnt;i++)add(i,1);
for(int i=1;i<=cntseg;i++){
int x=sa[i];
for(int j=ri[L[x]];j!=R[x];j=ri[j]){
point tmp=calc(x,A[j].se);
res++;
tot+=getdist(p,tmp);
}
add(L[x],-1),add(R[x],-1);
del(L[x]),del(R[x]);
}
return tot+(m-res)*mid;
}
intmain(){
ios::sync_with_stdio(false);
cin.tie(0),cout.tie(0);
cin>>n>>p.x>>p.y>>m;
p.x/=1000,p.y/=1000;
//fixedfor(int i=1;i<=n;i++){
cin>>seg[i].k>>seg[i].b;
seg[i].k/=1000,seg[i].b/=1000;
}
db l=0,r=4e9;
for(int i=1;i<=100;i++){
db mid=(l+r)/2;
if(check(mid)<=m)l=mid;
else r=mid;
}
//fixed cout.precision(20);
cout<<getans(l);
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 全程不用写代码,我用AI程序员写了一个飞机大战
· MongoDB 8.0这个新功能碉堡了,比商业数据库还牛
· 记一次.NET内存居高不下排查解决与启示
· 白话解读 Dapr 1.15:你的「微服务管家」又秀新绝活了
· DeepSeek 开源周回顾「GitHub 热点速览」