神奇的求圆周率代码

今天是圆周率日,我们就来讲一讲关于圆周率的故事吧。
很久以前,博主看到了一段惊为天人的代码,可以求出圆周率的前800位(修改参数可以实现任意精度):

1
2
3
int a=10000,b,c=2800,d,e,f[2801],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);}

经过考证,这段代码的作者是Dik T. Winter。下面我们来看看它是如何工作的。

在C语言中,for循环和while循环可以互相代替。这里涉及到for循环的执行顺序。例如

1
2
3
for (statement1; statement2; statement3) {
statements;
}

上面的for语句可以用下面的while语句来代替:

1
2
3
4
5
statement1;
while (statement2) {
statements;
statement3;
}

而逗号运算符和--运算符的功能也是需要注意的。要写出这样怪异的C程序,逗号运算符无疑是一个好的助手,它的作用是:从左到右依次计算各个表达式的值,并且返回最右边表达式的值。

了解这些内容后,我们将这段代码中的for循环展开,等价地写为while的形式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#include <stdio.h>
int a = 10000, b, c = 2800, d, e, f[2801], g;
main() {
for (int i = 0; i < c; i++) {
f[i] = a / 5;
}
while (c != 0) {
d = 0;
g = c * 2;
b = c;
while (1) {
d += f[b] * a;
g--;
f[b] = d % g;
d /= g;
g--;
b--;
if (b == 0) break;
d *= b;
}
c -= 14;
printf("%.4d", e + d / a);
e = d % a;
}
}

进一步观察,不难发现这里的a是常数,而g可以方便的代换掉。继续改写代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <stdio.h>
int b, c = 2800, d, e, f[2801];
main() {
for (int i = 0; i < c; i++) {
f[i] = 2000;
}
while (c != 0) {
b = c;
d = f[b] * 10000;
f[b] = d % (b * 2 - 1);
d /= b * 2 - 1;
b--;
while (1) {
d = d * b + f[b] * 10000;
f[b] = d % (b * 2 - 1);
d = d / (b * 2 - 1);
b--;
if (b == 0) break;
}
c -= 14;
printf("%.4d", e + d / 10000);
e = d % 10000;
}
}

不过,即使是这样,这段代码还是让人看得很迷糊。我们不妨开个挂,进行倒推吧。
我们来看看这段代码的数学对应。这是一个计算$\pi$的公式:
$$\frac{\pi}{2}=\sum_{k=0}^{\infty}\frac{k!}{(2k+1)!!}=\sum_{k=0}^{\infty}\frac{2^{k}k!^{2}}{(2k+1)!}=1+\frac{1}{3}\left(1+\frac{2}{5}\left(1+\frac{3}{7}(1+\cdots)\right)\right)$$
可以进一步写为
$$\pi=2\left(1+\frac{1}{3}\left(1+\frac{2}{5}\left(1+…\left(1+\frac{2799}{2\cdot2799+1}(1+0)\right)…\right)\right)\right)$$
根据维基百科的相关内容,该公式最早似乎由牛顿推出,而Wolfram收录的一个变式则写明作者是Beeler等人。参见PiFormulas
在这个迭代中,我们需要不断地生成形如$\frac{k}{2k+1}$的分式。对应于展开后的代码,不难发现g就是分母,而b则是分子。
当然,这还不是足够的。一个重要的问题是如何在逐次迭代中将误差进行传递,以实现任意精度的计算。整个程序中的唯一一个关系到输出的printf函数,每次会输出四位整数。那么c -= 14又是怎么来的呢?我们注意到$2^{14}>10^4=10000$,而$14\times 200=2800$,$200\times 4=800$,这些『神秘数字』就都说得通了。每次循环用14项获得4位精确的$\pi$值,它们是一个八位整数的前四位,因此需要模10000;然后将后四位乘以10000,代入下一次循环。

读者可以自行体会这个算法的精妙,特别是它是如何在不损失精度的情况下,完成整个计算过程。


参考文章:
Computing Pi in C
求圆周率Pi一万位程序分析

拓展阅读:
编程计算Pi的N种方法
Tiny programs for constants computation
圓周率 - 维基百科
圆周率近似值 - 维基百科
梅欽類公式 - 维基百科
The Pi-Search Page

🍭支持一根棒棒糖!
0%