<style> .reveal .slides { text-align: left; font-size:28px; } </style> # Computational Geometry ## 計算幾何 Introduction to Competitive Programming 2023/05/24 ---- - 浮點數誤差 - 向量與內外積 - 面積 - 各種計幾問題 - 點與線段相交 - 點與線段距離 - 線段相交 - 點與多邊形 - 凸包 - 旋轉卡尺 - 最遠點對 - 極角排序 - 計算幾何與三分搜 - 最小包覆圓 --- ## 浮點數誤差 浮點數誤差在計算幾何中是非常頭痛的問題， 如果在比賽中沒有遇到，那是代表你有好的計算幾何技巧或者是非常幸運 而如果遇到了則代表你接下來到比賽結束或者 AC 之前，都會非常的痛苦煩躁。 ---- 本周的課程有一半的時間都在教你怎麼避免浮點數 ~~可能不會有幾何的內容~~ ---- ## 浮點數誤差的例子 ```cpp= double x = 0.1; double y = 0.2; double z = x+y; cout << fixed << setprecision(18) << z << endl; ```  ---- ## 例子 2 可以嘗試跑以下程式碼 ```cpp= double x = 0.0; while(x != 1.0){ x += 0.1; } ``` 理論上在迴圈裡跑加 10 次就會結束了 但實際上執行會因為精確度誤差造成程式會停不下來 為了解決浮點數誤差的問題，通常會設一個很小的數字例如$10^{-9}$ 在比較大小關係時，如果在誤差範圍內則視為相等 ---- 改成當前的值與目標若差值小於 EPS 則視為相同 ```cpp= #define EPS 1e-9 double x=0.0; while(abs(x-1.0) > EPS){ x += 0.1; } ``` ---- 而 EPS 大小要設為多少呢 ? ---- EPS 大小通常會根據視不同題目而定  不同題目都有不同的精度限制 如果題目需要精確到小數點至少 9 位數，保險起見會多計算一位 ```cpp // 控制輸出到小數點後 10 位 cout << fixed << setprecision(10) << x << endl; ``` ---- 浮點數誤差而也不是越細越好，如果太計算的太細可能會變成... <span> TLE <!-- .element: class="fragment" data-fragment-index="1" --></span> ---- 選擇適當的型態是最重要的 | 型態 | 精度 | | -------- | -------- | | float | $10^{-7}$ | | double | $10^{-16}$ | | long double | $10^{-19}$ | 通常題目要求的精度範圍為 $10^{-6}\sim 10^{-9}$ 一般來說都會用 double，如果題目時限允許或精度要求更大會使用 long double ---- 浮點數比較大小 比較兩個浮點數 $a, b$ 的大小 | 判斷式 | 寫法 | | -------- | -------- | | a == b | abs(a-b) < EPS | | a != b | abs(a-b) > EPS | | a < b | a - b < -EPS | | a <= b | a - b < EPS | | a > b | a - b > EPS | | a >= b | a - b > -EPS | ---- ## 避免浮點數 - 距離遠近 由於浮點數誤差是個很麻煩的問題 因此在寫計算幾何時，如果能避免使用浮點數則會盡量使用整數型態 像是判斷距離，當前點到點 x, y 哪個比較近 距離公式 $\sqrt{x^2+y^2}$ 會開根號 但實際上只需求出 $x^2+y^2$ 即可判斷兩者距離大小 只需使用 long long 型態即可 ---- ## 避免浮點數 - 判斷兩個分數的大小 給你兩個分數 $\frac{x_1}{y_1}$ 和 $\frac{x_2}{y_2}$ 問誰比較大 ? 如果直接除會有機會遇到浮點數問題， ```cpp= if((double)x1/y1 > (double) x2/y2) ``` 判斷 $\frac{x_1}{y_1}$ 和 $\frac{x_2}{y_2}$ 等價於判斷 $x_1\times y_2$ 和 $x_2\times y_1$ 的大小 ```cpp= if(x1*y2 > x2*y1) ``` ---- ## 避免浮點數 - 浮點數運算 給你 n 個科目的 score，問你平均的 GPA 是多少?  如果題目為小數點固定後幾位，如本題為固定小數點後一位 做法為全部數字成 10 倍 ( 4.3 -> 43, 3.0 -> 30) 即可變成整數運算 ----  ```cpp long long getGPA(ll x){ if(x >= 90) return 43; if(x >= 85 ) return 40; if(x >= 80) return 37; if(x >= 77) return 33; if(x >= 73) return 30; if(x >= 70) return 27; if(x >= 67) return 23; if(x >= 63) return 20; if(x >= 60) return 17; return 0; } ``` 改成整數運算，能避免浮點數則避免 ---- ## 避免浮點數 - 浮點數輸出入 輸入為 $n$ 固定小數後 2 位的小數，問總和? input ``` 4 122.31 87.96 2222.12 23.20 ``` 改成整數輸入，整數輸出 ```cpp= void solve(){ cin >> n; int x, y; char ch; int sum = 0; for(int i = 0; i < n; i++){ cin >> x >> ch >> y; sum += x*100+y; } cout << sum/100 << "." << sum%100 << endl; ``` ---- 使用整數輸入除了避免浮點數誤差外，同時也會讓輸出入常數小很多 浮點數的輸出入非常的花時間 如果輸入為整數點但過程中需要用到浮點數建議以整數輸入後再轉型 ```cpp= struct Pt{ double x, y; Pt(){} Pt(double _x, double _y):x(_x), y(_y){} } p; int main(){ int a, b; cin >> a >> b; p = Pt(a, b); } ``` Test Problem: [CF 1496C](https://codeforces.com/contest/1496/problem/C) Relate Blog: https://codeforces.com/blog/entry/88549 --- ## 座標與向量 ---- ### 歐基里德距離 $d(p_1, p_2) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$ ### 曼哈頓距離 $d(p_1, p_2) = |x_1-x_2|+|y_1-y_2|$ ---- ## 向量 $\vec{v} = (a, b)$，代表往 +x 軸走單位 $a$，往 +y 軸走單位 $b$ 長度 $|\vec{v}|$ 為 $\sqrt{a^2+b^2}$ ### 向量四則運算 假設有兩向量 $v_{1} = (x_1, y_1), v_{2} = (x_2, y_2)$ - $v_1 + v_2 = (x_1 + x_2, y_1 + y_2)$ - $v_1 - v_2 = (x_1 - x_2, y_1 - y_2)$ - $v_1 × k = (x_1 × k, y_1 × k)$ (向量伸長成 $k$ 倍) - $v_1 ÷ k = (x_1 ÷ k, y_1 ÷ k)$ (向量縮短成 $k$ 倍) ---- ### 內積 $\vec{v_{1}} =(x_1, y_1), \vec{v_{2}}=(x_2, y_2)$ $\vec{v_{1}} \cdot \vec{v_{2}} = x_1x_2 + y_1y_2 = |\vec{v_1}||\vec{v_2}| cos \theta$ $(\theta$ 為兩向量夾角) 兩向量同向，則內積為正 兩向量反向，則內積為負 兩向量垂直，則內積為 0 ---- ### 外積 $\vec{v_{1}} =(x_1, y_1), \vec{v_{2}}=(x_2, y_2)$ $\vec{v_{1}} \times \vec{v_{2}} = x_1y_2 -x_2y_1=|\vec{v_1}||\vec{v_2}| sin \theta$ $(\theta$ 為兩向量夾角) $\vec{v_{1}} \times \vec{v_{2}} = - \vec{v_{2}} \times \vec{v_{1}}$ 具有負交換律 如果 $\vec{v_{1}}$ 到 $\vec{v_{2}}$ 為逆時鐘方向，則外積為正 如果 $\vec{v_{1}}$ 到 $\vec{v_{2}}$ 為順時鐘方向，則外積為負 兩向量平行則外積為 0 同時外積取絕對值(| $\vec{v_{1}} \times \vec{v_{2}}$ |)等同於兩向量夾起來的平行四邊形面積  ---- ### 儲存方法 #### 點 (儲存座標 x, y 表示) ```cpp #define ld long double struct Pt{ ll x, y; //ld x, y; }; ``` #### 線段 (儲存線段兩端點儲存) ```cpp struct Line{ Pt st, ed; }; ``` #### 圓 (儲存圓心與半徑) ```cpp= struct Circle{ Pt o; // 圓心 ld r; // 半徑 }; ``` #### 多邊形 ( 用 vector 儲存所有點) ```cpp= struct poly{ int n; // n 邊形 vector<Pt> pts; }; ``` ---- ### 運算子多載 為了方便實作，通常會進行運算子重載以實作向量的四則運算以及內外積運算 ```cpp= Pt plus(const Pt& a, const Pt& b){ return Pt(a.x+b.x, a.y+b.y); } int main(){ Pt a(1, 3), b(3, 19); Pt c = plus(a, b); } ``` 使用運算子重載 ```cpp= struct Pt { ll x, y; Pt(ll _x=0, ll _y=0):x(_x), y(_y) {} Pt operator+(const Pt &a){ return Pt(x+a.x, y+a.y); } }; int main(){ Pt a(1, 3), b(3, 19); Pt c = a+b; } ``` ---- ## default code ```cpp= #define ld long double struct Pt { ld x, y; Pt(ld _x=0, ld _y=0):x(_x), y(_y) {} Pt operator+(const Pt &a){ return Pt(x+a.x, y+a.y); } Pt operator-(const Pt &a){ return Pt(x-a.x, y-a.y); } Pt operator*(const ld &a){ return Pt(x*a, y*a); } Pt operator/(const ld &a){ return Pt(x/a, y/a); } ld operator*(const Pt &a){ //計算幾何程式碼中內積通常用*表示 return x*a.x + y*a.y; } ld operator^(const Pt &a){ //計算幾何程式碼中外積通常用^表示 return x*a.y - y*a.x; } bool operator<(const Pt &a) const { // 判斷兩點座標 先比 x 再比 y return x < a.x || (x == a.x && y < a.y); } }; ``` --- ### 三角形面積 - 海龍公式 $s = \frac{a+b+c}{2}$ (a,b,c 為三角形三邊長) $\triangle abc = \sqrt{s(s-a)(s-b)(s-c)}$ - 外積求法 |$\vec{v_1} \times \vec{v_2}$| 為所夾的平行四邊形面積 而平行四邊形的一半為三角形面積 因此 $\frac{|\vec{v_1} \times \vec{v_2}|}{2}$ 為三角形面積 ---- ### 多邊形面積 (鞋帶公式) 一個 $n$ 個點的多邊形($P_1, P_2,... P_n$) 可以拆成很多三角形，再分別套用外積公式 ($x_1y_2 -x_2y_1$) $\frac{1}{2}|\sum\limits_{i=1}^{n} \overrightarrow {OP_{i}} \times \overrightarrow {OP_{i+1}}|$ 複雜度 $O(N)$ ---- 使用外積公式 $(x_1y_2 -x_2y_1)$ 求面積會發現 $\frac{1}{2}|\sum\limits_{i=1}^{n} \overrightarrow {OP_{i}} \times \overrightarrow {OP_{i+1}}|$ 如果所有點都為整數點，則除了最後乘 $\frac{1}{2}$ 使得面積有可能為小數，否則在用外積公式時都為整數運算 因此整數點的面積不是整數就是 .5 結尾 過程中可以完全不用用到小數可以避免浮點數誤差 ---- $\frac{1}{2}|\sum\limits_{i=1}^{n} \overrightarrow {OP_{i}} \times \overrightarrow {OP_{i+1}}|$ 當全部點都是整數點時，只有外積總合為奇數面積才會有 .5，否則為整數 ```cpp= Pt arr[N]; void solve(){ cin >> n; for(int i = 0; i < n; i++) cin >> arr[i].x >> arr[i].y; for(int i = 0; i < n; i++){ area += arr[i].x * arr[(i+1)%n].y - arr[(i+1)%n].x*arr[i].y; } area = abs(area); if(area&1) cout << area/2 << ".5" << endl; else cout << area/2 << endl; } ``` --- ### 是否三點共線 如果 3 個點連成一點線，則圍起來的三角形面積為 0 則直接判斷外積是否為 0 即可 整數點判斷 ```cpp= bool collinearity(const Pt& a, const Pt& b, const Pt& c){ return (b-a)^(c-a) == 0; } ``` 小數點判斷 ```cpp= bool collinearity(const Pt& a, const Pt& b, const Pt& c){ return (b-a)^(c-a) < EPS; } ``` ---- ### 判斷點是否在線段上 必要條件為點與線段共線(外積為 0) 且 點往兩端點的方向會相反(內積為負)  整數點判斷 ```cpp= bool inLine(const Pt& p, const Line& li){ return collinearity(li.st, li.ed, p) && (li.st-p)*(li.ed-p) < 0; } ``` 小數點判斷 ```cpp= bool inLine(const Pt& p, const Line& li){ return collinearity(li.st, li.ed, p) && (li.st-p)*(li.ed-p) < EPS; } ``` ---- ## 點與線段最近距離 兩種 case 1. 點 $C$ 到線段 $\overline{AB}$ 的其中一個端點最近 2. 點 $C$ 到線段 $\overline{AB}$ 的垂直距離   ---- case 1. 點 $C$ 到線段 $\overline{AB}$ 的其中一個端點最近  可以分別判斷 - $\overrightarrow{AB}$ 與 $\overrightarrow{AC}$ - $\overrightarrow{BA}$ 與 $\overrightarrow{BC}$ 使用內積是否為反向，如其中一個為反向則其中一個端點最近 ---- case 2. 點 $C$ 到線段 $\overline{AB}$ 的垂直距離  垂直距離等價於 $\triangle ABC$ 以 $\overline{AB}$ 為底時的高 直接使用外積公式算出面積後除以 $\overline{AB}$ 的長度即可 ---- ### 判斷兩線段是否相交 判斷兩線段相交相信大家在高中都有學過，解方程式即可 但過程中會需要處理浮點數 而如果全部點都是整數點，會使用幾何(外積)的做法來作 ---- 如果相交則分成兩種 case 1. 交點在其中一個端點上 2. 交點不在端點上   ---- #### 第一種 case 交點在其中一個端點上  - $\texttt{inLine(line1.st, line2)}$ - $\texttt{inLine(line1.ed, line2)}$ - $\texttt{inLine(line2.st, line1)}$ - $\texttt{inLine(line2.ed, line1)}$ 判斷以上四種情況皆符合才是相交 ---- #### 第二種 case :交點不在端點上  呈現交叉情況 則點 A、B 會在線段 $\overline{CD}$相異側 則點 C、D 會在線段 $\overline{AB}$相異側 如果在相異側則外積相乘為負 $(\overrightarrow{AB}\times\overrightarrow{AC})(\overrightarrow{AB}\times\overrightarrow{AD}) <0$ $(\overrightarrow{CD}\times\overrightarrow{CA})(\overrightarrow{CD}\times\overrightarrow{CB}) <0$ ---- ## 點在多邊形內部 給一個多邊形 $P$，與點 $A$，問 $A$ 是否在 $P$ 內 且保證點會照逆時針或順時針方向給定 ---- 兩種判法 1. 射線法 2. 外積判斷 - 僅限於凸多邊形 ---- ### 1. 射線法 若點在多邊形內，則隨機選一個方向的射線出現會碰到**奇數**次邊 而如果碰到多邊形的點，**如果射線碰到多邊形的點則重選** (需要特判點是否在多邊形的邊或頂點上)   ---- ### 2. 外積判斷(僅限凸多邊形) 若點在凸多邊形內部，則點與多邊形上所有頂點 $P_1,P_2,...P_n$ 依序對於所有 $\overrightarrow{AP_i} \times \overrightarrow{AP_{i+1}}$ 皆為同號(全部都是正或者全部都是負的) 若有異號的情況則代表點在凸包外 若出現一個 0 則代表點在邊上，出現兩個 0 則代表點在頂點上 ---- 如果在多邊形內，則沿著多邊形每個點順時針 $\overrightarrow{AP_i} \times \overrightarrow{AP_{i+1}}$ 皆為同號(全部都是正或者全部都是負的)  ---- 若外積有異號的情況則代表點在凸包外 (順時針變成逆時針 或者相反)  --- ## 凸包 你有 $n$ 個點，然後你要找出一個凸多邊形可以圍住這 $n$ 個點且面積最小， 這個凸多邊形稱為凸包 你也可以想成凸包是用一條橡皮筋把牆壁上的圖釘圈起來的範圍  ----  ---- ## Andrew's Monotone Chain (競賽常用的做法) 作法為分別求出上下凸包，求完之後再和再一起， 上半凸包+下半凸包=完整的凸包  ---- 先求下凸包 1. 先將全部點照 $x$ 大小排序,再排 $y$ 大小 2. 依序嘗試將點加入凸包判斷是否合法 3. 不合法刪掉前面最後一個點，合法則將點加進凸包 ---- 如何判斷合法 首先先觀察凸包，若我們從左到右依序看點，會發現下一個向量都是在前一個向量的逆時鐘方向 像是 $\overrightarrow{P_2 P_3}$ 會在 $\overrightarrow{P_2 P_4}$ 的逆時鐘方向  ----  而 $\overrightarrow{P_{i} P_{i+1}}$ 若是在$\overrightarrow{P_{i} P_{i+2}}$ 的逆時鐘方向，則$\overrightarrow{P_{i} P_{i+1}} \times \overrightarrow{P_{i} P_{i+2}} > 0$ 因此如果 $\overrightarrow{P_{i} P_{i+1}} \times \overrightarrow{P_{i} P_{i+2}} \le 0$ 代表前一個點會在凸包裡面， 而不會在凸包上，就可以把前一個點刪掉，再把新的點加進去凸包裡面 ---- 而實作可以用一個stack去維護， 一開始先把前兩個點直接加進去stack之後每一個點$P_i$都去判斷， 向量是不是往逆時針方向 ( $\overrightarrow{P_{i-2} P_{i-1}} \times \overrightarrow{P_{i-2} P_{i}} > 0$ ) 若非法則把stack的top刪掉，一直刪到合法為止， 再把新的點加進去凸包 ---- ```cpp= struct Pt{ int x,y; Pt(){} Pt(int _x,int _y){ x=_x,y=_y; } Pt operator-(const Pt &a){ return Pt(x-a.x, y-a.y); } bool operator<(const Pt &a) const { return x < a.x || (x == a.x && y < a.y); } friend int cross(const Pt& o,const Pt& a,const Pt& b){ Pt lhs = o-a, rhs = o-b; return lhs.x*rhs.y - lhs.y*rhs.x; } }; vector<Pt> convex_hull(vector<Pt> hull){ sort(hull.begin(),hull.end()); int top=0; vector<Pt> stk; for(int i=0;i<hull.size();i++){ while(top>=2&&cross(stk[top-2],stk[top-1],hull[i])<=0) stk.pop_back(),top--; stk.push_back(hull[i]); top++; } for(int i=hull.size()-2,t=top+1;i>=0;i--){ while(top>=t&&cross(stk[top-2],stk[top-1],hull[i])<=0) stk.pop_back(),top--; stk.push_back(hull[i]); top++; } stk.pop_back(); return stk; } ``` 排序 $O(NlgN)$ + 找凸包 $O(N)$ 複雜度 $O(NlgN)$ --- ## 旋轉卡尺 * 最遠點對 * 最大三角形 * 最大四角形 ---- ### 最遠點對 一群點當中，距離最遠的兩個點叫作「最遠點對」 觀察一下可以發現，因為凸包上所有的點是包圍所有點的多邊形，因此最遠的兩點一定是凸包的其中兩點。 因此我們在做最遠點對的時候要先做凸包，接著可以用旋轉卡尺的概念去找最遠點對 ---- ## 旋轉卡尺  以兩條水平線夾緊多邊形，旋轉 360 度， 能讓兩條水平線最遠的距離為最遠點對距離 ----  實作上會使用雙指針來維護，每次把指針 i 移動凸包的下一個點， 指針 j 則移動到距離指針 i 最遠的點 ---- ## 找最遠的點 距離指針 i 最遠的點 j 可以用距離來判斷 點 j 往兩邊的點距離都會比到點 i 更短，因此實作上我們只需要判斷當前點 j 與相鄰點到點 i 的距離即可  ---- 如果下一個點 j+1 距離點 i 能更遠或相同，則指針 j 往下移動 ```cpp= double FarthestPair(vector<Pt> arr){ double ret=0; for(int i = 0, j = i+1; i<arr.size(); i++){ while(distance(arr[i], arr[j]) <= distance(arr[i], arr[(j+1)%arr.size())] ){ j = (j+1) % arr.size(); } ret = max(ret, distance(arr[i],arr[j])); } return ret; } ``` 旋轉卡尺 指針 $i$ 繞一圈 $O(N)$， 窮舉最遠點的指針 $j$ 最多跟著指針 $i$ 繞一圈， 均攤複雜度為 $O(N)$ ---- 總複雜度 $O(N\log N)$ 先求出凸包 $O(N\log N)$ + 旋轉卡尺 $O(N)$ ---- ### 最大三角形 給一堆點，求其中三個點圍成的三角形面積最大 ---- 一樣用旋轉卡尺，窮舉每個點當定點，第二個點每次向下轉一格， 而第三個點就跟著轉到找到最大為止 複雜度 $O(N^2)$ 窮舉點對 ---- ### 最大四角形 給一堆點，求其中四個點圍成的四角形面積最大 作法：把四角形拆成兩個三角形，兩邊分別旋轉卡尺 ----  複雜度 $O(N^2)$ --- ## 極角排序 給定 $N$ 個點，沿著原點 順/逆時針排序  ---- ## 做法1 --- 照角度排序 可以使用內建函式 atan2(y, x) 會回傳向量的弧度 弧度 $\times \frac{180}{π} =$ 角度 ```cpp= bool cmp(const Pt& lhs, const Pt rhs){ return atan2(lhs.y, lhs.x) < atan2(rhs.y, rhs.x); } sort(P.begin(), P.end(), cmp); ``` 但浮點數運算，超級慢 ---- ## 作法2 --- 用外積排序 把所有點依照 (0, 0) 分成左右兩半邊，分別對於兩半邊 分別依照外積排序，比較兩個向量的先後順序 ```cpp= bool cmp(const Pt& lhs, const Pt rhs){ if((lhs < Pt(0, 0)) ^ (rhs < Pt(0, 0))) return (lhs < Pt(0, 0)) < (rhs < Pt(0, 0)); return (lhs ^ rhs) > 0; } // 從 270 度開始逆時針排序 sort(P.begin(), P.end(), cmp); ```  ---- ## 經典例題 給平面上 $N$ 個點，保證不會有三點共線 每個點有點權 $w_i$ 選其中兩個點連成直線把平面切成兩半，選哪兩個點切能讓兩邊的點權差最少? - $N\le 2000$ - $|x_i|, |y_i|\le 10^9$ - $w_i\le 10^9$ 會發現 $N$ 才 2000 ---- 如果窮舉一個點，當要切的其中一個 另一個點再照某種順序窮舉計算兩邊的權重和  會發現我們可以用極角排序處理順序，接著再把點分成兩群 ---- 在過程中記錄，當前點所對應的分界線 用另一個指針紀錄  ---- 而由於要記錄兩個對角線相應的指針很麻煩 我們可以把其中一半的點映射到另外一邊 $(x, y)\to (-x, -y)$ $\forall (x, y) > (0, 0)$  $\to$  如此一來就不用另一半的雙指針了，只需根據顏色判斷為哪一側的點 ---- 窮舉每個點為中心 $O(N)$ 對於窮舉的中心點進行極角排序 $O(N\log N)$ 複雜度 $O(N^2\log N)$ --- ## 計算幾何與三分搜 計算幾何中很多求最小值的題目都有三分搜的性質，通常函數成凸/凹性函式 當遇到幾何的題目，如果沒有直接的想法，不妨想想看有沒有三分搜的性質 或許或讓題目變簡單很多 ---- ### 點與線段的最短距離 線段的每個位置與點之間的距離，會從最近點開始呈現一個凸函式 因此我們可以用三分搜來找出函式的最低點，為點 $A$ 與線段 $\overline{P_0P_1}$ 最近距離 線段 $S = \overline{P_0P_1}$ 我們可以用以下方程式來表示 $S(t) = t\times P_0 + (1-t)\times P_1$, $t\in [0, 1]$ 在 [0, 1] 之間找到的最小值 $S(t)$ 為最近距離 ---- 實作上，左右界分別為 0, 1， 搜出來的值代表用了幾 t% 的 $P_0$, 以及 (1-t)% 的 $P_1$ ```cpp= #define EPS 1e-8 Pt P0, P1, A; double check(double t){ // 回傳點 A, B 的距離 double B = P0 * t + P1 * (1.0-t); return dis(A, B); } double PointAndSegmentDistance(){ double l = 0.0, r = 1.0; while(r - l > EPS){ double ml = (l+l+r)/3.0, mr = (l+r+r)/3.0; if(check(ml) < check(mr)) r = mr; else l = ml; } return check(l); } ``` ---- ### 兩線段的最近距離 兩條線段如果相交，則最近距離為 0 否則最近距離為其中一條線段 A 的一個點到線段 B 的最近距離 變回原本的問題 $\to$ 點與線段的最短距離 ---- ## 最小包覆圓 平面上給 n 個點，在平面上選一個位置為圓心畫一個圓，圓要包住所有的點，半徑最小是多少 求出圓心位置與與最小半徑 - $1\le n\le 10^5$ - $1\le |x|, |y|\le 10^9$ ---- ### 簡單化題目 平面上給 n 個點，在平面上 x = 0 上畫一個圓，圓要包住所有的點， 問 y 要選在哪裡能使得包覆圓半徑最小 ---- ### 再簡單一點 平面上給 1 個點，在平面上 x = 0 上畫一個圓，圓要包住所有的點， 問 y 要選在哪裡能使得包覆圓半徑最小 ---- 可以發現，對大圓來說，要覆蓋住一個點，半徑一定是點與線的最短距離，圓如果向上或是向下移動時，圓的半徑都要增加才能覆蓋住整個點，就會變成一個很典型的單峰函數。  ---- ### 回到前一題 (平面上給 n 個點，在平面上 x = 0 上畫一個圓) 可以發現在這題就變成多個單峰函數的疊合，而疊合後的單峰函數還會是單峰函數，因此我們就可以直接在 x = 0 上三分搜，搜出讓半徑最小的 y。 ---- 搜出能讓半徑最小的 y，每次找出對於當前搜尋值 (0, y) 的最小半徑是多少 搜尋左右界則設成圓心可能出現的位置，也就是題目範圍的左右界 ```cpp= Pt arr[MXN]; double check(double y) { double cmax = 0; for(int i = 0; i < n; i++) { cmax = max(cmax,(arr[i].x - 0) * (arr[i].x - 0) + (arr[i].y - y) * (arr[i].y - y)); }// 過程中回傳距離^2 避免不必要的根號運算 return cmax; } double l = -1e9, r = 1e9; while(r - l > EPS) { double ml = (l+l+r) / 3, mr = (l+r+r) / 3; if (check(ml) < check(mr)) r = mr; else l = ml; } ``` ---- 回到原問題 平面上給 n 個點，在平面上選一個位置為圓心畫一個圓，圓要包住所有的點，半徑最小是多少 由於兩個維度都具有單調性，因此我們可以先搜其中一個維度 x 的最低點， 每次搜到一個 x 後，再去搜當前 x 所對應到的 y 值最小半徑是多少 ---- ```cpp [|18|19-23|10-17|2-9] Pt arr[MXN]; double checky(double x, double y) { double cmax = 0; for(int i = 0; i < n; i++) { cmax = max(cmax,(arr[i].x - x) * (arr[i].x - x) + (arr[i].y - y) * (arr[i].y - y)); }// 過程中回傳距離^2 避免不必要的根號運算 return cmax; } double checkx(double x){ double yl = -1e9, yr = 1e9; while(yr - yl > EPS) { double ml = (yl+yl+yr) / 3, mr = (yl+yr+yr) / 3; if (checky(x, ml) < checky(x, mr)) yr = mr; else yl = ml; } } double xl = -1e9, xr = 1e9; while(xr - xl > EPS) { double ml = (xl+xl+xr) / 3, mr = (xl+xr+xr) / 3; if (checkx(ml) < checkx(mr)) xr = mr; else xl = ml; } ``` 先搜 x $O(\log X)$, 再搜 y $O(\log Y)$, 再跑過每個點找此圓心對應的半徑 $O(N)$ 總複雜度 $O(N\log^2 N)$ ---- 而如果題目拓展到 3 維，求最小覆蓋球，求出圓心 $(x, y, z)$ 則改成搜三次，先搜 x, 再搜 y, 再搜 z。 複雜度 $O(N\log^3 N)$ --- ### 最小包覆圓 --- O(N) 作法 透過觀察圓的性質，使用機率分析複雜度 達到均攤 $O(N)$ ---- 要形成 N 個點的最小包覆圓，則圓上至少有兩個點  ---- 而也有可能三點或以上，但三個點可以決定一個圓 因此我們只需要窮舉最多三個點  ---- 而做法非常暴力，一開始把第一個點當成圓心，每次逐一加入一個新的點 $i$ 看是否新的點 $i$，在最小包覆圓外 1. 如果在圓內，則無視當前點 2. 如果在圓外，則此點必在前 i 個點的最小包覆圓上 會發現我們只需要處理第二個 case 就好，更新最小包覆圓 ---- 在已知點 $i$ 必在前 $i$ 個點的最小包覆圓上 使用點 $i$ 組成新的最小包覆圓，窮舉前面的 $i-1$ 個點，找出圓上的第二個點 $j$ 找到後在檢查是否剩下的點都在 $i, j$ 所組成的圓上 如果有點 $k$ 沒有符合，則將 $k$ 加入最小包覆圓，形成圓 $i, j, k$ ```cpp [|3-8|9-12|13-16|] Pt p[MXN]; pair<Pt,double> solve(){ double r2=0; Pt cen = p[0]; for (int i=0; i<n; i++){ if (dis(cen-p[i]) <= r2) continue; // 如果當前點在園內，則無視 cen = p[i]; // 否則當前點 i 必在圓上 r2 = 0; // 使用點 i 構成一個新的圓 for (int j=0; j<i; j++){ // 窮舉第二個點 j if (dis(cen-p[j]) <= r2) continue; cen=Pt((p[i].X+p[j].X)/2,(p[i].Y+p[j].Y)/2); // 找到第二個點 j 形成圓 r2 = dis(cen-p[j]); for (int k=0; k<j; k++){ if (dis(cen-p[k]) <= r2) continue; cen = center(p[i],p[j],p[k]); // 找出 ijk 形成的圓 r2 = dis(cen-p[k]); } } } return {cen,sqrt(r2)}; } ``` ---- 看起來複雜度為 $O(n^3)$ 但如果將點的順序打亂後，平均每次加入新點 i 時 新點不在前 i-1 個點組成的圓內的機率為 $\frac{2-1}{i}$, $\frac{3-1}{i}$ 也就是第 i 個點構成前 i 個點的最小包覆圓機率 當 $n$ 到很大 如 $10^5$ 時，添加一個新的點的複雜度為 O(1) 總共 $n$ 個點加入，平均複雜度為 $O(n)$ ---- 因此要求給定的 $n$ 個點最小包覆圓時，通常會先將點的順序打亂後再加入 ```cpp [|1,4|] mt19937 gen(chrono::steady_clock::now().time_since_epoch().count()); Pt p[MXN]; pair<Pt,double> solve(){ shuffle(p,p+n, gen); double r2=0; Pt cen = p[0]; for (int i=0; i<n; i++){ if (dis(cen-p[i]) <= r2) continue; // 如果當前點在園內，則無視 cen = p[i]; // 否則當前點 i 必在圓上 r2 = 0; // 使用點 i 構成一個新的圓 for (int j=0; j<i; j++){ // 窮舉第二個點 j if (dis(cen-p[j]) <= r2) continue; cen=Pt((p[i].X+p[j].X)/2,(p[i].Y+p[j].Y)/2); // 找到第二個點 j 形成圓 r2 = dis(cen-p[j]); for (int k=0; k<j; k++){ if (dis(cen-p[k]) <= r2) continue; cen = center(p[i],p[j],p[k]); // 找出 ijk 形成的圓 r2 = dis(cen-p[k]); } } } return {cen,sqrt(r2)}; } ``` --- ### Question Time and Practice [Homework Link](https://vjudge.net/contest/559868) 提交期限到下下星期六 6/10 23:59 截止 ---- ### 膜拜時間 交大有位幾何超人在去年 World Final 上拿到兩題幾何首殺，其中一題只有他們跟第一名過