[SPOJ8073] CIRU [自适应Simpson/格林公式][计算几何]

Link
Luogu - https://www.luogu.org/problemnew/show/SP8073
BZOJ - https://www.lydsy.com/JudgeOnline/problem.php?id=2178

---

我就不讲Simpson了

---

Green Formula
对平面单联通区域 $\mathrm{D}$ 及其正向轮廓线 $\mathrm{L}$ （分段光滑），且 $P,Q$ 在 $\mathrm{D}$ 上有一阶连续偏导
$$\underset{{}_{\mathrm{D}}}{\iint }\left(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}\right)dxdy=\underset{{}_{\mathrm{L}}}{\oint }Pdx+Qdy$$

---

取 $Q=x,P=-y$ 得到 $\underset{{}_{\mathrm{D}}}{\iint }2dxdy=\underset{{}_{\mathrm{L}}}{\oint }(-\text{−}-ydx+xdy)$
那么 $\underset{\mathrm{D}}{\iint }dxdy=\frac{1}{2}\underset{\mathrm{L}}{\oint }(-\text{−}-ydx+xdy)$

---

圆 $C$ 上弧 $L$ 参数方程$\{\begin{array}{l}\begin{array}{rclc}x& =& x_0+r\cos t& \\ y& =& y_0+r\sin t& \end{array},t\in [{\theta }_1,{\theta }_2]\end{array}$
有
$\begin{array}{rcl}S=\underset{L}{\oint }(-\text{−}-ydx+xdy)& =& \underset{{\theta }_1}{\overset{{\theta }_2}{\int }}-\text{−}-(y_0+r\sin t)d(x_0+r\cos t)+(x_0+r\cos t)d(y_0+r\sin t)\\ & =& \cdots \\ & =& \underset{{\theta }_1}{\overset{{\theta }_2}{\int }}(r^2+x_{0}r\cos t+y_{0}r\sin t)dt\end{array}$
（中间公式被吞过，懒得再打一遍了……反正就那样嘛，不懂微分怎么算的话百度微分运算法则）
然后最后对每个圆算一下产生贡献的弧，这个差分一下就可以辣

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int MAXN = 1005;
const double eps = 1e-12;
const double tpi = acos(-1.0);
int n;
double Ans;
struct Circle {
	double x, y, r;
	inline double f(const double& t) {
		return (t * r + sin(t) * x - cos(t) * y) * r;
	}
	inline double oint(const double& theta_l, const double theta_r) {
		return r * (r * (theta_r - theta_l) + x * (sin(theta_r) - sin(theta_l)) - y * (cos(theta_r) - cos(theta_l)));
//		return f(theta_r) - f(theta_l);
	}
	inline bool operator == (const Circle& o) const {
		return (x==o.x)&&(y==o.y)&&(r==o.r);
	}
	inline bool operator < (const Circle& o) const {
		return (x==o.x)?((y==o.y)?(r<o.r):(y<o.y)):(x<o.x);
	}
}C[MAXN];
inline double defabs(const double&x) {
	return (x<0)?(-x):x;
}
int tot;
pair<double, bool> Pt[(MAXN<<1)+5];
inline void Calc(const int& x) {
	tot = 0;
	register int Sm = 0;
	register double DistX, DistY, Dist, Angle, Cosine, Theta, Thetb;
	for (register int i = 1; i <= n; ++i) {
		if (i == x) continue;
		DistX = C[i].x - C[x].x;
		DistY = C[i].y - C[x].y;
		Dist = sqrt(DistX * DistX + DistY * DistY);
		if (C[x].r + Dist <= C[i].r + eps) return;
		if (C[i].r + Dist <= C[x].r + eps || Dist >= C[i].r + C[x].r - eps) continue;
		Angle = atan2(DistY, DistX),
		Cosine = acos((Dist*Dist+C[x].r*C[x].r-C[i].r*C[i].r)/(2.0*Dist*C[x].r));
		Theta = Angle - Cosine;
		Thetb = Angle + Cosine;
		if (Theta < -tpi) Theta += tpi * 2;
		if (Thetb >= tpi) Thetb -= tpi * 2;
		if (Theta > Thetb) ++Sm;
		Pt[++tot] = make_pair(Theta, 0);
		Pt[++tot] = make_pair(Thetb, 1);
	}
	Pt[0] = make_pair(-tpi, 0);
	Pt[++tot] = make_pair(tpi, 0);
	sort(Pt + 1,Pt + tot);
	for (register int i = 1; i <= tot; ++i) {
		if (!Sm) Ans += C[x].oint(Pt[i-1].first, Pt[i].first);
		(Pt[i].second)?(--Sm):(++Sm);
	}
}
int main() {
	scanf("%d", &n);
	for (register int i = 0, j = 1; j <= n; ++j) {
		++i;
		scanf("%lf%lf%lf", &C[i].x , &C[i].y, &C[i].r);
	}
	sort(C + 1, C + 1 + n);
	n = unique(C + 1, C + 1 + n) - C - 1;
	for (register int i = 1; i <= n; ++i) {
		Calc(i);
	}
	printf("%.3f", Ans / 2);
	return 0;
}

---

没有名字的一般计算几何方法
好！

---

当然也可以求面积交
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int MAXN = 1005;
const double eps = 1e-12;
const double tpi = acos(-1.0);
int n;
double Ans;
struct Circle {
	double x, y, r;
	inline double f(const double& t) {
		return (t * r + sin(t) * x - cos(t) * y) * r;
	}
	inline double oint(const double& theta_l, const double theta_r) {
		return r * (r * (theta_r - theta_l) + x * (sin(theta_r) - sin(theta_l)) - y * (cos(theta_r) - cos(theta_l)));
//		return f(theta_r) - f(theta_l);
	}
	inline bool operator == (const Circle& o) const {
		return (x==o.x)&&(y==o.y)&&(r==o.r);
	}
	inline bool operator < (const Circle& o) const {
		return (x==o.x)?((y==o.y)?(r<o.r):(y<o.y)):(x<o.x);
	}
}C[MAXN];
inline double defabs(const double&x) {
	return (x<0)?(-x):x;
}
int tot;
pair<double, bool> Pt[(MAXN<<1)+5];
inline void Calc(const int& x) {
	tot = 0;
	register int Sm = 0;
	register double DistX, DistY, Dist, Angle, Cosine, Theta, Thetb;
	for (register int i = 1; i <= n; ++i) {
		if (i == x) continue;
		DistX = C[i].x - C[x].x;
		DistY = C[i].y - C[x].y;
		Dist = sqrt(DistX * DistX + DistY * DistY);
		if (C[x].r + Dist <= C[i].r + eps) return;
		if (C[i].r + Dist <= C[x].r + eps || Dist >= C[i].r + C[x].r - eps) continue;
		Angle = atan2(DistY, DistX),
		Cosine = acos((Dist*Dist+C[x].r*C[x].r-C[i].r*C[i].r)/(2.0*Dist*C[x].r));
		Theta = Angle - Cosine;
		Thetb = Angle + Cosine;
		if (Theta < -tpi) Theta += tpi * 2;
		if (Thetb >= tpi) Thetb -= tpi * 2;
		if (Theta > Thetb) ++Sm;
		Pt[++tot] = make_pair(Theta, 0);
		Pt[++tot] = make_pair(Thetb, 1);
	}
	Pt[0] = make_pair(-tpi, 0);
	Pt[++tot] = make_pair(tpi, 0);
	sort(Pt + 1,Pt + tot);
	for (register int i = 1; i <= tot; ++i) {
		if (!Sm) Ans += C[x].oint(Pt[i-1].first, Pt[i].first);
		(Pt[i].second)?(--Sm):(++Sm);
	}
}
int main() {
	double Tem = 0;
	scanf("%d", &n);
	for (register int i = 1; i <= n; ++i) {
		scanf("%lf%lf%lf", &C[i].x , &C[i].y, &C[i].r);
		Tem += C[i].r * C[i].r * tpi;
	}
	sort(C + 1, C + 1 + n);
	n = unique(C + 1, C + 1 + n) - C - 1;
	for (register int i = 1; i <= n; ++i) {
		Calc(i);
	}
	printf("%.3f", Tem - Ans / 2);
	return 0;
}