圆的反演变换（*）

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4773

题意：给定两个圆，告诉半径和圆心，它们是相离的，在这两个圆外给定一个点p，求符合条件：过点p的圆且与已知的两个

圆外切的所有圆的总数和它们的圆心坐标和半径。

分析：根据题意，我们设已知两个圆的半径分别为和，它们的圆心分别为和，设点p的坐标为

并设要求的圆的圆心为，半径为，那么根据它们的关系我们可以很快就列出方程组：

然后，现在的问题就是来解，但是你会发现这个方程是很难解的，因此为了简化问题，我们利用反演变换来做。

那么，怎么用反演变换来做？ 首先，得知道什么是反演变换以及它有什么性质。

反演的定义：

已知一圆C，圆心为O，半径为r，如果P与P’在过圆心O的直线上，且，则称P与P'关于O互为反演。

反演的性质：

(1)除反演中心外，平面上的每一个点都只有唯一的反演点，且这种关系是对称的，位于反演圆上的点，保持在原处，位于

反演圆外部的点，变为圆内部的点，位于反演圆内部的点，变为圆外部的点。 举个最简单的例子，区间以1为反演

半径，那么反演后的区间就是，这就是一维反演，而圆的反演是二维反演。

(2)任意一条不过反演中心的直线，它的反形是经过反演中心的圆，反之亦然，特别地，过反演中心相交的圆，变为不过反

演中心的相交直线。

定理：不过反演中心的圆，它的反形是一个圆，反演中心是这两个互为反形的圆的一个位似中心，任意一对反演点是逆对应

点。用图形来解释，如下图：

那么，对于一个不过反演中心的圆，怎样求它的反形圆？

很容易知道我们只需要求出反形圆的圆心和半径就可以了。

对于上图我们设圆C1的半径为，C2的半径为，反演半径为

那么根据反演的定义有：

那么，消去得到：

这样我们就得到了反形圆的半径，那么还要求反形圆的圆心。

由于C1和O两点的坐标已知，而且我们知道O,C1，C2位于同一直线上，那么很明显对于C2的坐标，我们可以这样计算：

设O的坐标为，C1的坐标为，C2的坐标为

那么有：

至于由上面解处可以很容易得到，这样我们就完成了圆的反演变换。

由于本题的做法是这样的，先以点P为为反演中心，反演半径随便设置都可以，为了计算方便就设为1，把圆C1和圆C2反演后再求这两个圆的

公切线，再把这个公切线反演回去，那么就是一个过点P的圆，且与原来的C1和C2相切。

那么接下来就是如何计算两个圆的公切线了。这里只需要考虑公切线在同一侧的情况。那么，这个自己画图就能很容易计算了。

找到公切线后还要把它反演成圆，这个圆还经过P点，那么很容易得到了。

[cpp]  view plain  copy

1. #include <iostream>  
2. #include <string.h>  
3. #include <stdio.h>  
4. #include <math.h>  
5.   
6. using namespace std;  
7. double const eps = 1e-8;  
8.   
9. struct Point  
10. {  
11.     double x,y;  
12.     Point(double a = 1.0,double b = 1.0):x(a),y(b){}  
13.     Point operator + (const Point &a)  
14.     {  
15.         return Point(x+a.x,y+a.y);  
16.     }  
17.     Point operator - (const Point &a)  
18.     {  
19.         return Point(x-a.x,y-a.y);  
20.     }  
21.     Point operator * (const double a)  
22.     {  
23.         return Point(a*x,a*y);  
24.     }  
25.     Point Trans()  
26.     {  
27.         return Point(-y,x);  
28.     }  
29.     void Input()  
30.     {  
31.         scanf("%lf%lf",&x,&y);  
32.     }  
33. } ;  
34.   
35. struct Circle  
36. {  
37.     Point o;  
38.     double r;  
39.     Circle(Point a = Point(),double b = 1.0):o(a),r(b) {}  
40.     Point getPoint(double alpha)  
41.     {  
42.         return o + Point(r*cos(alpha),r*sin(alpha));  
43.     }  
44.     void Input()  
45.     {  
46.         o.Input();  
47.         scanf("%lf",&r);  
48.     }  
49.     void Output()  
50.     {  
51.         printf("%.8lf %.8lf %.8lf\n",o.x,o.y,r);  
52.     }  
53. } ;  
54.   
55. Point p;  
56. Circle c[15];  
57.   
58. double dist(Point A,Point B)  
59. {  
60.     return sqrt((A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y));  
61. }  
62.   
63. double cross(Point A,Point B,Point C)  
64. {  
65.     return (B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x);  
66. }  
67.   
68. int sign(double x)  
69. {  
70.     return (x > eps) - (x < -eps);  
71. }  
72.   
73. Circle Inverse(Circle C)  
74. {  
75.     Circle T;  
76.     double t = dist(C.o,p);  
77.     double x = 1.0 / (t - C.r);  
78.     double y = 1.0 / (t + C.r);  
79.     T.r = (x - y) / 2.0;  
80.     double s = (x + y) / 2.0;  
81.     T.o = p + (C.o - p) * (s / t);  
82.     return T;  
83. }  
84.   
85. void add(Point a,Point b,int &k)  
86. {  
87.     double t = cross(a,p,b);  
88.     if(t < 0) t = -t;  
89.     double d = dist(a,b);  
90.     t /= d;  
91.     if(t > eps)  
92.     {  
93.         double w = 0.5 / t;  
94.         Point dir = (b-a).Trans();  
95.         Point a1 = p + dir * (w / d);  
96.         Point b1 = p - dir * (w / d);  
97.         if(fabs(cross(a,b,a1)) < fabs(cross(a,b,b1)))  
98.            c[k++] = Circle(a1,w);  
99.         else  
100.            c[k++] = Circle(b1,w);  
101.     }  
102. }  
103.   
104. int Work()  
105. {  
106.     c[0] = Inverse(c[0]);  
107.     c[1] = Inverse(c[1]);  
108.     if(c[1].r > c[0].r) swap(c[1],c[0]);  
109.     Point v = c[1].o - c[0].o;  
110.     double alpha = atan2(v.y,v.x);  
111.     double d = dist(c[0].o,c[1].o);  
112.     double beta  = acos((c[0].r - c[1].r) / d);  
113.     int k = 2;  
114.     Point a = c[0].getPoint(alpha + beta);  
115.     Point b = c[1].getPoint(alpha + beta);  
116.     if(sign(cross(a,c[0].o,b)) == sign(cross(a,p,b)) &&  
117.        sign(cross(a,c[1].o,b)) == sign(cross(a,p,b)))  
118.         add(a,b,k);  
119.     a = c[0].getPoint(alpha - beta);  
120.     b = c[1].getPoint(alpha - beta);  
121.     if(sign(cross(a,c[0].o,b)) == sign(cross(a,p,b)) &&  
122.        sign(cross(a,c[1].o,b)) == sign(cross(a,p,b)))  
123.         add(a,b,k);  
124.     return k - 2;  
125. }  
126.   
127. int main()  
128. {  
129.     int T;  
130.     scanf("%d",&T);  
131.     while(T--)  
132.     {  
133.         c[0].Input();  
134.         c[1].Input();  
135.         p.Input();  
136.         int num = Work();  
137.         printf("%d\n",num);  
138.         for(int i=0;i<num;i++)  
139.             c[i+2].Output();  
140.     }  
141.     return 0;  
142. }