How to calculate the area of two circles' intersection?

Question

The topic link: https://codeforces.com/problemset/problem/600/D

For the question, I'm wrong answer on test28, which could look like this:

correct answer:119256.95877838134765625000

my answer: 120502.639190673828125

I guess it is caused by calculation accuracy, but I don't have evidence. Maybe algorithm itself is faulty, please point it out.

Algorithm ideas:

For any given two circles, in order to simplify the calculation, we can translate the origin of the coordinates to the center of one of the circles, and then rotate the other circle to the x-axis by rotating. For calculating the intersection area of the circles, before and after are equivalent, and finally a purple circle and a red circle are formed. In fact, the final intersection area is equal to the sum of the areas of the two sectors minus the area of the diamond in the middle(the figure below, Horizontal axis x, vertical axis y). However, before that, we must first calculate the intersection point of the two circles.

The coordinates of the center of the first circle at the beginning: .

The coordinates of the center of the second circle at the beginning: .

The coordinates of the center of the first circle after a series of transformations: .

The coordinates of the center of the second circle after a series of transformations: ,

.

The equations of two circles are combined:

are the radius of the first and second circles respectively,so:

we can use the sector area formula : ,

, .

In this place, there will be problems with the positive and negative values of the radian(the two figures below), but it can be proved that they can be absorbed in the final result.

The final result is the sum of the areas of the two arcs minus the area of the middle diamond.

mycode:

# 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

I used Wolfram Alpha to check your formula, so I think you're seeing catastrophic cancellation rather than an accumulation of small errors. I'd use William Kahan's formula to compute the triangle area, like so:

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