結果
問題 | No.622 点と三角柱の内外判定 |
ユーザー |
![]() |
提出日時 | 2017-12-03 11:16:43 |
言語 | C++11 (gcc 13.3.0) |
結果 |
AC
|
実行時間 | 2 ms / 1,500 ms |
コード長 | 8,505 bytes |
コンパイル時間 | 2,498 ms |
コンパイル使用メモリ | 178,600 KB |
実行使用メモリ | 5,248 KB |
最終ジャッジ日時 | 2024-11-28 08:46:54 |
合計ジャッジ時間 | 3,868 ms |
ジャッジサーバーID (参考情報) |
judge4 / judge1 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
other | AC * 32 |
ソースコード
#include <bits/stdc++.h>using namespace std;#define pb push_back#define rep(i,n) for(long long i = 0; i < (long long)(n); i++)typedef double number;const number eps = 1e-8;using arr = vector<number>;using matrix = vector<arr>;ostream &operator<<(ostream &o, const arr &v) {for (int i = 0; i < v.size(); i++) { cout << v[i] << " "; }return o;}ostream &operator<<(ostream &o, const matrix &v) {for (int i = 0; i < v.size(); i++) { cout << v[i] << endl;}return o;}// O( n^2 )matrix zero(int n) {matrix A(n, arr(n, 0));return A;}// O( n^2 )matrix transpose(matrix a) {int n = a.size();matrix ret = zero(n);rep(i, n) rep(j, n) {ret[i][j] = a[j][i];}return ret;}// O( n^2 )matrix identity(int n) {matrix A(n, arr(n, 0));for (int i = 0; i < n; ++i) A[i][i] = 1; // 積の単位元(和の単位元は?)return A;}// O( n^2 )arr mul(const matrix &A, const arr &x) {arr y(A.size(), 0);for (int i = 0; i < A.size(); ++i)for (int j = 0; j < A[0].size(); ++j)y[i] += A[i][j] * x[j]; // 加群の積と和の演算子return y;}// O( n^3 )matrix mul(const matrix &A, const matrix &B) {matrix C(A.size(), arr(B[0].size(), 0));for (int i = 0; i < C.size(); ++i)for (int j = 0; j < C[i].size(); ++j)for (int k = 0; k < A[i].size(); ++k)C[i][j] += A[i][k] * B[k][j]; // 加群の積と和の演算子return C;}// O( n^3 log e )matrix pow(const matrix &A, int e) {return e == 0 ? identity(A.size()) :e % 2 == 0 ? pow(mul(A, A), e/2) : mul(A, pow(A, e-1));}// O( n )number inner_product(const arr &a, const arr &b) {number ans = 0;for (int i = 0; i < a.size(); ++i)ans += a[i] * b[i];return ans;}// O( n )arr outer_product(const arr &a, const arr &b) {arr ret(3);ret[0] = a[1] * b[2] - a[2] * b[1];ret[1] = a[2] * b[0] - a[0] * b[2];ret[2] = a[0] * b[1] - a[1] * b[0];return ret;}// O( n^3 )number det(matrix A) {const int n = A.size();number D = 1;for (int i = 0; i < n; ++i) {int pivot = i;for (int j = i+1; j < n; ++j)if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;swap(A[pivot], A[i]);D *= A[i][i] * (i != pivot ? -1 : 1);if (abs(A[i][i]) < eps) break;for(int j = i+1; j < n; ++j)for(int k = n-1; k >= i; --k)A[j][k] -= A[i][k] * A[j][i] / A[i][i];}return D;}// O(n)number tr(const matrix &A) {number ans = 0;for (int i = 0; i < A.size(); ++i)ans += A[i][i];return ans;}// O( n^3 ).int rank(matrix A) {const int n = A.size(), m = A[0].size();int r = 0;for (int i = 0; r < n && i < m; ++i) {int pivot = r;for (int j = r+1; j < n; ++j)if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;swap(A[pivot], A[r]);if (abs(A[r][i]) < eps) continue;for (int k = m-1; k >= i; --k)A[r][k] /= A[r][i];for(int j = r+1; j < n; ++j)for(int k = i; k < m; ++k)A[j][k] -= A[r][k] * A[j][i];++r;}return r;}struct LUinfo {vector<number> value;vector<int> index;};// O( n^3 ), Gaussian forward eliminationLUinfo LU_decomposition(matrix A) {const int n = A.size();LUinfo data;for (int i = 0; i < n; ++i) {int pivot = i;for (int j = i+1; j < n; ++j)if (abs(A[j][i]) > abs(A[pivot][i])) pivot = j;swap(A[pivot], A[i]);data.index.push_back(pivot);// if A[i][i] == 0, LU decomposition failed.for(int j = i+1; j < n; ++j) {A[j][i] /= A[i][i];for(int k = i+1; k < n; ++k)A[j][k] -= A[i][k] * A[j][i];data.value.push_back(A[j][i]);}}for(int i = n-1; i >= 0; --i) {for(int j = i+1; j < n; ++j)data.value.push_back(A[i][j]);data.value.push_back(A[i][i]);}return data;}// O( n^2 ) Gaussian backward substitutionarr LU_backsubstitution(const LUinfo &data, arr b) {const int n = b.size();int k = 0;for (int i = 0; i < n; ++i){swap(b[data.index[i]], b[i]);for(int j = i+1; j < n; ++j)b[j] -= b[i] * data.value[k++];}for (int i = n-1; i >= 0; --i) {for (int j = i+1; j < n; ++j)b[i] -= b[j] * data.value[k++];b[i] /= data.value[k++];}return b;}// reduce Hessenberg form (inplace).// O ( n^3 )void hessenberg(matrix &A) {const int n = A.size();for (int k = 1; k <= n-2; ++k) {arr u(n);for (int i = k; i < n; ++i) u[i] = A[i][k-1];number ss = 0;for (int i = k+1; i < n; ++i) ss += u[i] * u[i];if (abs(ss) <= 0.0) continue;number s = sqrt( ss + u[k]*u[k] );if (u[k] > 0.0) s = -s;u[k] -= s;number uu = sqrt( ss + u[k]*u[k] );for (int i = k; i < n; ++i) u[i] /= uu;arr f(n), g(n);for (int i = 0; i < n; ++i)for (int j = k; j < n; ++j)f[i] += A[i][j] * u[j],g[i] += A[j][i] * u[j];number gamma = inner_product(u, g);for (int i = 0; i < n; ++i)f[i] -= gamma * u[i],g[i] -= gamma * u[i];for (int i = 0; i < n; ++i)for (int j = 0; j < n; ++j)A[i][j] = A[i][j] - 2*u[i]*g[j] - 2*f[i]*u[j];}}// find all eigenvalues using Hessenberg-QR Method// O( n^3 + M n^2 ) where M is the number of iterations.vector<number> eigenvalues(matrix A) {const int n = A.size();hessenberg(A);vector<number> s(n), c(n);for (int m = n; m >= 2; ) {if (abs(A[m-1][m-2]) < eps) { --m; continue; }number shift = A[m-1][m-1];for (int i = 0; i < m; ++i) A[i][i] -= shift;for (int k = 0; k < m-1; ++k) {number a = A[k][k], b = A[k+1][k], r = sqrt(a*a + b*b);s[k] = r == 0.0 ? 0.0 : b/r,c[k] = r == 0.0 ? 0.0 : a/r;for (int j = k; j < m; ++j) {number x = A[k][j], y = A[k+1][j];A[ k ][j] = c[k] * x + s[k] * y;A[k+1][j] = -s[k] * x + c[k] * y;}}for (int k = 0; k < m-1; ++k) {for (int i = 0; i <= k+1; ++i) {number x = A[i][k], y = A[i][k+1];A[i][ k ] = c[k] * x + s[k] * y;A[i][k+1] = -s[k] * x + c[k] * y;}}for (int i = 0; i < m; ++i) A[i][i] += shift;}vector<number> lambda;for (int i = 0; i < n; ++i)lambda.push_back( A[i][i] );return lambda;}// find the corresponding eigenvector from the eigenvalue.// O ( n^3 + M n^2 ) where M is the number of iterations.arr eigenvector(matrix A, number lambda) {const int n = A.size();arr y(n); y[0] = 1;for (int i = 0; i < n; ++i) A[i][i] -= lambda;LUinfo data = LU_decomposition(A);number mu, v2, v2s;do {arr v = LU_backsubstitution(data, y); // A v = ymu = inner_product(v, y);v2 = inner_product(v, v);v2s = sqrt(v2);for (int j = 0; j < n; ++j) y[j] = v[j] / v2s;} while (abs(1.0-mu*mu/v2) > eps);return y;}int main(void){vector<arr> s;rep(i, 3) {arr tmp(3); cin >> tmp[0] >> tmp[1] >> tmp[2];s.pb(tmp);}vector<arr> xyz(3, arr(3));rep(i, 3) {rep(j, 3) xyz[0][j] = s[1][j] - s[0][j];rep(j, 3) xyz[1][j] = s[2][j] - s[0][j];}xyz[2] = outer_product(xyz[0], xyz[1]);xyz = transpose(xyz);arr a(3); cin >> a[0] >> a[1] >> a[2];rep(i, 3) a[i] -= s[0][i];// cout << xyz << endl;LUinfo lu = LU_decomposition(xyz);arr ret = LU_backsubstitution(lu, a);// cout << mul(xyz, ret) <<"toi3je"<< endl;// cout << a << "#ijo" << endl;// cout << ret << " #ret" << endl;arr check(3);// cout << xyz << endl;rep(i, 3) {check[i] =ret[0] * xyz[i][0] +ret[1] * xyz[i][1] +ret[2] * xyz[i][2] +s[0][i];}// cout << check << "#check"<< endl;if (ret[0] >= 0 && ret[1] >= 0 && ret[0] + ret[1] <= 1) {cout << "YES" << endl;} else {cout << "NO" << endl;}return 0;}