<style> .reveal .slides { text-align: left; font-size:28px; } </style> # Computational Geometry ## 計算幾何 Introduction to Competitive Programming 2025/05/29 ---- - 浮點數誤差 - What is Computational Geometry - 向量與內外積 - 面積 - 各種計幾問題 - 點與線段相交 - 點與線段距離 - 線段相交 - 點與多邊形 - 凸包 - 旋轉卡尺 - 最遠點對 - 最大三角形 - 極角排序 - 皮克定理 - 計算幾何與三分搜 - 最小包覆圓 --- ## 浮點數誤差 浮點數誤差在計算幾何中是非常頭痛的問題， 如果在比賽中沒有遇到，那是代表你有好的計算幾何技巧或者是非常幸運 而如果遇到了則代表你接下來到比賽結束或者 AC 之前，都會非常的痛苦煩躁。  ---- ## 浮點數誤差 EX1 ```cpp= double x = 0.1; double y = 0.2; double z = x+y; cout << fixed << setprecision(18) << z << '\n'; ``` 結果如下，很明顯看到會有誤差  ---- ## 浮點數誤差 EX2 你原本預計只需要跑十次，但其實他永遠不會停。 ```cpp! for (double i = 0.0; i != 1.0; i += 0.1) cerr << "runing....\n"; ``` ---- ## 解決辦法 小改一下浮點數的大小判斷方式 我們設定一個很小的值（eps），若是兩個數字的差小於eps，就當做他們相等。 ```cpp! for (double i = 0.0; abs(i - 1.0) > eps ; i += 0.1) cerr << "runing....\n"; ``` ---- ## eps如何決定  eps通常根據題目的精度限制來決定，通常要設定的比題目要求的精度還要再精 ## 為什麼不設定超級小？ 可能會導致TLE，或者剛剛那個判不了相等的問題發生 ---- 選擇適當的型態是最重要的 | 型態 | 精度 | | -------- | -------- | | float | $10^{-7}$ | | double | $10^{-16}$ | | long double | $10^{-19}$ | 通常題目要求的精度範圍為 $10^{-6}\sim 10^{-9}$ 一般來說都會用 double，如果題目時限允許或精度要求更大會使用 long double ---- ## 浮點數比較大小 ```cpp! int sgn(double x) { return (x > -eps) - (x < eps); } // retrun -1 -> 負數 // retrun 0 -> 0 // retrun 1 -> 正數 ``` 觀察一下可以發現這個函數會回傳 x 的正負關係 要比較兩個浮點數的時候只需要相減看正負就好了 ---- ## 避免浮點數 - 距離 由於浮點數誤差是個很麻煩的問題 因此在寫程式時，如果能避免使用浮點數則會盡量使用整數型態 像是比較平面上的距離的時候 距離公式 $\sqrt{x^2+y^2}$ 會開根號 但實際上只需求出 $x^2+y^2$ 即可判斷兩者距離大小 使用整數型態即可 ---- ## 避免浮點數 - 分數比較 給你兩個分數 $\frac{x_1}{y_1}$ 和 $\frac{x_2}{y_2}$ 問誰比較大 ? 我們也要注意避免使用除法除完再做比較。 可以發現只需要同乘 $y_1 \cdot y_2$ 就可以比較大小關係 判斷 $x_1\times y_2$ 和 $x_2\times y_1$ 的大小 ```cpp= if (x1 * y2 > x2 * y1) ``` ---- ## 避免浮點數 - 浮點數運算 給 n 個浮點數，且都是5位小數，求總和 ```cpp! 1.12345 2.23456 3.45678 ``` 像是這樣 如果題目為小數點固定後幾位，如固定小數點後一位 做法為全部數字成 100000 倍 ( 1.12345 -> 112345 ) 轉換成整數運算 ---- ## 避免浮點數 - 浮點數輸入 如果可以的話也盡量避免使用浮點數來輸入 浮點數cin很慢，而且浮點數的輸入會有誤差 就算輸入的數字是整數也不要使用浮點數來直接輸入，會TLE！ ``` 4 1234 2345 3456 4567 ``` 請不要像下面這樣子！ 你有極大的可能會TLE ```cpp! int n; cin >> n; for (i32 i = 0; i < n; i++) { double x; cin >> x; } ``` ----   ---- ## EX 輸入為 $n$ 固定小數後 2 位的小數，問總和? input ``` 4 122.31 87.96 2222.12 23.20 ``` 改成整數輸入，整數輸出 ```cpp= int n; 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 << "." << setfill('0') << setw(2) << sum % 100 << '\n'; ``` ---- 講完以上非常重要的事情之後，可以來開始講計算幾何了 --- ## Computational Geometry ---- ## 向量 $\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}}$ |)等同於兩向量夾起來的平行四邊形面積  所以兩向量外積 / 2 等於求夾起來的三角形面積 ---- 點（向量）的儲存方式 利用重載運算子來實現剛剛所說的計算 利用template，因為有時候計算只需要整數，有時候需要浮點數 ```cpp! using ld = long double; template<class T> struct pt{ T x,y; pt(T _x,T _y):x(_x),y(_y){} pt():x(0),y(0){} pt operator * (T c){ return pt(x * c, y * c); } pt operator / (T c){ return pt(x / c, y / c); } pt operator + (pt a){ return pt(x + a.x, y + a.y); } pt operator - (pt a){ return pt(x - a.x, y - a.y); } T operator * (pt a){ return x * a.x + y * a.y; } // 內積 T operator ^ (pt a){ return x * a.y - y * a.x; } // 外積 bool operator < (pt a) const { return x < a.x || (x == a.x && y < a.y); } }; // using Pt = pt<long long>; // using Pt = pt<long double>; ``` ---- ### 多邊形面積 (鞋帶公式) 一個 $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)$ ---- ## 小技巧 如果你發現你最後需要 / 2 ，像是鞋帶公式（整數點的情況）。 你其實可以完全避免掉浮點數 be like： ```cpp! int area; if(area & 1) cout << area / 2 << ".5" << '\n'; else cout << area / 2 << '\n'; ``` ---- ## 三點共線 如果 3 個點連成一點線，則圍起來的三角形面積為 0 則直接判斷外積是否為 0 即可 ```cpp= bool collinearity(const Pt& a, const Pt& b, const Pt& c){ return sgn((b - a) ^ (c - a)) == 0; } ``` ---- ## 判斷點是否在線段上 必要條件為點與線段共線 (外積為 0) 且 兩端點往點的方向會相反 (內積為負)  ```cpp! auto ori(Pt a, Pt b, Pt c) { return (b - a) ^ (c - a); } struct Line { Pt a, b; Pt dir() const { return b - a; } }; bool PtOnSeg(Pt p, Line L) { return sgn(ori(L.a, L.b, p)) == 0 and sgn((p - L.a) * (p - L.b)) <= 0; } ``` ---- ## 點與線段最近距離 兩種 case 1. 點 $C$ 到線段 $\overline{AB}$ 的其中一個端點最近 2. 點 $C$ 到線段 $\overline{AB}$ 的垂直距離   ---- case 1. 點 $C$ 到線段 $\overline{AB}$ 的其中一個端點最近  可以分別判斷 - $\overrightarrow{AB}$ 與 $\overrightarrow{AC}$ - $\overrightarrow{BA}$ 與 $\overrightarrow{BC}$ 判斷是否為反向(內積 < 0)，若反向代表是case 1 ---- case 2. 點 $C$ 到線段 $\overline{AB}$ 的垂直距離  垂直距離等價於 $\triangle ABC$ 以 $\overline{AB}$ 為底時的高 $\triangle ABC$ 面積 * 2 除以 $\overline{AB}$ 的長度即可 ---- code ```cpp! ld PtSegDist(Pt p, Line l) { ld ans = min(abs(p - l.a), abs(p - l.b)); if (sgn(abs(l.a - l.b)) == 0) return ans; if (sgn((l.a - l.b) * (p - l.b)) < 0) return ans; // case1 if (sgn((l.b - l.a) * (p - l.a)) < 0) return ans; // case1 return min(ans, abs(ori(p, l.a, l.b)) / abs(l.a - l.b)); // case2 } ``` ---- ## 兩線段相交與否 如果相交則分成兩種 case 1. 交點在其中一個端點上 2. 交點不在端點上   ---- ## case1 交點在其中一個端點上  判斷這四種情況 - $\texttt{PtOnSeg(m.a, l)}$ - $\texttt{PtOnSeg(m.b, l)}$ - $\texttt{PtOnSeg(l.a, m)}$ - $\texttt{PtOnSeg(l.b, m)}$ 只要其中一個符合就代表有相交 ---- ## case2 交點不在端點上  呈現交叉情況 則點 A、B 會在線段 $\overline{CD}$相異側 則點 C、D 會在線段 $\overline{AB}$相異側 如果在相異側則外積相乘為負 ---- code ```cpp! auto ori(Pt a, Pt b, Pt c) { return (b - a) ^ (c - a); } bool PtOnSeg(Pt p, Line L) { return sgn(ori(L.a, L.b, p)) == 0 and sgn((p - L.a) * (p - L.b)) <= 0; } int PtSide(Pt p, Line L) { return sgn(ori(L.a, L.b, p)); } bool isInter(Line l, Line m) { if (PtOnSeg(m.a, l) or PtOnSeg(m.b, l) or PtOnSeg(l.a, m) or PtOnSeg(l.b, m)) return true; // case1 return PtSide(m.a, l) * PtSide(m.b, l) < 0 and PtSide(l.a, m) * PtSide(l.b, m) < 0; // case2 } ``` ---- ## 點在多邊形內部 逆時針枚舉多邊形的邊，判斷點有沒有在邊上 同時從點沿著 水平方向 向左向右發射一條射線，如果有一個 向量$\vec{AB}$(逆時針的邊向量) 和他有接觸，那就要判斷這個點在向量$\vec{AB}$的哪一側，左側的話就要將計數+1，右側的話就要將計數-1， ---- 我們觀察可以發現，不管是凹多邊形還是凸多邊形，計數為 0 一定在外部，否則在內部。   ---- 詳細看代碼 ```cpp! int inPoly(Pt p, const vector<Pt> &P) { const int n = P.size(); int cnt = 0; for (int i = 0; i < n; i++) { Pt a = P[i], b = P[(i + 1) % n]; if (PtOnSeg(p, {a, b})) return 1; // on edge if ((sgn(a.y - p.y) == 1) ^ (sgn(b.y - p.y) == 1)) cnt += sgn(ori(a, b, p)); } return cnt == 0 ? 0 : 2; // out, in } ``` --- ## 凸包 你有 $n$ 個點，然後你要找出一個凸多邊形可以圍住這 $n$ 個點且面積最小， 這個凸多邊形稱為凸包 你也可以想成凸包是用一條橡皮筋把牆壁上的圖釘圈起來的範圍   ---- ## Andrew's Monotone Chain (我們常用的做法) 作法為分別求出上下凸包，求完之後再和再一起， 上半凸包+下半凸包=完整的凸包  ---- ## 過程 首先sort（先按照x，再按照y） 可以發現 __最小的點__ 和 __最大的點__ 都一定在凸包上面 觀察到凸包上如果逆時針走，每次都只會是左轉的，如果有右轉就代表他不在凸包上面。 我們先升序找下凸包，再降序找上凸包 ---- 我們用一個單調棧來維護。 一旦發現即將進棧的點（P）和棧頂的兩個點（S₁、S₂，其中 S₁ 為top）行進方向向右旋轉，即外積小於 0，那S₁就不會在凸包上面 $\overrightarrow{S_2S_1}\times\overrightarrow{S_2P}<0$ 看圖很明顯可以發現，若是有p這樣的點，S₁不可能會在下凸包上面 所以S₁會被出棧，反復操作到沒有外積小於0的情況發生，再把 點（P） 進棧，這樣可以保持整個棧的單調性（都只會左轉）  ---- ## 模版 ```cpp! vector<Pt> Hull(vector<Pt> P) { sort(all(P)); P.erase(unique(all(P)), P.end()); P.insert(P.end(), P.rbegin() + 1, P.rend()); vector<Pt> stk; for (auto p : P) { auto it = stk.rbegin(); while (stk.rend() - it >= 2 and \ ori(*next(it), *it, p) <= 0 and \ (*next(it) < *it) == (*it < p)) { it++; } stk.resize(stk.rend() - it); stk.push_back(p); } stk.pop_back(); return stk; } ``` --- ## 旋轉卡尺 ---- ## 最遠點對 在平面上有 $N$ 個點，找出相距最遠的兩個點。 - $N\leq10^5$ 可以發現最遠點對會發生在凸包上面。 ---- 可以發現對於一條邊來說，其他的點和邊的距離具有凸性。  因為具有凸性，當我們逆時針枚舉每一條邊的時候，對應這條邊的最遠點也只會逆時針移動。 最遠點$P$ 離這條邊越遠，那這個點 $P$ 距離這條邊的兩個端點也會越遠。 ---- 利用雙指針的技巧，可以實現 $O(NlogN + N)$ 找到最遠點對。 - 凸包 $O(NlogN)$ - 旋轉卡尺 $O(N)$ ```cpp! auto abs2(Pt a) { return a * a; } auto RotatingCalipers(const vector<Pt> &hull) { // 最遠點對 回傳距離平方 int n = hull.size(); auto ret = abs2(hull[0]); ret = 0; if (hull.size() <= 2) return abs2(hull[0] - hull[1]); for (int i = 0, j = 2; i < n; i++) { Pt a = hull[i], b = hull[(i + 1) % n]; while(ori(hull[j], a, b) < (ori(hull[(j + 1) % n], a, b))) j = (j + 1) % n; chmax(ret, abs2(a - hull[j])); chmax(ret, abs2(b - hull[j])); } return ret; } ``` ---- ## 極角排序 給定 $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); ``` $atan2(y, x)$ 的回傳值 界在 $(-\pi, \pi]$ 排序完的順序看圖  但這個函數其實有點慢。 ---- ## 更快的做法 --- 外積判斷 用 X軸 把點分成上下兩部分， 如果兩個點屬於同一個部分，並且 $\vec{a} \times \vec{b} > 0$ 也就是 $\vec{a}$ 到 $\vec{b}$ 是逆時針，就代表 $arg(a) < arg(b)$ 細節請看code ```cpp! ld arg(Pt x) { return atan2(x.y, x.x); } bool argcmp(Pt a, Pt b) { // arg(a) < arg(b) int f = (Pt{a.y, -a.x} > Pt{} ? 1 : -1) * (a != Pt{}); int g = (Pt{b.y, -b.x} > Pt{} ? 1 : -1) * (b != Pt{}); return f == g ? (a ^ b) > 0 : f < g; } ``` ---- ## 經典例題 給平面上 $N$ 個點，每個點有點權 $w_i$ 選其中兩個點連成直線把平面切成兩半，選哪兩個點切能讓兩邊的點權差最少? - $N\le 2000$ - $|x_i|, |y_i|\le 10^9$ - $w_i\le 10^9$ - 保證不會有三點共線 ---- 很明顯我們可以找一個點當做原點極角排序一下 然後就可以照著排序好的順序枚舉另一個點和原點連成直線，轉一圈找到最小點權差 這樣是$O(NlogN + N)$ 但是每一個點都要做一次 $\rightarrow$ $O(N^2logN + N^2)$ ----   旋轉的時候每個節點最多只會變成另一半邊一次而已 旋轉的時候可以維護每個節點在哪一半 ---- ## 實作技巧 把其中一半的點旋轉$180^\circ$，並且打好標記，你會比較容易實作  $\to$  --- ## 皮克定理 (Pick Theorem) ### **$A = i + b/2 − 1$** - $A$: 多邊形的面積 - $i$: 多邊形內部的格子點 - $b$: 多邊形邊上的格子點 --- ## 三分搜 找連續凸函數的極值 ---- ## 複習 如果要找函數的最小值 ```cpp! using ld = long double; ld lo = 0, hi = mx; for (i32 t = 0; t < 200; t++) { ld ml = lo + (hi - lo) / 3.0; ld mr = hi - (hi - lo) / 3.0; if (f(ml) < f(mr)) hi = mr; else lo = ml; } ```  ---- ## [傳送帶](https://www.luogu.com.cn/problem/P2571) 在一個 2 維平面上有兩條傳送帶分別為線段 $\overline{AB}$ 和線段 $\overline{CD}$。 如果在$\overline{AB}$ 上的移動速度為 $P$，在 $\overline{CD}$ 上的移動速度為 $Q$，在平面上的移動速度 $R$。 現在想從 $A$ 點走到 $D$ 點，最少需要走多久? - 第一行包含 4 個整數，分別表示點 A 和點 B 的座標： - $A_x, A_y, B_x, B_y$ - 第二行包含 4 個整數，分別表示點 C 和點 D 的座標： - $C_x, C_y, D_x, D_y$ - 第三行包含 3 個整數： - $P, Q, R$ ---- 其實就是要決定在$\overline{AB}$ 和 $\overline{CD}$ 上走多久  ---- $x$ 是在 $\overline{AB}$ 上走的時間 $y$ 是在 $\overline{CD}$ 上走的時間 $f_1(x)$ 計算過了$x$單位時間他在$\overline{AB}$上的哪裡 $f_2(y)$ 計算過了$y$單位時間他在$\overline{CD}$上的哪裡 可以列出式子 $f(x, y) = \frac{\mathrm{dis}(f_1(x), f_2(y))}{R} + x + y$ 對這個式子做三分搜套三分搜 ---- 細節看code ```cpp! i64 p, q, r; Line AB, CD; i64 x1, y1, x2, y2; xin >> x1 >> y1 >> x2 >> y2; AB.a = Pt(x1, y1); AB.b = Pt(x2, y2); xin >> x1 >> y1 >> x2 >> y2; CD.b = Pt(x1, y1); CD.a = Pt(x2, y2); xin >> p >> q >> r; auto fff = [&](ld x, ld y) -> ld { Pt abpt = AB.a + unit(AB.dir()) * (p * x); Pt cdpt = CD.a + unit(CD.dir()) * (q * y); return (abs(abpt - cdpt) / (ld)r) + x + y; }; ld mxx = abs(AB.dir()) / (ld)p; ld mxy = abs(CD.dir()) / (ld)q; auto besty = [&](ld x) -> ld { ld lo = 0, hi = mxy; for (i32 ii = 0; ii < 200; ii++) { ld ml = lo + (hi - lo) / 3.0; ld mr = hi - (hi - lo) / 3.0; if (fff(x, ml) < fff(x, mr)) hi = mr; else lo = ml; } return fff(x, lo); }; ld lo = 0, hi = mxx; for (i32 ii = 0; ii < 200; ii++) { ld ml = lo + (hi - lo) / 3.0; ld mr = hi - (hi - lo) / 3.0; if (besty(ml) < besty(mr)) hi = mr; else lo = ml; } cout << besty(lo) << '\n'; ``` ---- ## 最小圓覆蓋 平面上給 n 個點，求出半徑最小的圓要包住所有的點。 求出圓心位置與與最小半徑 - $1\le n\le 10^5$ - $1\le |x|, |y|\le 10^9$ <img src="https://hackmd.io/_uploads/Skf2QDEzel.png" width="450" style="background:white; padding:4px;"> ---- ## 三分套三分 我們把最小圓覆蓋問題轉成這個式子。 如果能找到這個式子的最小值，那就找到了這個圓。 $$ f(x,y) \;=\; \max\limits_{1 \le i \le n} \sqrt{(x - x_i)^2 + (y - y_i)^2} $$ 並且這個函數是一個單峰函數 ---- 我們一樣用上一題的方式三分套三分 ---- 細節看code ```cpp! i32 n; cin >> n; vector<Pt> vec; for (i32 i = 0; i < n; i++) { double x, y; cin >> x >> y; vec.emplace_back(x, y); } auto f = [&](ld x, ld y) -> ld { ld ret = 0; for (i32 i = 0; i < n; i++) chmax(ret, abs(Pt(x, y) - vec[i])); return ret; }; ld cirx, ciry; auto besty = [&](ld x) -> ld { ld lo = -1e5, hi = 1e5; while ( sgn(lo - hi) < 0 ) { ld ml = lo + (hi - lo) / 3.0; ld mr = hi - (hi - lo) / 3.0; if (f(x, ml) < f(x, mr)) hi = mr; else lo = ml; } ciry = lo; return f(x, lo); }; ld lo = -1e5, hi = 1e5; while ( sgn(lo - hi) < 0 ) { ld ml = lo + (hi - lo) / 3.0; ld mr = hi - (hi - lo) / 3.0; if (besty(ml) < besty(mr)) hi = mr; else lo = ml; } cirx = lo; cout << besty(cirx) << '\n'; cout << cirx << ' ' << ciry << '\n'; ``` 但是我們很快發現這個方法不太好 因為精度很低，又有一點慢 ---- 底下這題精度要求 $1e-9$，很明顯三分搜精度不夠  所以我們講另一個辦法 ---- ## $Welzl's 算法$ ---- ## 原理 - 如果第 $i + 1$ 個點不在前 $i$ 個點的最小包覆圓上，那第 $i + 1$ 個點一定在新的圓上。 - 三個點可以確定一個圓 ---- ## 算法步驟 我們從所有未處理的點裡面隨機挑一個點加入答案集合中，然後維護這個答案集合，然後反復操作直到所有點都加進答案集合。 1. 如果我們挑的點在目前的最小包覆圓裡面，那很明顯可以直接加入，不需要再維護這個集合。 2. 如果我們挑的點不在目前的最小包覆圓裡面，那我們就需要找到新的最小包覆圓。 ---- ## 怎麼找到新的最小包覆圓？ 我們可以知道新加入的那個點一定在新的最小包覆圓的上面 那這樣我們可以枚舉答案集合裡面所有其他的點，找到另外兩個點可以構造出一個圓包住所有答案集合內的點。（三點定圓） (兩層迴圈枚舉另外兩個點) ---- 詳細看code ```cpp! Pt Center(Pt a, Pt b, Pt c) { Pt x = (a + b) / 2; Pt y = (b + c) / 2; return LineInter({x, x + rotate(b - a)}, {y, y + rotate(c - b)}); } Cir MEC(vector<Pt> P) { mt19937 rng(time(0)); shuffle(all(P), rng); Cir C = {P[0], 0.0}; for (int i = 0; i < P.size(); i++) { if (C.inside(P[i])) continue; C = {P[i], 0}; for (int j = 0; j < i; j++) { if (C.inside(P[j])) continue; C = {(P[i] + P[j]) / 2, abs(P[i] - P[j]) / 2}; for (int k = 0; k < j; k++) { if (C.inside(P[k])) continue; C.o = Center(P[i], P[j], P[k]); C.r = abs(C.o - P[i]); } } } return C; } ``` 有三圈 for 迴圈 $\cancel{O(N ^ 3)} \rightarrow {O(N)}$ 因為隨機挑選的話其實只會有$\frac 3 n$的幾率會挑中圓的支撐點（三點定圓） $O(N)\times\frac{3}{N}\times O(N)\times\frac{3}{N}\times O(N)\rightarrow O(N)$