二つの円の交点を求める

ABC157の F - Yakiniku Optimization Problem では二つの円の交点を求める必要がありましたが、その具体的な方法がeditorialでは触れられていなかったので、残しておくことにします。

方針

円の方程式を連立させて解くのは厳しいので、ベクトルを用います。
角度や点を図のように定義します。円の中心 と 半径 () が与えられているときに、点Aの座標を を利用して求めるのがこの記事の目的です。

最初に各辺の長さを求めていきます。
まず円の中心間の距離は です。次に三角形 で余弦定理を用いることで
　
がわかり、さらにこれを利用して
　
　
となります。

次に と を求めていきます。
まず に平行な単位ベクトル を求めます。 これは をその長さで割ることで
　
と求められます。

は単純に を 倍すればよく、
　
です。

は を反時計回りに 90 度回転させたもの を 倍すればよいです。ベクトルの回転の公式
　
で とすることで が求められ、これを用いて
　
となります。

以上のことから にて点A の座標を求めることができます。
Bについては最後に の代わりに を -90度回転させた を用いればよく、上記のベクトルの回転の公式で の代わりに とすればよいです。

ソースコード

2つの円が交点を持たない場合には対応していませんので、良しなに処理してください。

#include <bits/stdc++.h>
using namespace std;

struct Vector2D {
    // 平面ベクトル
    double x, y;

    Vector2D(double x, double y) : x(x), y(y) {}

    // ベクトルの和
    Vector2D operator+(const Vector2D p) { return Vector2D(x + p.x, y + p.y); }
    // ベクトルの差
    Vector2D operator-(const Vector2D p) { return Vector2D(x - p.x, y - p.y); }

    // ベクトルのスカラー倍
    Vector2D operator*(const double d) { return Vector2D(d * x, d * y); }
    Vector2D operator/(const double d) { return Vector2D(x / d, y / d); }

    // ベクトルの長さ
    double norm() { return sqrt(x * x + y * y); }

    // ベクトルの回転
    Vector2D rotate(double theta) {
        return Vector2D(x * cos(theta) - y * sin(theta), x * sin(theta) + y * cos(theta));
    }
};

struct Circle {
    // 中心
    Vector2D center;
    // 半径
    double radius;

    Circle(Vector2D center, double radius) : center(center), radius(radius) {}
};

vector<Vector2D> getIntersections(Circle C1, Circle C2) {
    double d = (C2.center - C1.center).norm();
    double rc = (C1.radius * C1.radius + d * d - C2.radius * C2.radius) / (2 * d);
    double rs = sqrt(C1.radius * C1.radius - rc * rc);

    Vector2D e1 = (C2.center - C1.center) / d;
    Vector2D e2 = e1.rotate(M_PI / 2);
    Vector2D e3 = e1.rotate(-M_PI / 2);
    Vector2D A = C1.center + e1 * rc + e2 * rs;
    Vector2D B = C1.center + e1 * rc + e3 * rs;

    return {A, B};
}