题面
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;
}