HDU 5733 tetrahedron（空间计算几何）

Problem Description

Given four points ABCD, if ABCD is a tetrahedron, calculate the inscribed sphere of ABCD.

 

Input

Multiple test cases (test cases  ≤100 ).

Each test cases contains a line of 12 integers  [−1e6,1e6]  indicate the coordinates of four vertices of ABCD.

Input ends by EOF.

 

Output

Print the coordinate of the center of the sphere and the radius, rounded to 4 decimal places.

If there is no such sphere, output "O O O O".

 

Sample Input

  
  

   
   0 0 0 2 0 0 0 0 2 0 2 0 
0 0 0 2 0 0 3 0 0 4 0 0 
  
  

 

Sample Output

  
  

   
   0.4226 0.4226 0.4226 0.4226
O O O O
  
  

 

2016多校赛第一场最后一题，据说是上一年的陈题，当时看模板有一个求外接圆的圆心和半径，就忙呼呼地照着敲上了，敲完之后发现样例不过，然后就gg了。。。发现是求内切圆，不会求于是就在傻乎乎地跑随机化算法，结果到结束都没有写完。赛后看题解说是公式（神公式。。。），几何功底实在是太差了，看了好些时间都没看明白这个公式啥意思，只能照着公式写了一发代码（这是我见过的最难代码化的公式 = =），以下是题解截图。。。

#include <bits/stdc++.h>

using namespace std;

const double eps = 1e-8;

struct Point3
{
	double x,y,z;
	Point3(){}
	Point3(double a,double b,double c):x(a),y(b),z(c){}
	Point3 operator - (const Point3 a)const{
		return Point3(x - a.x,y - a.y,z - a.z);
	}
	Point3 operator + (const Point3 a)const{
		return Point3(x + a.x,y + a.y,z + a.z);
	}
	Point3 operator * (const Point3 a)const{
		return Point3(	y * a.z - z * a.y,
						z * a.x - x * a.z,
						x * a.y - y * a.x);
	}
	Point3 operator * (const double a)const{
		return Point3(x * a,y * a,z * a);
	}
	double operator ~ (){
		return sqrt(x * x + y * y + z * z);
	}
}p[4];

struct Determinant
{
    double *arr;
    int n;
    Determinant(int a,double *tp):n(a){
		arr = new double[a * a];
		memcpy(arr,tp,sizeof(double) * a * a);
    }
	bool flag;
	int Dtransform(double D[][100],int dia,int n)
	{
		int r;
		for(r = dia + 1 ; r < n ; r++)
		{
			if(D[r][dia] != 0)
			{
				for(int c = dia ; c < n ; c++)
					swap(D[r][c],D[dia][c]);
				flag = !flag;
				break;
			}
		}
		return r;
	}
	double Dcalcular(double D[][100],int n)
	{
		int start;
		double temp,result = 1;
		for(int dia = 0 ; dia < n ; dia++)
		{
			if(fabs(D[dia][dia]) < eps)
				start = Dtransform(D,dia,n);
			else
				start = dia;
			if(start == n)
				continue;
			for(int r = start + 1 ; r < n ; r++)
			{
				if(fabs(D[r][dia]) < eps)
					continue;
				temp = D[r][dia] / D[dia][dia];
				for(int c = dia ; c < n ; c++)
					D[r][c] -= D[dia][c] * temp;
			}
		}
		for(int dia = 0 ; dia < n ; dia++)
		{
			if(fabs(D[dia][dia]) < eps)
				return 0;
			result *= D[dia][dia];
		}
		if(flag)
			result = -result;
		return result;
	}
	double value()
	{
		double tmp[100][100];
		for(int i = 0 ; i < n ; i++)
			for(int j = 0 ; j < n ; j++)
				tmp[i][j] = arr[n * i + j];
		return Dcalcular(tmp,n);
	}
};

int main()
{
	while(~scanf("%lf %lf %lf",&p[0].x,&p[0].y,&p[0].z))
	{
		for(int i = 1 ; i < 4 ; i++)
			scanf("%lf %lf %lf",&p[i].x,&p[i].y,&p[i].z);
		Point3 a,b,c,d;
		a = p[0],b = p[1],c = p[2],d = p[3];
		Point3 tp  =  d * ~((b - a) * (c - a))
					+ c * ~((b - a) * (d - a))
					+ b * ~((c - a) * (d - a))
					+ a * ~((c - b) * (d - b));
		double xx  =  ~((b - a) * (c - a))
					+ ~((b - a) * (d - a))
					+ ~((c - a) * (d - a))
					+ ~((c - b) * (d - b));
		tp = tp * (1.0 / xx);
		Determinant dt(4,new double[16]{a.x,a.y,a.z,1,
										b.x,b.y,b.z,1,
										c.x,c.y,c.z,1,
										d.x,d.y,d.z,1});
		double vol = dt.value();
		if(fabs(xx) < eps || fabs(vol) < eps)
		{
			printf("O O O O\n");
			continue;
		}
		printf("%.4f %.4f %.4f %.4f\n",tp.x,tp.y,tp.z,fabs(vol) / xx);
	}

    return 0;
}