洛谷题单指南-分治与倍增-P3517 [POI2011] WYK-Plot
原题链接：https://www.luogu.com.cn/problem/P3517
题意解读：有n个连续的点p1,p2,...,pn，将这n个点分成不超过m堆，每堆点连续，每一堆都缩成一个点qi，要使得原来的点p1~ps距离qi的最大值最小（最相似），求这个相似度，并计算一共分成几堆，以及每堆缩到的点qi的坐标。
解题思路：
要使得若干点缩到一个点，并且原来的点与新点距离的最大值最小，其实就是求以新点为圆心，能覆盖所有点的圆的最小半径。
1、最小圆覆盖
作为前置知识，先来了解一下最小圆覆盖问题：https://www.luogu.com.cn/problem/P1742
最小覆盖圆用到两个重要的性质：
性质一：最小覆盖圆唯一
性质二：若点 P不在点集 S的最小覆盖圆内部，则点 P必在P∪S的最小覆盖圆边上
由于过三点可以确定一个圆，可以枚举每三个点来确定圆是否包含其他点
但是直接枚举三个点会超时，为了规避掉最坏的情况，需要借助随机增量算法！
随机增量算法：
先将点的分布随机打散，可以借助random_shuffle函数，然后
取第一个点独立构成圆，圆心是第一个点，半径是0
for i:2->n
判断点i是否在当前圆内
如果不在圆内，以i点作为圆
for j:1->i-1
判断点j是否在当前圆内
如果不在圆内，以i，j两点构成圆
for k:1->j-1
判断点k是否在圆内
如果不在圆内，以i，j，k三点构成圆
两点构成圆：
点p1、p2构成圆的圆心是两点连线的中间点，半径是两点距离的一半
圆心坐标x = (p1.x + p2.x) / 2, y = (p1.y + p2.y) / 2
半径 = p1到圆心的距离
三点构成圆：
设三个点为(x1,y1)、(x2,y2)、(x3,y3)，圆心为(x,y)，半径为r，可得到方程组：
(x1-x)2 + (y1-y)2 = r2
(x2-x)2 + (y2-y)2 = r2
(x3-x)2 + (y3-y)2 = r2
展开后：
x12 - 2x1x + x2 + y12 - 2y1y + y2 = r2 (1)
x22 - 2x2x + x2 + y22 - 2y2y + y2 = r2 (2)
x32 - 2x3x + x2 + y32 - 2y3y + y2 = r2 (3)
用(1)-(2)，(2)-(3)，再移向后得到
2(x1-x2)x + 2(y1-y2)y = x12 - x22 + y12 - y22
2(x2-x3)x + 2(y2-y3)y = x22 - x32 + y22 - y32
设
a1 = (x1-x2)，b1 = (y1-y2)，c1 = x12 - x22 + y12 - y22
a2 = (x2-x3)，b2 = (y2-y3)，c2 = x22 - x32 + y22 - y32
解方程得到
x = (c1b2-c2b1)/(2a1b2-2a2b1)
y = (c1a2-c2a1)/(2a2b1-2a1b2)
半径是任意一个点到圆心的距离
100分代码（P1742）：
#include <bits/stdc++.h>
using namespace std;
const int N = 100005;
const double eps = 1e-9;
struct Point
{
double x, y;
} a[N];
int n;
//计算两点之间距离
double distance(Point p1, Point p2)
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y));
}
//两个点形成圆的圆心
Point circle(Point p1, Point p2)
{
return {(p1.x+p2.x) / 2, (p1.y+p2.y) / 2};
}
//三个点形成圆的圆心
Point circle(Point p1, Point p2, Point p3)
{
double x1 = p1.x, y1 = p1.y;
double x2 = p2.x, y2 = p2.y;
double x3 = p3.x, y3 = p3.y;
double a1 = x1 - x2, b1 = y1 - y2, c1 = x1*x1 - x2*x2 + y1*y1 - y2*y2;
double a2 = x2 - x3, b2 = y2 - y3, c2 = x2*x2 - x3*x3 + y2*y2 - y3*y3;
double x = (c1*b2 - c2*b1) / (2*a1*b2 - 2*a2*b1);
double y = (c1*a2 - c2*a1) / (2*a2*b1 - 2*a1*b2);
return {x, y};
}
int main()
{
cin >> n;
for(int i = 1; i <= n; i++) cin >> a[i].x >> a[i].y;
random_shuffle(a + 1, a + n + 1);
Point o = a[1]; //圆心
double r = 0; //半径
for(int i = 2; i <= n; i++)
{
if(distance(a[i], o) - r < eps) continue; //i点在圆内
o = a[i]; //更新圆为点i
for(int j = 1; j < i; j++)
{
if(distance(a[j], o) - r < eps) continue; //j点在圆内
o = circle(a[i], a[j]); //更新圆为i，j两点确定的圆
r = distance(a[i], o);
for(int k = 1; k < j; k++)
{
if(distance(a[k], o) - r < eps) continue; //k点在圆内
o = circle(a[i], a[j], a[k]); //更新圆为i，j，k三点确定的圆
r = distance(a[i], o);
}
}
}
cout << fixed << setprecision(10) << r << endl << o.x << " " << o.y;
}2、二分答案求结果
回到本题，要将n个点分成不超过m堆，使得每堆的最小圆覆盖半径最小，可以借助二分。
二分答案mid，找尽可能多的点最小圆覆盖半径不超过mid，这些点即可分为一堆
如何找尽可能多的点分到一堆，可以一个一个枚举，直到最小圆覆盖半径超过mid，则分为一堆，最后统计堆数，不超过m则返回true，否则false
这样枚举+求最小圆覆盖总体复杂度就在O(n^2)，再算上二分，复杂度一共是O(logn*n*n)
因此，还要继续优化！
3、倍增+二分找一堆的右端点
在找尽可能多的点分到一堆时，给定左端点l，可以通过倍增，先找到右端点r，使得最小圆覆盖半径超过mid
然后在l~r之间再二分准确的右端点R，这样l~R分为一堆
再更新l=R+1,继续找下一堆，最后记录堆数，看是否超过m
这样先倍增找一堆的范围，再在范围内二分找右端点，再计算最小圆覆盖半径，算法二分答案，总体复杂度在O(logn*logn*n)
100分代码：
#include <bits/stdc++.h>
using namespace std;
const int N = 1000005;
const double eps = 1e-9;
struct Point
{
double x, y;
} a[N], b[N], anso[N];
int cnt;
int n, m;
Point o; //圆心
double radius; //半径
//计算两点之间距离
double distance(Point p1, Point p2)
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y));
}
//两个点形成圆的圆心
Point circle(Point p1, Point p2)
{
return {(p1.x+p2.x) / 2, (p1.y+p2.y) / 2};
}
//三个点形成圆的圆心
Point circle(Point p1, Point p2, Point p3)
{
double x1 = p1.x, y1 = p1.y;
double x2 = p2.x, y2 = p2.y;
double x3 = p3.x, y3 = p3.y;
double a1 = x1 - x2, b1 = y1 - y2, c1 = x1*x1 - x2*x2 + y1*y1 - y2*y2;
double a2 = x2 - x3, b2 = y2 - y3, c2 = x2*x2 - x3*x3 + y2*y2 - y3*y3;
double x = (c1*b2 - c2*b1) / (2*a1*b2 - 2*a2*b1);
double y = (c1*a2 - c2*a1) / (2*a2*b1 - 2*a1*b2);
return {x, y};
}
//找l~r之间点的最小覆盖圆圆心和半径
void find(int l, int r)
{
int idx = 0;
for(int i = l; i <= r; i++) b[++idx] = a[i];
random_shuffle(b + 1, b + idx + 1);
o = b[1]; //圆心
radius = 0; //半径
for(int i = 2; i <= idx; i++)
{
if(distance(b[i], o) - radius < eps) continue; //i点在圆内
o = b[i]; //更新圆为点i
for(int j = 1; j < i; j++)
{
if(distance(b[j], o) - radius < eps) continue; //j点在圆内
o = circle(b[i], b[j]); //更新圆为i，j两点确定的圆
radius = distance(b[i], o);
for(int k = 1; k < j; k++)
{
if(distance(b[k], o) - radius < eps) continue; //k点在圆内
o = circle(b[i], b[j], b[k]); //更新圆为i，j，k三点确定的圆
radius = distance(b[i], o);
}
}
}
}
bool check(double maxx)
{
cnt = 0;
for(int i = 1; i <= n; i++) //枚举左边界
{
int j;
for(j = 1; i + (1 << j) - 1 <= n; j++) //倍增找右边界
{
find(i, i + (1 << j) - 1);
if(radius - maxx > eps) break; //一段点的最小圆覆盖半径超过限定值
}
int l = i + (1 << j - 1) - 1, r = min(n, i + (1 << j) - 1), ans = i;
//二分找到准确的右端点，使得找到的一段最小圆覆盖半径都在maxx以内
while(l <= r)
{
int mid = (l + r) >> 1;
find(i, mid);
if(radius - maxx < eps) ans = mid, l = mid + 1;
else r = mid - 1;
}
find(i, ans); //更新圆心、半径
i = ans;
cnt++;
anso[cnt] = o;
if(cnt > m) return false;
}
return true;
}
int main()
{
scanf("%d%d", &n, &m);
for(int i = 1; i <= n; i++) scanf("%lf%lf", &a[i].x, &a[i].y);
find(1, n); //在所有点中计算最小圆覆盖的圆心、半径，得到的半径就是最大半径
double l = 0, r = radius;
//二分找答案
while(r - l > eps)
{
double mid = (l + r) / 2;
if(check(mid)) r = mid;
else l = mid;
}
check(r); //更新一次正确答案对应的圆心、半径等
printf("%.8lf\n", r);
printf("%d\n", cnt);
for(int i = 1; i <= cnt; i++) printf("%.8lf %.8lf\n", anso[i].x, anso[i].y);
return 0;
}
