题意：
求覆盖平面$n$个点的最小半径圆。

思路：
Welzl’s算法：

• 最直接的方法就是模拟退火，不过这里不用。
• 这里用的是Welzl’s算法求最小覆盖圆，期望复杂度为$O(n)$。基于的结论是：只要得出了覆盖前$i$个点的最小覆盖圆，那么第$i+1$个点如果在圆外，那么第$i+1$个点就应该在圆上，因此可以根据第$i+1$个点与圆的位置，迭代的更新圆上的三个点（或两个点）。
• 对于最小覆盖球也有同样的结论，但是因为要枚举球上4个点，且3点有共面的情况，所以实现比较复杂，推荐使用模拟退火。

三点求圆心问题：
上课无聊证明了一下

$(x0-x1{)}^2+(y0-y1{)}^2=r^2(1)$
$(x0-x2{)}^2+(y0-y2{)}^2=r^2(2)$
$(x0-x3{)}^2+(y0-y3{)}^2=r^2(3)$

$(x0-x1{)}^2-(x0-x2{)}^2+(y0-y1{)}^2-(y0-y2{)}^2=0;(4)$
$(x0-x1{)}^2-(x0-x3{)}^2+(y0-y1{)}^2-(y0-y3{)}^2=0;(5)$

$-2x0x1+2x0x2-2y0y1+2y0y2=x{2}^2-x{1}^2+y{2}^2-y{1}^2;=>x0(x2-x1)+y0(y2-y1)=(x{2}^2-x{1}^2+y{2}^2-y{1}^2)/2;(6)$

$-2x0x1+2x0x3-2y0y1+2y0y3=x{3}^2-x{1}^2+y{3}^2-y{1}^2=>x0(x3-x1)+y0(y3-y1)=(x{3}^2-x{1}^2+y{3}^2-y{1}^2)/2(7)$

$a1=x2-x1,b1=y2-y1,c1=(x{2}^2-x{1}^2+y{2}^2-y{1}^2)/2(8)$
$a2=x3-x1,b2=y3-y1,c2=(x{3}^2-x{1}^2+y{3}^2-y{1}^2)/2(9)$

{$a1\ast x0+b1\ast y0=c1(10)$
{$a2\ast x0+b2\ast y0=c2(11)$

根据行列式求解规则，结合(10),(11)得到结果。

#include<cstdio>
#include<cstring>
#include<string>
#include<set>
#include<iostream>
#include<stack>
#include<queue>
#include<vector>
#include<algorithm>
#include<assert.h>
#include<cmath>
#include<map>
typedef long long ll;
using namespace std;
const int maxn = 1e6 + 10;
const double eps = 1e-8;
int cmp(double x) {
    if (fabs(x) < eps) return 0;
    if (x > 0) return 1;
    return -1;
}

struct point {
    double x, y;
    point() {}
    point(double a, double b) : x(a), y(b) {}
    friend point operator-(const point& a, const point& b) {
        return point(a.x - b.x, a.y - b.y);
    }
    double norm() {
        return sqrt(pow(x, 2) + pow(y, 2));
    }
}p[maxn];
double dist(const point& a, const point& b) {
    return (a - b).norm();
}

void circle_center(point p0, point p1, point& cp) {
    cp.x = (p0.x + p1.x) / 2;
    cp.y = (p0.y + p1.y) / 2;
}

void circle_center(point p0, point p1, point p2, point& cp) {
    double a1 = p1.x - p0.x,b1 = p1.y - p0.y,c1 = (p1.x * p1.x - p0.x * p0.x + p1.y * p1.y - p0.y * p0.y) / 2;
    double a2 = p2.x - p0.x,b2 = p2.y - p0.y,c2 = (p2.x * p2.x - p0.x * p0.x + p2.y * p2.y - p0.y * p0.y) / 2;
    cp.x = (c1 * b2 - b1 * c2) / (a1 * b2 - b1 * a2);
    cp.y = (a1 * c2 - c1 * a2) / (a1 * b2 - b1 * a2);
}

double radius; point center;
bool point_in(const point& p) {
    return cmp(dist(p, center) - radius) < 0;
}

void min_circle_cover(point a[], int n)
{
    radius = 0;
    center = a[0];
    for (int i = 1; i < n; i++) if (!point_in(a[i])) {
        center = a[i], radius = 0;
        for (int j = 0; j < i; j++) if (!point_in(a[j])) {
            circle_center(a[i], a[j], center);
            radius = dist(a[j], center);
            for (int k = 0; k < j; k++) if (!point_in(a[k])) {
                circle_center(a[i], a[j], a[k], center);
                radius = dist(a[k], center);
            }
        }
    }
}

int main()
{
    int n;
    scanf("%d", &n);
    for (int i = 0; i < n; i++)
    {
        double x, y;
        scanf("%lf%lf", &x, &y);
        p[i] = point(x, y);
    }
    random_shuffle(p, p + n);
    min_circle_cover(p, n);
    printf("%.10lf\n%.10lf %.10lf",radius,center.x, center.y);
    
    return 0;
}