【計算幾何模板】求兩個圓的交點

問題描述：

給兩個相交的圓，第一個圓的圓心爲\((x_1, \, y_1)\)，半徑爲\(r_1\)，第二個圓的圓心爲\((x_2, \, y_2)\)，半徑爲\(r_2\)，求兩個圓的交點。

問題分析：

《訓練指南》上求兩圓交點的模板用了atan2，acos等庫函數，精度損失比較嚴重。
下面介紹一種精度損失較小的作法：
原文地址

首先回顧一下圓的兩種表示方法：

• 圓的標準方程：\((x-x_0)^2+(y-y_0)^2=r^2\)
• 圓的參數方程：\(\left\{\begin{matrix} x=x_0+r\cdot cos\theta \\ y=y_0+r\cdot sin\theta \end{matrix}\right.\)

將第一個圓的參數方程帶入第二個圓的標準方程：
\((x_1+r_1cos\theta-x_2)^2+(y_1+r_1sin\theta-y_2)^2=r_2^2\)
展開後獲得：
\(2r_1(x_1-x_2)cos\theta+2r_1(y_1-y_2)sin\theta=r_2^2-r_1^2-(x_1-x_2)^2-(y_1-y_2)^2\)
令：
\(a=2r_1(x_1-x_2)\)
\(b=2r_1(y_1-y_2)\)
\(c=r_2^2-r_1^2-(x_1-x_2)^2-(y_1-y_2)^2\)
原式變爲：
\(a cos\theta+b sin\theta=c\)
令\(cos\theta=x, \, sin\theta=\sqrt{1-x^2}\)，關於\(sin\theta\)的正負後面再判斷。
代入方程，獲得，\(ax+b\sqrt{1-x^2}=c\)
移項再兩邊平方，\((ax-c)^2=b^2(1-x^2)\)
整理得，\((a^2+b^2)x^2-2acx+c^2-b^2=0\)
下面就是解一元二次方程了。

將\(sin\theta\)和\(cos\theta\)代回到第一個圓的參數方程，能獲得交點的座標。
若是該點不在第二個圓上，說明\(sin\theta\)是個負數，還須要對這個交點稍做調整。
還有一種特殊狀況就是，若是已經肯定有兩個不一樣的交點，但解出來的\(cos\theta\)值只有一個。
說明對應的\(sin\theta\)值必然一正一負。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <vector>

typedef long double LD;

const LD eps = 1e-10;

int dcmp(LD x) {
	if(fabs(x) < eps) return 0;
	return x < 0 ? -1 : 1;
}

LD sqr(LD x) { return x * x; }

struct Point
{
	LD x, y;
	Point(LD x = 0, LD y = 0):x(x), y(y) {}
	void read() { cin >> x >> y; }
};

Point operator - (const Point& A, const Point& B) {
	return Point(A.x - B.x, A.y - B.y);
}

bool operator == (const Point& A, const Point& B) {
	return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.x) == 0;
}

LD Dot(const Point& A, const Point& B) {
	return A.x * B.x + A.y * B.y;
}

LD Length(const Point& A) { return sqrt(Dot(A, A)); }

struct Circle
{
	Point c;
	LD r;
	Circle() {}
	Circle(Point c, LD r):c(c), r(r) {}
};

int getCircleIntersection(Circle C1, Circle C2) {
	LD &r1 = C1.r, &r2 = C2.r;
	LD &x1 = C1.c.x, &x2 = C2.c.x, &y1 = C1.c.y, &y2 = C2.c.y;
	LD d = Length(C1.c - C2.c);
	if(dcmp(fabs(r1-r2) - d) > 0) return -1;
	if(dcmp(r1 + r2 - d) < 0) return 0;
	LD d2 = Dot(C1.c - C2.c, C1.c - C2.c);
	LD a = r1*(x1-x2)*2, b = r1*(y1-y2)*2, c = r2*r2-r1*r1-d*d;
	LD p = a*a+b*b, q = -a*c*2, r = c*c-b*b;

	LD cosa, sina, cosb, sinb;
	//One Intersection
	if(dcmp(d - (r1 + r2)) == 0 || dcmp(d - fabs(r1 - r2)) == 0) {
		cosa = -q / p / 2;
		sina = sqrt(1 - sqr(cosa));
		Point p(x1 + C1.r * cosa, y1 + C1.r * sina);
		if(!OnCircle(p, C2)) p.y = y1 - C1.r * sina;
		inter.push_back(p);
		return 1;
	}
	//Two Intersections
	LD delta = sqrt(q * q - p * r * 4);
	cosa = (delta - q) / p / 2;
	cosb = (-delta - q) / p / 2;
	sina = sqrt(1 - sqr(cosa));
	sinb = sqrt(1 - sqr(cosb));
	Point p1(x1 + C1.r * cosa, y1 + C1.r * sina);
	Point p2(x1 + C1.r * cosb, y1 + C1.r * sinb);
	if(!OnCircle(p1, C2)) p1.y = y1 - C1.r * sina;
	if(!OnCircle(p2, C2)) p2.y = y1 - C1.r * sinb;
	if(p1 == p2)  p1.y = y1 - C1.r * sina;
	inter.push_back(p1);
	inter.push_back(p2);
	return 2;
}