计算圆周率的外星人程序

问题的来源

徐少侠朋友在我的上一篇随笔的3楼的评论中给出了一个只有四行的计算圆周率的 C 程序。我在网络上查找了一下来源，找到这么一个计算圆周率的外星人程序：

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

这个只有 158 个字符的 C 程序将输出 2,400 位圆周率的值。上述程序为了追求短小，格式很不规范。稍微整理一下，得到如下程序：

1:  #include <stdio.h>
2:  
3:  int main()
4:  {
5:    int f[8401], a, b, c = sizeof(f) / sizeof(f[0]) - 1, d, e, g;
6:    for (a = 10000, b = e = 0; b != c; ) f[b++] = a / 5;
7:    for (; g = c * 2; c -= 14, printf("%.4d", e + d / a), e = d % a)
8:      for (d = 0, b = c; d += f[b] * a, f[b] = d % --g, --b; d *= b) d /= g--;
9:  }

但是从这个程序中还是看不什么名堂来。

计算圆周率的一个公式

现在让我们来看看计算圆周率以下公式：

下面的 C 程序用于验证这个公式：

1:  #include <stdio.h>
2:  
3:  void main()
4:  {
5:    double pi = 2;
6:    for (int k = 49; k > 0; k--) pi = pi * k / (2 * k + 1) + 2;
7:    printf("pi = %.15lf...\n", pi);
8:  }

这个程序是从上述公式的最里面一层开始计算的。运行结果为：

pi = 3.141592653589793...

可见所得到的圆周率的小数点后全部 15 位数字都是正确的，而这只需要迭代 49 次。但是 double 的精度就这么十五六位，如果我们要计算圆周率到小数点后二千多位，就必须使用大整数进行运算。下面就是实现这个公式的 C 程序：

01:  #include <stdio.h>
02:  
03:  const int DIGITS = 2400;  // must: DIGITS % LEN == 0
04:  const int BASE   = 10000; // BASE == 10 ** LEN
05:  const int LEN    = 4;
06:  const int TIMES  = 14;
07:  
08:  int get_item(int pi[], int n)
09:  {
10:    int item = 0;
11:    for (int k = n - 1; k >= 0; k--) {
12:      item += pi[k] * BASE;
13:      pi[k] = item % (k * 2 + 1);
14:      item /= (k * 2 + 1);
15:      if (k > 0) item *= k;
16:    }
17:    return item;
18:  }
19:  
20:  int main()
21:  {
22:    int pi[DIGITS / LEN * TIMES], n = sizeof(pi) / sizeof(pi[0]);
23:    for (int i = 0; i < n; i++) pi[i] = BASE / 5; // BASE / 5 == 2.000
24:    for (int remainder = 0; n > 0; n -= TIMES) {
25:      int item = get_item(pi, n);
26:      printf("%.*d", LEN, remainder + item / BASE);
27:      remainder = item % BASE;
28:    }
29:  }

这个程序的运行结果和本文开头的外星人程序的运行结果一模一样，也是输出 2,400 位圆周率的值。在上述程序中：

• 第 05 行，指定每次循环计算的位数 LEN = 4。
• 第 04 行，指定计算的基数 BASE ，这个基数是 10 的 LEN 次方。
• 第 03 行，指定要计算的圆周率的位数 DIGITS，必须是 LEN 倍数。
• 第 06 行，指定为了达到所需的精度，必须将上述公式迭代的次数是循环次数的多少倍: TIMES 。
• 第 22 行，定义数组 pi ，用来存放进位等中间结果。以及数组 pi 的大小 n 。
• 第 23 行，将数组 pi 的值初始化为 2 ，因为上述公式中每次迭代都要加 2 。
• 第 24 到 27 行，主循环，每次循环输出 LEN 位圆周率的值。
• 第 08 到 18 行，获取每一项的值，供主循环使用。
• 第 11 到 16 行，按照上述公式执行迭代过程。从最里面一层开始迭代。
• 第 12 行，执行每次迭代中的加 2 。
• 第 13 行，在每次迭代中设置进位，以便进行下一次迭代。
• 第 14 行，执行每次迭代中的除以 2k + 1 。
• 第 15 行，执行每次迭代中的乘以 k 。

翻译为外星人程序

将上述 C 程序的变量作以下替换：

BASE	a
k	b
n	c
item	d
remainder	e
pi	f
k * 2 + 1	g

再作一些小调整，基本上就翻译为外星人程序了。外星人程序追求的是源程序的最短，用了很多巧妙的方法，但是非常难于理解。如果不是参加源程序短小的竞赛的话，不建议这样写程序。

进一步的话题

我把上面最后一个 C 程序第 03 行 DIGITS 的值改为 10,000 后，发现输出的圆周率的值还是正确的。所以说实际上这个外星人程序至少可以计算圆周率到小数点后一万位。网络上流传的 2,400 位，甚至 800 位，太保守了。