【学习笔记】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;
            //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;
}
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;
    //fixed
    for(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);
}