计算圆周率

在 Sphere Online Judge (SPOJ) 网站上有这么一道题目：

Digits of Pi: In this problem you have to find as many digits of PI as possible.

Output: Output must contain as many digits of PI as possible (not more than 1,000,000).

Score: The score awarded to your program will be the first position of the digit where the first difference occured.

Time limit: 25s

Source limit: 4,096 Bytes

这道题目要求你在 25 秒之内尽可能地输出圆周率 π 的有效数字。注意题目中将源程序的大小限制在 4KB 以内。下面的程序来源于在2005年9月的随笔“计算圆周率的C程序”，作了少许改动，并增加了一个 pi2 函数：

01:  #include <stdio.h>
02:  #include <stdlib.h>
03:  
04:  const int DIGITS = 21987;
05:  
06:  void pi2()
07:  {
08:    static char s[] = ")[c?jEia#44#34R......f%Zf"; // 2,762 chars => 5,000 digits
09:    for (int i = 0; i < sizeof(s) - 1; i++) {
10:      if (s[i] == '#') printf("0%c", s[++i]);
11:      else {
12:        int v = s[i];
13:        if (v == '$') v = '\\';
14:        printf("%02d", v - '%' + 10);
15:      }
16:    }
17:  }
18:  
19:  int main()
20:  {
21:    int t0[] = {176, 28, 48, 96}, k0[] = {1, 1, 0, 1}, n0[] = {57, 239, 682, 12943};
22:    int m, n, r, s, i, j, k, p, d = DIGITS, z = sizeof(t0) / sizeof(t0[0]);
23:    int* t = (int *)calloc((d += 5) + 1, sizeof(int));
24:    int* pi = (int *)calloc(d + 1, sizeof(int));
25:    for (i = d; i >= 0; i--) pi[i] = 0;
26:    for (p = 0; p < z; p++) {
27:      for (k=k0[p], n=n0[p], t[i=j=d]=t0[p], i--; i >= 0; i--) t[i] = 0;
28:      for (r = 0, i = j; i >= 0; i--) {
29:        r = (m = 10 * r + t[i]) % n;
30:        t[i] = m / n;
31:        k ? (pi[i] += t[i]) : (pi[i] -= t[i]);
32:      }
33:      while (j > 0 && t[j] == 0) j--;
34:      for (k = !k, s = 3, n *= n; j > 0; k = !k, s += 2) {
35:        for (r = 0, i = j; i >= 0; i--) {
36:          r = (m = 10 * r + t[i]) % n;
37:          t[i] = m / n;
38:        }
39:        while (j > 0 && t[j] == 0) j--;
40:        for (r = 0, i = j; i >= 0; i--) {
41:          r = (m = 10 * r + t[i]) % s;
42:          m /= s;
43:          k ? (pi[i] += m) : (pi[i] -= m);
44:        }
45:      }
46:    }
47:    for (n = i = 0; i <= d; pi[i++] = r) {
48:      n = (m = pi[i] + n) / 10;
49:      if ((r = m % 10) < 0) r += 10, n--;
50:    }
51:    printf("3."); 
52:    for (i = d - 1; i >= 5; i--) putchar((int)pi[i] + '0');
53:    pi2();
54:    return 0;
55:  }

在该网站提交，运行时间 24.04 秒，内存占用 1.8 MB，成绩是“26,989”，目前在 C99 strict 语言中排名第一位。但是这道题目的最好成绩是“720,002”，使用 Haskell 语言，运行时间 24.35 秒，内存占用 41 MB。不知道他们使用什么算法，成绩这么好。 :)

在上述 C 程序中，算法是用泰勒公式计算反正切值，计算到小数点后 21,987 位(源程序第 04 行，第 21 到 52 行)。使用的公式是(第 21 行指定公式的参数)：

π = 176 * arctan(1/57) + 28 * arctan(1/239) - 48 * arctan(1/682) + 96 * arctan(1/12943)  [Stormer]

上述公式中的每一项使用泰勒公式展开如下：

c * arctan(1/x) = c/x - c/(3*x³) + c/(5*x⁵) - c/(7*x⁷) + ...

然后使用 pi2 函数接着输出 5,000 位(第 53 行，第 06 到 17  行)。总共输出 26,989 位(包括“3.”这两位，第 51 行)。

在 pi2 函数中首先使用一个长度为 2,762 的静态字符串来表示事先计算出来的 5,000 位数字(第 08 行)。由于前 128 个 ASCII 码中只有 95 个可打印字符(从 0x32 到 0x7E)，所以使用 0x37(%) 到 0x7E 这 90 个连续的 ASCII 字符表示 10 到 99 (第 12 到 14 行)。为了避免‘\’的特殊含义，把其中的 0x5C(\) 用 0x36($) 代替(第 13 行)。使用 0x35(#) 后面跟一位数字来表示 00 到 09 (第 10 行)。

此外，根据“计算圆周率的C#程序”这篇随笔稍做改动的 C# 程序，在该网站提交后的成绩是“16,278”，运行时间 24.40 秒，内存占用 12 MB，目前在 C# 语言中排名第一位。

参考资料

1. 维基百科: 圆周率
2. Wikipedia: Pi
3. Wikipedia: Machin-like formaula
4. Wikipedia: Taylor series