#include <bits/stdc++.h>
#pragma comment(linker, "/STACK:102400000,102400000")
using namespace std;
#define LL long long
#define pii pair<int,int>
#define MP make_pair
#define ls i << 1
#define rs ls | 1
#define md (ll + rr >> 1)
#define lson ll, md, ls
#define rson md + 1, rr, rs
#define Pi acos(-1.0)
#define mod 1000000007
#define eps 1e-12
#define inf 0x3f3f3f3f
#define N 200010
#define M 2000020
int dcmp(double x){
if(fabs(x) < eps) return 0;
return x < 0 ? -1 : 1;
}
struct point{
double x, y, z;
point(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
point operator - (const point &b) const {
return point(x - b.x, y - b.y, z - b.z);
}
point operator + (const point &b) const {
return point(x + b.x, y + b.y, z + b.z);
}
point operator * (const double &k) const {
return point(x * k, y * k, z * k);
}
point operator / (const double &k) const {
return point(x / k, y / k, z / k);
}
void input(){
scanf("%lf%lf%lf", &x, &y, &z);
}
double len(){
return sqrt(x * x + y * y + z * z);
}
void output(){
printf("%.9f %.9f %.9f\n", x, y, z);
}
};
double dot(point a, point b){
return a.x * b.x + a.y * b.y + a.z * b.z;
}
point cross(point a, point b){
return 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 area(point a, point b, point c){
return cross(b - a, c - a).len();
}
double volume(point a, point b, point c, point d){
point tmp = cross(b - a, c - a);
return dot(d - a, tmp);
}
point A, B, C, D;
double a[20][20], x[20];
void Gauss(){
int i, j, k, col, max_r;
for(k = 0, col = 0; k < 3; ++k, col++){
max_r = k;
for(i = k + 1; i < 3; ++i){
if(fabs(a[i][col]) > fabs(a[max_r][col]))
max_r = i;
}
if(k != max_r){
for(j = col; j < 3; ++j)
swap(a[k][j], a[max_r][j]);
swap(x[k], x[max_r]);
}
x[k] /= a[k][col];
for(j = col + 1; j < 3; ++j) a[k][j] /= a[k][col];
a[k][col] = 1;
for(i = 0; i < 3; ++i){
if(i != k){
x[i] -= x[k] * a[i][k];
for(j = col + 1; j < 3; ++j) a[i][j] -= a[k][j] * a[i][col];
a[i][col] = 0;
}
}
}
}
int main(){
while(scanf("%lf%lf%lf", &A.x, &A.y, &A.z) != EOF){
B.input();
C.input();
D.input();
double V = fabs(volume(A, B, C, D));
if(dcmp(V) == 0){
puts("O O O O"); continue;
}
double r = V / (area(A, B, C) + area(B, C, D) + area(A, C, D) + area(A, B, D));
point p1 = cross(B - A, C - A);
point p2 = cross(D - A, B - A);
point p3 = cross(C - A, D - A);
point p4 = cross(D - B, C - B);
if(dcmp(dot(p1, D - A)) < 0) p1 = p1 * -1;
if(dcmp(dot(p2, C - A)) < 0) p2 = p2 * -1;
if(dcmp(dot(p3, B - A)) < 0) p3 = p3 * -1;
if(dcmp(dot(p4, A - B)) < 0) p4 = p4 * -1;
point d1 = p1 / p1.len() * r;
point d2 = p2 / p2.len() * r;
point d3 = p3 / p3.len() * r;
point d4 = p4 / p4.len() * r;
a[0][0] = p1.x, a[0][1] = p1.y, a[0][2] = p1.z;
a[1][0] = p2.x, a[1][1] = p2.y, a[1][2] = p2.z;
a[2][0] = p3.x, a[2][1] = p3.y, a[2][2] = p3.z;
x[0] = p1.x * (A.x + d1.x) + p1.y * (A.y + d1.y) + p1.z * (A.z + d1.z);
x[1] = p2.x * (A.x + d2.x) + p2.y * (A.y + d2.y) + p2.z * (A.z + d2.z);
x[2] = p3.x * (A.x + d3.x) + p3.y * (A.y + d3.y) + p3.z * (A.z + d3.z);
Gauss();
printf("%.4f %.4f %.4f %.4f\n", x[0], x[1], x[2], r);
/*
a[0][0] = p1.x, a[0][1] = p1.y, a[0][2] = p1.z;
a[1][0] = p2.x, a[1][1] = p2.y, a[1][2] = p2.z;
a[2][0] = p4.x, a[2][1] = p4.y, a[2][2] = p4.z;
x[0] = p1.x * (B.x + d1.x) + p1.y * (B.y + d1.y) + p1.z * (B.z + d1.z);
x[1] = p2.x * (B.x + d2.x) + p2.y * (B.y + d2.y) + p2.z * (B.z + d2.z);
x[2] = p4.x * (B.x + d4.x) + p4.y * (B.y + d4.y) + p4.z * (B.z + d4.z);
Gauss();
printf("%.4f %.4f %.4f %.4f\n", x[0], x[1], x[2], r);
a[0][0] = p1.x, a[0][1] = p1.y, a[0][2] = p1.z;
a[1][0] = p3.x, a[1][1] = p3.y, a[1][2] = p3.z;
a[2][0] = p4.x, a[2][1] = p4.y, a[2][2] = p4.z;
x[0] = p1.x * (C.x + d1.x) + p1.y * (C.y + d1.y) + p1.z * (C.z + d1.z);
x[1] = p3.x * (C.x + d3.x) + p3.y * (C.y + d3.y) + p3.z * (C.z + d3.z);
x[2] = p4.x * (C.x + d4.x) + p4.y * (C.y + d4.y) + p4.z * (C.z + d4.z);
Gauss();
printf("%.4f %.4f %.4f %.4f\n", x[0], x[1], x[2], r);
*/
}
return 0;
}
hdu 5733(几何)
最新推荐文章于 2020-03-17 17:07:52 发布