#include <bits/stdc++.h>

using namespace std;

#define mp make_pair
#define pb push_back
#define sz(s) ((int) ((s).size()))

#ifdef DEBUG
#define eprintf(...) fprintf(stderr, __VA_ARGS__), fflush(stderr)
#else
#define eprintf(...) ;
#endif

typedef long double ld;
typedef long long ll;

const ld eps = 1e-11;
const ld pi = acos((ld) -1.0);
const ld min_x_center = 0;
const ld max_x_center = 100;
const ld max_r = 100;
const ld min_x_point = min_x_center - max_r;
const ld max_x_point = max_x_center + max_r;

mt19937 mrand(0);

int rnd(int x) {
  ll val = ((ll) 1 << 32) % x;
  while (true) {
    ll cur = mrand();
    if (cur >= ((ll) 1 << 32) - val) {
      continue;
    }
    return cur % x;
  }
  assert(0);
}

ld rnd_double() {
  static const ld to_divide = pow(2, 32);
  ld res = 0;
  for (int i = 0; i < 3; i++) {
    res += mrand();
    res /= to_divide;
  }
  return res;
}

ld sqr(ld a) {
  return a * a;
}

struct point {
  ld x, y, z;

  point() : x(0), y(0), z(0) {}

  point(ld x, ld y, ld z) : x(x), y(y), z(z) {}

  point operator+(const point & b) const {
    return point(x + b.x, y + b.y, z + b.z);
  }

  point operator-(const point & b) const {
    return point(x - b.x, y - b.y, z - b.z);
  }

  ld operator*(const point & b) const {
    return x * b.x + y * b.y + z * b.z;
  }

  point operator-() const{
    return point(-x, -y, -z);
  }

  point operator*(const ld & b) const {
    return point(x * b, y * b, z * b);
  }

  point operator/(const ld & b) const {
    return point(x / b, y / b, z / b);
  }

  point operator^(const point & b) const {
    return point(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);
  }

  ld slen() const {
    return x * x + y * y + z * z;
  }

  ld len() const {
    return sqrt(slen());
  }

  point norm() const {
    ld l = len();
    return (*this) / l;
  }
};

point rnd_point(ld mn, ld mx) {
  return point(mn + (mx - mn) * rnd_double(),
               mn + (mx - mn) * rnd_double(), 
               mn + (mx - mn) * rnd_double());
}

ld get_angle(point a, point b) {
  return atan2((a ^ b).len(), a * b);
}

//solves vs[i] * x = cs[i], x * x = 1
vector<point> solve_eq(vector<point> vs, vector<ld> cs) {
  point w = vs[0] ^ vs[1];
  ld det = w.slen();
  if (det < eps) {
    return {};
  }
  ld x = (cs[0] * vs[1].slen() - cs[1] * (vs[0] * vs[1])) / det;
  ld y = (cs[1] * vs[0].slen() - cs[0] * (vs[0] * vs[1])) / det;
  point o = vs[0] * x + vs[1] * y;
  if (o.slen() > 1 - eps) {
    return {};
  }
  ld toadd = sqrt(max((ld) 0, 1 - o.slen()));
  w = w.norm();
  /*assert(abs((o + w * toadd) * vs[0] - cs[0]) < eps);
  assert(abs((o - w * toadd) * vs[0] - cs[0]) < eps);
  assert(abs((o - w * toadd) * (o - w * toadd) - 1) < eps);
  assert(abs((o + w * toadd) * (o + w * toadd) - 1) < eps);*/
  return {o + w * toadd, o - w * toadd};
}

bool is_inside(const vector<point> & cs, const vector<ld> & rs, const point & cur) {
  for (int j = 0; j < 3; j++) {
    if ((cur - cs[j]).slen() < sqr(rs[j])) {
      return true;
    }
  }
  vector<point> vs;
  vector<point> tocheck;
  for (int j = 0; j < 3; j++) {
    vs.pb(cs[j] - cur);
    tocheck.pb(vs.back().norm());
  }
  for (int j = 0; j < 3; j++) {
    for (int k = j + 1; k < 3; k++) {
      vector<point> cur = solve_eq({vs[j], vs[k]}, {rs[j], rs[k]});
      for (auto r : cur) {
        tocheck.pb(r);
        if (j != 0 || k != 1) {
          break;
        }
      }
    }
  }
  for (auto r : tocheck) {
    bool ok = true;
    for (int j = 0; j < 3; j++) {
      if (r * vs[j] < rs[j] - eps) {
        ok = false;
        break;
      }
    }
    if (ok) {
      return false;
    }
  }
  return true;
}

ld solve(const vector<point> & cs, const vector<ld> & rs) {
  ld res = 0;
  vector<point> centers; //centers of cubes
  centers.pb(point((max_x_point + min_x_point) / 2,
                   (max_x_point + min_x_point) / 2,
                   (max_x_point + min_x_point) / 2));
  ld cura = (max_x_point - min_x_point);
  ld currad = (max_x_point - min_x_point) / 2 * sqrt((ld) 3);
  while (sz(centers) < 5e6) {
    vector<point> ncenters;
    ld ncura = cura / 2;
    ld ncurrad = currad / 2;
    for (auto p : centers) {
      for (int i = -1; i <= 1; i += 2) {
        for (int j = -1; j <= 1; j += 2) {
          for (int k = -1; k <= 1; k += 2) {
            ncenters.pb(point(p.x + i * ncura / 2, p.y + j * ncura / 2, p.z + k * ncura / 2));
          }
        }
      }
    }
    centers.clear();
    vector<ld> nrs;
    for (auto r : rs) {
      nrs.pb(r + ncurrad);
    }
    ld curvol = pow(ncura, 3);
    for (auto p : ncenters) {
      if (is_inside(cs, nrs, p)) {
        bool ok = true;
        for (int i = -1; i <= 1 && ok; i += 2) {
          for (int j = -1; j <= 1 && ok; j += 2) {
            for (int k = -1; k <= 1 && ok; k += 2) {
              ok &= is_inside(cs, rs, point(p.x + i * ncura / 2, p.y + j * ncura / 2, p.z + k * ncura / 2));
            }
          }
        }
        if (ok) {
          res += curvol;
          continue;
        }
        centers.pb(p);
      }
    }
    cura = ncura;
    currad = ncurrad;
    eprintf("remains = %.3f\n", (double) (pow(cura, 3) * sz(centers)));
    eprintf("current result = %.3f\n", (double) res);
    eprintf("number of cubes = %d\n", sz(centers));
  }
  eprintf("!\n");
  static const int steps = 1e7;
  int toadd = 0;
  for (int i = 0; i < steps; i++) {
    int j = rnd(sz(centers));
    point cur = centers[j] + point(-cura / 2 + cura * rnd_double(),
                                   -cura / 2 + cura * rnd_double(),
                                   -cura / 2 + cura * rnd_double());
    if (is_inside(cs, rs, cur)) {
      toadd++;
    }
  }
  return res + pow(cura, 3) * sz(centers) * toadd / steps;
}

vector<point> cs;
vector<ld> rs;

bool read() {
  cs.clear();
  rs.clear();
  for (int i = 0; i < 3; i++) {
    int rad, x, y, z;
    if (scanf("%d%d%d%d", &x, &y, &z, &rad) < 4) {
      return false;
    }
    cs.pb(point(x, y, z));
    rs.pb(rad);
  }
  return true;
}

void solve() {
  printf("%.18f\n", (double) solve(cs, rs));
}

int main() {
  while (read()) {
    solve();
  }
  return 0;
}