银河

SKYIV STUDIO

  博客园 :: 首页 :: 博问 :: 闪存 :: :: :: 订阅 订阅 :: 管理 ::

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*x3) + c/(5*x5) - c/(7*x7) + ...

然后使用 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
posted on 2010-07-13 00:04  银河  阅读(4132)  评论(17编辑  收藏  举报