BZOJ2280 [Poi2011]Plot 二分+倍增+最小圆覆盖
题目传送门
https://lydsy.com/JudgeOnline/problem.php?id=2280
题解
显然对于一段的 \(q_i\)，就是这一段的最小圆覆盖的圆心。
考虑二分一个 \(mid\)，表示每一段的最小圆覆盖的半径的最大值。
然后我们可以在点序列上从开头一直往后走，直到半径大于 \(mid\)，就算上一段并重新开始记录最小圆覆盖，然后把总的段数与 \(m\) 比较。
看上去没什么问题啊，时间复杂度 \(O(n\log eps^{-1})\)，但是怎么 BZOJ 上时间限制有 300s，比紫荆花之恋还长啊。
等等，好像最小圆覆盖的时间复杂度的正确性需要把整个序列都打乱啊。不好，这样就不能保证了题目中要的连续的一段了。
那么我们想一想怎么更正这个算法。
每次从当前的起点 \(i\) 向右扩展的时候，可以二分一个 \(mid\)，然后把 \(i..mid\) 这一段求一下最小圆覆盖？不对不对，这样二分的话，时间复杂度是 \(O(nm\log n\log eps^{-1})\)。
但是我们想要的是，对于长度为 \(l\) 的一段，能够通过 \(l\log l\) 把它二分出来，而不是 \(n\log n\)。
所以这里有一个小 trick，可以一个一个测试 \(i\) 右边 \(1, 2, 4, 8, \cdots\) 位之间的每一段满不满足条件。这样倍增下去，最后我们就可以得到一个信息 \(2^k-1 \leq l < 2^k-1\)。然后我们就可以在这个范围内二分了。这样因为每一次二分的 \(mid\) 都是和 \(l\) 在同一级别的，所以复杂度为 \(O(l\log l)\)。
这样就可以在 \(O(n\log l\log eps^{-1})\) 的时间复杂度内 AC 了。
下面是代码。如上文所属，这个代码的时间复杂度为 \(O(n\log l \log eps^{-1})\)。
#include<bits/stdc++.h>
#define fec(i, x, y) (int i = head[x], y = g[i].to; i; i = g[i].ne, y = g[i].to)
#define dbg(...) fprintf(stderr, __VA_ARGS__)
#define File(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define fi first
#define se second
#define pb push_back
template<typename A, typename B> inline char smax(A &a, const B &b) {return a < b ? a = b , 1 : 0;}
template<typename A, typename B> inline char smin(A &a, const B &b) {return b < a ? a = b , 1 : 0;}
typedef long long ll; typedef unsigned long long ull; typedef std::pair<int, int> pii;
template<typename I>
inline void read(I &x) {
int f = 0, c;
while (!isdigit(c = getchar())) c == '-' ? f = 1 : 0;
x = c & 15;
while (isdigit(c = getchar())) x = (x << 1) + (x << 3) + (c & 15);
f ? x = -x : 0;
}
const int N = 100000 + 7;
const double eps = 1e-8;
int n, m, cnt;
inline int dcmp(const double &x) { return fabs(x) <= eps ? 0 : (x < 0 ? -1 : 1); }
struct Point {
double x, y;
inline Point(const double &x = 0, const double &y = 0) : x(x), y(y) {}
} a[N], b[N], ans[N];
inline double dist(const Point &a, const Point &b) { return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); }
inline Point get_C(const Point &p1, const Point &p2, const Point &p3) {
double c1 = (p1.x * p1.x + p1.y * p1.y) - (p2.x * p2.x + p2.y * p2.y), a1 = 2 * (p1.x - p2.x), b1 = 2 * (p1.y - p2.y);
double c2 = (p2.x * p2.x + p2.y * p2.y) - (p3.x * p3.x + p3.y * p3.y), a2 = 2 * (p2.x - p3.x), b2 = 2 * (p2.y - p3.y);
double x = (b1 * c2 - b2 * c1) / (a2 * b1 - a1 * b2), y = (a1 * c2 - a2 * c1) / (a1 * b2 - a2 * b1);
return Point(x, y);
}
inline std::pair<Point, double> calc(int L, int R) {
int n = 0;
for (int i = L; i <= R; ++i) b[++n] = a[i];
std::random_shuffle(b + 1, b + n + 1);
Point O = b[1];
double r = 0;
for (int i = 2; i <= n; ++i) if (dcmp(dist(O, b[i]) - r) > 0) {
O = b[i], r = 0;
for (int j = 1; j < i; ++j) if (dcmp(dist(O, b[j]) - r) > 0) {
O = Point((b[i].x + b[j].x) / 2, (b[i].y + b[j].y) / 2), r = dist(b[i], O);
for (int k = 1; k < j; ++k) if (dcmp(dist(O, b[k]) - r) > 0) {
O = get_C(b[i], b[j], b[k]);
r = dist(b[i], O);
}
}
}
return std::make_pair(O, r);
}
inline bool check(const double &mid) {
cnt = 0;
int now = 1;
while (now <= n) {
int i = 1;
for (; now + (1 << i) - 1 <= n; ++i)
if (calc(now, now + (1 << i) - 1).se - mid > eps) break;
int l = now + (1 << (i - 1)) - 1, r = std::min(now + (1 << i) - 2, n);
while (l < r) {
int mid2 = (l + r + 1) >> 1;
if (calc(now, mid2).se - mid < eps) l = mid2;
else r = mid2 - 1;
}
ans[++cnt] = calc(now, l).fi, now = l + 1;
}
return cnt <= m;
}
inline void work() {
double l = 0, r = 2000000;
while (r - l >= eps) {
double mid = (l + r) / 2;
if (check(mid)) r = mid;
else l = mid;
}
printf("%.8lf\n", r);
check(r);
printf("%d\n", cnt);
for (int i = 1; i <= cnt; ++i) printf("%.8lf %.8lf\n", ans[i].x, ans[i].y);
}
inline void init() {
srand(time(0) + (ull)new char);
scanf("%d%d", &n, &m);
for (int i = 1; i <= n; ++i) scanf("%lf%lf", &a[i].x, &a[i].y);
}
int main() {
#ifdef hzhkk
freopen("hkk.in", "r", stdin);
#endif
init();
work();
fclose(stdin), fclose(stdout);
return 0;
}
