第一次在Virtual Judge上刷到UVa 11355 Cool Points這道題時我一度以為是個腦筋急轉(zhuǎn)彎點還能分“酷”和“不酷”等看清題面才發(fā)現(xiàn)這是一道非常典型的計算幾何應(yīng)用題——給定一個凸多邊形作為隨機撒點的范圍再給若干個圓問隨機選一個點落在所有圓之外的概率。說白了就是“多邊形與圓的面積交”加一個簡單的概率換算。這道題對準備區(qū)域賽、或者想系統(tǒng)補計算幾何模板的人來說性價比很高。它不像那些動輒后綴自動機、網(wǎng)絡(luò)流的難題卻能把幾何里最常用的基礎(chǔ)工具串起來有向面積、極角差、扇形面積、線段與圓求交、浮點誤差處理。我自己當(dāng)年在這題上被精度問題和符號問題折磨了大半天所以想著把完整的推導(dǎo)和實現(xiàn)整理出來給后來的人省點時間。1. 題意還原從隨機撒點到面積占比的概率建模1.1 輸入輸出和題面到底在說什么先說清楚UVa 11355的原題是英文題面網(wǎng)上能找到的版本細節(jié)有點出入我按自己記憶和常見解題報告的描述還原一下核心設(shè)定——一個凸多邊形區(qū)域內(nèi)部有一些圓問隨機選點落在所有圓外的概率也就是“Cool Point”的比例。具體輸入格式請以你手里的原始題面為準這里更關(guān)鍵的是理解模型。典型的輸入結(jié)構(gòu)大概是這樣的先給出多邊形頂點數(shù)n和n個頂點坐標然后給出圓的個數(shù)m和每個圓的圓心坐標、半徑。多邊形的頂點按什么順序給不一定但一定構(gòu)成凸多邊形。輸出通常是保留若干位小數(shù)的概率。我第一次做的時候犯過一個低級錯誤盯著“概率”兩個字想用蒙特卡洛模擬去隨機撒點測頻率。幸虧看了眼數(shù)據(jù)范圍果斷放棄這種題要的是精確解不是統(tǒng)計近似。隨機模擬是驗證答案的手段不是解法本身。1.2 均勻隨機下的概率就是面積比這里需要把“概率”翻譯成“面積”。如果點是在多邊形內(nèi)均勻隨機選取的那么落在某個區(qū)域的概率就等于該區(qū)域面積占多邊形總面積的比例。這是幾何概型最樸素的定義。所以Cool Points的答案其實就是Cool概率 1 - (所有圓覆蓋區(qū)域面積 / 多邊形面積)關(guān)鍵點在于“所有圓覆蓋區(qū)域面積”這六個字。如果多個圓之間存在重疊直接對每個圓與多邊形的交面積求和重復(fù)的部分會被算兩次最后算出的概率會偏小。多數(shù)情況下這道題的數(shù)據(jù)會保證圓與圓不相交或者你只需要按題面給定的約束去處理。如果你的版本沒有明確說明那就要去求“圓的面積并”而不是簡單累加這一點我在第5節(jié)會專門展開講也是很多人WA得莫名其妙的地方。1.3 這道題真正在考什么剝掉概率這層外殼本質(zhì)是計算幾何里最基礎(chǔ)也最?？嫉膯栴}一個凸多邊形和一個圓相交交集的面積怎么算。聽起來簡單做起來全是細節(jié)。圓在內(nèi)部、圓在外部、圓心在多邊形內(nèi)、圓心在多邊形外、圓邊界穿過某條邊、圓剛好包裹一個頂點……每種情況都可能是合法的輸入。直接寫一堆if-else判斷點與圓的位置關(guān)系代碼很快會爛成一鍋粥。好在這類問題有一條經(jīng)典的捷徑把多邊形拆成若干個以圓心為頂點的有向三角形逐個計算每個三角形與圓的交面積再相加。這個思路能統(tǒng)一處理所有case不需要為“圓心是否在多邊形內(nèi)”單獨討論也是我接下來要重點拆解的部分。2. 核心幾何多邊形面積的有向三角形分解2.1 任意簡單多邊形的有向面積公式先回顧一個多邊形最基本的公式給定按順序排列的頂點(P_0, P_1, \dots, P_{n-1})它的有向面積是[ S \frac{1}{2} \sum_{i0}^{n-1} \left( P_i.x \times P_{i1}.x? \right) ]寫成代碼就是double cross(const Point a, const Point b) { return a.x * b.y - a.y * b.x; } double polygonArea(const vectorPoint poly) { double res 0; int n poly.size(); for (int i 0; i n; i) { res cross(poly[i], poly[(i 1) % n]); } return fabs(res) * 0.5; }這個公式對凸多邊形和凹多邊形都成立前提是頂點按順時針或逆時針順序給出。算出來是帶符號的逆時針為正順時針為負取絕對值就是普通面積。這個“有符號”特性不是毛病反而是后面算法的基石。2.2 以圓心為頂點拆三角形為什么這么拆現(xiàn)在假設(shè)圓心是C我們想把多邊形(P_0P_1\dots P_{n-1})與圓C的交集面積拆成一堆小塊的代數(shù)和。考慮圓心C和任意一條邊(P_iP_{i1})組成的三角形(C-P_i-P_{i1})。把所有這些有向三角形的面積加起來剛好等于整個多邊形的有向面積。這是個非常重要的性質(zhì)不管圓心在多邊形內(nèi)還是多邊形外這個等式都成立因為三角形帶符號后外部區(qū)域會被正負抵消。既然多邊形的面積能這么拆那么多邊形與圓的交面積也能這么拆對每一條邊計算三角形(C-P_i-P_{i1})與圓C的交集面積帶符號全部累加最后取絕對值。為什么這樣拆是有效的因為“和”滿足線性疊加而每個三角形與圓的交面積都可以在一個統(tǒng)一框架下計算。這個框架的關(guān)鍵是把圓心C平移到原點這樣圓就變成了以原點為圓心、半徑為r的標準圓所有計算都圍繞原點到兩個頂點的向量展開各種幾何量都變得非常干凈。2.3 符號問題必須先想清楚我自己第一次寫時在這里翻過車累加每個三角形與圓的交面積時忘了它是“有符號”的畫個圖以為自己推錯了。正確的認識是如果多邊形頂點是逆時針順序每條邊(P_iP_{i1})與圓心C形成的三角形其有向面積符號由叉積(P_iP_{i1} \times CP_i)決定圓心在多邊形內(nèi)部時每個三角形的符號通常相同圓心在多邊形外部時一部分三角形是正的一部分是負的加起來才等于多邊形面積。所以核心函數(shù)triCircleIntersect必須保留符號最終累加完再取fabs。中間任何一步取絕對值都會破壞這個“符號抵消”機制這也是網(wǎng)上很多人照抄模板卻一直WA的隱藏原因之一。3. 圓與三角形求交的完整推導(dǎo)五種情形補全3.1 基礎(chǔ)工具極角差與扇形面積圓與三角形求交繞不開“扇形面積”。設(shè)圓心在原點兩個向量a和b從原點出發(fā)那么以原點為圓心、半徑為r的扇形O-a-b的面積是[ S_{\text{sector}} \frac{1}{2} r^2 \theta ]其中(\theta)是向量a轉(zhuǎn)到向量b的有向夾角。求這個夾角最穩(wěn)的方式不是用acos(dot / (len*len))而是用atan2double sectorArea(const Point a, const Point b, double r) { double ang atan2(cross(a, b), dot(a, b)); return ang * r * r * 0.5; }atan2(cross, dot)返回的是向量a到向量b的有向夾角范圍在([-\pi, \pi])天然帶符號。用acos的問題有兩個一是值域只有([0, \pi])正負信息會丟二是當(dāng)夾角接近0或(\pi)時浮點誤差會被放大。實測下來atan2方案不僅代碼統(tǒng)一精度也穩(wěn)得多。3.2 情形A兩個端點都在圓內(nèi)這是最簡單的情況。三角形三個點圓心O、端點A、端點B中A和B都在圓內(nèi)或圓上而O是圓心當(dāng)然也在圓內(nèi)。由于三角形是凸組合整個三角形都在圓內(nèi)交集面積就是三角形本身。double triCircleIntersect(Point a, Point b, double r) { double s cross(a, b); if (fabs(s) EPS) return 0.0; // 退化三角形 double da len(a), db len(b); double res 0.0; if (da r EPS db r EPS) { res s * 0.5; } // ... 其余情形 return res; }注意保留叉積的符號如果a在b的順時針方向s會是負的這樣累加才能正確反映“三角形在圓心外側(cè)”的情況。3.3 情形B一個端點在圓內(nèi)一個在圓外假設(shè)a在圓內(nèi)b在圓外。線段ab必然與圓相交且只有一個交點p。此時三角形與圓的交集由一個鈍角三角形a-O-p和一個扇形O-p-b拼成Point lineCircleIntersect(Point a, Point b, double r, bool ok) { // 僅用于 a 在內(nèi)、b 在外 或 b 在內(nèi)、a 在外 的情況 ok false; Point d b - a; double A dot(d, d); double B 2 * dot(a, d); double C dot(a, a) - r * r; double delta B * B - 4 * A * C; if (delta -EPS) return Point(); delta max(0.0, delta); double t1 (-B - sqrt(delta)) / (2 * A); double t2 (-B sqrt(delta)) / (2 * A); double t t1; if (t -EPS || t 1 EPS) t t2; if (t -EPS || t 1 EPS) { ok false; return Point(); } ok true; return a d * t; }在triCircleIntersect里對應(yīng)分支是else if (da r EPS db r EPS) { bool ok; Point p lineCircleIntersect(a, b, r, ok); res cross(a, p) * 0.5 sectorArea(p, b, r); } else if (da r EPS db r EPS) { bool ok; Point p lineCircleIntersect(b, a, r, ok); res sectorArea(a, p, r) cross(p, b) * 0.5; }要注意lineCircleIntersect(a, b, r)求的是線段ab上靠近a的交點所以當(dāng)b在內(nèi)、a在外時要先傳(b, a)拿到離b近的交點再讓扇形從a切到p。3.4 情形C兩個端點都在圓外且線段ab與圓相交這是最有推導(dǎo)價值的一種情形。a和b都在圓外線段ab穿過圓有兩個交點p、q。此時三角形O-a-b與圓的交集其實是“完整扇形O-a-b”減去“弓形區(qū)域弦pq對應(yīng)的圓外部分”。完整扇形面積是sectorArea(a, b, r)弓形面積是扇形O-p-q減去三角形O-p-q的面積。整理一下交集面積剛好等于vectorPoint segmentCircleIntersect(Point a, Point b, double r) { vectorPoint res; Point d b - a; double A dot(d, d); double B 2 * dot(a, d); double C dot(a, a) - r * r; double delta B * B - 4 * A * C; if (delta -EPS) return res; double root sqrt(max(0.0, delta)); double t1 (-B - root) / (2 * A); if (t1 -EPS t1 1 EPS) res.push_back(a d * t1); double t2 (-B root) / (2 * A); if (t2 -EPS t2 1 EPS) { Point p a d * t2; if (res.empty() || len(p - res.back()) EPS) res.push_back(p); } return res; }在triCircleIntersect中else { vectorPoint pts segmentCircleIntersect(a, b, r); if (pts.size() 2) { Point p pts[0], q pts[1]; if (dot(p - a, q - a) 0) swap(p, q); res sectorArea(a, p, r) cross(p, q) * 0.5 sectorArea(q, b, r); } else { res sectorArea(a, b, r); } }這里cross(p, q) * 0.5就是三角形O-p-q的有向面積。整個式子的幾何含義從a到p是扇形p到q是三角形q到b是扇形三段拼接正好構(gòu)成了圓內(nèi)部、三角形覆蓋的那塊區(qū)域。這是最容易被畫圖誤解的地方建議自己動手畫幾個不同位置的圓和邊驗證一下公式。3.5 情形D兩個端點都在圓外且線段ab與圓不相交如果線段ab與圓沒有交點要么整條邊離圓很遠要么垂足落在a或b的外側(cè)。此時三角形O-a-b與圓的交集就是夾角(\angle aOb)對應(yīng)的整個扇形前提是這個角小于(\pi)。而三角形的內(nèi)角天然小于(\pi)所以直接返回res sectorArea(a, b, r);有人會擔(dān)心如果圓心O到線段ab的垂足落在線段上但距離剛好大于r此時ab雖與圓無交點但三角形里包含的是不是整個圓不是。因為三角形以O(shè)為頂點只有夾角范圍內(nèi)的扇形屬于三角形而該扇形本身不可能把整個圓包進去除非夾角是(2\pi)那不可能。再補充一種邊界線段ab與圓相切pts.size()為1同樣走這個分支。相切時圓弧沒有缺口完整的扇形正好就是交集所以結(jié)果也是對的。3.6 五情形匯總條件交集面積備注兩端點都在圓內(nèi)三角形面積直接用叉積除以2a內(nèi)b外三角形a-O-p 扇形O-p-bp為線段ab與圓的交點a外b內(nèi)扇形O-a-p 三角形p-O-bp為交點兩端點都在圓外ab與圓交于兩點扇形a-O-p 三角形p-O-q 扇形q-O-b兩交點p、q按方向排序兩端點都在圓外ab與圓不相交或相切完整扇形a-O-b相切也歸入此類這張表配合代碼看基本就能覆蓋所有正常情況。真正在寫題時還需要警惕的是浮點誤差和退化情況這些我放在第5節(jié)詳細講。4. 可以直接抄的C實現(xiàn)與復(fù)雜度說明4.1 結(jié)構(gòu)體與基礎(chǔ)函數(shù)把上面所有片段拼起來就是一個比較完整的計算幾何模板。先定義Point和基礎(chǔ)運算#include bits/stdc.h using namespace std; const double EPS 1e-10; const double PI acos(-1.0); struct Point { double x, y; Point(double x_ 0, double y_ 0) : x(x_), y(y_) {} Point operator - (const Point rhs) const { return Point(x - rhs.x, y - rhs.y); } Point operator (const Point rhs) const { return Point(x rhs.x, y rhs.y); } Point operator * (double k) const { return Point(x * k, y * k); } }; double cross(const Point a, const Point b) { return a.x * b.y - a.y * b.x; } double dot(const Point a, const Point b) { return a.x * b.x a.y * b.y; } double len(const Point a) { return sqrt(dot(a, a)); }4.2 核心函數(shù)三角形與圓求交下面這個函數(shù)是整個算法的心臟輸入a、b是相對于圓心的向量r是半徑返回帶符號的交面積double sectorArea(const Point a, const Point b, double r) { double ang atan2(cross(a, b), dot(a, b)); return ang * r * r * 0.5; } double triCircleIntersect(Point a, Point b, double r) { double s cross(a, b); if (fabs(s) EPS) return 0.0; double da len(a), db len(b); double res 0.0; if (da r EPS db r EPS) { res s * 0.5; } else if (da r EPS db r EPS) { bool ok; Point p lineCircleIntersect(a, b, r, ok); res cross(a, p) * 0.5 sectorArea(p, b, r); } else if (da r EPS db r EPS) { bool ok; Point p lineCircleIntersect(b, a, r, ok); res sectorArea(a, p, r) cross(p, b) * 0.5; } else { vectorPoint pts segmentCircleIntersect(a, b, r); if (pts.size() 2) { Point p pts[0], q pts[1]; if (dot(p - a, q - a) 0) swap(p, q); res sectorArea(a, p, r) cross(p, q) * 0.5 sectorArea(q, b, r); } else { res sectorArea(a, b, r); } } return res; }4.3 主流程與答案計算有了核心函數(shù)主流程非常簡單struct Circle { Point c; double r; }; double polygonCircleIntersect(const vectorPoint poly, const Circle cir) { double res 0; int n poly.size(); for (int i 0; i n; i) { Point a poly[i] - cir.c; Point b poly[(i 1) % n] - cir.c; res triCircleIntersect(a, b, cir.r); } return fabs(res); } int main() { int n; while (scanf(%d, n) 1 n) { vectorPoint poly(n); for (int i 0; i n; i) { scanf(%lf%lf, poly[i].x, poly[i].y); } double polyArea polygonArea(poly); int m; scanf(%d, m); double covered 0.0; for (int i 0; i m; i) { Circle cir; scanf(%lf%lf%lf, cir.c.x, cir.c.y, cir.r); covered polygonCircleIntersect(poly, cir); } double ans 1.0 - covered / polyArea; printf(%.5lf\n, ans); } return 0; }復(fù)雜度是(O(n \times m))每個圓都要掃一遍多邊形的每條邊每條邊只進行常數(shù)次求交和三角函數(shù)運算。n和m通常都在幾十到幾百的量級跑起來毫無壓力。即使n和m都到1000也只是百萬級運算不存在性能瓶頸。4.4 輸出精度和EPS選擇輸出保留幾位小數(shù)以題目要求為準。常見的幾何題精度要求是1e-5或1e-6所以我習(xí)慣用%.5lf起步。EPS設(shè)成1e-10主要用來處理“點在圓上”“相切”這類臨界情況。注意EPS不是越大越好設(shè)太大容易把真正的交點過濾掉設(shè)太小又沒法覆蓋浮點誤差累積1e-10在多數(shù)double運算場景下是安全的。5. 我在調(diào)試中踩過的四個坑5.1 有符號面積忘記取絕對值這個坑最隱蔽。多邊形頂點逆時針給時三角形與圓的交面積累加結(jié)果為正但如果不小心把輸入當(dāng)成順時針累加結(jié)果就是負的最后概率會變成一個大于1或者負數(shù)的奇怪值。解決辦法是polygonCircleIntersect最后統(tǒng)一fabs(res)。同時polygonArea也要fabs否則面積可能是負的概率就亂套了。我見過有人的寫法是每個三角形單獨取fabs再累加這在圓心位于多邊形外時會算出完全錯誤的結(jié)果因為外部區(qū)域的三角形面積符號需要相互抵消。5.2 夾角用acos導(dǎo)致正負丟失最開始我用acos(dot(a, b) / (len(a) * len(b)))算扇形角結(jié)果在夾角接近(\pi)時角度會跳變累加結(jié)果對不上。后來全部換成atan2(cross, dot)一次修改解決所有問題。核心原因是acos只能返回([0, \pi])的非負角而我們需要的是帶方向的角。比如從向量a順時針轉(zhuǎn)到b角度應(yīng)當(dāng)是負的atan2能正確返回acos做不到。5.3 圓與圓重疊問題如果題面沒有明確“圓互不相交”直接累加每個圓與多邊形的交面積重疊部分會被重復(fù)計算。我印象里UVa 11355的常見數(shù)據(jù)假設(shè)是互不相交但我吃過一次虧后養(yǎng)成了習(xí)慣先看約束條件不確定就做并集處理。圓的并集面積算法比這題本身復(fù)雜得多需要處理圓與圓相交的弧段。這里給出一個簡單思路如果把所有圓限制在同一個多邊形內(nèi)且圓的個數(shù)不多可以先對所有圓兩兩求交點把每個圓被其他圓覆蓋的弧段剔除再用三角剖分求并集面積。但這不是UVa 11355的常規(guī)解法更像一個擴展話題。真正考試或訓(xùn)練時優(yōu)先確認題面假設(shè)別在沒必要的復(fù)雜度上浪費時間。5.4 相切時的交點個數(shù)問題線段與圓相切時segmentCircleIntersect理論上有兩個相等的交點但由于浮點誤差可能返回1個或2個。我在代碼里用len(p - res.back()) EPS去重保證pts面積里不會混入同一個點兩次。如果去重沒寫好pts.size() 2分支里的cross(p, q) * 0.5會變成一個近似零的量看起來無害但累積起來可能讓答案差幾個單位的精度。這種問題不會讓你WA得明顯只會在邊界數(shù)據(jù)上卡你一手調(diào)試時非常難受。6. 一點個人體會寫計算幾何題最大的敵人不是數(shù)學(xué)而是浮點、符號、邊界這三種“臟活”。UVa 11355恰好把這三樣都湊齊了值得反復(fù)做幾遍。我后來把它整理成手寫模板的起點凡是涉及多邊形與圓面積交、多邊形與圓并集、圓與圓相交面積的題目基本都是從這份代碼擴展出去的。如果你正在準備比賽建議不要直接復(fù)制這份模板了事而是自己照著推一遍、敲一遍。尤其是triCircleIntersect里那五個分支只有親手畫過圖、親手被樣例卡過才能真正理解為什么每個分叉都要這樣寫。之后你再看類似題目就不會再覺得幾何是“玄學(xué)”了。