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