（注意，题解中的数学推导可能并不严谨。）

随机游走的结论是，两维的坐标是独立的，均值为 $0$，方差为 $\frac{T}{2}$ 的正态分布。即 $(x,y)$ 点的概率密度函数为 $$ f(x,y)=\frac{1}{\pi T}e^{-\frac{x^2+y^2}{T}} $$ 所求的概率即为 $f(x,y)$ 在一个半径为 $R$，圆心离原点距离为 $d$ 的圆内的积分。用格林公式得到上述积分等于 $\frac{1}{2\pi}\frac{1-e^{-\frac{x^2+y^2}{T}}}{x^2+y^2}(x\mathrm d y-y\mathrm dx)$ 在圆周上的积分。直接使用自适应 Simpson 积分在 $[0,2\pi]$ 上积分即可。注意当 $x^2+y^2\to 0$ 时，需要用到 $\frac{1-e^{-\frac{x^2+y^2}{T}}}{x^2+y^2}$ 的极限等于 $\frac{1}{T}$。

第二项则是上述的结果关于 $T$ 再做一个积分。如果再套一层自适应积分，效率可能太低，所以我们交换顺序，先对 $T$ 进行积分，然后进行相同的自适应 Simpson 积分。即要求 $$ \int_0^T \frac{1-e^{-\frac{r^2}{t}}}{r^2}\mathrm dt $$ 这个积分不是初等函数，结果是 $$ \frac{T}{r^2}\left(1-e^{-\frac{r^2}{T}}\right)-\mathrm{Ei}\left(-\frac{r^2}{T}\right) $$ 其中 $\mathrm{Ei}$ 是指数积分，可以用 C++ 标准库 `std::expint` 来计算。注意这个函数在 $r\to 0$ 时是发散的，因为当 $x\to 0$ 时 $$ \mathrm{Ei}(x)\sim \gamma+\ln|x|+O(x^2) $$ 但是其在 $x\in(-\epsilon,\epsilon)$ 处的积分是 $O(\epsilon |\ln\epsilon|)$ 的，我们可以合理地排除掉 $r^2$ 很接近于 $0$ 的部分。事实上，只要保证自适应 Simpson 积分时不会恰好命中 $r=0$ 导致算出 $+\infty$ 从而无限迭代，就能够计算出正确的结果。而考虑到浮点数自带的误差，不需要任何特殊处理也能够通过。