【CF607E】Cross Sum（二分+计算几何）

点此看题面

• 给定一个点\(p\)和\(n\)条直线，求到点\(p\)最近的\(m\)个交点到它的距离之和（若某个点多次作为交点则将重复计算）。
• \(n\le5\times10^4,m\le3\times10^7\)

---

一道恶心的卡精度题？样例都和答案相差\(1\)，调了半个多小时没调出来？

心态爆炸之后仔细一看题目限制，并不是绝对误差不超过\(10^{-6}\)，而是相对或绝对误差不超过\(10^{-6}\)！

所以终究还是我瞎。。。

二分半径+树状数组/\(set\)

考虑我们二分一下最近的\(m\)个交点中到点\(p\)最远的点距\(x\)。

那么现在就是要判断在一个半径为\(x\)的圆内是否存在大于等于\(m\)个交点。

在有一个圆的情况下判两直线是否存在交点是很简单的，就是看它们与圆的两个交点在圆周上对应区间是否相交。

我们按左端点从小到大枚举每个区间，则区间\(j\)与\(i\)（\(j<i\)，即\(j\)的左端点小于\(i\)的右端点）相交意味着\(j\)的右端点在\(i\)的左右端点之间。

因此，只要用树状数组就可以维护了。

而最终要求出最近的\(m\)个交点，因为已经二分出了答案，那么交点总数是\(O(m)\)的。

此时线段相交仍然可以利用上面的转化判定，我们只要把树状数组改成\(set\)，枚举先前出现过的右端点在\(i\)的左右端点之间的所有区间，就能求出所有交点了。

代码：\(O(nlog^2n+mlogn)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 50000
#define DB long double
#define eps 1e-20
using namespace std;
int n,m;
struct P {DB x,y;I P(Con DB& a=0,Con DB& b=0):x(a),y(b){}}o;
struct S {DB k,b;I S(Con DB& x=0,Con DB& y=0):k(x),b(y){}}s[N+5];
struct Event
{
	int p;DB A;I Event(CI x=0,Con DB& a=0):p(x),A(a){}
	I bool operator < (Con Event& o) Con {return fabs(A-o.A)>eps?A<o.A:p<o.p;}
}g[2*N+5];
int ct,R[N+5];I void Get(Con DB& x)//求出线段与圆的交点，表示成区间
{
	RI i;DB d,X1,X2;for(ct=0,i=1;i<=n;++i)
	{
		if((d=4*s[i].k*s[i].k*s[i].b*s[i].b-4*(s[i].k*s[i].k+1)*(s[i].b*s[i].b-x*x))<-eps) continue;//解方程X^2+(kX+b)^2=x^2，delta<0无解
		X1=(-2*s[i].k*s[i].b+sqrt(d))/(2*(s[i].k*s[i].k+1)),X2=(-2*s[i].k*s[i].b-sqrt(d))/(2*(s[i].k*s[i].k+1));//x=(-b±sqrt(delta))/(2a)
		R[i]=0,g[++ct]=Event(i,atan2(s[i].k*X1+s[i].b,X1)),g[++ct]=Event(i,atan2(s[i].k*X2+s[i].b,X2));//根据角度转化为对应区间
	}
	for(sort(g+1,g+ct+1),i=ct;i;--i) !R[g[i].p]&&(R[g[i].p]=i);//记录每个区间的右端点位置
}
struct TreeArray
{
	int a[2*N+5];I void Cl() {for(RI i=1;i<=ct;++i) a[i]=0;}//清空
	I void U(RI x) {W(x<=ct) ++a[x],x+=x&-x;}I int Q(RI x,RI t=0) {W(x) t+=a[x],x-=x&-x;return t;}//单点修改；前缀查询
}T;
I bool Check(Con DB& x)//验证答案
{
	RI i,t=0;for(Get(x),T.Cl(),i=1;i<=ct;++i) if(i^R[g[i].p])//对于每个左端点
		{if((t+=T.Q(R[g[i].p])-T.Q(i-1))>=m) return true;T.U(R[g[i].p]);}return false;//加上先前右端点在当前区间内的个数
}
I DB IP_Dis(CI i,CI j) {DB x=(s[j].b-s[i].b)/(s[i].k-s[j].k),y=s[i].k*x+s[i].b;return sqrt(x*x+y*y);}//计算交点到原点距离
set<pair<int,int> > G;set<pair<int,int> >::iterator it;I DB Calc(Con DB& x)//计算答案
{ 
	RI i,c=0;DB t=0;for(Get(x),i=1;i<=ct;++i) if(i^R[g[i].p])
	{
		for(it=G.upper_bound(make_pair(i,0));it!=G.end()&&it->first<R[g[i].p];++it) ++c,t+=IP_Dis(g[i].p,it->second);//枚举先前右端点在区间内的所有线段
		G.insert(make_pair(R[g[i].p],g[i].p));
	}return t+(m-c)*x;//c记录计入答案的总点数，加上(m-c)*x补正贡献
}
int main()
{
	RI i;ios::sync_with_stdio(false),cin>>n>>o.x>>o.y>>m,o.x/=1000,o.y/=1000;
	for(i=1;i<=n;++i) cin>>s[i].k>>s[i].b,s[i].k/=1000,(s[i].b/=1000)+=s[i].k*o.x-o.y;//方便起见，修改解析式使题目中给定点变成原点
	RI t=0;for(i=1;i<=n;++i) fabs(s[i].b)<eps&&++t;if(1LL*t*(t-1)/2>=m) return puts("0"),0;//特判给定点作为交点次数超过m
	DB l=0,r=4e9,mid;for(RI i=1;i<=100;++i) Check(mid=(l+r)/2)?r=mid:l=mid;//二分半径
	return cout<<fixed<<setprecision(8)<<Calc(r)<<endl,0;//计算答案
}