Poj Solution 2069

http://poj.org/problem?id=2069

#include<iostream>
#include"stdio.h"
#include"math.h"
using namespace std;
#define sq(a) ((a.x)*(a.x)+(a.y)*(a.y)+(a.z)*(a.z))

int n;
typedef double  det[3][3];
struct point
{
    double x,y,z;
};

double hls(det a)
{
    return a[0][0]*a[1][1]*a[2][2]+a[0][1]*a[1][2]*a[2][0]+a[0][2]*a[1][0]*a[2][1]
          -a[0][2]*a[1][1]*a[2][0]-a[0][1]*a[1][0]*a[2][2]-a[0][0]*a[1][2]*a[2][1];
}

bool qiujie(det s,double s0,double s1,double s2,point &o)
{
    det t={{s0,s[0][1],s[0][2]},{s1,s[1][1],s[1][2]},{s2,s[2][1],s[2][2]}};
    double h=hls(s);
    
    if(h==0)return 0;
    o.x=hls(t)/h;
    t[0][0]=s[0][0];t[1][0]=s[1][0];t[2][0]=s[2][0];
    t[0][1]=s0;t[1][1]=s1;t[2][1]=s2;
    o.y=hls(t)/h;
    
    t[0][1]=s[0][1];t[1][1]=s[1][1];t[2][1]=s[2][1];
    t[0][2]=s0;t[1][2]=s1;t[2][2]=s2;
    o.z=hls(t)/h;
    return 1;
}

bool qiuxin(point a,point b,point c,point d,point &o)
{
    det s={{a.x-b.x,a.y-b.y,a.z-b.z},{b.x-c.x,b.y-c.y,b.z-c.z},{c.x-d.x,c.y-d.y,c.z-d.z}};
    double s0=(sq(a)-sq(b))/2,s1=(sq(b)-sq(c))/2,s2=(sq(c)-sq(d))/2;
    
    return qiujie(s,s0,s1,s2,o);
}

bool qiuxin(point a,point b,point c,point &o)
{    
    det s={{a.x-b.x,a.y-b.y,a.z-b.z},{b.x-c.x,b.y-c.y,b.z-c.z}};
    double s0=(sq(a)-sq(b))/2,s1=(sq(b)-sq(c))/2,s2=0,t;
    t=(b.y-a.y)*(c.z-a.z)-(c.y-a.y)*(b.z-a.z);
    s[2][0]=t;s2+=a.x*t;
    t=(b.x-a.x)*(c.z-a.z)-(c.x-a.x)*(b.z-a.z);
    s[2][1]=-t;s2+=-t*a.y;
    t=(b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
    s[2][2]=t;s2+=t*a.z;
    
    return qiujie(s,s0,s1,s2,o);
}
double jl(point &a,point &b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
point p[31];

bool is_ok(point o,double r)
{
    int i;
    for(i=0;i<n;i++)
    if(jl(p[i],o)>r+0.00001)break;
    
    return i==n;
}

int main()
{
    int i,j,k,l;point o;
    double answer,r;
    while(1)
    {
        cin>>n;
        if(n==0)break;
        
        for(i=0;i<n;i++)
        cin>>p[i].x>>p[i].y>>p[i].z;
        
        answer=999999;
        
        for(i=0;i<n;i++)
        for(j=i+1;j<n;j++)
        {
            o.x=(p[i].x+p[j].x)/2;
            o.y=(p[i].y+p[j].y)/2;
            o.z=(p[i].z+p[j].z)/2;
            
            r=jl(o,p[i]);
            if(is_ok(o,r)&&r<answer)answer=r;
            
            for(k=j+1;k<n;k++)
            {
                if(qiuxin(p[i],p[j],p[k],o))
                {
                    r=jl(o,p[i]);
                    if(is_ok(o,r)&&r<answer)answer=r;
                }
                
                for(l=k+1;l<n;l++)
                if(qiuxin(p[i],p[j],p[k],p[l],o))
                {
                    r=jl(o,p[i]);
                    if(is_ok(o,r)&&r<answer)answer=r;
                }    
            }
        }
        
        printf("%.5lfn",answer);
    }
    return 0;
}