今天逛知乎的时候看到了一份计算 $\pi$ 的 代码。我试图跑了一下，发现其输出的 2400 位居然全部是对的。我试图学习了一下其中的姿势，感觉还是很神奇的。

int a=10000,b,c=8400,d,e,f[8401],g;main(){
    for(;b-c;)f[b++]=a/5;
    for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
        for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

然后我就试图分析这段代码到底是在干啥。

首先我们发现这代码中只有一个 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 的值确实是一个好近似。