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

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";
  }
}