hdu 5733(几何)

#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;
}