fork download
copy 
  1. #include <bits/stdc++.h>
  2.  
  3. using namespace std;
  4.  
  5. #define mp make_pair
  6. #define pb push_back
  7. #define sz(s) ((int) ((s).size()))
  8.  
  9. #ifdef DEBUG
  10. #define eprintf(...) fprintf(stderr, __VA_ARGS__), fflush(stderr)
  11. #else
  12. #define eprintf(...) ;
  13. #endif
  14.  
  15. typedef long double ld;
  16. typedef long long ll;
  17.  
  18. const ld eps = 1e-11;
  19. const ld pi = acos((ld) -1.0);
  20. const ld min_x_center = 0;
  21. const ld max_x_center = 100;
  22. const ld max_r = 100;
  23. const ld min_x_point = min_x_center - max_r;
  24. const ld max_x_point = max_x_center + max_r;
  25.  
  26. mt19937 mrand(0);
  27.  
  28. int rnd(int x) {
  29.   ll val = ((ll) 1 << 32) % x;
  30.   while (true) {
  31.     ll cur = mrand();
  32.     if (cur >= ((ll) 1 << 32) - val) {
  33.       continue;
  34.     }
  35.     return cur % x;
  36.   }
  37.   assert(0);
  38. }
  39.  
  40. ld rnd_double() {
  41.   static const ld to_divide = pow(2, 32);
  42.   ld res = 0;
  43.   for (int i = 0; i < 3; i++) {
  44.     res += mrand();
  45.     res /= to_divide;
  46.   }
  47.   return res;
  48. }
  49.  
  50. ld sqr(ld a) {
  51.   return a * a;
  52. }
  53.  
  54. struct point {
  55.   ld x, y, z;
  56.  
  57.   point() : x(0), y(0), z(0) {}
  58.  
  59.   point(ld x, ld y, ld z) : x(x), y(y), z(z) {}
  60.  
  61.   point operator+(const point & b) const {
  62.     return point(x + b.x, y + b.y, z + b.z);
  63.   }
  64.  
  65.   point operator-(const point & b) const {
  66.     return point(x - b.x, y - b.y, z - b.z);
  67.   }
  68.  
  69.   ld operator*(const point & b) const {
  70.     return x * b.x + y * b.y + z * b.z;
  71.   }
  72.  
  73.   point operator-() const{
  74.     return point(-x, -y, -z);
  75.   }
  76.  
  77.   point operator*(const ld & b) const {
  78.     return point(x * b, y * b, z * b);
  79.   }
  80.  
  81.   point operator/(const ld & b) const {
  82.     return point(x / b, y / b, z / b);
  83.   }
  84.  
  85.   point operator^(const point & b) const {
  86.     return point(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);
  87.   }
  88.  
  89.   ld slen() const {
  90.     return x * x + y * y + z * z;
  91.   }
  92.  
  93.   ld len() const {
  94.     return sqrt(slen());
  95.   }
  96.  
  97.   point norm() const {
  98.     ld l = len();
  99.     return (*this) / l;
  100.   }
  101. };
  102.  
  103. point rnd_point(ld mn, ld mx) {
  104.   return point(mn + (mx - mn) * rnd_double(),
  105.                mn + (mx - mn) * rnd_double(), 
  106.                mn + (mx - mn) * rnd_double());
  107. }
  108.  
  109. ld get_angle(point a, point b) {
  110.   return atan2((a ^ b).len(), a * b);
  111. }
  112.  
  113. //solves vs[i] * x = cs[i], x * x = 1
  114. vector<point> solve_eq(vector<point> vs, vector<ld> cs) {
  115.   point w = vs[0] ^ vs[1];
  116.   ld det = w.slen();
  117.   if (det < eps) {
  118.     return {};
  119.   }
  120.   ld x = (cs[0] * vs[1].slen() - cs[1] * (vs[0] * vs[1])) / det;
  121.   ld y = (cs[1] * vs[0].slen() - cs[0] * (vs[0] * vs[1])) / det;
  122.   point o = vs[0] * x + vs[1] * y;
  123.   if (o.slen() > 1 - eps) {
  124.     return {};
  125.   }
  126.   ld toadd = sqrt(max((ld) 0, 1 - o.slen()));
  127.   w = w.norm();
  128.   /*assert(abs((o + w * toadd) * vs[0] - cs[0]) < eps);
  129.   assert(abs((o - w * toadd) * vs[0] - cs[0]) < eps);
  130.   assert(abs((o - w * toadd) * (o - w * toadd) - 1) < eps);
  131.   assert(abs((o + w * toadd) * (o + w * toadd) - 1) < eps);*/
  132.   return {o + w * toadd, o - w * toadd};
  133. }
  134.  
  135. bool is_inside(const vector<point> & cs, const vector<ld> & rs, const point & cur) {
  136.   for (int j = 0; j < 3; j++) {
  137.     if ((cur - cs[j]).slen() < sqr(rs[j])) {
  138.       return true;
  139.     }
  140.   }
  141.   vector<point> vs;
  142.   vector<point> tocheck;
  143.   for (int j = 0; j < 3; j++) {
  144.     vs.pb(cs[j] - cur);
  145.     tocheck.pb(vs.back().norm());
  146.   }
  147.   for (int j = 0; j < 3; j++) {
  148.     for (int k = j + 1; k < 3; k++) {
  149.       vector<point> cur = solve_eq({vs[j], vs[k]}, {rs[j], rs[k]});
  150.       for (auto r : cur) {
  151.         tocheck.pb(r);
  152.         if (j != 0 || k != 1) {
  153.           break;
  154.         }
  155.       }
  156.     }
  157.   }
  158.   for (auto r : tocheck) {
  159.     bool ok = true;
  160.     for (int j = 0; j < 3; j++) {
  161.       if (r * vs[j] < rs[j] - eps) {
  162.         ok = false;
  163.         break;
  164.       }
  165.     }
  166.     if (ok) {
  167.       return false;
  168.     }
  169.   }
  170.   return true;
  171. }
  172.  
  173. ld solve(const vector<point> & cs, const vector<ld> & rs) {
  174.   ld res = 0;
  175.   vector<point> centers; //centers of cubes
  176.   centers.pb(point((max_x_point + min_x_point) / 2,
  177.                    (max_x_point + min_x_point) / 2,
  178.                    (max_x_point + min_x_point) / 2));
  179.   ld cura = (max_x_point - min_x_point);
  180.   ld currad = (max_x_point - min_x_point) / 2 * sqrt((ld) 3);
  181.   while (sz(centers) < 5e6) {
  182.     vector<point> ncenters;
  183.     ld ncura = cura / 2;
  184.     ld ncurrad = currad / 2;
  185.     for (auto p : centers) {
  186.       for (int i = -1; i <= 1; i += 2) {
  187.         for (int j = -1; j <= 1; j += 2) {
  188.           for (int k = -1; k <= 1; k += 2) {
  189.             ncenters.pb(point(p.x + i * ncura / 2, p.y + j * ncura / 2, p.z + k * ncura / 2));
  190.           }
  191.         }
  192.       }
  193.     }
  194.     centers.clear();
  195.     vector<ld> nrs;
  196.     for (auto r : rs) {
  197.       nrs.pb(r + ncurrad);
  198.     }
  199.     ld curvol = pow(ncura, 3);
  200.     for (auto p : ncenters) {
  201.       if (is_inside(cs, nrs, p)) {
  202.         bool ok = true;
  203.         for (int i = -1; i <= 1 && ok; i += 2) {
  204.           for (int j = -1; j <= 1 && ok; j += 2) {
  205.             for (int k = -1; k <= 1 && ok; k += 2) {
  206.               ok &= is_inside(cs, rs, point(p.x + i * ncura / 2, p.y + j * ncura / 2, p.z + k * ncura / 2));
  207.             }
  208.           }
  209.         }
  210.         if (ok) {
  211.           res += curvol;
  212.           continue;
  213.         }
  214.         centers.pb(p);
  215.       }
  216.     }
  217.     cura = ncura;
  218.     currad = ncurrad;
  219.     eprintf("remains = %.3f\n", (double) (pow(cura, 3) * sz(centers)));
  220.     eprintf("current result = %.3f\n", (double) res);
  221.     eprintf("number of cubes = %d\n", sz(centers));
  222.   }
  223.   eprintf("!\n");
  224.   static const int steps = 1e7;
  225.   int toadd = 0;
  226.   for (int i = 0; i < steps; i++) {
  227.     int j = rnd(sz(centers));
  228.     point cur = centers[j] + point(-cura / 2 + cura * rnd_double(),
  229.                                    -cura / 2 + cura * rnd_double(),
  230.                                    -cura / 2 + cura * rnd_double());
  231.     if (is_inside(cs, rs, cur)) {
  232.       toadd++;
  233.     }
  234.   }
  235.   return res + pow(cura, 3) * sz(centers) * toadd / steps;
  236. }
  237.  
  238. vector<point> cs;
  239. vector<ld> rs;
  240.  
  241. bool read() {
  242.   cs.clear();
  243.   rs.clear();
  244.   for (int i = 0; i < 3; i++) {
  245.     int rad, x, y, z;
  246.     if (scanf("%d%d%d%d", &x, &y, &z, &rad) < 4) {
  247.       return false;
  248.     }
  249.     cs.pb(point(x, y, z));
  250.     rs.pb(rad);
  251.   }
  252.   return true;
  253. }
  254.  
  255. void solve() {
  256.   printf("%.18f\n", (double) solve(cs, rs));
  257. }
  258.  
  259. int main() {
  260.   while (read()) {
  261.     solve();
  262.   }
  263.   return 0;
  264. }
  265.  
Success #stdin #stdout 0s 15248KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/ozbLzf
language:
C++14 (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!