一段关于 $\pi$ 的代码
今天逛知乎的时候看到了一份计算 \(\pi\) 的 代码。我试图跑了一下，发现其输出的 2400 位居然全部是对的。我试图学习了一下其中的姿势，感觉还是很神奇的。
1 | int a=10000,b,c=8400,d,e,f[8401],g;main(){ |
然后我就试图分析这段代码到底是在干啥。
首先我们发现这代码中只有一个 printf ，其输出内容是 e+d/a 。观察可知，e 只在一个地方变过，就是 e=d%a 。
再观察 f 的变化。首先可以看到 d+=f[b]*a ，然后我们改变了 f[b] 。这初看起来很像高精操作，但是发现每次 f[b] 模的数是不同的，这让代码并不是那么好分析。但是注意到每个 f[b] 都是乘了 a 。如果我们把 f 看成是 \(\pi\) 的某种近似的表现形式的话，这段代码就可以理解了：第四行的 for 实际上是在计算 f*a ，其结果取下整后存入 d，小数部分继续放在 f 里面。
接下来就要猜 f 是怎样一个存储方式。我们注意到，第四行里面的 g 就是 2b-1，而 f[b] 一旦减少了 g，f[b-1] 就会增加 b-1 （而不是 1，因为第四行的这个 for 最后面会令 d*=b）。于是我们可以大胆猜测：令 \(F\) 为 f 表示的数，则有
\[F=\sum_{b=1}^c f_b w_b,\]
其中 \(w_i\) 为这个程序的一组参数。由 \(f_b\) 与 \(f_{b-1}\) 的换算规则可以推出：
\[(2b-1)w_b = (b-1) w_{b-1},\]
稍稍一算即可知
\[w_b = \frac{(b-1)!}{(2b-1)!!} w_1.\]
这样所有的 \(w\) 就全部依赖于 \(w_1\) 了。最后注意到，初始时的 \(f_1\) 每增加 1，输出的数增加 \(1\) （用 \(b=1\) 模拟一遍第四行即可），故有 \[w_1 = 1.\]
由于这里的 c 只是一定程度上控制的循环次数，我们可以继续探究 c 无穷大时 \(F\) 的值。
为了研究 \(F\) ，我们需要一个更好的 \(w_b\) 的表现形式。注意到：
\[
\int_0^{\pi/2} \sin^{2n-1} x\ \text{d}x = \frac{(2n-2)!!}{(2n-1)!!} = 2^{n-1} \frac{w_n}{w_1},
\]
故当 c 趋于 \(\infty\) 时我们有：
\[
F = \sum_{b=1}^\infty f_b w_b = \sum_{b=1}^\infty \frac{a}{5} 2^{1-b} \int_0^{\pi/2} \sin^{2b-1} x\ \text{d}x = \frac{a}{5} \sum_{b=1}^\infty \int_0^{\pi/2} 2^{1-b} \sin^{2b-1} x\ \text{d}x,
\]
我懒得讨论这里的积分号与求和号是否可交换了，让我们先试试吧！
\[
\begin{aligned}
F &= \frac{a}{5} \int_0^{\pi/2}\sum_{b=1}^\infty 2^{1-b} \sin^{2b-1} x\ \text{d}x \\
&= \frac{a}{5} \int_0^{\pi/2} \sin x \sum_{b=0}^\infty \left( \frac{\sin^2 x}{2}\right)^b \text{d}x\\
&= \frac{a}{5} \int_0^{\pi/2} \frac{\sin x}{1 - \frac{\sin^2 x}{2} } \text{d}x \\
& = \frac{a}{5} \int_0^{\pi/2} \frac{2}{1 + \cos^2 x} \text{d} (-\cos x) \\
& = \frac{2 a}{5} \left(\arctan \left(-\cos \frac{\pi}{2} \right) - \arctan (-\cos 0)\right) \\
& = \frac{a}{10} \pi.
\end{aligned}
\]
于是我们发现，\(F = 1000 \pi\)！这也和程序的输出吻合，因为它从 3141 开始输出。我们还可以发现，\(w_b\) 为指数级衰减（因为有 \(2^{n-1} \frac{w_n}{w_1} = \frac{(2n-2)!!}{(2n-1)!!} < 1\)），故上限为 c 的值确实是一个好近似。