bzoj 1336 最小圆覆盖

最小圆覆盖

问题：给定平面上的一个点集，求半径最小的一个圆，使得点集中的点都在其内部或上面。

随机增量算法：

　　定义：点集A的最小圆覆盖是Circle(A)

　　定理：如果Circle(A)=C₁，且a不被C₁覆盖，那么a在Circle(AU{a})的边界上。

　　证明：换一种找最小圆覆盖的思路，我们初始化一些圆，圆心为A中的点，半径为0，并且让半径慢慢变大，必定存在一个时刻，所有圆的交集由空变为非空，那个最开始的非空交集是一个点，并且就是我们最小圆覆盖的圆心位置。当A中的所有点代表的圆有交集时，点a代表的圆还没有到达那个点（否则点a就被C₁覆盖掉了），我们让半径继续增大，必然会有一个时刻点a代表的圆与A代表的圆的公共区域相交，这个点就是AU{a}的最小圆覆盖的圆心，它到点a的距离就是半径。

　　算法：

c = ( p[1] )
for i = 2 to n
    if p[i] in c then continue
    c = ( p[i] )
    for j = 1 to i-1
        if p[j] in c then continue
        c = ( p[i], p[j] )
        for k = 1 to j-1
            if p[k] in c then continue
            c = ( p[i], p[j], p[k] )

　　（(p[i],p[j])代表包含这两个点的最小的圆）

　　第一层循环的循环不变量是：c是p[1],p[2],...,p[i-1]的最小圆覆盖。

　　第二层循环的循环不变量是：c是p[1],p[2],...,p[j-1]和p[i]的最小圆覆盖。

　　第一层循环的循环不变量是：c是p[1],p[2],...,p[k-1],p[i]和p[j]的最小圆覆盖。

　　转移用上面的定理证明。

　　有个性质：上面伪代码的第10行中的三个点不可能共线（只需分别证明三个点中的一个不会在另外两个代表的线段上就行了）。

/**************************************************************
    Problem: 1336
    User: idy002
    Language: C++
    Result: Accepted
    Time:372 ms
    Memory:2372 kb
****************************************************************/
 
#include <cstdio>
#include <cmath>
#include <algorithm>
#define line(a,b) ((b)-(a))
#define N 100010
#define eps 1e-10
using namespace std;
 
int sg( double x ) { return (x>-eps)-(x<eps); }
struct Vector {
    double x, y;
    Vector(){}
    Vector( double x, double y ):x(x),y(y){}
    Vector operator+( const Vector &b ) const { return Vector(x+b.x,y+b.y); }
    Vector operator-( const Vector &b ) const { return Vector(x-b.x,y-b.y); }
    Vector operator*( double b ) const { return Vector(x*b,y*b); }
    Vector operator/( double b ) const { return Vector(x/b,y/b); }
    double operator^( const Vector &b ) const { return x*b.y-y*b.x; }
    double len() { return sqrt(x*x+y*y); }
    Vector normal() { return Vector(-y,x); }
};
typedef Vector Point;
Point inter( Point P, Vector u, Point Q, Vector v ) {
    return P+u*((line(P,Q)^v)/(u^v));
}
struct Circle {
    Point o;
    double r;
    Circle(){}
    Circle( Point &a ) {
        o = a;
        r = 0;
    }
    Circle( Point &a, Point &b ) {
        o = (a+b)/2;
        r = (a-b).len()/2;
    }
    Circle( Point &a, Point &b, Point &c ) {
        Point P=(a+b)/2, Q=(b+c)/2;
        Vector u=(a-b).normal(), v=(b-c).normal();
        o = inter(P,u,Q,v);
        r = (a-o).len();
    }
    bool contain( Point &a ) {
        return sg( line(a,o).len()-r ) <= 0;
    }
};
 
int n;
Point pts[N];
Circle cir;
 
int main() {
    scanf( "%d", &n );
    for( int i=1; i<=n; i++ ) {
        double x, y;
        scanf( "%lf%lf", &x, &y );
        pts[i] = Point(x,y);
    }
    random_shuffle( pts+1, pts+1+n );
    cir = Circle(pts[1]);
    for( int i=2; i<=n; i++ ) {
        if( cir.contain(pts[i]) ) continue;
        cir = Circle(pts[i]);
        for( int j=1; j<i; j++ ) {
            if( cir.contain(pts[j]) ) continue;
            cir = Circle(pts[i],pts[j]);
            for( int k=1; k<j; k++ ) {
                if( cir.contain(pts[k]) ) continue;
                cir = Circle(pts[i],pts[j],pts[k]);
            }
        }
    }
    printf( "%.10lf
", cir.r );
    printf( "%.10lf %.10lf
", cir.o.x, cir.o.y );
}