自认为推导的式子没有问题

设两圆 $O_1, O_2$ 的圆心分别为 $(x_1, y_1), (x_2, y_2)$，半径分别为 $r_1, r_2$，则有：

交点所在直线 $2(x_2 - x_1)x + 2(y_2 - y_1)y = r_1^2 - r_2^2 - x_1^2 + x_2^2 - y_1^2 + y_2^2$

关于交点横坐标的一元二次方程（$y_1 ≠ y_2$，设 $D = r_1^2 - r_2^2 - x_1^2 + x_2^2 - y_1^2 + y_2^2$） $\Big[\Big(\frac{x_2 - x_1}{y_2 - y_1}\Big)^2 + 1\Big]x^2 + \Big[\frac{y_1(x_2 - x_1)}{y_2 - y_1} - \frac{D(x_2 - x_1)}{(y_2 - y_1)^2} - 2x_1\Big]x + \Big[\frac{D}{2(y_2 - y_1)}\Big]^2 - \frac{Dy_1}{y_2 - y_1} + x_1^2 + y_1^2 - r_1^2 = 0$

$y_1 = y_2$ 的情况见代码。

$\colorbox{RED}{\textcolor{WHITE}{\texttt{WA}}} \texttt{ Code:}$

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

#define EPS 1e-6

inline bool eq(long double a, long double b) {
	return a > b ? (a - b < EPS) : (b - a < EPS);
}

struct coord {
	long double x, y;
	
	bool operator <(const coord &ref) const {
		if (!eq(x, ref.x)) {
			return x < ref.x;
		}
		return y < ref.y;
	}
	
	bool operator ==(const coord &ref) const {
		return eq(x, ref.x) && eq(y, ref.y);
	}
	
	double dist(const coord &ref) {
		return sqrt((x - ref.x) * (x - ref.x) + (y - ref.y) * (y - ref.y));
	}
};

inline long double sqr(long double a) {
	return a * a;
}

struct circle {
	coord o;
	long double r;
	
	bool operator ==(const circle &ref) const {
		return o == ref.o && eq(r, ref.r);
	}
	
	vector<coord> intersect(const circle &ref) {
		vector<coord> result;
		long double dc = o.dist(ref.o), d = r + ref.r - dc;
//		printf("dc = %.2Lf, d = %.2Lf\n", dc, d);
		if (d < 0 || dc < abs(r - ref.r)) {
			return result;
		}
		if (eq(d, 0)) {
			result.push_back({o.x + (ref.o.x - o.x) * r / (r + ref.r), o.y + (ref.o.y - o.y) * r / (r + ref.r)});
			return result;
		}
		if (eq(o.y, ref.o.y)) {
			long double x = (sqr(r) - sqr(ref.r) + sqr(ref.o.x) - sqr(o.x)) / (ref.o.x - o.x) / 2;
			long double a = 1, b = -2 * ref.o.y, c = sqr(x) - 2 * o.x * x + sqr(o.x) + sqr(o.y) - sqr(r);
			long double y = (-b - sqrt(sqr(b) - 4 * a * c)) / a / 2;
			result.push_back({x, y});
			y = (-b + sqrt(sqr(b) - 4 * a * c)) / a / 2;
			result.push_back({x, y});
		} else {
			long double a = sqr((ref.o.x - o.x) / (ref.o.y - o.y)) + 1;
			d = sqr(r) - sqr(ref.r) - sqr(o.x) + sqr(ref.o.x) - sqr(o.y) + sqr(ref.o.y);
			long double b = o.y * (ref.o.x - o.x) / (ref.o.y - o.y) - d * (ref.o.x - o.x) / sqr(ref.o.y - o.y) - 2 * o.x;
			long double c = sqr(d / (ref.o.y - o.y) / 2) - d * o.y / (ref.o.y - o.y) + sqr(o.x) + sqr(o.y) - sqr(r);
//			printf("x = (-%.2Lf - sqrt(sqr(%.2Lf) - 4 * %.2Lf * %.2Lf)) / %.2Lf / 2\n", b, b, a, c, a);
			long double x = (-b - sqrt(sqr(b) - 4 * a * c)) / a / 2;
			long double y = (sqr(r) - sqr(ref.r) + sqr(ref.o.x) - sqr(o.x) + sqr(ref.o.y) - sqr(o.y)) / (ref.o.y - o.y) / 2 - (o.x - ref.o.x) * x / (o.y - ref.o.y);
			result.push_back({x, y});
			x = (-b + sqrt(sqr(b) - 4 * a * c)) / a / 2;
			y = (sqr(r) - sqr(ref.r) + sqr(ref.o.x) - sqr(o.x) + sqr(ref.o.y) - sqr(o.y)) / (ref.o.y - o.y) / 2 - (o.x - ref.o.x) * x / (o.y - ref.o.y);
			result.push_back({x, y});
		}
		return result;
	}
} c1, c2;

int main() {
//	freopen("test.in", "r", stdin);
//	freopen("test.out", "w", stdout);
	while (scanf("%Lf%Lf%Lf%Lf%Lf%Lf", &c1.o.x, &c1.o.y, &c1.r, &c2.o.x, &c2.o.y, &c2.r) != EOF) {
		if (c1 == c2) {
			puts("THE CIRCLES ARE THE SAME");
			continue;
		}
		auto ans = c1.intersect(c2);
		if (!ans.size()) {
			puts("NO INTERSECTION");
			continue;
		}
		for (auto i : ans) {
			printf("(%.3Lf,%.3Lf)", i.x, i.y);
		}
		putchar('\n');
	}
	return 0;
}