1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162 | #include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
using namespace std;
const double EPS = 1e-8; // 精度系数
const double PI = acos(-1.0); // π
const int N = 4;
struct Point {
double x, y;
Point(double x = 0, double y = 0) : x(x), y(y) {}
const bool operator<(Point A) const { return x == A.x ? y < A.y : x < A.x; }
}; // 点的定义
typedef Point Vector; // 向量的定义
Vector operator+(Vector A, Vector B) {
return Vector(A.x + B.x, A.y + B.y);
} // 向量加法
Vector operator-(Vector A, Vector B) {
return Vector(A.x - B.x, A.y - B.y);
} // 向量减法
Vector operator*(Vector A, double p) {
return Vector(A.x * p, A.y * p);
} // 向量数乘
Vector operator/(Vector A, double p) {
return Vector(A.x / p, A.y / p);
} // 向量数除
int dcmp(double x) {
if (fabs(x) < EPS)
return 0;
else
return x < 0 ? -1 : 1;
} // 与0的关系
double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y; } // 向量点乘
double Length(Vector A) { return sqrt(Dot(A, A)); } // 向量长度
double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x; } // 向量叉乘
Point GetLineProjection(Point P, Point A, Point B) {
Vector v = B - A;
return A + v * (Dot(v, P - A) / Dot(v, v));
} // 点在直线上投影
struct Circle {
Point c;
double r;
Circle() : c(Point(0, 0)), r(0) {}
Circle(Point c, double r = 0) : c(c), r(r) {}
Point point(double a) {
return Point(c.x + cos(a) * r, c.y + sin(a) * r);
} // 输入极角返回点坐标
}; // 圆
// a[i] 和 b[i] 分别是第i条切线在圆A和圆B上的切点
int getTangents(Circle A, Circle B, Point* a, Point* b) {
int cnt = 0;
if (A.r < B.r) {
swap(A, B);
swap(a, b);
}
double d2 =
(A.c.x - B.c.x) * (A.c.x - B.c.x) + (A.c.y - B.c.y) * (A.c.y - B.c.y);
double rdiff = A.r - B.r;
double rsum = A.r + B.r;
if (dcmp(d2 - rdiff * rdiff) < 0) return 0; // 内含
double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
if (dcmp(d2) == 0 && dcmp(A.r - B.r) == 0) return -1; // 无限多条切线
if (dcmp(d2 - rdiff * rdiff) == 0) { // 内切，一条切线
a[cnt] = A.point(base);
b[cnt] = B.point(base);
++cnt;
return 1;
}
// 有外公切线
double ang = acos(rdiff / sqrt(d2));
a[cnt] = A.point(base + ang);
b[cnt] = B.point(base + ang);
++cnt;
a[cnt] = A.point(base - ang);
b[cnt] = B.point(base - ang);
++cnt;
if (dcmp(d2 - rsum * rsum) == 0) { // 一条内公切线
a[cnt] = A.point(base);
b[cnt] = B.point(PI + base);
++cnt;
} else if (dcmp(d2 - rsum * rsum) > 0) { // 两条内公切线
double ang = acos(rsum / sqrt(d2));
a[cnt] = A.point(base + ang);
b[cnt] = B.point(PI + base + ang);
++cnt;
a[cnt] = A.point(base - ang);
b[cnt] = B.point(PI + base - ang);
++cnt;
}
return cnt;
} // 两圆公切线 返回切线的条数，-1表示无穷多条切线
Circle Inversion_C2C(Point O, double R, Circle A) {
double OA = Length(A.c - O);
double RB = 0.5 * ((1 / (OA - A.r)) - (1 / (OA + A.r))) * R * R;
double OB = OA * RB / A.r;
double Bx = O.x + (A.c.x - O.x) * OB / OA;
double By = O.y + (A.c.y - O.y) * OB / OA;
return Circle(Point(Bx, By), RB);
} // 点 O 在圆 A 外，求圆 A 的反演圆 B，R 是反演半径
Circle Inversion_L2C(Point O, double R, Point A, Vector v) {
Point P = GetLineProjection(O, A, A + v);
double d = Length(O - P);
double RB = R * R / (2 * d);
Vector VB = (P - O) / d * RB;
return Circle(O + VB, RB);
} // 直线反演为过 O 点的圆 B，R 是反演半径
bool theSameSideOfLine(Point A, Point B, Point S, Vector v) {
return dcmp(Cross(A - S, v)) * dcmp(Cross(B - S, v)) > 0;
} // 返回 true 如果 A B 两点在直线同侧
int main() {
int T;
scanf("%d", &T);
while (T--) {
Circle A, B;
Point P;
scanf("%lf%lf%lf", &A.c.x, &A.c.y, &A.r);
scanf("%lf%lf%lf", &B.c.x, &B.c.y, &B.r);
scanf("%lf%lf", &P.x, &P.y);
Circle NA = Inversion_C2C(P, 10, A);
Circle NB = Inversion_C2C(P, 10, B);
Point LA[N], LB[N];
Circle ansC[N];
int q = getTangents(NA, NB, LA, LB), ans = 0;
for (int i = 0; i < q; ++i)
if (theSameSideOfLine(NA.c, NB.c, LA[i], LB[i] - LA[i])) {
if (!theSameSideOfLine(P, NA.c, LA[i], LB[i] - LA[i])) continue;
ansC[ans++] = Inversion_L2C(P, 10, LA[i], LB[i] - LA[i]);
}
printf("%d\n", ans);
for (int i = 0; i < ans; ++i) {
printf("%.8f %.8f %.8f\n", ansC[i].c.x, ansC[i].c.y, ansC[i].r);
}
}
return 0;
}
|
