0%

一段关于 $\pi$ 的代码

今天逛知乎的时候看到了一份计算 \(\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] 一旦减少了 gf[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,输出的数增加 \(\frac{10}{a}\),故有 \[w_1 = \frac{10}{a}.\]

由于这里的 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 2^{2-b} \int_0^{\pi/2} \sin^{2b-1} x\ \text{d}x, \] 我懒得讨论这里的积分号与求和号是否可交换了,让我们先试试吧! \[ \begin{array}{rl} F &= 2 \int_{0}^{\pi/2} \sum_{b=1}^\infty 2^{1-b} \sin^{2b-1} x\ \text{d}x \\ &= 2 \int_0^{\pi/2} \sin x \sum_{b=0}^\infty \left( \frac{\sin^2 x}{2}\right)^b \text{d}x\\ &= \int_0^{\pi/2} \frac{2 \sin x}{1 - \frac{\sin^2 x}{2} } \text{d}x \\ & = \int_0^{\pi/2} \frac{4}{1 + \cos^2 x} \text{d} (-\cos x) \\ & = 4\left(\arctan (-\cos \frac{\pi}{2}) - \arctan (-\cos 0)\right) \\ & = \pi. \end{array} \]

于是我们发现,\(F = \pi\)!而且观察发现,\(w_b\) 衰减的非常快,为指数级衰减,故上限为 c 的值确实是一个好近似。