如何计算两个圆相交的面积？

Question

话题链接：https://codeforces.com/problemset/problem/600/D

对于这个问题，我在 test28 上的答案是错误的，可能如下所示：

正确答案：119256.95877838134765625000

我的回答：120502.639190673828125

我猜是计算准确性造成的，但我没有证据。可能算法本身有问题，请指出。

算法思路：

对于任意给定的两个圆，为了简化计算，我们可以将坐标原点平移到其中一个圆的中心，然后通过旋转将另一个圆旋转到x轴。计算圆的相交面积，前后等价，最后形成紫色圆和红色圆。实际上，最终相交的面积等于两个扇区的面积之和减去中间菱形的面积（下图，横轴x，纵轴y）。不过，在这之前，我们必须先计算出两个圆的交点。

开始处第一个圆的中心坐标： .

开始第二个圆的圆心坐标： .

经过一系列变换后的第一个圆的中心坐标： .

经过一系列变换后的第二个圆的中心坐标： ,

.

两个圆的方程组合起来：

分别是第一个和第二个圆的半径，所以：

我们可以使用扇形面积公式：，

, .

在这个地方，弧度的正负值（下图两个）会有问题，但可以证明可以吸收到最终结果中。

最后的结果是两条弧的面积之和减去中间菱形的面积。

我的代码：

# define MY_PI 3.14159265358979323846
using namespace std;

typedef long double ld;

ld pow2(ld x){return x * x;}
int main()
{
    cout.precision(25);
    ld x1, y1, r1, x2, y2, r2;   // (x1, y1) is the center point of the circle, r1 is the radius of the circle, the other is similar.
    cin >> x1 >> y1 >> r1 >> x2 >> y2 >> r2;

    ld c = sqrt(pow2(x2 - x1) + pow2(y2 - y1));
    if(r1 + r2 < c)         // when r1+r2 < c, there is no intersection
    {
        cout << 0 << endl;
        return 0;
    }
    if((max(r1, r2) - min(r1, r2)) >= c)  // one circle is inside the other circle.
    {
        cout << MY_PI * pow2(min(r1, r2)) << endl;
        return 0;
    }

    ld x = (pow2(r1) - pow2(r2)) / (2 * c) + c * 0.5;
    ld y = sqrt(pow2(r1) - pow2(x));
    ld angle1 = 2.0 * acos(x / r1);
    ld angle2 = 2.0 * acos((c - x) / r2);
    ld ar1 = angle1 * pow2(r1) * 0.5;
    ld ar2 = angle2 * pow2(r2) * 0.5;
    ld res = ar1 + ar2 - y * c;

    cout << res << endl;

    return 0;
}

Answer 1

我使用 Wolfram Alpha 来检查您的公式，所以我认为您看到的是 catastrophic cancellation，而不是小错误的累积。我会使用William Kahan's formula 来计算三角形面积，如下所示：

#include <algorithm>
#include <cmath>
#include <iostream>
#include <limits>
typedef long double ld;

struct Circle {
  ld x;
  ld y;
  ld r;
};

static std::istream& operator>>(std::istream& i, Circle& c) {
  return i >> c.x >> c.y >> c.r;
}

static ld Pi() { return std::acos(-1); }

static ld Square(ld x) { return x * x; }

static void SortDescending2(ld& a, ld& b) {
  if (a < b) {
    using std::swap;
    swap(a, b);
  }
}

static void SortDescending3(ld& a, ld& b, ld& c) {
  SortDescending2(a, b);
  SortDescending2(b, c);
  SortDescending2(a, b);
}

static ld KahanAreaOfTriangle(ld a, ld b, ld c) {
  SortDescending3(a, b, c);
  return 0.25 * std::sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) *
                          (a + (b - c)));
}

static ld AreaOfIntersection(const Circle& c1, const Circle& c2) {
  ld R = c1.r;
  ld r = c2.r;
  ld d = std::hypot(c2.x - c1.x, c2.y - c1.y);
  if (R + r <= d) {
    return 0;
  }
  SortDescending2(R, r);
  if (d <= R - r) {
    return Pi() * Square(r);
  }
  ld area = 2 * KahanAreaOfTriangle(R, r, d);
  ld y = area / d;
  ld A = std::asin(y / R) * Square(R);
  ld a = std::asin(y / r);
  if (std::hypot(r, d) <= R) {
    a = Pi() - a;
  }
  a *= Square(r);
  SortDescending2(A, a);
  return (A - area) + a;
}

int main() {
  Circle c1;
  Circle c2;
  if (std::cin >> c1 >> c2) {
    std::cout.precision(std::numeric_limits<ld>::max_digits10);
    std::cout << AreaOfIntersection(c1, c2) << "\n";
  }
}

Answer 2

不要使用太多中间浮动变量来获得最终答案。小的不准确加起来会导致巨大的不准确，您可以在问题的预期和实际输出中清楚地看到。您可以做的是最小化这些中间浮动变量。例如

ld c = sqrt(pow2(x2 - x1) + pow2(y2 - y1));
if(r1 + r2 < c)         // when r1+r2 < c, there is no intersection
{
    cout << 0 << endl;
    return 0;
}

您可以通过这样做消除 c 作为浮动变量

int c = pow(x2-x1) + pow(y2-y1);
if (pow(r1+r2) > pow(c)) {
   // do your stuff
}