2069.txt

北大ACM题目例程 详细的解答过程 程序实现 算法分析

#include"iostream.h"
#include"stdio.h"
#include"math.h"
#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("%.5lf\n",answer);
	}
	return 0;
}