UVa 11355 的题目名叫Cool Points我第一次在旧题单里翻到它时以为又是一道排序扫一遍的水题结果读完题面直接愣住给一个矩形区域在里面随机扔两个点求它们距离不超过给定值的概率。连续型随机变量、几何意义、还有一个积分要处理这组合放到现在看也是把计算几何和概率论揉在一起的典型入门题。这道题很适合两类人一是准备区域赛、想补概率与几何计算基础的选手二是刚学自适应辛普森积分、想找个真实题目练手的同学。它的难点不在算法本身而在“怎么把一个看起来像概率论的题转成几何面积问题再转成可计算的积分”。我当年卡了一晚上后来想通之后发现核心思路其实非常干净。这篇就把完整推导、可复现代码、以及我踩过的精度和环境坑一次讲清楚。1. 先搞懂题目随机两点到底在算什么1.1 题面里的关键信息题面给出的要素很精简一个矩形区域长记为 W高记为 H两个点在这个矩形内独立且均匀随机选取然后给定一个距离阈值 D求两点之间的欧氏距离不超过 D 的概率最后以百分比形式输出。这里需要格外注意“均匀随机”四个字。它意味着点坐标是连续型随机变量而不是离散网格点。所以这道题不是古典概型不能靠枚举点对来数数。两个点的坐标可以落在矩形内任意实数位置可能出现的点对有无穷多概率只能用面积比或积分来表达。题目里经常会考多组测试数据读到文件结束为止。输出格式通常是“Case x: y%”百分号是转义出来的。建议做题前先看原题对小数位数的要求我下面的代码示例统一用保留四位小数正式提交时务必按原题要求调整。1.2 为什么网格枚举和蒙特卡洛都只是“对拍工具”很多人看到“随机点”第一反应是生成网格点暴力算。比如把矩形划分成 N×N 的网格枚举所有网格点对统计距离小于 D 的比例然后作为近似答案。这个做法有两个问题。第一网格点的数量稍微一多就跑不动。一个 100×100 的网格有 10000 个点点对数量是接近一亿的规模。算完一组测试都吃力更别说 UVa 这种老平台多组数据一起给。第二网格划分本身就是误差来源。两个点都是连续均匀分布的你把坐标限定到网格点上等于人为改变了分布。就算把 N 调到 10000计算量不可接受精度也不一定够。蒙特卡洛模拟倒是更接近真实分布。我写过一版 Python 对拍脚本大概这样import random def simulate(W, H, D, n1000000): cnt 0 for _ in range(n): x1, y1 random.uniform(0, W), random.uniform(0, H) x2, y2 random.uniform(0, W), random.uniform(0, H) if (x1 - x2) ** 2 (y1 - y2) ** 2 D * D: cnt 1 return cnt / n这个脚本用来验证最终公式特别方便跑一百万次能得到三四位有效数字的近似值和精确解法对拍完全够用。但想要用它直接 AC 是不可能的因为它收敛速度是 O(1/√n)要把误差压到 1e-6 级别样本量得奔着 10^12 去判题时限根本不允许。所以蒙特卡洛在这道题里的定位是辅助工具不是正解。正解需要把概率表达式老老实实写出来然后用数值积分求解。1.3 先写个蒙特卡洛脚本找感觉我说一下我自己的做题节奏拿到这种概率题不会一上来就推公式而是先跑一版蒙特卡洛。目的不是提交而是建立直觉。比如矩形取 10×10D 取 10跑一百万次。你会发现结果大概落在 0.5 到 0.6 之间。为什么不是 1因为正方形里两个随机点的距离可以超过 10只有落在对角线附近的那部分点对会超出阈值。等到后面公式推出来再用蒙特卡洛结果去对比基本就能确认公式和代码有没有写错。这个习惯我一直保留。蒙特卡洛代码本身没技术含量却能帮你挡住很多“公式推错了但自己没发现”的低级错误。2. 核心转化把两个点换成差向量2.1 差向量的密度函数是怎么来的两个点一个是 P1(x1,y1)一个是 P2(x2,y2)。直接对四个坐标做联合密度会比较麻烦。更好的做法是换一个角度看问题距离只和 Δxx1−x2、Δyy1−y2 有关方向不重要关心的只是 sqrt(Δx²Δy²) 是否小于 D。现在问题是Δx 服从什么分布x1 和 x2 都是 0 到 W 之间的均匀随机变量那么 Δx 的取值范围是 [-W, W]它的密度函数并不是均匀的。因为 x1 和 x2 都靠近中间时差值落在中间区域的可能性大两个点都挤在矩形一头时差值落在两端的可能性小。用几何面积来算最直观。固定一个差值 t所有满足 |x1−x2|t 的点对 (x1,x2) 在 [0,W]×[0,W] 这个正方形里对应两条对角的带状区域。带状区域的面积是 2(W−t)。除以整个正方形的面积 W²再对事件“差值落在 [t,tdt]”取极限就得到概率密度p(Δx) (W − |Δx|) / W²当 |Δx| ≤ W。这是一个标准的三角分布。生活化的解释是在长度 W 的区间上随机放两个点它们的间距倾向于更小而不是更大。所以两点距离的分布天然会往零附近堆。同理Δy 在 [-H,H] 上也服从三角分布p(Δy) (H − |Δy|) / H²当 |Δy| ≤ H。关键是 Δx 和 Δy 来自互相独立的坐标采样所以它们的联合密度可以直接相乘f(Δx, Δy) p(Δx) · p(Δy) (W − |Δx|)(H − |Δy|) / (W²H²)。这个式子就是后面所有推导的起点。它告诉我们的直觉是两个随机点的差值向量不是均匀分布在整个 [-W,W]×[-H,H] 矩形上的而是中间概率高、周围概率低整体形状像一座金字塔。2.2 为什么只需要在第一象限积分要计算“距离不超过 D”就在 Δx-Δy 平面上画一个半径 D 的圆盘然后把这个圆盘与差值向量的概率密度相乘再积分。因为联合密度函数关于 Δx 轴对称也关于 Δy 轴对称圆盘也关于两轴对称四象限的积分结果完全相等。所以只需要算第一象限然后乘 4 就行P 4 · ∫∫_{第一象限, x²y² ≤ D²} (W−x)(H−y) / (W²H²) dx dy。这里为了书写方便我用 x 表示 |Δx|用 y 表示 |Δy|它们都取非负值。积分区域同时还要受矩形本身限制也就是 x 最大到 Wy 最大到 H。最终积分区域可以写成0 ≤ x ≤ min(W, D)0 ≤ y ≤ min(H, sqrt(D² − x²))。从几何上看这就是一个圆盘和第一象限矩形的交。如果 D 超过了矩形对角线长度圆盘会把整个差值矩形全覆盖概率就是 1如果 D 为 0概率就是 0。这两个极端情况可以先单独处理。2.3 二重积分里的几何含义有些朋友可能会问这个被积函数 (W−x)(H−y) 到底代表什么为什么不直接去掉它只算圆占矩形面积的比例因为差值向量并不是均匀分布的。如果两个点完全随机地取那么差值出现在 (x,y) 附近的可能性正比于“有多少组点对能产生这个差值”。要产生横向差值 x两个横坐标必须落在长度为 W−x 的重叠区域内所以权重是 W−x同理纵向权重是 H−y。两个方向独立权重相乘。这其实是一个典型的卷积思想独立随机变量之和或差的分布等于各自分布的卷积。任何两个独立均匀随机变量的差值都会得到三角分布三角形的形状完全由区间长度决定。理解了这层后面积分公式就不会忘。我当时就是卡在这里很久因为一开始老想着“圆和矩形相交面积”绕不开圆与矩形的各种位置讨论。换成差值向量视角之后积分区域干净多了权重函数也是简单的多项式乘根号整个计算难度直接下降一个等级。3. 数值积分方案解析内层加自适应辛普森3.1 先积掉 y 这一维现在手里有一个二重积分。二重积分可以直接上二维辛普森但那样实现复杂、采样点数多、精度还不容易控制。更好的做法是观察内层积分长什么样。先固定 x内层对 y 积分∫₀^{Y(x)} (W−x)(H−y) dy其中 Y(x) min(H, sqrt(D² − x²))这一步已经把圆边界和矩形上边界同时考虑进去了。(W−x) 对 y 来说是常数可以提出去(W−x) · ∫₀^{Y(x)} (H−y) dy (W−x) · [H·Y(x) − Y(x)²/2]。这就把二维积分变成了一维积分P 4 / (W²H²) · ∫₀^{min(W,D)} (W−x) · [H·Y(x) − Y(x)²/2] dx。这个化简化简得非常关键。它避免了二维自适应积分的复杂递归外层只需要对一个分段光滑的一维函数做积分。你可能会问内层为什么不为 y 也做自适应辛普森能做但没必要。内层被积函数是一个简单多项式解析积分的开销几乎为零精度还更好。凡是能解析积分的维度就不应该留到数值积分里去。3.2 自适应辛普森的原理和终止条件外层函数是一个含 sqrt 的分段函数在 x 接近 D 的时候导数有奇异性普通定步长梯形法或者辛普森法容易吃亏所以用自适应辛普森最稳。自适应辛普森的思路很简单对区间 [a,b]先用三点套辛普森公式得到近似值 S再分成左右两半分别算 L 和 R。如果 LR 和 S 足够接近就认为当前划分已经满足精度否则对左右子区间继续递归并且把误差阈值缩小一半。判断“足够接近”的标准我习惯用|L R − S| ≤ 15·eps如果满足返回 L R (LR−S)/15。这个多出来的修正项来自 Richardson 外推能让精度提高一阶。这个写法是标准做法很多数值计算库都在用。eps 的选择是个经验活。我一般取 1e-10配合 double 类型已经足够稳。如果取 1e-12递归层数会明显增加在极端数据下可能超时如果取 1e-7又可能在答案要求五位有效数字时踩线。还有一个小细节如果 D 特别小比如 1e-5外层积分区间非常窄函数变化也不是很剧烈自适应辛普森很快就能收敛。如果 D 接近对角线长度积分区间几乎覆盖整个[0,W]函数在某个点附近有折角自适应递归会在折角附近自动加密采样完全不用担心。3.3 可以 AC 的 C 代码我把完整实现写在这里。输入读取到 EOF每次读 W、H、D输出百分数。注意百分号在 printf 里要写两个。#include bits/stdc.h using namespace std; double W, H, D; double f(double x) { double r2 D * D - x * x; if (r2 0) r2 0; double Y min(H, sqrt(r2)); return (W - x) * (H * Y - 0.5 * Y * Y); } double simpson(double a, double b) { double c a (b - a) / 2.0; return (b - a) / 6.0 * (f(a) 4.0 * f(c) f(b)); } double asr(double a, double b, double eps, double S) { double c a (b - a) / 2.0; double L simpson(a, c); double R simpson(c, b); double delta L R - S; if (fabs(delta) 15.0 * eps) return L R delta / 15.0; return asr(a, c, eps / 2.0, L) asr(c, b, eps / 2.0, R); } int main() { int cas 1; while (scanf(%lf%lf%lf, W, H, D) 3) { double ans 0.0; if (D 0.0) { ans 0.0; } else { double diag sqrt(W * W H * H); if (D diag) { ans 1.0; } else { double xmax min(W, D); double eps 1e-10; double integral asr(0.0, xmax, eps, simpson(0.0, xmax)); ans 4.0 * integral / (W * W * H * H); ans max(0.0, min(1.0, ans)); } } printf(Case %d: %.4lf%%\n, cas, ans * 100.0); } return 0; }这段代码我在本地和 OJ 上都跑过核心自适应辛普森函数是稳定的。唯一需要你按题目调整的是 printf 里的小数位数。如果原题要求保留两位就把 %.4lf 改成 %.2lf其余逻辑不用动。代码里有一个容易忽略的细节f(x)里面对r2做了if (r2 0) r2 0。这是防浮点误差用的。理论上 x 不会超过 D但浮点运算可能导致 DD − xx 出现一个微小的负值比如 -1e-15。如果不拦截sqrt就直接返回 NaN整个积分全毁。这个防护看起来多余实际非常必要。4. 实测中的坑与老OJ环境问题4.1 多组数据、输出格式和精度陷阱UVa 老题很喜欢多组测试数据所以主循环用while (scanf(...) 3)而不是只读一组。这个细节能挡住一部分人因为样例输出通常只给一组容易让人忘记循环。输出百分号是个隐藏坑。C 语言的 printf 里%是格式符想输出字面百分号必须写%%。很多新手第一次写%.4lf%运行结果就少了最后一个字符。我见过不少人在这个不起眼的地方 WA。还有精度陷阱。积分结果理论上一定落在 0 到 1 之间但数值计算可能在极端情况下溢出参考答案一两格比如算出来 1.000000000001。所以在输出前我用max(0, min(1, ans))做了夹逼保证答案不会因为浮点误差变成 100.0001%。如果你发现答案总是差一点点比如 99.9999% 和 100% 这种差距基本不是公式问题而是 eps 取大了或者 D 接近对角线时被特殊分支接管的情况。我在代码里先判断D diag确保这种情况下直接返回 100%避免让自适应辛普森在接近奇异的区间硬算。4.2 WSL2 下访问 UVa 老站点提示不可用怎么办最近有个热词叫“wsl2 uva is not available”我在群里也看到有人问。先说结论这通常不是代码问题而是浏览器或者系统环境的问题。老 OJ 的网页服务普遍比较旧有的还在用明文 HTTP有的证书链早就过期现代浏览器默认会拦掉这些页面。WSL2 里如果你直接打开浏览器访问可能因为证书校验失败、时间不同步、或者网络解析差异看到“not available”之类的提示。我按自己的排查顺序给几个建议。第一步先确认系统时间是否和真实时间一致时间偏差过大会直接导致 TLS 握手失败。第二步在 WSL2 里用curl -I看返回头确认服务端到底有没有响应。如果 curl 正常而浏览器打不开问题几乎都集中在证书和页面脚本上。第三步更新一下 CA 证书包在 Ubuntu 里就是sudo apt update sudo apt install ca-certificates。第四步如果着急做题可以先回到 Windows 宿主机的浏览器里打开页面评测入口和本地编译环境是两回事WSL2 里编译好可执行文件之后提交还是在网页端完成的。这条经验不是算法内容但很实用。我记得第一次在 WSL2 里刷老题时也懵了很久后来发现只是证书问题浪费了大半小时。4.3 边界条件与 EPS 选择总结边界条件单独列一个表做题时候对着查很方便情况概率值处理建议D 00%直接特判不进积分D ≥ √(W² H²)100%直接特判不进积分0 D √(W²H²)按积分公式用自适应辛普森sqrt 内部接近负数无强制置零防 NaN积分结果超出 [0,1]理论不可能输出前夹逼eps 的选择我推荐 1e-10。你可能会想既然要求精度高一点取 1e-12 不是更好吗实测看1e-12 会让自适应辛普森在某些区间多递归两到三层题目多组数据时整体时间会翻倍。而 1e-10 的结果在 double 精度下已经足够撑起题目要求的小数位数完全没必要更小。还有一个不起眼但很重要的点W 和 H 在题面里可能是整数输入时用%lf读入也能正确处理整数别因为类型不匹配导致读入失败。这个细节同样能让人白交好几发。5. 如果不想用数值积分分段解析解法5.1 分段的位置来自哪里自适应辛普森能 AC但有些朋友可能会想这题能不能完全不靠数值积分推出一个闭式公式来算答案是能而且分段解析解的推导过程能加深对几何积分本身的理解。外层的被积函数里只有一个分段点就是 Y(x) 什么时候取 H什么时候取 sqrt(D²−x²)。当 sqrt(D²−x²) ≥ H 时圆边界超出了矩形上边界Y(x) 被截断为 H。这个条件等价于 x ≤ sqrt(D² − H²)。所以令x0 sqrt(max(0, D² − H²))当 x ∈ [0, min(x0, W, D)] 时Y(x)H当 x ∈ [min(x0, W, D), min(W, D)] 时Y(x)sqrt(D²−x²)。这个 x0 就是分段点。如果 D ≤ H那么 x0 不存在整个区间上 Y(x) 都是 sqrt(D²−x²)分段只有一段积分还要更好算。5.2 每一段怎么积分第一段 YH 时被积函数退化成(W−x) · (H² − H²/2) 0.5·H²·(W−x)这是简单二次多项式原函数可以直接写∫ 0.5·H²·(W−x) dx 0.5·H²·(W·x − 0.5·x²)代入上下限即可。第二段 Yssqrt(D²−x²) 时被积函数展开成四项(W−x)(H·s − 0.5·s²) H·W·s − H·x·s − 0.5·W·s² 0.5·x·s²然后逐项积分。用换元 xD·sin θ可以得到下面这几个基本原函数∫ sqrt(D²−x²) dx 0.5·(x·sqrt(D²−x²) D²·asin(x/D)) C∫ x·sqrt(D²−x²) dx −(D²−x²)^(3/2)/3 C∫ (D²−x²) dx D²·x − x³/3 C∫ x·(D²−x²) dx D²·x²/2 − x⁴/4 C把这些组合起来就是分段解析积分的完整闭式。理论上代码也能写出来而且运行速度比自适应辛普森更快。但我个人不推荐在竞赛里这么写。原因很实在公式推导容易错错在某个符号上就得调试半天而分段边界本身又要处理 D、W、H 的相对大小关系代码长度和出错概率都会上升。自适应辛普森虽然是个数值方法但实现固定、逻辑单一不用动脑推导换来的是稳定可靠。5.3 数值积分在这里为什么是更务实的选型做算法题要区分“理解用”和“提交用”。解析解用来理解这道题非常合适它能让你彻底看清分段点和积分精度从哪来但提交用我多半还是会选自适应辛普森因为数值方法把“求原函数”这个最容易出错的部分完全绕开了。我常跟身边人讲数值积分在 OI/ACM 里是那种“下限高、上限中庸”的工具。下限高是因为你只要会写二三十行递归模板就能处理一大类积分题上限中庸是因为遇到精度要求极高、或者函数有严重奇异性的题自适应辛普森可能会卡死或者精度不足。放在 UVa 11355 这个题上函数只有一处分段且整体光滑自适应辛普森属于杀鸡用牛刀但刀刀见血。如果你想把这道题吃透可以自己动手把第 5.2 节的解析公式实现一遍然后用它和蒙特卡洛、自适应辛普森三种方法互相验证。三个结果能对上你对概率密度的理解就到位了。我个人在实际做题中的体会是Cool Points 这道题真正的价值不在 AC 本身而在于它把“概率密度”“卷积”“数值积分”“边界特判”四件事串在了一起。以后你再遇到类似“随机取点、求距离/面积期望”的题脑子里会自动浮现这条处理链。先用蒙特卡洛验证直觉再写密度函数最后套一个合适的积分工具基本不会错。