Problem
给定 \(N\) 条不平行的直线 \(A_ix+B_iy+C=0\)，\(N\) 条直线总共会有 \(\frac{N(N-1)}{2}\) 个交点（包含在同一个位置的点，即相同位置算不同的点），找出距离原点第 \(K\) 近的交点的距离。
$2\le N\le 5\times 10^4 $，\(1\le K\le \frac{N(N-1)}{2}\)，$-1000\le|A_i|,|B_i|,|C_i|\le1000(1\le i\le N) $，保证 \(A_i\)，\(B_i\) 均非零。
Sol
发现题目点的个数其实是有 \(\mathcal{O}(n^2)\) 个的，显然不好直接统计，考虑二分答案。
具体地，二分出距离 \(r\)，变为求原点与交点的距离 \(\le r\) 的个数。思考直线 \(L_i\) 与直线 \(L_j\) 的交点会被统计的条件。可以画出一个 \(x^2 +y^2 = r^2\) 的圆，即变为交点在圆内。由于不好直接考虑 \(L_i\)，\(L_j\) 的交点，可以算出所有直线与圆的交点 \((x_{i1}, y_{i1}), (x_{i2}, y_{i2})\)。将两个交点都转化为一个连向原点的与 \(x\) 轴的夹角，这样可以得到一个度数区间 \([l_i, r_i](0\le l_i\le r_i<2\pi)\)，然后两条直线满足交点的条件就是这两个度数区间相交但不包含。统计相交但不包含的区间对数可以通过按左端点从小到大排序后用树状数组统计。
时间复杂度：\(\mathcal{O(n\log n\log \frac{V}{\text{eps}})}\)。这里的 \(\text{eps}\) 是程序中设置的精度。
\(x^2 + y^2 = R^2\) 与 \(Ax + By + C = 0\) 的交点。
- \(B \neq 0\) 时，有 \(y = -\frac{Ax + C}{B}\)：
$ y $ 可以用 $\frac{Ax + C}{B} $ 计算。
- \(B = 0\) 时，\(x = -\frac{C}{A}\)。
Code
#include<bits/stdc++.h>
#define ll long long
#define sz(a) ((int) (a).size())
#define vi vector < int >
#define pb emplace_back
using namespace std;
#define fi first
#define se second
#define double long double
const double eps = 1e-10;
int n, k;
double a[50010], b[50010], c[50010], lsh[100010];
pair < double, double > tmp[50010];
struct node {
int x, y, c;
} pos[50010];
struct BIT {
int n;
vi v;
BIT(int _n = 100000) : v(_n + 10) { n = _n; }
void add(int x, int y) {
for(; x <= n; x += x & -x)
v[x] += y;
}
int ask(int x) {
int ret = 0;
for(; x; x -= x & -x)
ret += v[x];
return ret;
}
int ask(int x, int y) { return ask(y) - ask(x - 1); }
} ;
double sqr(double x) { return x * x; }
double deg(double a, double b) {
return atan2(a, b);
}
pair < double, double > calc(double a, double b, double c, double R) {
double sa = sqr(a), sb = sqr(b), sc = sqr(c), sr = sqr(R);
pair < double, double > p1, p2;
if(fabs(b) > eps) {
double delta = (sa + sb) * sb * sr - sb * sc;
if(delta < 0)
return make_pair(1011451423, 1011451423);
p1.fi = (-a * c + sqrt(delta)) / (sa + sb);
p2.fi = (-a * c - sqrt(delta)) / (sa + sb);
p1.se = -(a * p1.fi + c) / b;
p2.se = -(a * p2.fi + c) / b;
}
else {
p1.fi = p2.fi = -c / a;
p1.se = sqrt(sr - sc / sa);
p2.se = -p1.se;
}
double l = deg(p1.fi, p1.se), r = deg(p2.fi, p2.se);
if(l > r)
swap(l, r);
return make_pair(l, r);
}
map < pair < int, int >, int > mp;
bool check(double mid) {
int num = 0;
for(int i = 1; i <= n; ++i) {
pair < double, double > ret = calc(a[i], b[i], c[i], mid);
if(ret.fi != 1011451423 || ret.se != 1011451423)
++num, tmp[num] = ret, lsh[2 * num - 1] = tmp[num].fi, lsh[2 * num] = tmp[num].se;
}
sort(lsh + 1, lsh + 2 * num + 1);
int len = unique(lsh + 1, lsh + 2 * num + 1) - lsh - 1;
mp.clear();
for(int i = 1; i <= num; ++i) {
pos[i].x = lower_bound(lsh + 1, lsh + len + 1, tmp[i].fi) - lsh;
pos[i].y = lower_bound(lsh + 1, lsh + len + 1, tmp[i].se) - lsh;
++mp[make_pair(pos[i].x, pos[i].y)];
}
int m = 0;
for(auto it : mp)
++m, pos[m].x = it.fi.fi, pos[m].y = it.fi.se, pos[m].c = it.se;
sort(pos + 1, pos + m + 1, [](node a, node b) { return a.x != b.x ? a.x < b.x : a.y < b.y; });
BIT T(len);
int cnt = 0;
for(int i = 1; i <= m; ++i) {
cnt += T.ask(pos[i].x, pos[i].y);
T.add(pos[i].y, pos[i].c);
}
return cnt >= k;
}
int main() {
ios :: sync_with_stdio(false);
cin.tie(0); cout.tie(0);
cin >> n >> k;
for(int i = 1; i <= n; ++i)
cin >> a[i] >> b[i] >> c[i];
double L = eps, R = 5e6;
while(R - L > eps) {
double mid = (L + R) / 2;
if(check(mid))
R = mid;
else
L = mid;
}
cout << fixed << setprecision(10) << L << "\n";
return 0;
}