第一次在Virtual Judge上刷到UVa 11355 Cool Points这道题时我一度以为是个脑筋急转弯点还能分“酷”和“不酷”等看清题面才发现这是一道非常典型的计算几何应用题——给定一个凸多边形作为随机撒点的范围再给若干个圆问随机选一个点落在所有圆之外的概率。说白了就是“多边形与圆的面积交”加一个简单的概率换算。这道题对准备区域赛、或者想系统补计算几何模板的人来说性价比很高。它不像那些动辄后缀自动机、网络流的难题却能把几何里最常用的基础工具串起来有向面积、极角差、扇形面积、线段与圆求交、浮点误差处理。我自己当年在这题上被精度问题和符号问题折磨了大半天所以想着把完整的推导和实现整理出来给后来的人省点时间。1. 题意还原从随机撒点到面积占比的概率建模1.1 输入输出和题面到底在说什么先说清楚UVa 11355的原题是英文题面网上能找到的版本细节有点出入我按自己记忆和常见解题报告的描述还原一下核心设定——一个凸多边形区域内部有一些圆问随机选点落在所有圆外的概率也就是“Cool Point”的比例。具体输入格式请以你手里的原始题面为准这里更关键的是理解模型。典型的输入结构大概是这样的先给出多边形顶点数n和n个顶点坐标然后给出圆的个数m和每个圆的圆心坐标、半径。多边形的顶点按什么顺序给不一定但一定构成凸多边形。输出通常是保留若干位小数的概率。我第一次做的时候犯过一个低级错误盯着“概率”两个字想用蒙特卡洛模拟去随机撒点测频率。幸亏看了眼数据范围果断放弃这种题要的是精确解不是统计近似。随机模拟是验证答案的手段不是解法本身。1.2 均匀随机下的概率就是面积比这里需要把“概率”翻译成“面积”。如果点是在多边形内均匀随机选取的那么落在某个区域的概率就等于该区域面积占多边形总面积的比例。这是几何概型最朴素的定义。所以Cool Points的答案其实就是Cool概率 1 - (所有圆覆盖区域面积 / 多边形面积)关键点在于“所有圆覆盖区域面积”这六个字。如果多个圆之间存在重叠直接对每个圆与多边形的交面积求和重复的部分会被算两次最后算出的概率会偏小。多数情况下这道题的数据会保证圆与圆不相交或者你只需要按题面给定的约束去处理。如果你的版本没有明确说明那就要去求“圆的面积并”而不是简单累加这一点我在第5节会专门展开讲也是很多人WA得莫名其妙的地方。1.3 这道题真正在考什么剥掉概率这层外壳本质是计算几何里最基础也最常考的问题一个凸多边形和一个圆相交交集的面积怎么算。听起来简单做起来全是细节。圆在内部、圆在外部、圆心在多边形内、圆心在多边形外、圆边界穿过某条边、圆刚好包裹一个顶点……每种情况都可能是合法的输入。直接写一堆if-else判断点与圆的位置关系代码很快会烂成一锅粥。好在这类问题有一条经典的捷径把多边形拆成若干个以圆心为顶点的有向三角形逐个计算每个三角形与圆的交面积再相加。这个思路能统一处理所有case不需要为“圆心是否在多边形内”单独讨论也是我接下来要重点拆解的部分。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 以圆心为顶点拆三角形为什么这么拆现在假设圆心是C我们想把多边形(P_0P_1\dots P_{n-1})与圆C的交集面积拆成一堆小块的代数和。考虑圆心C和任意一条边(P_iP_{i1})组成的三角形(C-P_i-P_{i1})。把所有这些有向三角形的面积加起来刚好等于整个多边形的有向面积。这是个非常重要的性质不管圆心在多边形内还是多边形外这个等式都成立因为三角形带符号后外部区域会被正负抵消。既然多边形的面积能这么拆那么多边形与圆的交面积也能这么拆对每一条边计算三角形(C-P_i-P_{i1})与圆C的交集面积带符号全部累加最后取绝对值。为什么这样拆是有效的因为“和”满足线性叠加而每个三角形与圆的交面积都可以在一个统一框架下计算。这个框架的关键是把圆心C平移到原点这样圆就变成了以原点为圆心、半径为r的标准圆所有计算都围绕原点到两个顶点的向量展开各种几何量都变得非常干净。2.3 符号问题必须先想清楚我自己第一次写时在这里翻过车累加每个三角形与圆的交面积时忘了它是“有符号”的画个图以为自己推错了。正确的认识是如果多边形顶点是逆时针顺序每条边(P_iP_{i1})与圆心C形成的三角形其有向面积符号由叉积(P_iP_{i1} \times CP_i)决定圆心在多边形内部时每个三角形的符号通常相同圆心在多边形外部时一部分三角形是正的一部分是负的加起来才等于多边形面积。所以核心函数triCircleIntersect必须保留符号最终累加完再取fabs。中间任何一步取绝对值都会破坏这个“符号抵消”机制这也是网上很多人照抄模板却一直WA的隐藏原因之一。3. 圆与三角形求交的完整推导五种情形补全3.1 基础工具极角差与扇形面积圆与三角形求交绕不开“扇形面积”。设圆心在原点两个向量a和b从原点出发那么以原点为圆心、半径为r的扇形O-a-b的面积是[ S_{\text{sector}} \frac{1}{2} r^2 \theta ]其中(\theta)是向量a转到向量b的有向夹角。求这个夹角最稳的方式不是用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])正负信息会丢二是当夹角接近0或(\pi)时浮点误差会被放大。实测下来atan2方案不仅代码统一精度也稳得多。3.2 情形A两个端点都在圆内这是最简单的情况。三角形三个点圆心O、端点A、端点B中A和B都在圆内或圆上而O是圆心当然也在圆内。由于三角形是凸组合整个三角形都在圆内交集面积就是三角形本身。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会是负的这样累加才能正确反映“三角形在圆心外侧”的情况。3.3 情形B一个端点在圆内一个在圆外假设a在圆内b在圆外。线段ab必然与圆相交且只有一个交点p。此时三角形与圆的交集由一个钝角三角形a-O-p和一个扇形O-p-b拼成Point lineCircleIntersect(Point a, Point b, double r, bool ok) { // 仅用于 a 在内、b 在外 或 b 在内、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里对应分支是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的交点所以当b在内、a在外时要先传(b, a)拿到离b近的交点再让扇形从a切到p。3.4 情形C两个端点都在圆外且线段ab与圆相交这是最有推导价值的一种情形。a和b都在圆外线段ab穿过圆有两个交点p、q。此时三角形O-a-b与圆的交集其实是“完整扇形O-a-b”减去“弓形区域弦pq对应的圆外部分”。完整扇形面积是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是扇形三段拼接正好构成了圆内部、三角形覆盖的那块区域。这是最容易被画图误解的地方建议自己动手画几个不同位置的圆和边验证一下公式。3.5 情形D两个端点都在圆外且线段ab与圆不相交如果线段ab与圆没有交点要么整条边离圆很远要么垂足落在a或b的外侧。此时三角形O-a-b与圆的交集就是夹角(\angle aOb)对应的整个扇形前提是这个角小于(\pi)。而三角形的内角天然小于(\pi)所以直接返回res sectorArea(a, b, r);有人会担心如果圆心O到线段ab的垂足落在线段上但距离刚好大于r此时ab虽与圆无交点但三角形里包含的是不是整个圆不是。因为三角形以O为顶点只有夹角范围内的扇形属于三角形而该扇形本身不可能把整个圆包进去除非夹角是(2\pi)那不可能。再补充一种边界线段ab与圆相切pts.size()为1同样走这个分支。相切时圆弧没有缺口完整的扇形正好就是交集所以结果也是对的。3.6 五情形汇总条件交集面积备注两端点都在圆内三角形面积直接用叉积除以2a内b外三角形a-O-p 扇形O-p-bp为线段ab与圆的交点a外b内扇形O-a-p 三角形p-O-bp为交点两端点都在圆外ab与圆交于两点扇形a-O-p 三角形p-O-q 扇形q-O-b两交点p、q按方向排序两端点都在圆外ab与圆不相交或相切完整扇形a-O-b相切也归入此类这张表配合代码看基本就能覆盖所有正常情况。真正在写题时还需要警惕的是浮点误差和退化情况这些我放在第5节详细讲。4. 可以直接抄的C实现与复杂度说明4.1 结构体与基础函数把上面所有片段拼起来就是一个比较完整的计算几何模板。先定义Point和基础运算#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 核心函数三角形与圆求交下面这个函数是整个算法的心脏输入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 主流程与答案计算有了核心函数主流程非常简单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; }复杂度是(O(n \times m))每个圆都要扫一遍多边形的每条边每条边只进行常数次求交和三角函数运算。n和m通常都在几十到几百的量级跑起来毫无压力。即使n和m都到1000也只是百万级运算不存在性能瓶颈。4.4 输出精度和EPS选择输出保留几位小数以题目要求为准。常见的几何题精度要求是1e-5或1e-6所以我习惯用%.5lf起步。EPS设成1e-10主要用来处理“点在圆上”“相切”这类临界情况。注意EPS不是越大越好设太大容易把真正的交点过滤掉设太小又没法覆盖浮点误差累积1e-10在多数double运算场景下是安全的。5. 我在调试中踩过的四个坑5.1 有符号面积忘记取绝对值这个坑最隐蔽。多边形顶点逆时针给时三角形与圆的交面积累加结果为正但如果不小心把输入当成顺时针累加结果就是负的最后概率会变成一个大于1或者负数的奇怪值。解决办法是polygonCircleIntersect最后统一fabs(res)。同时polygonArea也要fabs否则面积可能是负的概率就乱套了。我见过有人的写法是每个三角形单独取fabs再累加这在圆心位于多边形外时会算出完全错误的结果因为外部区域的三角形面积符号需要相互抵消。5.2 夹角用acos导致正负丢失最开始我用acos(dot(a, b) / (len(a) * len(b)))算扇形角结果在夹角接近(\pi)时角度会跳变累加结果对不上。后来全部换成atan2(cross, dot)一次修改解决所有问题。核心原因是acos只能返回([0, \pi])的非负角而我们需要的是带方向的角。比如从向量a顺时针转到b角度应当是负的atan2能正确返回acos做不到。5.3 圆与圆重叠问题如果题面没有明确“圆互不相交”直接累加每个圆与多边形的交面积重叠部分会被重复计算。我印象里UVa 11355的常见数据假设是互不相交但我吃过一次亏后养成了习惯先看约束条件不确定就做并集处理。圆的并集面积算法比这题本身复杂得多需要处理圆与圆相交的弧段。这里给出一个简单思路如果把所有圆限制在同一个多边形内且圆的个数不多可以先对所有圆两两求交点把每个圆被其他圆覆盖的弧段剔除再用三角剖分求并集面积。但这不是UVa 11355的常规解法更像一个扩展话题。真正考试或训练时优先确认题面假设别在没必要的复杂度上浪费时间。5.4 相切时的交点个数问题线段与圆相切时segmentCircleIntersect理论上有两个相等的交点但由于浮点误差可能返回1个或2个。我在代码里用len(p - res.back()) EPS去重保证pts面积里不会混入同一个点两次。如果去重没写好pts.size() 2分支里的cross(p, q) * 0.5会变成一个近似零的量看起来无害但累积起来可能让答案差几个单位的精度。这种问题不会让你WA得明显只会在边界数据上卡你一手调试时非常难受。6. 一点个人体会写计算几何题最大的敌人不是数学而是浮点、符号、边界这三种“脏活”。UVa 11355恰好把这三样都凑齐了值得反复做几遍。我后来把它整理成手写模板的起点凡是涉及多边形与圆面积交、多边形与圆并集、圆与圆相交面积的题目基本都是从这份代码扩展出去的。如果你正在准备比赛建议不要直接复制这份模板了事而是自己照着推一遍、敲一遍。尤其是triCircleIntersect里那五个分支只有亲手画过图、亲手被样例卡过才能真正理解为什么每个分叉都要这样写。之后你再看类似题目就不会再觉得几何是“玄学”了。