#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);
ld sqr(ld a) {
return a * a;
}
ld cube(ld a) {
return a * 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;
}
};
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();
return {o + w * toadd, o - w * toadd};
}
//solves if there is no tangent plane
ld solve_easy(vector<point> cs, vector<ld> rs) {
ld cos[2], val[2], toadd[2];
point w[2];
for (int i = 0; i < 2; i++) {
cos[i] = (rs[2] - rs[i]) / (cs[i] - cs[2]).len();
w[i] = (cs[i] - cs[2]).norm();
val[i] = (1 - cos[i]) * 2 * pi;
toadd[i] = val[i] * (-cube(rs[2]) + cube(rs[i])) / 3 +
(1 - sqr(cos[i])) * (cs[i] - cs[2]).len() * (sqr(rs[2]) + rs[2] * rs[i] + sqr(rs[i])) / 3 * pi;
}
ld res = 4 * pi / 3 * cube(rs[2]);
if (max(acos(cos[0]), acos(cos[1])) > get_angle(w[0], w[1])) {
if (cos[0] < cos[1]) {
res += toadd[0];
} else {
res += toadd[1];
ld c = (cs[1] - cs[2]).slen() + rs[1] * rs[2];
if ((cs[0] - cs[2]).slen() > c) {
ld ncos = (rs[1] - rs[0]) / (cs[0] - cs[1]).len();
ld nval = (1 - ncos) * 2 * pi;
ld ntoadd = nval * (-cube(rs[1]) + cube(rs[0])) / 3 +
(1 - sqr(ncos)) * (cs[1] - cs[0]).len() * (sqr(rs[1]) + rs[1] * rs[0] + sqr(rs[0])) / 3 * pi;
res += ntoadd;
}
}
} else {
res += toadd[0] + toadd[1];
}
return res;
}
ld solve(vector<point> cs, vector<ld> rs) {
//normals of tangent planes
point p1, p2;
{
//find normals
vector<point> ans = solve_eq({cs[1] - cs[0], cs[2] - cs[0]}, {rs[0] - rs[1], rs[0] - rs[2]});
if (sz(ans) == 0) {
return solve_easy(cs, rs);
} else {
p1 = ans[0], p2 = ans[1];
}
}
ld vol = ((cs[1] - cs[0]) ^ (cs[2] - cs[0])) * p1;
if (vol < 0) {
vol *= -1;
swap(p1, p2);
}
ld res = (rs[0] + rs[1] + rs[2]) / 3 * vol; //volume of "prism"
//in the next loop we calculate volume of 3 cones and 3 parts of balls
vector<ld> vols;
for (int i = 0; i < 3; i++) {
point v = cs[1] - cs[0];
point w = (-v).norm();
point a1 = p1 ^ w, a2 = p2 ^ w;
ld alpha = get_angle(a1, a2); //"angle" of the part of the cone
if ((p1 ^ w) * p2 < 0) {
alpha = 2 * pi - alpha;
}
//volume of the part of the cone
res += (sqr(rs[0]) + rs[0] * rs[1] + sqr(rs[1])) / 6 * v.len() * sqr((w ^ p1).len()) * alpha;
{
//now we calculate volume of the part of the ball
ld e = (rs[1] - rs[0]) / v.len();
ld cur = (1 - e) * alpha;
ld gamma = get_angle(w ^ p1, p2 ^ p1);
if ((p1 ^ p2) * w > 0) {
gamma *= -1;
}
cur -= alpha + 2 * gamma - pi;
vols.pb(cur);
}
rotate(cs.begin(), cs.begin() + 1, cs.end());
rotate(rs.begin(), rs.begin() + 1, rs.end());
}
ld sum = 0;
for (int i = 0; i < 3; i++) {
ld val = vols[i] - vols[(i + 2) % 3];
while (val < 0) {
val += 4 * pi;
}
while (val > 4 * pi) {
val -= 4 * pi;
}
sum += val;
res += rs[i] * rs[i] * rs[i] * val / 3;
}
return res;
}
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;
}