fork(1) 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.  
  21. ld sqr(ld a) {
  22.   return a * a;
  23. }
  24.  
  25. ld cube(ld a) {
  26.   return a * a * a;
  27. }
  28.  
  29. struct point {
  30.   ld x, y, z;
  31.  
  32.   point() : x(0), y(0), z(0) {}
  33.  
  34.   point(ld x, ld y, ld z) : x(x), y(y), z(z) {}
  35.  
  36.   point operator+(const point & b) const {
  37.     return point(x + b.x, y + b.y, z + b.z);
  38.   }
  39.  
  40.   point operator-(const point & b) const {
  41.     return point(x - b.x, y - b.y, z - b.z);
  42.   }
  43.  
  44.   ld operator*(const point & b) const {
  45.     return x * b.x + y * b.y + z * b.z;
  46.   }
  47.  
  48.   point operator-() const{
  49.     return point(-x, -y, -z);
  50.   }
  51.  
  52.   point operator*(const ld & b) const {
  53.     return point(x * b, y * b, z * b);
  54.   }
  55.  
  56.   point operator/(const ld & b) const {
  57.     return point(x / b, y / b, z / b);
  58.   }
  59.  
  60.   point operator^(const point & b) const {
  61.     return point(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);
  62.   }
  63.  
  64.   ld slen() const {
  65.     return x * x + y * y + z * z;
  66.   }
  67.  
  68.   ld len() const {
  69.     return sqrt(slen());
  70.   }
  71.  
  72.   point norm() const {
  73.     ld l = len();
  74.     return (*this) / l;
  75.   }
  76. };
  77.  
  78. ld get_angle(point a, point b) {
  79.   return atan2((a ^ b).len(), a * b);
  80. }
  81.  
  82. //solves vs[i] * x = cs[i], x * x = 1
  83. vector<point> solve_eq(vector<point> vs, vector<ld> cs) {
  84.   point w = vs[0] ^ vs[1];
  85.   ld det = w.slen();
  86.   if (det < eps) {
  87.     return {};
  88.   }
  89.   ld x = (cs[0] * vs[1].slen() - cs[1] * (vs[0] * vs[1])) / det;
  90.   ld y = (cs[1] * vs[0].slen() - cs[0] * (vs[0] * vs[1])) / det;
  91.   point o = vs[0] * x + vs[1] * y;
  92.   if (o.slen() > 1 - eps) {
  93.     return {};
  94.   }
  95.   ld toadd = sqrt(max((ld) 0, 1 - o.slen()));
  96.   w = w.norm();
  97.   return {o + w * toadd, o - w * toadd};
  98. }
  99.  
  100. //solves if there is no tangent plane
  101. ld solve_easy(vector<point> cs, vector<ld> rs) {
  102.   ld cos[2], val[2], toadd[2];
  103.   point w[2];
  104.   for (int i = 0; i < 2; i++) {
  105.     cos[i] = (rs[2] - rs[i]) / (cs[i] - cs[2]).len();
  106.     w[i] = (cs[i] - cs[2]).norm();
  107.     val[i] = (1 - cos[i]) * 2 * pi;
  108.     toadd[i] = val[i] * (-cube(rs[2]) + cube(rs[i])) / 3 +
  109.                (1 - sqr(cos[i])) * (cs[i] - cs[2]).len() * (sqr(rs[2]) + rs[2] * rs[i] + sqr(rs[i])) / 3 * pi;
  110.   }
  111.   ld res = 4 * pi / 3 * cube(rs[2]);
  112.   if (max(acos(cos[0]), acos(cos[1])) > get_angle(w[0], w[1])) {
  113.     if (cos[0] < cos[1]) {
  114.       res += toadd[0];
  115.     } else {
  116.       res += toadd[1];
  117.       ld c = (cs[1] - cs[2]).slen() + rs[1] * rs[2];
  118.       if ((cs[0] - cs[2]).slen() > c) {
  119.         ld ncos = (rs[1] - rs[0]) / (cs[0] - cs[1]).len();
  120.         ld nval = (1 - ncos) * 2 * pi;
  121.         ld ntoadd = nval * (-cube(rs[1]) + cube(rs[0])) / 3 +
  122.                    (1 - sqr(ncos)) * (cs[1] - cs[0]).len() * (sqr(rs[1]) + rs[1] * rs[0] + sqr(rs[0])) / 3 * pi;
  123.         res += ntoadd;
  124.       }
  125.     }
  126.   } else {
  127.     res += toadd[0] + toadd[1];
  128.   }
  129.   return res;
  130. }
  131.  
  132. ld solve(vector<point> cs, vector<ld> rs) {
  133.   //normals of tangent planes
  134.   point p1, p2;
  135.   {
  136.     //find normals
  137.     vector<point> ans = solve_eq({cs[1] - cs[0], cs[2] - cs[0]}, {rs[0] - rs[1], rs[0] - rs[2]});
  138.     if (sz(ans) == 0) {
  139.       return solve_easy(cs, rs);
  140.     } else {
  141.       p1 = ans[0], p2 = ans[1];
  142.     }
  143.   }
  144.   ld vol = ((cs[1] - cs[0]) ^ (cs[2] - cs[0])) * p1;
  145.   if (vol < 0) {
  146.     vol *= -1;
  147.     swap(p1, p2);
  148.   }
  149.   ld res = (rs[0] + rs[1] + rs[2]) / 3 * vol; //volume of "prism"
  150.   //in the next loop we calculate volume of 3 cones and 3 parts of balls
  151.   vector<ld> vols;
  152.   for (int i = 0; i < 3; i++) {
  153.     point v = cs[1] - cs[0];
  154.     point w = (-v).norm();
  155.     point a1 = p1 ^ w, a2 = p2 ^ w;
  156.     ld alpha = get_angle(a1, a2); //"angle" of the part of the cone
  157.     if ((p1 ^ w) * p2 < 0) {
  158.       alpha = 2 * pi - alpha;
  159.     }
  160.     //volume of the part of the cone
  161.     res += (sqr(rs[0]) + rs[0] * rs[1] + sqr(rs[1])) / 6 * v.len() * sqr((w ^ p1).len()) * alpha;
  162.     {
  163.       //now we calculate volume of the part of the ball
  164.       ld e = (rs[1] - rs[0]) / v.len();
  165.       ld cur = (1 - e) * alpha;
  166.       ld gamma = get_angle(w ^ p1, p2 ^ p1);
  167.       if ((p1 ^ p2) * w > 0) {
  168.         gamma *= -1;
  169.       }
  170.       cur -= alpha + 2 * gamma - pi;
  171.       vols.pb(cur);
  172.     }
  173.     rotate(cs.begin(), cs.begin() + 1, cs.end());
  174.     rotate(rs.begin(), rs.begin() + 1, rs.end());
  175.   }
  176.   ld sum = 0;
  177.   for (int i = 0; i < 3; i++) {
  178.     ld val = vols[i] - vols[(i + 2) % 3];
  179.     while (val < 0) {
  180.       val += 4 * pi;
  181.     }
  182.     while (val > 4 * pi) {
  183.       val -= 4 * pi;
  184.     }
  185.     sum += val;
  186.     res += rs[i] * rs[i] * rs[i] * val / 3;
  187.   }
  188.   return res;
  189. }
  190.  
  191. vector<point> cs;
  192. vector<ld> rs;
  193.  
  194. bool read() {
  195.   cs.clear();
  196.   rs.clear();
  197.   for (int i = 0; i < 3; i++) {
  198.     int rad, x, y, z;
  199.     if (scanf("%d%d%d%d", &x, &y, &z, &rad) < 4) {
  200.       return false;
  201.     }
  202.     cs.pb(point(x, y, z));
  203.     rs.pb(rad);
  204.   }
  205.   return true;
  206. }
  207.  
  208. void solve() {
  209.   printf("%.18f\n", (double) solve(cs, rs));
  210. }
  211.  
  212. int main() {
  213.   while (read()) {
  214.     solve();
  215.   }
  216.   return 0;
  217. }
  218.  
Success #stdin #stdout 0s 15248KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/jvLDAR
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!