在 Sphere Online Judge (SPOJ) 网站上有这么一道题目:
这道题目要求你在 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 行指定公式的参数):
上述公式中的每一项使用泰勒公式展开如下:
然后使用 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# 语言中排名第一位。