[BZOJ2178]圆的面积并(格林公式)

题面

http://darkbzoj.tk/problem/2178

题解

前置知识

• 格林公式：《高等数学》11.3节

因为本题对精度要求不高（并不是，测试点13警告），也可以用自适应Simpson近似（题解）；不过那毕竟是近似，格林公式相较于Simpson近似，速度既快，还可以求出精确值。

我们有格林公式

[{iintlimits_{D}}({frac{partial Q}{partial x}}-{frac{partial P}{partial y}})dxdy={ointlimits_{L^+}}Pdx+Qdy ]

其中D是我们要求面积的图形，(L)是此图形边界曲线，(+)保证了其方向：当观察者在(L)上按此方向行走时，D内在他近处的那一部分总在他的左边。

取(P(x,y)=-y，Q(x,y)=x)，那么等式左边就是两倍面积，所以

[S={frac{1}{2}}{ointlimits_{L^+}}xdy-ydx ]

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

对于此题，把(L^+)分解成每一个圆上的没有被其他圆覆盖部分，而这个部分一定由一些圆弧组成。所以我们要做的就是对于圆弧(I)，求出

[{ointlimits_{I逆时针方向}}xdy-ydx ]

设这条圆弧是在圆((x-x_0)^2+(y-y_0)^2=r^2)上、对应的极角是([{ heta_l},{ heta_r}])。（如果跨过x轴负半轴，那么拆成两段）换元，上式就是

[{int_{ heta_l}^{ heta_r}}((x_0+r cos heta)frac{d(y_0+rsin heta)}{d heta}-(y_0+rsin heta){frac{d(x_0+rcos heta)}{d heta}})d heta ]

[={int_{ heta_l}^{ heta_r}}((x_0+r cos heta)rcos heta+(y_0+rsin heta)rsin heta)d heta ]

[={int_{ heta_l}^{ heta_r}}(x_0rcos heta+y_0rsin heta+r^2)d heta ]

[=r^2( heta_r- heta_l)^2+x_0r(sin( heta_r)-sin( heta_l))-y_0r(cos( heta_r)-cos( heta_l)) ]

这就解决了本题。

时间复杂度(O(n^2log 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
",(double)ans);
	return 0;
}