ACM比賽計算幾何就比較重要了，高中只學了個凸包，今兒從圓的反演學起吧。
先來看一道需要用反演解決的題：HDU4773 Problem of Apollonius
題意：給定兩個圓(x1,y1,r1)、(x2,y2,r2)，它們是相離的，在這兩個圓外給定一個點p(x0,y0)。求符合條件：過點p的圓且與已知的兩個圓外切的所有圓的總數和它們的圓心坐標和半徑。
SOL：如果不知道反演，這道題的做法顯然是列出兩個方程：
然后發現這個東西根本就沒有辦法搞。所以就需要用到反演來解答此題。
反演的定義：已知一圓C，圓心為O，半徑為r，如果P與P’在過圓心O的直線上，且
，則稱P與P'關於O互為反演點。
幾條反演的重要性質：
（1）過反演中心的圓，反形(經過反演之后的圖形)為不過反演中心的直線。
（2）不過反演中心的直線，反形為過反演中心的圓。
證明的話比較顯然，對於上面這個圖，就是證明$|OM|*|OM^{'}| = |OC|*|OC^{'}|$就好了，也就是$\Delta OM^{'}C^{'} \sim \Delta OCM$。這個的話根據那個直角很容易得出。
（3）不過反演中心的圓，反形也為不過反演中心的圓，並且反演中心為這兩個互為反形的圓的位似中心。
比如上圖，O為反演中心，A和A'、B和B'都是互為反演點。證明從略。
這樣的話顯而易見的就得到了定理(4)：相切兩圓的反象仍相切，若切點恰是反演中心，則其反象為兩平行線。
另外可以推出，相離的兩圓反演(反演中心不在圓上)后仍然相離；兩圓相切，若反演中心在某圓上，則為反形為相切的直線與圓。
有了這個定理，我們就可以這樣解這個題：
首先我們把給定的p點看做反演中心，然后自定義反演半徑的大小（為保證精度，不宜過小），接着，求出題目給出的兩個圓c1，c2各自的反形c1’，c2’。然后，作出c1’和c2’的外公切線，再將公切線反演回去就得到答案所需的圓了。
為什么這樣做是對的呢？我們倒着考慮，最終的圖形為3個圓，答案圓與另外兩圓外切，如果將這三個圓反演回去，答案圓因為過反演中心，所以反形是一條直線，另外兩個題目給出的圓因為沒有過反演中心，所以反形是兩個圓，由於相切性不變，這條直線與這兩個反形圓均相切，即為公切線。
然而我們不能隨意求公切線，不僅要求的是公切線，為了防止出現內切的情況，我們要求的是“外公切線”，並且要使得P點和反形圓的圓心在外公切線的同一側。
（這里可以這樣想，反演之前如果圓A在圓B左側，反演之后圓A就會在在圓B右側，對反演成直線的圓也是同理的。）
所以接下來的問題就是要求圓和直線的反形以及兩個圓的公切線了。
先說如何求圓的反形。對於上面那個圖，我們先列出式子：
$r_1$和$OC_1$是已知的，然后消去$OC_2$，得到：$r_2 = \frac{r_1*r^{2}}{OC_1^2-r_1^2}$。
然后對於反演圓的圓心，我們這樣求：
設O的坐標為(x0,y0)，C1的坐標為(x1,y1)，C2的坐標為(x2,y2)，然后就有：
至於$OC_2$的長度，我們可以直接根據反演的定義求出$OA^{'}$的長度，然后加上r2即可。
好了，我們完成了求一個圓的反形。那么直線的反形呢？
先求出p到直線的距離，這樣根據定義就求出了反形圓的半徑。至於圓心很煩，但在本題可以通過與已知圓圓心與切點的向量的像似向量求出，詳見代碼。
好了，最后我們只剩下求兩個圓的公切線了，轉化成求切點。
方法就是求出先求出兩圓圓心連線的傾角，在求出與切線垂直的半徑與兩圓圓心連線的夾角，這樣就求出了切點與圓心連線與x軸的夾角，再用向量解決之。
好了，很煩人啊，下面給出代碼：
1 #include <cmath> 2 #include <cstdio> 3 #include <cstring> 4 #include <iostream> 5 #include <algorithm> 6 using namespace std; 7 const int LEN = 1e3 + 5; 8 const double zero = 1e-8; 9 int i, j, k, n, m, s, t, tot; 10 double R; 11 struct node { 12 double x, y; 13 void in() { 14 scanf("%lf %lf", &x, &y); 15 } 16 } p; 17 node operator + (const node &x, const node &y) { 18 return (node){x.x + y.x, x.y + y.y}; 19 } 20 node operator - (const node &x, const node &y) { 21 return (node){x.x - y.x, x.y - y.y}; 22 } 23 node operator * (const node &x, const double &t) { 24 return (node){x.x * t, x.y * t}; 25 } 26 struct circle { 27 node o; 28 double r; 29 void in() { 30 o.in(); 31 scanf("%lf", &r); 32 } 33 void out() { 34 printf("%.8f %.8f %.8f\n", o.x, o.y, r); 35 } 36 } a, b, ans[LEN]; 37 double sqr(double x) { 38 return x * x; 39 } 40 double dist(const node &x, const node &y) { 41 return sqrt(sqr(x.x - y.x) + sqr(x.y - y.y)); 42 } 43 circle inversion(const circle &a) { 44 circle res; 45 double d1 = dist(p, a.o); 46 res.r = a.r / (sqr(d1) - sqr(a.r)) * R * R; 47 double d2 = R * R / (d1 - a.r) -res.r; 48 res.o = p + (a.o - p) * (d2 / d1); 49 return res; 50 } 51 node moveit(const node &x, const double &a, const double &d) { 52 return (node){x.x + d * cos(a), x.y + d * sin(a)}; 53 } 54 double cross(const node &x, const node &y) { 55 return y.y * x.x - y.x * x.y; 56 } 57 int sign(double x) { 58 if (fabs(x) < zero) { 59 return 0; 60 } else if (x > 0) { 61 return 1; 62 } else { 63 return -1; 64 } 65 } 66 void addans(const node &x, const node &y) { 67 tot++; 68 double d = fabs(cross(y - p, x - p)) / dist(x, y); 69 ans[tot].r = sqr(R) / (2.0 * d); 70 ans[tot].o = p + (x - a.o) * (ans[tot].r / a.r); 71 } 72 void solve() { 73 a = inversion(a); 74 b = inversion(b); 75 if (a.r < b.r) { 76 swap(a, b); 77 } 78 node c = b.o - a.o; 79 double a1 = atan2(c.y, c.x); 80 double a2 = acos((a.r - b.r) / dist(a.o, b.o)); 81 node t1 = moveit(a.o, a1 + a2, a.r); 82 node t2 = moveit(b.o, a1 + a2, b.r); 83 if (sign(cross(t2 - t1, a.o - t1)) == sign(cross(t2 - t1, p - t1))) { 84 addans(t1, t2); 85 } 86 t1 = moveit(a.o, a1 - a2, a.r); 87 t2 = moveit(b.o, a1 - a2, b.r); 88 if (sign(cross(t2 - t1, a.o - t1)) == sign(cross(t2 - t1, p - t1))) { 89 addans(t1, t2); 90 } 91 } 92 int main() { 93 R = 15.0; 94 int T; 95 scanf("%d", &T); 96 while (T--) { 97 tot = 0; 98 a.in(), b.in(); 99 p.in(); 100 solve(); 101 printf("%d\n", tot); 102 for (int i = 1; i <= tot; i++) { 103 ans[i].out(); 104 } 105 } 106 return 0; 107 }
以上內容有部分搬自https://blog.csdn.net/acdreamers/article/details/16966369#comments，如涉及版權問題，請聯系2265755689@qq.com，本人將及時刪除。
花了九牛二虎之力AC了反演的第一題，再來看看這一題：HDU6097Mindis
題意：有一個圓心在原點的圓，給定圓的半徑，給定P、Q兩點坐標（PO=QO，P、Q不在圓外），取圓上一點D，求PD+QD的最小值。
SOL:做P點關於圓的反演點P',$\Delta OPD \sim \Delta ODP'$，相似比是|OP| : r，Q點同理。
極小化PD+QD可以轉化為極小化P'D+Q'D。
當P'Q'與圓有交點時，答案為兩點距離，否則最優值在中垂線上取到。
時間復雜度 O(1)。
這個題就比較簡單了，比較坑的地方是兩個點重疊和P,Q為原點的情況。
code:
1 #include <cmath> 2 #include <cstdio> 3 #include <cstring> 4 #include <iostream> 5 #include <algorithm> 6 using namespace std; 7 int i, j, k, n, m, s, t; 8 const double zero = 1e-8; 9 double R; 10 struct node { 11 double x, y; 12 void in() { 13 scanf("%lf %lf", &x, &y); 14 } 15 } a, b; 16 double sqr(const double &x) { 17 return x * x; 18 } 19 double dist(const node &x, const node &y) { 20 return sqrt(sqr(x.x - y.x) + sqr(x.y - y.y)); 21 } 22 node moveit(const node &p, const double &a, const double &d) { 23 return (node){p.x + d * cos(a), p.y + d * sin(a)}; 24 } 25 node inversion(const node &x) { 26 double d = sqr(R) / dist((node){0, 0}, x); 27 double a = atan2(x.y, x.x); 28 return moveit((node){0, 0}, a, d); 29 } 30 double cross(const node &x, const node &y) { 31 return y.y * x.x - x.y * y.x; 32 } 33 int sign(const double &x) { 34 if (fabs(x) < zero) { 35 return 0; 36 } else if (x > 0) { 37 return 1; 38 } else { 39 return -1; 40 } 41 } 42 double getdist(const node &x, const node &y) { 43 double L = dist(x, y); 44 if (fabs(L) == 0) { 45 return 2.0 * (dist((node){0, 0}, x) - R); 46 } 47 double d = fabs(cross(x, y)) / L; 48 if (sign(R - d) != -1) { 49 return L; 50 } else { 51 return 2.0 * sqrt(sqr(L / 2) + sqr(d - R)); 52 } 53 } 54 int main() { 55 int T; 56 scanf("%d", &T); 57 while (T--) { 58 scanf("%lf", &R); 59 a.in(), b.in(); 60 double t = dist((node){0, 0}, a); 61 if (sign(t) != 1) { 62 printf("%.8f\n", 2.0 * R); 63 continue; 64 } 65 t = dist((node){0, 0}, a) / R; 66 a = inversion(a); 67 b = inversion(b); 68 printf("%.8f\n", getdist(a, b) * t); 69 } 70 return 0; 71 }
接下來又是一個練習題：HDU6158The Designer
題意：
給出最大的兩個圓，求帶標號的那些小圓的前N個圓的面積和。
SOL:
將兩個大圓的切點作為反演中心反演兩個大圓成兩條直線，根據反演前后相切性不變的道理，直接將兩直線間夾着的那些圓反演回去求個面積就好了。
另外此題N很大，所以當一定時候面積很小時可以直接忽略。
1 #include <cmath> 2 #include <cstdio> 3 #include <cstring> 4 #include <iostream> 5 #include <algorithm> 6 using namespace std; 7 const double pi = acos(-1); 8 const double zero = 1e-14; 9 int i, j, k, n, m, s, t; 10 double r1, r2, R, x1, x2, ans, r; 11 struct node { 12 double x, y; 13 } p; 14 struct circle { 15 node o; 16 double r; 17 } now; 18 double sqr(const double &x) { 19 return x * x; 20 } 21 int sign(const double &x) { 22 if (fabs(x) < zero) { 23 return 0; 24 } else if (x > 0) { 25 return 1; 26 } else { 27 return -1; 28 } 29 } 30 double dist(const node &x, const node &y) { 31 return sqrt(sqr(x.x - y.x) + sqr(x.y - y.y)); 32 } 33 circle inversion(const node &p, const double &r) { 34 circle res; 35 res.r = r / (sqr(dist(p, (node){0, 0})) - sqr(r)) * sqr(R); 36 res.o = (node){0, 0}; 37 return res; 38 } 39 int main() { 40 R = 20.0; 41 int T; 42 scanf("%d", &T); 43 while (T--) { 44 scanf("%lf %lf", &r1, &r2); 45 if (r1 < r2) { 46 swap(r1, r2); 47 } 48 scanf("%d", &n); 49 x1 = sqr(R) / (2.0 * r1); 50 x2 = sqr(R) / (2.0 * r2); 51 r = (x2 - x1) / 2.0; 52 p = (node){(x1 + x2) / 2.0 , 0}; 53 now = inversion(p, r); 54 ans = pi * sqr(now.r); 55 for (int i = 2; i <= n; i++) { 56 if ((i & 1) == 0) { 57 p.y += 2.0 * r; 58 now = inversion(p, r); 59 } 60 ans += pi * sqr(now.r); 61 if (sign(pi * sqr(now.r)) != 1) { 62 break; 63 } 64 } 65 printf("%.5f\n", ans); 66 } 67 return 0; 68 }
圖片來自https://blog.csdn.net/coder__yang/article/details/77489327，如有侵權請聯系2265755689@qq.com,本人將及時刪除。
好了，大概反演就是這樣子了。以后等遇到題目是在多練練吧。