神奇的求圆周率代码

发表于 2019-03-14 分类于程序设计 Disqus：本文字数： 788 阅读时长 ≈ 3 分钟

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

statement1;
while (statement2) {
    statements;
    statement3;
}

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

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

#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 可以方便的代换掉。继续改写代码。

#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 版本：

#!/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