#include<cassert>
#include<cstdio>
#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
#define REP(i,b,n) for(int i=b;i<n;i++)
#define rep(i,n)   REP(i,0,n)

const double eps = 1e-10;

class Point{
public:
  double x,y,z;
  double dist(){
    return sqrt(x*x+y*y+z*z);
  }
  Point operator+(Point a)const{
    return (Point){x+a.x,y+a.y,z+a.z};
  }
  Point operator-(Point a)const{
    return (Point){x-a.x,y-a.y,z-a.z};
  }
  Point operator*(double a)const{
    return (Point){a*x,a*y,a*z};
  }
  Point operator/(double a)const{
    return (Point){x/a,y/a,z/a};
  };
};

//Lib
double dot(Point A,Point B){
  return A.x*B.x + A.y*B.y + A.z*B.z;
}

//Lib
Point CircumCenter(Point A,Point B,Point C){
  Point ret;
  double a=(B-C).dist(),b=(C-A).dist(),c=(A-B).dist();
  double cosA=dot(B-A,C-A)/(b*c);
  double s=(c-b*cosA)/(2*c*(1-cosA*cosA));
  double t=(b-c*cosA)/(2*b*(1-cosA*cosA));
  ret=(B-A)*s+(C-A)*t+A;
  return ret;
}

