神奇的求圆周率代码

今天是圆周率日,我们就来讲一讲关于圆周率的故事吧。
很久以前,笔者看到了一段惊为天人的代码,可以求出圆周率的前 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;
}
}

不过,即使是这样,这段代码还是让人看得很迷糊。我们不妨开个挂,进行倒推吧。
我们来看看这段代码的数学对应。这是一个计算的公式:

可以进一步写为

根据维基百科的相关内容,该公式最早似乎由牛顿推出,而 Wolfram 收录的一个变式则写明作者是 Beeler 等人(参见 PiFormulas)。
在这个迭代中,我们需要不断地生成形如的分式。对应于展开后的代码,不难发现 g 就是分母,而 b 则是分子。
当然,这还不是足够的。一个重要的问题是如何在逐次迭代中将误差进行传递,以实现任意精度的计算。整个程序中的唯一一个关系到输出的 printf 函数,每次会输出四位整数。那么 c -= 14 又是怎么来的呢?我们注意到,而,这些「神秘数字」就都说得通了。每次循环用 14 项获得 4 位精确的值,它们是一个八位整数的前四位,因此需要模 10000;然后将后四位乘以 10000,代入下一次循环。

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

最后再附上一个 Python 版本:

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
#!/usr/bin/env python3

a = 10000
c = 2800
e = 0
f = [a // 5 for i in range(c + 1)]

while c != 0:
d = 0
g = c * 2
b = c
while True:
d += f[b] * a
g -= 1
f[b] = d % g
d //= g
g -= 1
b -= 1
if b == 0:
break
d *= b

c -= 14
print("%.4d" % (e + d // a), end = "")
e = d % a

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

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