题解：P17221 [ICPC 2017 Nanning R] Resonators
lailai0916
·
·
题解
题意简述
给定平面上的若干圆。
设它们覆盖区域的并集为 \Omega，求：
\iint_\Omega(x^2+y^2)\mathrm{d}x\mathrm{d}y
解题思路
面积并集不便直接分割，但圆并的边界只由若干圆弧组成。
可以用格林公式把面积分转化为边界上的线积分。
取：
\begin{aligned}
P & =-\frac{y(x^2+y^2)}{4} \\
Q & =\frac{x(x^2+y^2)}{4}
\end{aligned}
它们满足：
\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}=x^2+y^2
由格林公式：
\iint_\Omega(x^2+y^2)\mathrm{d}x\mathrm{d}y
=\frac{1}{4}\oint_{\partial\Omega}(x^2+y^2)(x\mathrm{d}y-y\mathrm{d}x)
沿每个圆逆时针遍历未被其他圆覆盖的边界弧。
此时圆内部始终位于行进方向左侧。
外边界会按逆时针统计，洞的边界会自然按顺时针拼接，
二者都符合格林公式对正向边界的要求。
设当前圆的圆心为 (a,b)，半径为 r，参数方程为：
\begin{aligned}
x & =a+r\cos t \\
y & =b+r\sin t
\end{aligned}
代入可得：
x\mathrm{d}y-y\mathrm{d}x=r(r+a\cos t+b\sin t)\mathrm{d}t
把它与 x^2+y^2 相乘并积分，得到原函数：
F(t)=\frac{r}{4}\left((2r(a^2+b^2)+r^3)t
+\frac{r(a^2-b^2)}{2}\sin 2t-rab\cos 2t
+(a^2+b^2+3r^2)(a\sin t-b\cos t)\right)
因此未覆盖圆弧 [\alpha,\beta] 的贡献为：
F(\beta)-F(\alpha)
接下来求每个圆上被其他圆覆盖的角区间。
若当前圆被另一个圆完整包含，
它的整条边界都不会出现在圆并边界上，可以直接跳过。
若当前圆完整包含另一个圆，后者不会覆盖当前圆边界。
两圆外离或外切时也没有正长度覆盖区间。
剩下的是两圆相交的情况。
设当前圆半径为 r，另一圆半径为 R，
圆心距为 d，圆心连线方向角为 \theta。
由余弦定理，两个交点相对 \theta 的偏角为：
\delta=\arccos\left(\frac{r^2+d^2-R^2}{2rd}\right)
另一个圆覆盖当前圆的角区间正是
若区间跨过 $0$，就把它拆成位于
$[0,2\pi]$ 内的两段。
把所有覆盖区间排序并求并集。
扫描过程中遇到的空隙就是未覆盖弧，
用原函数之差累加其贡献。
相切只涉及单个点，不影响积分。
完全相同的圆必须先去重。
否则它们会互相判定为完整覆盖，导致这条边界被重复删除。
每个圆产生 $O(n)$ 个角区间，排序需要 $O(n\log n)$。
总时间复杂度为 $O(n^2\log n)$，空间复杂度为 $O(n)$。
## 参考代码
```cpp
#include <bits/stdc++.h>
using namespace std;
using ld=long double;
const int N=1005;
const ld eps=1e-15L;
const ld pi=acosl(-1.0L);
struct circle
{
int x,y,r;
};
circle c[N];
bool operator<(const circle &a,const circle &b)
{
if(a.x!=b.x)return a.x<b.x;
if(a.y!=b.y)return a.y<b.y;
return a.r<b.r;
}
bool operator==(const circle &a,const circle &b)
{
return a.x==b.x&&a.y==b.y&&a.r==b.r;
}
ld calc(const circle &a,ld t)
{
ld x=a.x;
ld y=a.y;
ld r=a.r;
ld s=x*x+y*y;
return r/4*((2*r*s+r*r*r)*t+r*(x*x-y*y)/2*sinl(2*t)-r*x*y*cosl(2*t)+(s+3*r*r)*(x*sinl(t)-y*cosl(t)));
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T;
cin>>T;
cout<<fixed<<setprecision(3);
while(T--)
{
int n;
cin>>n;
for(int i=1;i<=n;i++)cin>>c[i].x>>c[i].y>>c[i].r;
sort(c+1,c+n+1);
n=unique(c+1,c+n+1)-c-1;
ld ans=0;
vector<pair<ld,ld>> seg;
seg.reserve(n*2);
for(int i=1;i<=n;i++)
{
seg.clear();
bool cover=0;
for(int j=1;j<=n;j++)
{
if(i==j)continue;
int x=c[j].x-c[i].x;
int y=c[j].y-c[i].y;
int d2=x*x+y*y;
int sub=c[j].r-c[i].r;
if(sub>=0&&d2<=sub*sub)
{
cover=1;
break;
}
if(sub<=0&&d2<=sub*sub)continue;
int sum=c[i].r+c[j].r;
if(d2>=sum*sum)continue;
ld mid=atan2l(y,x);
if(mid<0)mid+=2*pi;
ld d=sqrtl(d2);
ld v=((ld)c[i].r*c[i].r+d2-(ld)c[j].r*c[j].r)/(2*c[i].r*d);
ld len=acosl(max(-1.0L,min(1.0L,v)));
ld l=mid-len;
ld r=mid+len;
if(l<0)
{
seg.push_back({l+2*pi,2*pi});
l=0;
}
if(r>2*pi)
{
seg.push_back({0,r-2*pi});
r=2*pi;
}
seg.push_back({l,r});
}
if(cover)continue;
sort(seg.begin(),seg.end());
ld pos=0;
for(auto [l,r]:seg)
{
if(l>pos+eps)ans+=calc(c[i],l)-calc(c[i],pos);
pos=max(pos,r);
if(pos>=2*pi-eps)break;
}
if(pos<2*pi-eps)ans+=calc(c[i],2*pi)-calc(c[i],pos);
}
cout<<ans<<'\n';
}
return 0;
}
```