#define _USE_MATH_DEFINES #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; template class Operators { public: template const T1 operator+(const T2& right) const{ T1 ans = static_cast( *this ); ans += right; return ans; } template const T1 operator-(const T2& right) const{ T1 ans = static_cast( *this ); ans -= right; return ans; } template const T1 operator*(const T2& right) const{ T1 ans = static_cast( *this ); ans *= right; return ans; } template const T1 operator/(const T2& right) const{ T1 ans = static_cast( *this ); ans /= right; return ans; } bool operator!=(const T1& right) const{ const T1& left = static_cast( *this ); return !(left == right); } bool operator>(const T1& right) const{ const T1& left = static_cast( *this ); return right < left; } bool operator<=(const T1& right) const{ const T1& left = static_cast( *this ); return !(right < left); } bool operator>=(const T1& right) const{ const T1& left = static_cast( *this ); return !(left < right); } }; class Point : public Operators { public: double x, y, z; Point(){ x = y = z = 0.0; } Point(double x0, double y0, double z0){ x = x0; y = y0; z = z0; } Point& operator+=(const Point& p){ x += p.x; y += p.y; z += p.z; return *this; } Point& operator-=(const Point& p){ x -= p.x; y -= p.y; z -= p.z; return *this; } Point& operator*=(double a){ x *= a; y *= a; z *= a; return *this; } Point& operator/=(double a){ x /= a; y /= a; z /= a; return *this; } double length() const{ return sqrt(x * x + y * y + z * z); } double dist(const Point& p) const{ return sqrt(pow(x - p.x, 2) + pow(y - p.y, 2) + pow(z - p.z, 2)); } }; // 行列式を計算する double determinant(vector > mat) { const double EPS = 1.0e-10; const int n = mat.size(); double ret = 1.0; for(int i=0; i abs(mat[tmp][i])) tmp = j; } if(abs(mat[tmp][i]) < EPS) return 0.0; if(tmp != i){ mat[i].swap(mat[tmp]); ret *= -1.0; } for(int j=i+1; j=i; --k){ mat[j][k] -= mat[i][k] * (mat[j][i] / mat[i][i]); } } ret *= mat[i][i]; } return ret; } /***************************************************************************************************/ // 三角形の面積を求める(ヘロンの公式) // a, b, c : 3辺の長さ /***************************************************************************************************/ double triangleArea(double a, double b, double c) { double s = (a + b + c) / 2.0; return sqrt(s * (s - a) * (s - b) * (s - c)); } /***************************************************************************************************/ // 四面体の体積を求める // x, y, z : 4つの頂点の座標 /***************************************************************************************************/ double tetrahedronVolume(const vector& p) { vector > mat(4, vector(4, 1.0)); for(int i=0; i<4; ++i){ mat[i][0] = p[i].x; mat[i][1] = p[i].y; mat[i][2] = p[i].z; } return abs(determinant(mat)) / 6.0; } int main() { int n; Point p; cin >> n >> p.x >> p.y >> p.z; vector q(n); for(int i=0; i> q[i].x >> q[i].y >> q[i].z; double ans = 0.0; for(int i=0; i a = {q[i], q[j], q[k], p}; double volume = tetrahedronVolume(a); double area = triangleArea(q[i].dist(q[j]), q[j].dist(q[k]), q[k].dist(q[i])); ans += volume / area * 3.0; } } } printf("%.10f\n", ans); return 0; }