Advertisement
Guest User

Untitled

a guest
Jul 18th, 2019
88
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 7.83 KB | None | 0 0
  1. #include <bits/stdc++.h>
  2.  
  3. #define ll long long
  4. #define db long double
  5. #define mp make_pair
  6. #define pb push_back
  7.  
  8. using namespace std;
  9.  
  10. const db eps = 1e-9;
  11.  
  12. template<class Type>
  13. struct pt{
  14.     Type x, y, z;
  15.     pt() {}
  16.     pt<Type> (Type x, Type y, Type z): x(x), y(y), z(z) {}
  17.     pt<Type> operator- (const pt<Type> &nxt) const { return pt<Type>(x - nxt.x, y - nxt.y, z - nxt.z); }
  18.     Type operator* (const pt<Type> &nxt) const { return x * nxt.x + y * nxt.y + z * nxt.z; }
  19.     pt<Type> operator% (const pt<Type> &nxt) const { return pt<Type>(y * nxt.z - z * nxt.y, z * nxt.x - x * nxt.z, x * nxt.y - y * nxt.x); }
  20.     operator pt<db>() const { return pt<db>(x, y, z); }
  21.     db len() const {
  22.         return sqrt(x * x + y * y + z * z);
  23.     }
  24.     pt<db> resize(db new_len){
  25.         if (len() < eps) return pt<db>(0, 0, 0);
  26.         db cur_len = len();
  27.         new_len /= cur_len;
  28.         return pt<db>(x * new_len, y * new_len, z * new_len);
  29.     }
  30. };
  31.  
  32. template<class Type>
  33. struct plane{
  34.     pt<Type> a, b, c;
  35.     plane() {}
  36.     plane(pt<Type>& a, pt<Type>& b, pt<Type>& c): a(a), b(b), c(c) {}
  37. };
  38.  
  39. template<int SZ>
  40. struct matrix{
  41.     int a[SZ][SZ];
  42.  
  43.     ll det(){
  44.         vector<int> t(SZ);
  45.         for (int i = 0; i < SZ; i++) t[i] = i;
  46.         ll ans = 0;
  47.         do {
  48.             ll now = 1;
  49.             for (int i = 0; i < SZ; i++) now *= a[i][t[i]];
  50.             for (int i = 0; i < SZ; i++) for (int j = 0; j < i; j++) if (t[i] < t[j]) now *= -1;
  51.             ans += now;
  52.         } while(next_permutation(t.begin(), t.end()));
  53.         return ans;
  54.     }
  55. };
  56.  
  57. ll calcDirectedVolume(pt<ll>& a, pt<ll>& b, pt<ll>& c, pt<ll>& d){
  58.     pt<ll> w[3] = {b - a, c - a, d - a};
  59.     matrix<3> m;
  60.     for (int i = 0; i < 3; i++){
  61.         m.a[i][0] = w[i].x;
  62.         m.a[i][1] = w[i].y;
  63.         m.a[i][2] = w[i].z;
  64.     }
  65.     return m.det();
  66. }
  67.  
  68. vector<plane<ll>> slowBuild3DConvexHull(vector<pt<ll>> &a){
  69.     vector<plane<ll>> ans;
  70.     for (int i = 0; i < a.size(); i++) for (int j = i + 1; j < a.size(); j++) for (int k = j + 1; k < a.size(); k++){
  71.         int w[3] = {0, 0, 0};
  72.         for (int s = 0; s < a.size(); s++){
  73.             ll val = calcDirectedVolume(a[i], a[j], a[k], a[s]);
  74.             if (val > 0) val = 2;
  75.             else if (val == 0) val = 1;
  76.             else val = 0;
  77.             w[val]++;
  78.         }
  79.         if ((w[0] > 0) + (int)(w[2] > 0) <= 1) ans.push_back(plane<ll>(a[i], a[j], a[k]));
  80.     }
  81.     return ans;
  82. }
  83.  
  84. bool wasEdge[1001][1001];
  85.  
  86. vector<plane<ll>> build3DConvexHull(vector<pt<ll>> &a){
  87.     vector<tuple<int, int, int>> pl;
  88.    
  89.     for (int i = 0; i < 4; i++) for (int j = i + 1; j < 4; j++) for (int k = j + 1; k < 4; k++){
  90.         int last = (0 ^ 1 ^ 2 ^ 3) ^ (i ^ j ^ k);
  91.         if (calcDirectedVolume(a[i], a[j], a[k], a[last]) > 0){
  92.             pl.push_back(make_tuple(i, k, j));
  93.         } else {
  94.             pl.push_back(make_tuple(i, j, k));
  95.         }
  96.     }
  97.  
  98.     for (int i = 4; i < a.size(); i++){
  99.         vector<int> rem;
  100.         for (int j = (int)pl.size() - 1; j >= 0; j--){
  101.             int w[3] = { get<0>(pl[j]), get<1>(pl[j]), get<2>(pl[j]) };
  102.             if (calcDirectedVolume(a[w[0]], a[w[1]], a[w[2]], a[i]) > 0){
  103.                 rem.push_back(j);
  104.                 wasEdge[w[0]][w[1]] = 1;
  105.                 wasEdge[w[1]][w[2]] = 1;
  106.                 wasEdge[w[2]][w[0]] = 1;
  107.             }
  108.         }
  109.         if (rem.size() == 0) continue;
  110.         for (int v : rem){
  111.             int w[3] = { get<0>(pl[v]), get<1>(pl[v]), get<2>(pl[v]) };
  112.             for (int j = 0; j < 3; j++){
  113.                 int k = j == 2 ? 0 : (j + 1);
  114.                 if (wasEdge[w[j]][w[k]] + (int)wasEdge[w[k]][w[j]] == 1){
  115.                     pl.push_back(make_tuple(i, w[j], w[k]));
  116.                 }
  117.                 wasEdge[w[j]][w[k]] = 0;
  118.                 wasEdge[w[k]][w[j]] = 0;
  119.             }
  120.             swap(pl[v], pl.back());
  121.             pl.pop_back();
  122.         }
  123.     }
  124.  
  125.     vector<plane<ll>> ans;
  126.     for (const auto &c : pl) ans.push_back(plane<ll>(a[get<0>(c)], a[get<1>(c)], a[get<2>(c)]));
  127.     return ans;
  128. }
  129.  
  130. namespace SCP{ //Smallest Circle Problem
  131.     struct pt{
  132.         db x, y;
  133.         pt() {}
  134.         pt(db x, db y): x(x), y(y) {}
  135.         pt operator- (const pt &nxt) const { return pt(x - nxt.x, y - nxt.y); }
  136.         db len(){
  137.             return sqrt(x * x + y * y);
  138.         }
  139.     };
  140.  
  141.     struct line{
  142.         db a, b, c;
  143.     };
  144.  
  145.     db getSquare(db r){
  146.         return M_PI * r * r;
  147.     }
  148.  
  149.     pt getMedian(pt &a, pt &b){
  150.         return pt((a.x + b.x) / 2, (a.y + b.y) / 2);
  151.     }
  152.  
  153.     pair<pt, db> SCP(pt &a, pt &b){
  154.         return make_pair(getMedian(a, b), (a - b).len() / 2);
  155.     }
  156.  
  157.     pt intersectLines(line &l1, line &l2){
  158.         if (abs(l1.a * l2.b - l2.a * l1.b) < eps) throw 42;
  159.         db x = (l2.c * l1.b - l1.c * l2.b) / (l1.a * l2.b - l2.a * l1.b);
  160.         db y = (l2.c * l1.a - l1.c * l2.a) / (l1.b * l2.a - l2.b * l1.a);
  161.         return pt(x, y);
  162.     }
  163.  
  164.     pair<pt, db> SCP(pt &a, pt &b, pt &c){
  165.         pt o1 = getMedian(a, b);
  166.         pt o2 = getMedian(b, c);
  167.         line l1, l2;
  168.         l1.a = (b - a).x; l1.b = (b - a).y; l1.c = -(l1.a * o1.x + l1.b * o1.y);
  169.         l2.a = (b - c).x; l2.b = (b - c).y; l2.c = -(l2.a * o2.x + l2.b * o2.y);
  170.         try {
  171.             pt o = intersectLines(l1, l2);
  172.             return make_pair(o, (o - a).len());
  173.         } catch(...) {
  174.             throw;
  175.         }
  176.     }
  177.  
  178.     bool inCircle(pt &a, pt &O, db r){
  179.         return (O - a).len() <= r + eps;
  180.     }
  181.  
  182.     pair<pt, db> recSolve(vector<pt> &a, vector<pt> &b){
  183.         assert(b.size() <= 3);
  184.         if (b.size() == 3){
  185.             auto [O, r] = SCP(b[0], b[1], b[2]);
  186.             bool ok = 1;
  187.             for (auto p : a) if (!inCircle(p, O, r)){
  188.                 ok = 0;
  189.                 break;
  190.             }
  191.             if (ok) return make_pair(O, r);
  192.             else return make_pair(O, -2);
  193.         } else {
  194.             if (a.size() == 0){
  195.                 if (b.size() == 0) return make_pair(pt(0, 0), 0);
  196.                 if (b.size() == 1) return make_pair(b[0], 0);
  197.                 if (b.size() == 2) return SCP(b[0], b[1]);
  198.             } else {
  199.                 pt p = a.back(); a.pop_back();
  200.                 auto [O, r] = recSolve(a, b);
  201.                 a.push_back(p);
  202.                 if (inCircle(p, O, r)) return make_pair(O, r);
  203.                 a.pop_back(), b.push_back(p);
  204.                 auto res = recSolve(a, b);
  205.                 a.push_back(p), b.pop_back();
  206.                 return res;
  207.             }
  208.         }
  209.     }
  210.  
  211.     db solve(vector<pt> &a){
  212.         if (a.size() == 1) return 0;
  213.         random_shuffle(a.begin(), a.end());
  214.         vector<pt> b;
  215.         db ans = recSolve(a, b).second;
  216.         return getSquare(ans);
  217.     }
  218. }
  219.  
  220. db solve(const vector<pt<db>> &a, const pt<db> &O, const pt<db> &I, const pt<db> &J){
  221.     vector<SCP::pt> b;
  222.     for (int i = 0; i < a.size(); i++){
  223.         b.push_back(SCP::pt(I * (a[i] - O), J * (a[i] - O)));
  224.     }
  225.     return SCP::solve(b);
  226. }
  227.  
  228. int main(){
  229. #ifdef denisson
  230.     freopen("input.txt", "r", stdin);
  231. #endif
  232.     ios_base::sync_with_stdio(0); cin.tie(0);
  233.  
  234.     int n;
  235.     vector<pt<ll>> a;
  236.  
  237.     cin >> n;
  238.     a.resize(n);
  239.     for (int i = 0; i < n; i++) cin >> a[i].x >> a[i].y >> a[i].z;
  240.  
  241.     auto pl = build3DConvexHull(a);
  242.  
  243.     db ans = 1e18;
  244.  
  245.     for (const auto &p : pl){
  246.         vector<pt<db>> pr;
  247.         pt<db> norm = ((p.b - p.a) % (p.c - p.a)).resize(1);
  248.         db height = 0;
  249.         for (int i = 0; i < n; i++){
  250.             auto proj = norm.resize((pt<db>)(a[i] - p.a) * norm);
  251.             pr.push_back((pt<db>)a[i] - proj);
  252.             height = max(height, proj.len());
  253.         }
  254.         ans = min(ans, height * solve(pr, p.a, (p.b - p.a).resize(1), ((pt<db>)(p.b - p.a) % norm).resize(1)));
  255.     }
  256.  
  257.     cout.precision(10);
  258.     cout << fixed << ans;
  259. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement