UVa 11355 Cool Points：随机点距离概率与自适应辛普森积分实战

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, n=1000000): 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×10，D 取 10，跑一百万次。你会发现结果大概落在 0.5 到 0.6 之间。为什么不是 1？因为正方形里两个随机点的距离可以超过 10，只有落在对角线附近的那部分点对会超出阈值。等到后面公式推出来，再用蒙特卡洛结果去对比，基本就能确认公式和代码有没有写错。

这个习惯我一直保留。蒙特卡洛代码本身没技术含量，却能帮你挡住很多“公式推错了但自己没发现”的低级错误。

2. 核心转化：把两个点换成差向量

2.1 差向量的密度函数是怎么来的

两个点一个是 P1=(x1,y1)，一个是 P2=(x2,y2)。直接对四个坐标做联合密度会比较麻烦。更好的做法是换一个角度看问题：距离只和 Δx=x1−x2、Δy=y1−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,t+dt]”取极限，就得到概率密度：

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 最大到 W，y 最大到 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。如果 L+R 和 S 足够接近，就认为当前划分已经满足精度；否则对左右子区间继续递归，并且把误差阈值缩小一半。

判断“足够接近”的标准，我习惯用：

|L + R − S| ≤ 15·eps

如果满足，返回 L + R + (L+R−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 选择总结

边界条件单独列一个表，做题时候对着查很方便：

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 每一段怎么积分

第一段 Y=H 时，被积函数退化成：

(W−x) · (H² − H²/2) = 0.5·H²·(W−x)

这是简单二次多项式，原函数可以直接写：

∫ 0.5·H²·(W−x) dx = 0.5·H²·(W·x − 0.5·x²)

代入上下限即可。

第二段 Y=s=sqrt(D²−x²) 时，被积函数展开成四项：

(W−x)(H·s − 0.5·s²) = H·W·s − H·x·s − 0.5·W·s² + 0.5·x·s²

然后逐项积分。用换元 x=D·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 本身，而在于它把“概率密度”“卷积”“数值积分”“边界特判”四件事串在了一起。以后你再遇到类似“随机取点、求距离/面积期望”的题，脑子里会自动浮现这条处理链。先用蒙特卡洛验证直觉，再写密度函数，最后套一个合适的积分工具，基本不会错。