Description：
题解：
simpson积分，据说是用二次函数拟合，拟合圆这类图形有奇效。
即(g(x,y)=int_x^y f(x)~dx≈(y-x)/6*(f(x)+4f((x+y)/2)+f(y)))
自适应的话就是一直分治，直到(|g(a,b)-g(a,(a+b)/2)-g((a+b)/2,b)|<eps)即可。
用simpson积分可以搞掉三维空间的一维。
变成二维多个圆并的面积。
这个问题可以直接几何法解决：
1.去掉被包含的圆，多个圆相同保留一个
2.枚举每一个圆，求它没有被覆盖的圆弧，每段圆弧用几角区间表示
3.对每个几角区间，再分成([-pi,0],[pi,0])两部分，用上部分的积分-下部分的积分
4.一段圆弧的积分：
变成弦的积分+一段小地方的积分（扇形-三角形）
注意这段圆弧是可能经过x轴的，画图发现小地方的贡献系数永远是+1，但弦的积分贡献系数不是，所以分开算。
过程需要用到一些三角函数知识，还行。
simpson要卡精度：取那一维的相邻两个端点间的每段来积分。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
#define db double
const db pi = acos(-1);
const int N = 50;
int n;
struct nod {
db x, y, z, r;
} a[N];
struct P {
db x, y;
P(db _x = 0, db _y = 0) {
x = _x, y = _y;
}
};
P operator - (P a, P b) { return P(a.x - b.x, a.y - b.y);}
db len(P a) { return sqrt(a.x * a.x + a.y * a.y);}
db getang(P a) { return atan2(a.y, a.x);}
struct Q {
P x; db r;
} b[N]; int b0;
int cmpb(Q a, Q b) { return a.r < b.r;}
const db eps = 1e-7;
int pd(Q a, Q b) {
return len(a.x - b.x) <= b.r - a.r + eps;
}
int bz[N][N];
P c[N * 2]; int c0;
int cmpc(P a, P b) {
return a.x < b.x;
}
db cw(Q a, db x, db y) {
return (a.x.y + sin(x) * a.r + a.x.y + sin(y) * a.r) * abs((cos(x) - cos(y)) * a.r) / 2;
}
db cw2(Q a, db x, db y) {
return a.r * a.r * (y - x) / 2 - (a.r * a.r * sin(y - x) / 2);
}
db calc(Q a) {
db ans = 0;
c[++ c0] = P(-pi, 0);
c[++ c0] = P(0, 0);
c[++ c0] = P(pi, 0);
sort(c + 1, c + c0 + 1, cmpc);
int s = 0;
fo(i, 1, c0) {
if(i > 1 && s == 0 && c[i].x - c[i - 1].x > eps) {
ans += cw(a, c[i - 1].x, c[i].x) * (c[i].x <= eps ? -1 : 1);
ans += cw2(a, c[i - 1].x, c[i].x);
}
s += round(c[i].y);
}
return ans;
}
db solve() {
if(b0 == 0) return 0;
db ans = 0;
sort(b + 1, b + b0 + 1, cmpb);
fo(i, 1, b0) fo(j, 1, b0) bz[i][j] = pd(b[i], b[j]);
fo(i, 1, b0) {
int ok = 1;
fo(j, i + 1, b0) if(bz[i][j]) ok = 0;
if(!ok) continue;
c0 = 0;
fo(j, 1, b0) if(!bz[j][i]) {
if(len(b[j].x - b[i].x) > b[i].r + b[j].r - eps) continue;
db d = getang(b[j].x - b[i].x);
db r1 = b[i].r, r2 = b[j].r, d1 = len(b[j].x - b[i].x);
db e = acos((d1 * d1 + r1 * r1 - r2 * r2) / (2 * d1 * r1));
db x = d - e, y = d + e;
if(x < -pi) x += 2 * pi;
if(y > pi) y -= 2 * pi;
c[++ c0] = P(x, 1);
c[++ c0] = P(y, -1);
if(x > y) c[++ c0] = P(-pi, 1);
}
ans += calc(b[i]);
}
return ans;
}
db f(db z) {
b0 = 0;
fo(i, 1, n) {
db d = abs(a[i].z - z);
if(d >= a[i].r) continue;
b[++ b0].x = P(a[i].x, a[i].y);
b[b0].r = sqrt(a[i].r * a[i].r - d * d);
}
return solve();
}
db g(db a, db b) {
return (b - a) / 6 * (f(a) + f(b) + f((a + b) / 2) * 4);
}
db d[N]; int d0;
db dg(db a, db b) {
db mi = (a + b) / 2;
db z1 = g(a, mi), z2 = g(mi, b), z3 = g(a, b);
if(abs(z1 + z2 - z3) < 1e-7) return z1 + z2;
return dg(a, mi) + dg(mi, b);
}
int main() {
freopen("sphere.in", "r", stdin);
freopen("sphere.out", "w", stdout);
scanf("%d", &n);
db l = 1e9, r = -1e9;
fo(i, 1, n) {
scanf("%lf %lf %lf %lf", &a[i].x, &a[i].y, &a[i].z, &a[i].r);
l = min(l, a[i].z - a[i].r);
r = max(r, a[i].z + a[i].r);
d[++ d0] = a[i].z - a[i].r;
d[++ d0] = a[i].z + a[i].r;
}
sort(d + 1, d + d0 + 1);
db ans = 0;
fo(i, 2, d0) if(d[i] - d[i - 1] > eps)
ans += dg(d[i - 1], d[i]);
pp("%.3lf
", ans);
}