#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;
}