題面

http://darkbzoj.tk/problem/2178

題解

前置知識

• 格林公式：《高等數學》11.3節

因為本題對精度要求不高（並不是，測試點13警告），也可以用自適應Simpson近似（題解）；不過那畢竟是近似，格林公式相較於Simpson近似，速度既快，還可以求出精確值。

我們有格林公式

\[{\iint\limits_{D}}({\frac{\partial Q}{\partial x}}-{\frac{\partial P}{\partial y}})dxdy={\oint\limits_{L^+}}Pdx+Qdy \]

其中D是我們要求面積的圖形，\(L\)是此圖形邊界曲線，\(+\)

保證了其方向：當觀察者在\(L\)上按此方向行走時，D內在他近處的那一部分總在他的左邊。

取\(P(x,y)=-y，Q(x,y)=x\)，那麼等式左邊就是兩倍面積，所以

\[S={\frac{1}{2}}{\oint\limits_{L^+}}xdy-ydx \]

（其實如果不懂格林公式，等式右邊的式子也可以按照多邊形求和公式的方式理解）

對於此題，把\(L^+\)分解成每一個圓上的沒有被其他圓覆蓋部分，而這個部分一定由一些圓弧組成。所以我們要做的就是對於圓弧\(I\)，求出

\[{\oint\limits_{I逆時針方向}}xdy-ydx \]

設這條圓弧是在圓\((x-x_0)^2+(y-y_0)^2=r^2\)

上、對應的極角是\([{\theta_l},{\theta_r}]\)。（如果跨過x軸負半軸，那麼拆成兩段）換元，上式就是

\[{\int_{\theta_l}^{\theta_r}}((x_0+r \cos \theta)\frac{d(y_0+r\sin\theta)}{d\theta}-(y_0+r\sin\theta){\frac{d(x_0+r\cos\theta)}{d\theta}})d\theta \]

\[={\int_{\theta_l}^{\theta_r}}((x_0+r \cos \theta)r\cos\theta+(y_0+r\sin\theta)r\sin\theta)d\theta \]

\[={\int_{\theta_l}^{\theta_r}}(x_0r\cos\theta+y_0r\sin\theta+r^2)d\theta \]

\[=r^2(\theta_r-\theta_l)^2+x_0r(\sin(\theta_r)-\sin(\theta_l))-y_0r(\cos(\theta_r)-\cos(\theta_l)) \]

這就解決了本題。

時間複雜度\(O(n^2\log n)\)。

程式碼

• 實現有一些細節，比如重合、內含、外離、r=0等等，需要特判。

#include<bits/stdc++.h>

using namespace std;

#define ld long double
#define In inline
#define rg register

const int N = 1000;
const ld eps = 1e-9;
const ld pi = acosl(-1.0);

typedef pair<ld,ld>pll;

In int read(){
	int s = 0,ww = 1;
	char ch = getchar();
	while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
	while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
	return s * ww;
}

In int sgn(ld x){
	return x < -eps ? -1 : x > eps;
}

In ld sqr(ld x){
	return x * x;
}

struct vec{
	ld x,y;
	vec(){}
	vec(ld _x,ld _y){x = _x,y = _y;}
	In friend vec operator + (vec a,vec b){
		return vec(a.x + b.x,a.y + b.y);
	}
	In friend vec operator - (vec a,vec b){
		return vec(a.x - b.x,a.y - b.y);
	}
	In friend ld len(vec a){
		return sqrt(sqr(a.x) + sqr(a.y));
	}
};

struct cir{
	vec c;ld r;
	cir(){}
	cir(vec _c,ld _r){c = _c,r = _r;};
	In friend bool HaveIts(cir a,cir b){ //非外離或外切
		ld dis = len(a.c - b.c);
		return sgn(dis - (a.r+b.r)) < 0;
	}
	In friend bool include(cir a,cir b){ //b內含於或內切於a
		if(sgn(a.r - b.r) < 0)return 0;
		ld dis = len(a.c - b.c);
		return sgn(dis - fabs(a.r-b.r)) <= 0;
	}
}l[N+5]; //loop

In bool cmp1(cir a,cir b){
	int k;
	if((k=sgn(a.c.x-b.c.x)) != 0)return k < 0;
	if((k=sgn(a.c.y-b.c.y)) != 0)return k < 0;
	return sgn(a.r - b.r) < 0;
}

In bool cmp2(cir a,cir b){
	return sgn(a.c.x - b.c.x) == 0
		&& sgn(a.c.y - b.c.y) == 0
		&& sgn(a.r - b.r) == 0;
} //equal

pll intv[2*N+5]; //覆蓋部分

int n;

In ld calc(ld L,ld R,int id){
	return sqr(l[id].r) * (R - L) + l[id].r * (l[id].c.x*(sinl(R)-sinl(L)) - l[id].c.y*(cosl(R)-cosl(L)));
}

ld f(int id){
	int cnt = 0;
	if(sgn(l[id].r) <= 0)return 0;
	for(rg int i = 1;i <= n;i++){
		if(i == id)continue;
		if(include(l[i],l[id]))return 0; //被內含
		if(!HaveIts(l[i],l[id]) || include(l[id],l[i]))continue; //無交點
		ld a = atan2l(l[i].c.y - l[id].c.y,l[i].c.x - l[id].c.x);
		ld r1 = l[id].r,r2 = l[i].r,dis = len(l[i].c - l[id].c);
		ld t = acosl((sqr(r1)+sqr(dis)-sqr(r2)) / (2*r1*dis));
		ld l = a - t,r = a + t;
		while(sgn(l + pi) < 0)l += 2 * pi;
		while(sgn(r - pi) >= 0)r -= 2 * pi;
		if(sgn(l-r) <= 0)intv[++cnt] = make_pair(l,r);
		else{
			intv[++cnt] = make_pair(l,pi);
			intv[++cnt] = make_pair(-pi,r);
		}
	}
	if(cnt)sort(intv + 1,intv + cnt + 1);
	ld ans = 0;
	ld L = -pi,R = -pi;
	for(rg int i = 1;i <= cnt;i++){
		if(sgn(intv[i].first-R) >= 0)ans += calc(R,intv[i].first,id),L = R = intv[i].first;
		if(sgn(intv[i].second-R) >= 0)R = intv[i].second;
	}
	ans += calc(R,pi,id);
	return ans / 2;
}

int main(){
	n = read();
	for(rg int i = 1;i <= n;i++){
		int x = read(),y = read(),r = read();
		l[i] = cir(vec(x,y),r);
	}
	sort(l + 1,l + n + 1,cmp1);
	n = unique(l + 1,l + n + 1,cmp2) - l - 1; //去重
	ld ans = 0;
	for(rg int i = 1;i <= n;i++)ans += f(i);
	printf("%.3lf\n",(double)ans);
	return 0;
}