\(\mathcal{Description}\)
Link.
给定平面上 \(n\) 个点，求最小的能覆盖其中至少 \(m\) 个点的圆半径及一个可能的圆心。
\(n\le500\)，坐标值 \(X\in[0,10^4]\)。
\(\mathcal{Solution}\)
不难想到二分答案 \(r\)，以每个点为圆心，\(r\) 为半径作圆，若 \(r\) 合法则能找到一个被至少 \(m\) 个圆覆盖的点。
但是圆的交极难处理，结合数据范围，考虑通过一些枚举操作来简化问题——钦定最终圆心 \(O\) 在以 \(P_i\) 为圆心的圆内，同时考虑到 \(O\) 的最优性，这一限制可以进一步加强为：圆心 \(O\) 在以 \(P_i\) 为圆心的圆周上，于是可以在圆周上按极角做扫描线。到此，我们得到 \(\mathcal O(n^2\log n\log \frac{X}{\epsilon})\) 的算法，应该能过。
接下来是计算几何算法的常见套路：基于数据随机化降低复杂度。设钦定 \(O\) 在以 \(P_i\) 为圆心的圆周上时，最小的合法半径为 \(r_i\)。假设以 \(i=1..n\) 的顺序依次计算 \(r_i\)，此时我们可以用 \(\mathcal O(n\log n)\) 的时间判断某个 \(r_i\) 是否是前缀最小值。若 \(r_i\) 是前缀最小值，我们还需要 \(\mathcal O(n\log n\log\frac{X}{\epsilon})\) 的时间求出 \(r_i\) 的具体值。由于数据随机，所以期望复杂度为
\[\begin{aligned}
T(n) &= \sum_{i=1}^n\left(\frac{\mathcal O\left(n\log n\log\frac{X}{\epsilon}\right)}{i}+\mathcal O(n\log n)\right)\\
&= \mathcal O\left(n\log^2n\log\frac{X}{\epsilon}\right).
\end{aligned}
\]
相对于期望复杂度的优化，感觉第一步通过适当枚举加强限制，简化问题的思路更为重要。不要一心想着优化复杂度就不敢考虑暴力了。
\(\mathcal{Code}\)
/*+Rainybunny+*/
#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)
typedef std::pair<double, bool> PDB;
#define fi first
#define se second
const int MAXN = 500;
const double EPS = 1e-7, PI = acos(-1.);
int n, m;
inline double dabs(const double u) { return u < 0. ? -u : u; }
inline int dcmp(const double u) { return dabs(u) < EPS ? 0 : u < 0 ? -1 : 1; }
inline double norm(double u) {
u = fmod(u, 2. * PI);
return u < 0. ? u + 2. * PI : u;
}
struct Point {
double x, y;
inline Point operator * (const double k) const {
return { k * x, k * y };
}
inline Point operator + (const Point& u) const {
return { x + u.x, y + u.y };
}
inline bool operator == (const Point& u) const {
return !dcmp(x - u.x) && !dcmp(y - u.y);
}
friend inline double dist(const Point& u, const Point& v) {
return sqrt((u.x - v.x) * (u.x - v.x) + (u.y - v.y) * (u.y - v.y));
}
} pnt[MAXN + 5];
inline bool check(const Point& O, const double R, Point& resO, double& resR) {
static std::vector<PDB> evt; evt.clear();
int cur = 0;
rep (i, 1, n) if (!(pnt[i] == O)) {
double d = dist(O, pnt[i]);
if (d > 2. * R) continue;
double alp = acos(0.5 * d / R);
double the = norm(atan2(pnt[i].y - O.y, pnt[i].x - O.x));
double u = norm(the - alp), v = norm(the + alp);
if (u > v) ++cur;
evt.push_back({ u, 1 }), evt.push_back({ v, 0 });
}
std::sort(evt.begin(), evt.end());
if (cur >= m - 1) {
if (R < resR) resO = O + Point{ 0, resR = R };
return true;
}
for (auto [alp, opt]: evt) {
if ((cur += opt ? 1 : -1) >= m - 1) {
if (R < resR) resO = O + Point{ cos(alp), sin(alp) } * (resR = R);
return true;
}
}
return false;
}
int main() {
scanf("%d %d", &n, &m);
rep (i, 1, n) scanf("%lf %lf", &pnt[i].x, &pnt[i].y);
std::shuffle(pnt + 1, pnt + n + 1, std::mt19937(20120712));
double ans = 7.5e3; Point ansP = { -1., -1. };
rep (i, 1, n) {
if (!check(pnt[i], ans, ansP, ans)) continue;
double l = 0., r = ans;
while (l + EPS < r) {
double mid = 0.5 * (l + r);
if (check(pnt[i], mid, ansP, ans)) r = mid;
else l = mid;
}
}
printf("%f\n%f %f\n", ans, ansP.x, ansP.y);
return 0;
}