O(n)最小圆覆盖算法 （随机增量法）

为什么叫随机增量我也不知道
首先考虑一个点集A，设他的最小覆盖圆为g(A)。 g(A)是存在且唯一的。 （假设存在两个，则必定存在更小覆盖圆）
而且g(A)必定满足以下两者之一：
1. 圆上有三个（或以上）A中的点
2. 圆上有两个A中的点并且其连线为直径。
若上述条件不成立则显然有更小覆盖圆。（感性调整）

就像这样：

记K(A)为A中在g(A)上的点集合，这样的点下文叫做关键点。

下面给出一个关键的定理

若点b不在g(A)内，则b属于K(A+b)。

反证：
假设b不属于K(A+b)，
则b在g(A+b)内（不是关键点），那么g(A)与g(A+b)是同圆（有没有他都一样）。又因为b不在g(A)内，矛盾。

设g(A,x)为，点集A的最小覆盖圆，并且满足x在圆上。 （x,y不是A中的点）
设g(A,x,y)为，点集A的最小覆盖圆，并且满足x,y在圆上。
类似的，我们可以证明：
若点b不在g(A,x)内，则b是g(A,x)的关键点
若点b不在g(A,x,y)内，则b是g(A,x,y)的关键点

先对点随机排序。 这对我们的复杂度分析至关重要。
已知1..i-1的最小覆盖圆C的圆心与半径。
g(1)=圆心为点1，半径为0.

从2开始，我们顺序求g(1..i)。

FirstStep： 对于一个新的i

若i在C中则不需要额外操作。
否则，由上面证明的定理可得i在g(1..i) （点1..i的最小覆盖圆，下面都这么简写了）上。
只要求得g(1..i-1,i)即可得到g(1..i)。，相当于此时进入一个子过程变相求当前g(i)

转化1

初始圆=圆心为点i，半径为0圆。
顺序求g(1..i-1,i)

依次令j=从1到i-1，若点j不在g(1..j-1,i)中（若在则不需要额外操作），则同理得j在g(1..j,i)上。
只要求得g(1..j-1,i,j)即可得到g(1..j,i)，相当于此时进入一个子过程变相求当前g(j,i)
（点集为1..j并且i在圆上的最小覆盖圆，表示有点不严谨理解一下）

转化2

初始圆=点i,j为直径所成的圆。
顺序求g(1..j-1,i,j)

依次令k=从1到j-1，若k不在g(1..k-1,i,j)中，则同理得k在g(1..k,i,j)上。

这下我们就知道g(1..k,i,j)上有三个点i,j,k了，不需要继续找在圆上的点了。 因此将圆设为找到的定圆，继续直到求出g(1..j-1,i,j)为止。

时间复杂度的证明

前提是随机序列！
关键分析，第一层来说： 对于一个新的点i，i在不当前圆中的概率为3/i.
证明： 因为g(1..i)只由至多三个关键点决定。这三个关键点是从1..i中选三个。选中i时才会进入下一层 （蜜汁证明….怪怪的但随机情况好像正确）

第二三层类似，然后分析一下复杂度就可以得到期望O(n)的结论。

实现

g是当前圆心，r是当前半径。

    for (int i=1; i<=10*n; i++) { //随机打乱。
        int x=rand()%n+1,y=rand()%n+1;
        swap(p[x],p[y]);
    }
    g=(P) {p[1].x,p[1].y}; r=0;
    for (int i=2; i<=n; i++) if (!zai(p[i])) {
        g=(P){p[i].x,p[i].y}; r=0;
        for (int j=1; j<i; j++) if (!zai(p[j])) {
            g=mid(p[i],p[j]);
            r=sqrt(dis2(p[i],p[j]))*0.5;
            for (int k=1; k<j; k++) if (!zai(p[k])) {
                make(p[i],p[j],p[k]);
            }
        }
    }

make是求出三点所得的圆。 长这样

struct P{db x,y;};
struct L{db k,b;};
L zx(P g,db k,db kx) {
    if (kx==0) return (L){G,g.x}; return (L){k/kx,g.y-k/kx*g.x};
}
P jiaodian(L u,L v) {
    if (v.k==G) swap(u,v);
    db x=0;
    if (u.k==G) x=u.b; else x=(v.b-u.b)/(u.k-v.k);
    return (P) {x,v.b+v.k*x};
}
void make(P a,P b,P c) {
    P yx=jiaodian(zx(mid(a,b),-(b.x-a.x),(b.y-a.y)),zx(mid(a,c),-(c.x-a.x),(c.y-a.y)));
    g=yx; r=sqrt(dis2(yx,a));
}

就是求中垂线交点。 注意特判与y轴平行的情况。 （本来想打向量的但是就这么一丢丢特判一下算了。。。）