#include <vector>
using namespace std;
- const double EPS = 1e-8; //精度系数
- const double PI = acos(-1.0); //π
+ 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;
- }
+ 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; //向量的定义
+ 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); } //向量数除
+ 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的关系
+ 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; } //向量叉乘
+ 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));
- } //点在直线上投影
+ 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); } //输入极角返回点坐标
- }; //圆
+ 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表示无穷多条切线
+ 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 是反演半径
+ 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 是反演半径
+ 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 两点在直线同侧
+ 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);
- }
+ 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;
+ }
+
+ return 0;
}
```
OSDN Git Service
