题解：CF600D Area of Two Circles' Intersection
dangbowen1008
·
·
题解
记一些卡精度技巧，这下再也不怕 ymg 讲邪恶计算几何了。
题面很短。有点像 CSP - S 2021 阅读程序第一题（求两球的体积交）。
记圆心距为 d，两圆半径分别为 r_1,r_2。
好日子先过，先把平凡解判掉：
-
相离或外切：当 d^2 \ge (r_1+r_2)^2 时，两圆没有重叠面积，答案为 0。
-
内含或内切：当 d^2 \le (r_1-r_2)^2 时，小圆完全在大圆内部，交集就是小圆本身，答案为 \pi \min(r_1, r_2)^2。
对于相交的情况，先画图。
相交面积 S 可以如此拆解：
\begin{align*}
S &= (S_{\text{扇形}O_1AB} - S_{\triangle O_1AB}) + (S_{\text{扇形}O_2AB} - S_{\triangle O_2AB}) \\
&= S_{\text{扇形}O_1AB} + S_{\text{扇形}O_2AB} - S_{\text{筝形}O_1AO_2B} \\
&= S_{\text{扇形}O_1AB} + S_{\text{扇形}O_2AB} - 2S_{\triangle O_1O_2A}
\end{align*}
注意到 \triangle O_1O_2A 的三边都是已知的，说明它是一个可解三角形，按照最直观的想法，我们会对由两圆心和一个交点构成的 \triangle O_1O_2A 使用余弦定理：
\cos\angle AO_1O_2 = \frac{d^2 + r_1^2 - r_2^2}{2dr_1}
从而可以用 \arccos 算出一扇形半圆心角 \alpha=\angle AO_1O_2，对称地得另一扇形半圆心角，从而求得扇形面积。注意到我们求出的 \alpha 是半圆心角。完整的扇形圆心角为 2\alpha，代入扇形面积公式 \frac1 2 r^2 (2\alpha)，恰好抵消，即面积为 r^2 \alpha。
对于三角形的面积，考虑使用海伦公式。但这题卡精度，为了避免除以 2 的浮点运算，我们将其代数展开：
16 S_{\triangle}^2 = (d+r_1+r_2)(d+r_1-r_2)(d-r_1+r_2)(r_1+r_2-d)
记 Y=\sqrt{\text{RHS}}=\sqrt{(d+r_1+r_2)(d+r_1-r_2)(d-r_1+r_2)(r_1+r_2-d)}，则 S_{\triangle O_1O_2A}=\frac Y 4。
依据以上式子计算即可。
\text{Ans} = r_1^2 \arccos \left( \frac{d^2 + r_1^2 - r_2^2}{2dr_1} \right) + r_2^2 \arccos \left( \frac{d^2 + r_2^2 - r_1^2}{2dr_2} \right) - \frac{1}{2} Y
……？
WA on test 28。
精度爆炸的离谱（差出了 \approx 1000）。
我们来分析一下是哪里出了问题。这篇题解 声称，“精度问题来自海伦公式”。AI 则声称，这是错误的。精度问题实际上来自于 std::acos。将其换成 <cmath> 中的 std::atan2（其中 \operatorname{atan2}(Y,X) 在数学中可以理解为 \arctan(\frac Y X)）就不会爆精度。
:::info[为什么 \arccos 会炸掉？为什么 \arctan 不会？]
相反，$\arctan(x)$ 的导数是 $\frac 1 {1+x^2}$。这个导数是有界的，不会趋近于无穷大。
:::
:::info[为什么其他题解的 $\arccos$ 没有炸掉？]
记扇形圆心角为 $\theta$。现在的 [高赞题解](https://www.luogu.com.cn/article/if90wuuf) 是这样算 $S_1=S_{\text{扇形}O_1AB} - S_{\triangle O_1AB}$ 的：
$$
S_1 = \frac{1}{2} r^2 \theta - \frac{1}{2} r^2 \sin \theta = \frac{1}{2} r^2 (\theta - \sin \theta)
$$
我们定义函数 $f(x) = x - \sin x$。
代码里传入的不是真实的 $\theta$，而是带有误差的 $\theta + \Delta$。用泰勒展开估算误差：$f(\theta + \Delta) \approx f(\theta) + f'(\theta) \cdot \Delta, f'(x) = 1 - \cos x$。当 $\theta \to 0$ 时，$1 - \cos \theta \approx 0$。这意味着，尽管 `std::acos` 制造了一个极其庞大的角度误差 $\Delta$（比如导致面积本应偏离 $10^6$），但这个误差被乘以了 $1 - \cos \theta$（一个接近 $10^{-8}$ 的极小数），导致看上去没有误差。
:::
于是尝试凑一凑对边比邻边。设两圆交点构成的公共弦 $AB$ 与连心线 $O_1O_2$ 垂直交于点 $H$。在 $\text{Rt}\triangle O_1HA$ 中，顶角即为半圆心角 $\alpha$，斜边为 $r_1$，对边为半弦长 $\frac1 2|AB|$，记邻边为 $d_1=O_1H$，记 $d_2=d-d_1=O_2H$。
在 $\text{Rt}\triangle O_1HA$ 和 $\text{Rt}\triangle O_2HA$ 中，
$$
d_1^2 - d_2^2 = r_1^2 - r_2^2
$$
解得 $d_1 = \frac{d^2 + r_1^2 - r_2^2}{2d}$。
所以
$$
\alpha = \operatorname{atan2}(\frac1 2|AB|,\frac{d^2 + r_1^2 - r_2^2}{2d})=\operatorname{atan2}(d|AB|,d^2+r_1^2-r_2^2)
$$
注意到 $d|AB|$ 恰好是筝形的两条对角线之积，联想到筝形的另一个面积公式。
$$
\frac 1 2 d|AB|=\frac 1 2 |O_1O_2| \times |AB|=S_{\text{筝形}O_1AO_2B}=\frac 1 2Y \implies d|AB|=Y
$$
好妙！这里和前文的 $Y$ 联系上了。所以 $\alpha=\operatorname{atan2}(Y,d^2+r_1^2-r_2^2)$。
因此答案如果用 $\operatorname{atan2}$ 表述，就是
$$
\text{Ans} = r_1^2 \operatorname{atan2}( Y,d^2+r_1^2-r_2^2) + r_2^2 \operatorname{atan2}( Y,d^2+r_2^2-r_1^2) - \frac{1}{2} Y
$$
就可以过了。
不过好像还有一个坑没有填。
:::info[为什么要使用 `std::atan2` 而不是 `std::atan`？]
很多人误以为 `std::atan2(Y, X)` 只是 `std::atan(Y / X)` 的语法糖，但这个想法在计算几何中是致命的。
首先，上面我们讲的做法看似完美无瑕。
……但你做了几百道题的直觉提醒你：你有没有遍历整个实数空间？你有没有忽略边界情况？
一个我们可能忽略的情况。
当一个圆的圆心（不妨设为 $O_1$）在 $\odot O_2$ 内部时，四边形 $O_1AO_2B$ 可能会变为一个燕尾形，而不是一个筝形。
不过好在，在数学考试里，只要我们肯耐下心再分类讨论一遍，你就可以发现你**得到了完全相同的公式**。
可是在信息竞赛里面你要考虑得就多了。
这种情况中，$\alpha = \angle AO_1O_2 > \frac \pi 2$。可是无论是在数学里的 $\arctan(x)$ 还是信息里的 `std::atan` 的值域都是 $(-\frac \pi 2,\frac \pi 2)$ 啊？！
所以此时如果你使用 `std::atan` 是无法算出真正的 $\alpha$ 的，此时需要特判。
但是 `std::atan2` 就可以产生在 $(-\pi,\pi]$ 内的所有角度，因为它可以通过传入的 `Y` 和 `X` 的符号来判断角度的象限。所以我们为了避免额外特判，使用 `std::atan2`。
:::
:::success[Code]
使用了 AI 优化代码风格。
```cpp
#include <iostream>
#include <cmath>
#include <iomanip>
#include <algorithm>
using namespace std;
typedef long double ld;
typedef long long ll;
inline ll sqr(ll x) { return x * x; }
int main() {
ios::sync_with_stdio(false); cin.tie(nullptr);
ll x1, y1, r1, x2, y2, r2;
cin >> x1 >> y1 >> r1 >> x2 >> y2 >> r2;
ll d_2 = sqr(x1 - x2) + sqr(y1 - y2);
ld ans = 0;
if (d_2 >= sqr(r1 + r2)) {
ans = 0;
} else if (d_2 <= sqr(r1 - r2)) {
ans = sqr(min(r1, r2)) * acosl(-1.0L);
} else {
ll r1_2 = sqr(r1), r2_2 = sqr(r2);
ld Y = sqrtl((ld)(sqr(r1 + r2) - d_2) * (d_2 - sqr(r1 - r2)));
ans = r1_2 * atan2l(Y, d_2 + r1_2 - r2_2) + r2_2 * atan2l(Y, d_2 + r2_2 - r1_2) - 0.5 * Y;
}
cout << fixed << setprecision(10) << ans << "\n";
return 0;
}
```
:::
---
关于数学公式的问题，我找不到其他的方式来表示 $S_{\text{扇形}}$ 等等了。我专门翻了下人教版九年级上册数学书，他用的也是这个形式。求管理员通过……或者可以提供细致的建议。