【学习笔记】CF607E Cross Sum

最后一道数据结构，不能再多了。

而且需要一点计算几何的知识，有点难搞。

分为两个部分求解。

首先考虑找到距离$\le r$的交点数量。发现这等价于圆上两段圆弧相交，因此将圆上的点离散化后排序，用一个主席树来求就做完了。

然后是距离求和。这看起来非常棘手。事实上，只要把所有交点都找出来就做完了。首先可以放心的将圆环从一个位置断开。其次，考虑以某种顺序将所有直线依次删掉。发现当按照长度从大到小删时，假设这条线段是$[l_i,r_i]$，发现只要满足$l_j\in [l_i,r_i]$或者$r_j\in [l_i,r_i]$，那么直线$i,j$就一定相交。那么前面一个问题也得到了解决，只需用树状树组维护就做完了。

当然，事实上我们可以$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>
using namespace std;
const int N=1e5+5;
struct seg{db k,b;
}seg[N];
struct point{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];
bool cmp(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;
}
void add(int x,int y){for(;x<=cnt;x+=x&-x)bit[x]+=y;
}
db getdist(point x,point y){return sqrt((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};
}
void del(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;//fixedpos[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;//fixedpos[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;
}
int main(){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;}//fixedcout.precision(20);cout<<getans(l);
}