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.
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".
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;
}