JZOJ 6678. 【2020.05.01省选模拟】苏菲的世界 (simpson积分+几何法求多个圆的并的面积)

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