题解 P1742 【最小圆覆盖】
totorato
2018-12-14 08:52:11
# 最小圆覆盖 luoguP1742
给定$n$个点，求一个最小的圆包围所有的点。
### 随机增量法
**定理1：**如果点$p$不在集合$S$的最小覆盖圆内，则$p$一定在$S\cup\{p\}$的最小覆盖圆上。
根据这个定理，我们可以分三次确定前$i$个点的最小覆盖圆。
- 1.令前$i-1$个点的最小覆盖圆为$C$
- 2.如果第$i$个点在$C$内，则前$i$个点的最小覆盖圆也是$C$
- 3.如果不在，那么第$i$个点一定在前$i$个点的最小覆盖圆上，接着确定前$i-1$个点中还有哪两个在最小覆盖圆上。因此，设当前圆心为$P_i$，半径为0，做固定了第$i$个点的前$i$个点的最小圆覆盖。
- 4.固定了一个点：不停地在范围内找到第一个不在当前最小圆上的点$P_j$，设当前圆心为$(P_i+P_j)/2$，半径为$\frac{1}{2}|P_iP_j|$，做固定了两个点的，前$j$个点外加第$i$个点的最小圆覆盖。
- 5.固定了两个点：不停地在范围内找到第一个不在当前最小圆上的点$P_k$，设当前圆为$P_i,P_j,P_k$的外接圆。
写成伪代码的形式非常简单：
```cpp
圆 C;
for(i=1 to n)
{
if(P[i] 不在 C 内)
{
C = {P[i], 0};
for(j=1 to i-1)
{
if(P[j] 不在 C 内)
{
C = {0.5*(P[i]+P[j]), 0.5*dist(P[i], P[j])};
for(k=1 to j-1)
{
if(P[k] 不在 C 内)
{
C = 外接圆(P[i], P[j], P[k]);
}
}
}
}
}
}
```
只需要三个模式完全相同的for循环就可以搞定。
但是对于这个算法，还有三个地方需要明确：如何求外接圆，以及复杂度是多少。
### 外接圆
可以使用初中的中垂线法。设三个不共线点为$A,B,C$，那么我们可以求$AB,AC$中垂线的交点。
以下我们需要四个工具去完成这个任务：
- 求两个向量的中点
- 将一个向量旋转90度
- 用直线上某一点坐标和其方向向量表示一条直线
- 求两条直线的交点
第一个任务，求两个向量的中点，将横纵坐标相加除以二即可。
第二个任务，将一个向量旋转90度，如果是顺时针旋转，即$(x,y)$变成$(y,-x)$
第三个任务，已经说白了，两个向量就可以表示一条直线。
第四个任务稍微有一点复杂，但是不是很复杂。众所周知，二维向量的叉积可以表示两向量所形成的平行四边形的有向面积。我们可以利用这一点解决直线交点的问题。
用$p_1,v_1$表示$I$，即：$I=p_1+tv_1$
又$I$在绿色直线上，所以又向量$(I-p_2)$和$v_2$构成的平行四边形面积为$0$，即$(p_1+tv_1-p_2)\times v_2 = 0$
解方程得，$t=\frac{(p_2-p_1)\times v_2}{v_1\times v_2}$，交点$I$就是$p_1+tv_1$
利用上述四个工具，我们可以求出三个点的外接圆。
### 复杂度证明
由于一堆点最多只有$3$个点确定了最小覆盖圆，因此$n$个点中每个点参与确定最小覆盖圆的概率不大于$\frac{3}{n}$
所以，每一层循环在第$i$个点处调用下一层的概率不大于$\frac{3}{i}$
那么设算法的三个循环的复杂度分别为$T_1(n),T_2(n),T_3(n)$，则有：
$$\begin{aligned}T_1(n) & = O(n) + \sum_{i=1}^{n}{\frac{3}{i}T_2(i)}\\T_2(n) & = O(n) + \sum_{i=1}^{n}{\frac{3}{i}T_3(i)}\\T_3(n) & = O(n)\end{aligned}$$
不难解得，$T_1(n)=T_2(n)=T_3(n)=O(n)$
### 代码
```cpp
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
struct vec
{
double x, y;
vec (const double& x0 = 0, const double& y0 = 0) : x(x0), y(y0) {}
vec operator + (const vec& t) const {return vec(x+t.x, y+t.y);}
vec operator - (const vec& t) const {return vec(x-t.x, y-t.y);}
vec operator * (const double& t) const {return vec(x*t, y*t);}
vec operator / (const double& t) const {return vec(x/t, y/t);}
const double len2 () const {return x*x + y*y;}
const double len () const {return sqrt(len2());}
vec norm() const {return *this/len();}
vec rotate_90_c () {return vec(y, -x);}
};
double dot(const vec& a, const vec& b) {return a.x*b.x + a.y*b.y;}
double crs(const vec& a, const vec& b) {return a.x*b.y - a.y*b.x;}
vec lin_lin_int(const vec& p0, const vec& v0, const vec& p1, const vec& v1)
{
double t = crs(p1-p0, v1) / crs(v0, v1);
return p0 + v0 * t;
}
vec circle(const vec& a, const vec& b, const vec& c)
{
return lin_lin_int((a+b)/2, (b-a).rotate_90_c(), (a+c)/2, (c-a).rotate_90_c());
}
int n;
vec pot[100005];
int main()
{
scanf("%d", &n);
for(int i=1; i<=n; i++) scanf("%lf%lf", &pot[i].x, &pot[i].y);
random_shuffle(pot+1, pot+n+1);
vec o;
double r2 = 0;
for(int i=1; i<=n; i++)
{
if((pot[i]-o).len2() > r2)
{
o = pot[i], r2 = 0;
for(int j=1; j<i; j++)
{
if((pot[j]-o).len2() > r2)
{
o = (pot[i]+pot[j])/2, r2 = (pot[j]-o).len2();
for(int k=1; k<j; k++)
{
if((pot[k]-o).len2() > r2)
{
o = circle(pot[i], pot[j], pot[k]), r2 = (pot[k]-o).len2();
}
}
}
}
}
}
printf("%.10lf\n%.10lf %.10lf\n", sqrt(r2), o.x, o.y);
return 0;
}
```