//Lib
Point normal_vector(Point a,Point b){
  Point ret =(Point){a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
  double len=ret.dist();
  ret=ret/len;
  return ret;
}

double isp(Point a,Point b,Point c){
  return (a-c).dist()+(b-c).dist() < (a-b).dist() + eps;
}

bool is_contain_all(Point a,double r,vector<Point> &in){
  rep(i,in.size()){
    if ((in[i]-a).dist() <r+eps);
    else return false;
  }
  return true;
}

double solve(Point cir,Point norm,vector<Point> &in,Point check,
	     double d,double f){//d=(A+B)/2).dist(),f=(cir-(A+B)/2)/dist()
  double l = 0,r=40;
  double mid,ret=1e6;
  //  while(r-l >eps){
  rep(i,100){
    double mid = (l+r)/2;
    Point vec = cir+norm*mid;
    double radius = (vec-check).dist();
    if (isnan(radius))exit(1);
    if (is_contain_all(vec,radius,in))r=mid-eps,ret=min(ret,radius);
    else l=mid+eps;
  }
  return ret;
}

double solve_online(Point a,Point b,vector<Point> &in){
  Point mid = (a+b)/2;
  double radius=(mid-a).dist();
  if (is_contain_all(mid,radius,in))return radius;
  else return 1e10;
}

const int N = 4;
//n matrix size  a00+x0+...=c0
//k current step
//mat Matrix
bool select_pivot(int n,int k,double mat[N][N]){
  int s=k;//selected pivot
  double val=fabs(mat[k][k]);
  REP(i,k+1,n)if (fabs(mat[i][k]) >val)val=fabs(mat[i][k]),s=i;

  if ( val<eps)return false;//cannot solve this linear equations.
  REP(i,k,n+1)swap(mat[k][i],mat[s][i]);
  return true;//no problem for this pivot;
}


Point gauss(int n,double mat[N][N]){
  rep(j,n){
    bool isok=false;
    rep(i,n){
      if (mat[i][j] > eps)isok=true;
    }
    if (!isok)return (Point){-1000,-1000,-1000};
  }
  
  //step1
  rep(k,n){
    if ( select_pivot(n,k,mat)==false)return (Point){-1000,-1000,-1000};
    REP(i,k+1,n+1)mat[k][i]/=mat[k][k];//devided by m[k][k];
    mat[k][k]=1;//for looking the equation
    REP(i,k+1,n){//ツ甘クツ氾ゥ
      REP(j,k+1,n+1)mat[i][j]+=-mat[i][k]*mat[k][j];
      mat[i][k]=0;//For looking the equation
    }
  }
  //step2
  for(int i=n-1;i>=0;i--){//solve ith row
    REP(j,i+1,n){
      mat[i][n]-=mat[i][j]*mat[j][n];
    }
  }
  /*
  cout <<"begin"<<endl;
  rep(i,n){
    rep(j,n+1)cout << mat[i][j] <<" " ;
    cout << endl;
  }
  cout << "end"<<endl<<endl;
  */
  //return (Point){mat[0][1],mat[1][2],mat[2][3]};
  return (Point){mat[0][3],mat[1][3],mat[2][3]};
}


double solve(vector<Point> & in){
  double ret = 1e6;
  double mat[N][N];
  rep(i,in.size()){
    REP(j,i+1,in.size()){
      REP(k,j+1,in.size()){
	REP(l,k+1,in.size()){
	  mat[0][0]=2*(in[j].x-in[i].x);
	  mat[0][1]=2*(in[j].y-in[i].y);
	  mat[0][2]=2*(in[j].z-in[i].z);
	  mat[0][3]=(in[i].x*in[i].x-in[j].x*in[j].x+
		      in[i].y*in[i].y-in[j].y*in[j].y+
		      in[i].z*in[i].z-in[j].z*in[j].z);
	  mat[1][0]=2*(in[k].x-in[i].x);
	  mat[1][1]=2*(in[k].y-in[i].y);
	  mat[1][2]=2*(in[k].z-in[i].z);
	  mat[1][3]=(in[i].x*in[i].x-in[k].x*in[k].x+
		      in[i].y*in[i].y-in[k].y*in[k].y+
		      in[i].z*in[i].z-in[k].z*in[k].z);
	  mat[2][0]=2*(in[l].x-in[i].x);
	  mat[2][1]=2*(in[l].y-in[i].y);
	  mat[2][2]=2*(in[l].z-in[i].z);
	  mat[2][3]=(in[i].x*in[i].x-in[l].x*in[l].x+
		      in[i].y*in[i].y-in[l].y*in[l].y+
		      in[i].z*in[i].z-in[l].z*in[l].z);
	  Point center = gauss(3,mat);
	  double r[4]  ={
	    (in[i]-center).dist(),
	    (in[j]-center).dist(),
	    (in[k]-center).dist(),
	    (in[l]-center).dist(),
	  }; 
	  /*	  
	  cout << (in[i]-center).dist() <<" " <<
	    (in[j]-center).dist() <<" " << 
	    (in[k]-center).dist() << " " <<
	    (in[l]-center).dist() << endl;
	  */	  	  
	  //cout << center.x <<" " << center.y <<" " << center.z << endl;
	  
	  if (abs(r[0]-r[1])>eps || abs(r[0]-r[2])>eps||
	      abs(r[0]-r[3])>eps)continue;
	  //if (is_contain_all(center,r[0],in))ret=min(r[0],ret);
	}
      }
    }
  }


  rep(i,in.size()){
    REP(j,i+1,in.size()){
      REP(k,j+1,in.size()){
	if (isp(in[i],in[j],in[k])||isp(in[i],in[k],in[j])
	    ||isp(in[j],in[k],in[i])){
	  continue;
	}
	Point cir = CircumCenter(in[i],in[j],in[k]);
	double r =(cir-in[i]).dist();
	/*
	cout << (cir-in[i]).dist() <<" " 
	     << (cir-in[j]).dist() <<" " 
	     << (cir-in[k]).dist() << endl;
	*/
	if (is_contain_all(cir,r,in))ret=min(ret,r);
	
	Point norm= normal_vector(in[j]-in[i],in[k]-in[i]);
	Point tmp = (in[i]+in[j])/2-in[i];
	Point tmp2= (cir-(in[i]+in[j])/2);
	ret=min(ret,solve(cir,norm,in,in[i],tmp.dist(),tmp2.dist()));
	norm=norm*-1;
	ret=min(ret,solve(cir,norm,in,in[i],tmp.dist(),tmp2.dist()));
	
      }
    }
  }
  
  rep(i,in.size()){
    REP(j,i+1,in.size()){
      Point center=(in[i]+in[j])/2;
      double r = (center-in[i]).dist();
      if (is_contain_all(center,r,in))ret=min(ret,r);
    }
  }

  return ret;
}

int main(){
  int n;
  while(cin>>n && n){
    vector<Point> in;
    rep(i,n){
      double x,y,z;
      cin>>x>>y>>z;
      in.push_back((Point){x,y,z});
    }

    printf("%.8lf\n",solve(in));
  }
}
### Prompt

Develop a solution in Cpp to the problem described below:

During a voyage of the starship Hakodate-maru (see Problem A), researchers found strange synchronized movements of stars. Having heard these observations, Dr. Extreme proposed a theory of "super stars". Do not take this term as a description of actors or singers. It is a revolutionary theory in astronomy.

According to this theory, stars we are observing are not independent objects, but only small portions of larger objects called super stars. A super star is filled with invisible (or transparent) material, and only a number of points inside or on its surface shine. These points are observed as stars by us.

In order to verify this theory, Dr. Extreme wants to build motion equations of super stars and to compare the solutions of these equations with observed movements of stars. As the first step, he assumes that a super star is sphere-shaped, and has the smallest possible radius such that the sphere contains all given stars in or on it. This assumption makes it possible to estimate the volume of a super star, and thus its mass (the density of the invisible material is known).

You are asked to help Dr. Extreme by writing a program which, given the locations of a number of stars, finds the smallest sphere containing all of them in or on it. In this computation, you should ignore the sizes of stars. In other words, a star should be regarded as a point. You may assume the universe is a Euclidean space.



Input

The input consists of multiple data sets. Each data set is given in the following format.


n
x1 y1 z1
x2 y2 z2
...
xn yn zn


The first line of a data set contains an integer n, which is the number of points. It satisfies the condition 4 ≤ n ≤ 30.

The locations of n points are given by three-dimensional orthogonal coordinates: (xi, yi, zi) (i = 1,..., n). Three coordinates of a point appear in a line, separated by a space character.

Each value is given by a decimal fraction, and is between 0.0 and 100.0 (both ends inclusive). Points are at least 0.01 distant from each other.

The end of the input is indicated by a line containing a zero.

Output

For each data set, the radius ofthe smallest sphere containing all given points should be printed, each in a separate line. The printed values should have 5 digits after the decimal point. They may not have an error greater than 0.00001.

Example

Input

4
10.00000 10.00000 10.00000
20.00000 10.00000 10.00000
20.00000 20.00000 10.00000
10.00000 20.00000 10.00000
4
10.00000 10.00000 10.00000
10.00000 50.00000 50.00000
50.00000 10.00000 50.00000
50.00000 50.00000 10.00000
0


Output

7.07107
34.64102

### Response

```cpp
#include<cassert>
#include<cstdio>
#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
using namespace std;
#define REP(i,b,n) for(int i=b;i<n;i++)
#define rep(i,n)   REP(i,0,n)

const double eps = 1e-10;

class Point{
public:
  double x,y,z;
  double dist(){
    return sqrt(x*x+y*y+z*z);
  }
  Point operator+(Point a)const{
    return (Point){x+a.x,y+a.y,z+a.z};
  }
  Point operator-(Point a)const{
    return (Point){x-a.x,y-a.y,z-a.z};
  }
  Point operator*(double a)const{
    return (Point){a*x,a*y,a*z};
  }
  Point operator/(double a)const{
    return (Point){x/a,y/a,z/a};
  };
};

//Lib
double dot(Point A,Point B){
  return A.x*B.x + A.y*B.y + A.z*B.z;
}

//Lib
Point CircumCenter(Point A,Point B,Point C){
  Point ret;
  double a=(B-C).dist(),b=(C-A).dist(),c=(A-B).dist();
  double cosA=dot(B-A,C-A)/(b*c);
  double s=(c-b*cosA)/(2*c*(1-cosA*cosA));
  double t=(b-c*cosA)/(2*b*(1-cosA*cosA));
  ret=(B-A)*s+(C-A)*t+A;
  return ret;
}

//Lib
Point normal_vector(Point a,Point b){
  Point ret =(Point){a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};
  double len=ret.dist();
  ret=ret/len;
  return ret;
}

double isp(Point a,Point b,Point c){
  return (a-c).dist()+(b-c).dist() < (a-b).dist() + eps;
}

bool is_contain_all(Point a,double r,vector<Point> &in){
  rep(i,in.size()){
    if ((in[i]-a).dist() <r+eps);
    else return false;
  }
  return true;
}

double solve(Point cir,Point norm,vector<Point> &in,Point check,
	     double d,double f){//d=(A+B)/2).dist(),f=(cir-(A+B)/2)/dist()
  double l = 0,r=40;
  double mid,ret=1e6;
  //  while(r-l >eps){
  rep(i,100){
    double mid = (l+r)/2;
    Point vec = cir+norm*mid;
    double radius = (vec-check).dist();
    if (isnan(radius))exit(1);
    if (is_contain_all(vec,radius,in))r=mid-eps,ret=min(ret,radius);
    else l=mid+eps;
  }
  return ret;
}

double solve_online(Point a,Point b,vector<Point> &in){
  Point mid = (a+b)/2;
  double radius=(mid-a).dist();
  if (is_contain_all(mid,radius,in))return radius;
  else return 1e10;
}

const int N = 4;
//n matrix size  a00+x0+...=c0
//k current step
//mat Matrix
bool select_pivot(int n,int k,double mat[N][N]){
  int s=k;//selected pivot
  double val=fabs(mat[k][k]);
  REP(i,k+1,n)if (fabs(mat[i][k]) >val)val=fabs(mat[i][k]),s=i;

  if ( val<eps)return false;//cannot solve this linear equations.
  REP(i,k,n+1)swap(mat[k][i],mat[s][i]);
  return true;//no problem for this pivot;
}


Point gauss(int n,double mat[N][N]){
  rep(j,n){
    bool isok=false;
    rep(i,n){
      if (mat[i][j] > eps)isok=true;
    }
    if (!isok)return (Point){-1000,-1000,-1000};
  }
  
  //step1
  rep(k,n){
    if ( select_pivot(n,k,mat)==false)return (Point){-1000,-1000,-1000};
    REP(i,k+1,n+1)mat[k][i]/=mat[k][k];//devided by m[k][k];
    mat[k][k]=1;//for looking the equation
    REP(i,k+1,n){//ツ甘クツ氾ゥ
      REP(j,k+1,n+1)mat[i][j]+=-mat[i][k]*mat[k][j];
      mat[i][k]=0;//For looking the equation
    }
  }
  //step2
  for(int i=n-1;i>=0;i--){//solve ith row
    REP(j,i+1,n){
      mat[i][n]-=mat[i][j]*mat[j][n];
    }
  }
  /*
  cout <<"begin"<<endl;
  rep(i,n){
    rep(j,n+1)cout << mat[i][j] <<" " ;
    cout << endl;
  }
  cout << "end"<<endl<<endl;
  */
  //return (Point){mat[0][1],mat[1][2],mat[2][3]};
  return (Point){mat[0][3],mat[1][3],mat[2][3]};
}


double solve(vector<Point> & in){
  double ret = 1e6;
  double mat[N][N];
  rep(i,in.size()){
    REP(j,i+1,in.size()){
      REP(k,j+1,in.size()){
	REP(l,k+1,in.size()){
	  mat[0][0]=2*(in[j].x-in[i].x);
	  mat[0][1]=2*(in[j].y-in[i].y);
	  mat[0][2]=2*(in[j].z-in[i].z);
	  mat[0][3]=(in[i].x*in[i].x-in[j].x*in[j].x+
		      in[i].y*in[i].y-in[j].y*in[j].y+
		      in[i].z*in[i].z-in[j].z*in[j].z);
	  mat[1][0]=2*(in[k].x-in[i].x);
	  mat[1][1]=2*(in[k].y-in[i].y);
	  mat[1][2]=2*(in[k].z-in[i].z);
	  mat[1][3]=(in[i].x*in[i].x-in[k].x*in[k].x+
		      in[i].y*in[i].y-in[k].y*in[k].y+
		      in[i].z*in[i].z-in[k].z*in[k].z);
	  mat[2][0]=2*(in[l].x-in[i].x);
	  mat[2][1]=2*(in[l].y-in[i].y);
	  mat[2][2]=2*(in[l].z-in[i].z);
	  mat[2][3]=(in[i].x*in[i].x-in[l].x*in[l].x+
		      in[i].y*in[i].y-in[l].y*in[l].y+
		      in[i].z*in[i].z-in[l].z*in[l].z);
	  Point center = gauss(3,mat);
	  double r[4]  ={
	    (in[i]-center).dist(),
	    (in[j]-center).dist(),
	    (in[k]-center).dist(),
	    (in[l]-center).dist(),
	  }; 
	  /*	  
	  cout << (in[i]-center).dist() <<" " <<
	    (in[j]-center).dist() <<" " << 
	    (in[k]-center).dist() << " " <<
	    (in[l]-center).dist() << endl;
	  */	  	  
	  //cout << center.x <<" " << center.y <<" " << center.z << endl;
	  
	  if (abs(r[0]-r[1])>eps || abs(r[0]-r[2])>eps||
	      abs(r[0]-r[3])>eps)continue;
	  //if (is_contain_all(center,r[0],in))ret=min(r[0],ret);
	}
      }
    }
  }


  rep(i,in.size()){
    REP(j,i+1,in.size()){
      REP(k,j+1,in.size()){
	if (isp(in[i],in[j],in[k])||isp(in[i],in[k],in[j])
	    ||isp(in[j],in[k],in[i])){
	  continue;
	}
	Point cir = CircumCenter(in[i],in[j],in[k]);
	double r =(cir-in[i]).dist();
	/*
	cout << (cir-in[i]).dist() <<" " 
	     << (cir-in[j]).dist() <<" " 
	     << (cir-in[k]).dist() << endl;
	*/
	if (is_contain_all(cir,r,in))ret=min(ret,r);
	
	Point norm= normal_vector(in[j]-in[i],in[k]-in[i]);
	Point tmp = (in[i]+in[j])/2-in[i];
	Point tmp2= (cir-(in[i]+in[j])/2);
	ret=min(ret,solve(cir,norm,in,in[i],tmp.dist(),tmp2.dist()));
	norm=norm*-1;
	ret=min(ret,solve(cir,norm,in,in[i],tmp.dist(),tmp2.dist()));
	
      }
    }
  }
  
  rep(i,in.size()){
    REP(j,i+1,in.size()){
      Point center=(in[i]+in[j])/2;
      double r = (center-in[i]).dist();
      if (is_contain_all(center,r,in))ret=min(ret,r);
    }
  }

  return ret;
}

int main(){
  int n;
  while(cin>>n && n){
    vector<Point> in;
    rep(i,n){
      double x,y,z;
      cin>>x>>y>>z;
      in.push_back((Point){x,y,z});
    }

    printf("%.8lf\n",solve(in));
  }
}
```