[BZOJ1336]-[Balkan2002]Alien-最小圆覆盖

说在前面

这个算法还是挺简单的…
然而me并不能理解其正确性
网上的博客也都没有证明
于是手动模拟算法用心感知正确性hhhhh

---

题目

BZOJ1336传送门

题面

给出N个点，让你画一个最小的包含所有点的圆。
点数不超过100000

输入输出格式

输入格式：
第一行一个整数N，表示点的个数
接下来每行两个实数（x，y）表示一个点的坐标

输出格式：
输出共两行，第一行为半径，第二行为圆心坐标，有SPJ

---

解法

就是一个板子题，关于算法流程，最小圆覆盖算法 by commonc
me感性理解（口胡）的正确性
对一个点集作最小圆覆盖，等于对这个点集的凸包作最小圆覆盖
假设已经获得了包含前i个点的最小圆，现在要让它包含第i+1个点。如果第i+1个点本来就在圆内那么不用管。不然的话，这个点一定在原来的圆可以覆盖的凸包外，那么新的圆一定要过这个点。（这一个for是求包含前i+1个点的最小圆覆盖）
现在求包含第i+1个点和前j个点的最小圆覆盖，对于第二个点，按照for j = 1 to i的顺序依次加入，如果第j个点在圆外，同理，新圆一定要过第j个点。
现在求包含第i+1个点和第j个点和前k个点的最小圆覆盖，对于第三个点，按照for k = 1 to j-1的顺序依次加入，原理同上。并且一定不会出现有两个点u,v（i>j>u,v）不能被同时包含的情况（如果出现，那么圆i,u,v应该包含了点j，且先被枚举到，即在第二层循环时j已在圆内，根本就不会再进第三层循环）。
那么一定可以求出「包含第i+1个点和第j个点和前k个点的最小圆覆盖」，当k=j-1时求出「包含第i+1个点和前j个点的最小圆覆盖」，当j=i时求出「包含前i+1个点的最小圆覆盖」。所以最后能够求出所有点的最小圆覆盖

时间复杂度是期望$\Theta(n)$的，上面那篇博客有证明

---

下面是自带大长度的代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

int N ;
double R , PI = 3.1415926535897932384 , eps = 1e-10 ;
struct Vector{
    double x , y ;
    double len(){
        return sqrt( x * x + y * y ) ;
    }
    Vector rotate( double rad ){
        return ( Vector ){ x*cos(rad) - y*sin(rad) , x*sin(rad) + y*cos(rad) } ;
    }
}p[100005] , O ;

typedef Vector Point ;
Vector operator + ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x + B.x , A.y + B.y } ; }
Vector operator - ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x - B.x , A.y - B.y } ; }
Vector operator * ( const Vector &A , const double &p ){ return ( Vector ){ A.x * p , A.y * p } ;}
Vector operator / ( const Vector &A , const double &p ){ return ( Vector ){ A.x / p , A.y / p } ;}
double cross( const Vector &A , const Vector &B ){ return A.x * B.y - A.y * B.x ; }
double dot( const Vector &A , const Vector &B ){ return A.x * B.x + A.y * B.y ; } 
double dis ( const Point &A , const Point &B ){
    double t1 = A.x - B.x , t2 = A.y - B.y ;
    return sqrt( t1 * t1 + t2 * t2 ) ;
}
bool inside( const Point &A ){ return dis( A , O ) <= R ; }
double dcmp( double x ){
    if( x > -eps && x < eps ) return 0 ;
    return x > eps ? 1 : -1 ;
}

Point glt( const Point &P , const Vector &u , const Point &Q , const Vector &v ){
    double t = cross( Q - P , v ) / cross( u , v ) ;
    return P + u * t ;
}

Point cal( Point &A , Point &B , Point &C ){
    Vector a = A - B , b = B - C ;
    if( dcmp( cross( a , b ) ) == 0 ){
        double t1 = dis( A , B ) , t2 = dis( B , C ) , t3 = dis( A , C ) ;
        if( t1 >= t2 && t1 >= t3 ) return ( A + B ) / 2 ;
        if( t2 >= t1 && t2 >= t3 ) return ( B + C ) / 2 ;
        if( t3 >= t1 && t3 >= t2 ) return ( A + C ) / 2 ;
    } return glt( ( A + B ) / 2 , a.rotate( PI / 2 ) , ( B + C ) / 2 , b.rotate( PI / 2 ) ) ;
}

void solve(){
    for( int i = 1 ; i <= N ; i ++ )
    if( !inside( p[i] ) ){
        O = p[i] ; R = 0 ;
        for( int j = 1 ; j < i ; j ++ )
        if( !inside( p[j] ) ){
            O = ( p[i] + p[j] ) / 2 ;
            R = dis( p[i] , O ) ;
            for( int k = 1 ; k < j ; k ++ )
            if( !inside( p[k] ) ){
                O = cal( p[i] , p[j] , p[k] ) ;
                R = dis( p[i] , O ) ;
            }
        }
    }
    printf( "%f\n" , R ) ;
    printf( "%f %f\n" , O.x , O.y ) ;
}

int main(){
    srand( 20031109 ) ;
    scanf( "%d" , &N ) ;
    for( int i = 1 ; i <= N ; i ++ )
        scanf( "%lf%lf" , &p[i].x , &p[i].y ) ;
    random_shuffle( p + 1 , p + N + 1 ) ;
    solve() ;
}