1231-Super Star

問題概要

最小包含球を求めよ

解法

全ての候補点を求め，その候補点に対し二分探索を行なう．
候補点は
・4点からの距離が等しい点
・3点からの距離が等しく，その3点と同一平面上にある点
・2点からの距離が等しく，その2点による直線上にある点
の3種類であり，列挙はN^4で出来る．
よって，N^4*N*二分探索の計算回数で求めることが出来て，AOJだとギリギリ間に合う．

候補点を求める時に一緒に半径も求めると二分探索は必要ないかもしれない

実装(C++)

#include <complex>
#include <iostream>
#include <vector>
#include <iomanip>
#include <cstdio>
#include <string>
#include <cassert>
#include <algorithm>
#include <iterator>
#include <valarray>
const double EPS=1e-8;
using namespace std;
typedef valarray<double> Point;
typedef vector<double> Array;
typedef vector<Array> Matrix;
double size(const Point &a){
    return sqrt((a*a).sum());
}
double size2(const Point &a){
    return ((a*a).sum());
}
double dot(const Point &a,const Point &b){
    return (a*b).sum();
}
//連立方程式を解く
vector<double> solve(Matrix M){
    int n=M.size();
    for(int i=0;i<n;i++){
        if(abs(M[i][i])<EPS){
            for(int j=i+1;j<n;j++){
                if(abs(M[j][i])>EPS){
                    swap(M[i],M[j]);
                    break;
                }
            }
        }
        if(abs(M[i][i])<EPS)return vector<double>();//解無し
        double t=M[i][i];
        for(int j=0;j<n+1;j++){
            M[i][j]/=t;
        }
        for(int j=0;j<n;j++){
            if(i==j)continue;
            t=M[j][i];
            for(int k=0;k<n+1;k++){
                M[j][k]-=M[i][k]*t;
            }
        }
    }
    vector<double> res;
    for(int i=0;i<n;i++){
        res.push_back(M[i][n]);
    }
    return res;
}
int n;
vector<Point> points;
//半径を二分探索によって定める
double solve(const Point &a,double res_){
    double s=0,e=res_*1.5;
    double res;
    while(e-s>5e-6){
        res=(s+e)/2;
        int j;
        for(j=0;j<n;j++){
            if(size(points[j]-a)>res){
                break;
            }
        }
        if(j==n){
            e=res;
        }else{
            s=res;
        }
    }
    return res;
}

double solve(){
    double res=300;
    //4点から等距離にある点を全て求める
    for(int i=0;i<n;i++){
        for(int j=i+1;j<n;j++){
            for(int k=j+1;k<n;k++){
                for(int l=k+1;l<n;l++){
                    Matrix M(3,Array(4));
                    Point *p[3];
                    p[0]=&points[j];
                    p[1]=&points[k];
                    p[2]=&points[l];
                    for(int q=0;q<3;q++){
                        M[q][3]=0;
                        for(int r=0;r<3;r++){
                            M[q][r]=2*((*p[q])[r]-points[i][r]);
                            M[q][3]+=(*p[q])[r]*(*p[q])[r];
                        }
                        M[q][3]-=size2(points[i]);
                    }
                    vector<double> tmp=solve(M);
                    if(tmp.empty())continue;
                    Point k(3);
                    for(int q=0;q<3;q++){
                        k[q]=tmp[q];
                    }
                    res=min(res,solve(k,res));
                }
            }
        }
    }

    //3点から等しい距離にある点を全て求める
    for(int i=0;i<n;i++){
        for(int j=i+1;j<n;j++){
            for(int k=j+1;k<n;k++){
                Point a=points[j]-points[i];
                Point b=points[k]-points[i];
                double A=dot(a,a);
                double B=dot(a,b);
                double C=dot(b,b);
                Matrix M(2,Array(3));
                M[0][0]=2*B;
                M[0][1]=2*C;
                M[0][2]=C;
                M[1][0]=2*A;
                M[1][1]=2*B;
                M[1][2]=A;
                vector<double> tmp=solve(M);
                if(tmp.empty())continue;
                Point k=a*tmp[0]+b*tmp[1]+points[i];
                res=min(res,solve(k,res));
            }
        }
    }
    //2点から等しい距離にある点
    for(int i=0;i<n;i++){
        for(int j=i+1;j<n;j++){
            Point k=(points[i]+points[j])*.5;
            res=min(res,solve(k,res));
        }
    }
    return res;
}
int main(){
    while(1){
        cin>>n;
        if(!n)break;
        points=vector<Point>(n,Point(3));
        for(int i=0;i<n;i++){
            for(int j=0;j<3;j++){
                cin>>points[i][j];
            }
        }
        cout<<setprecision(5)<<setiosflags(ios::fixed)<<solve()<<endl;

    }

}