学步园

前几天，我发表了一篇随笔：“BigArithmetic - 提供任意精度的算术运算的静态类”。现在，让我们使用 BigArithmetic 类来计算圆周率。

我们需要一个计算 π 的分析算法。有用的算法是二次收敛的，即每一次迭代使有效位数增加一倍。计算 π 的二次收敛算法是基于 ACM 法(算术几何平均法)。首先设置初始值为：

然后，当 i = 0, 1, ... 重复迭代：

直到足够的精度。

下面就是源程序代码：

 1 using System;
 2 
 3 namespace Skyiv.Numeric
 4 {
 5   /// <summary>
 6   /// 圆周率
 7   /// </summary>
 8   static class Pi
 9   {
10     //= pp.780-781, mppi, 20.6 任意精度的运算
11     /// <summary>
12     /// 计算圆周率到小数点后 digits 位数字。
13     /// 计算结果保存在返回的字节数组中，每个字节存放两个十进制数字，字节数组的第一个元素是“03”。
14     /// 字节数组的长度可能大于 (digits + 1) / 2 + 1，只保证小数点后前 digits 个十进制数字是准确的。
15     /// </summary>
16     /// <param name="digits">小数点后的十进制数字个数</param>
17     /// <returns>存放圆周率的字节数组</returns>
18     public static byte[] Compute(int digits)
19     {
20       if (digits < 0) throw new ArgumentOutOfRangeException("digits", "can't less than zero");
21       int n = Math.Max(5, (digits + 1) / 2 + 2);
22       byte[] pi = new byte[n + 1];
23       byte[] x = new byte[n + 1], y = new byte[n << 1];
24       byte[] sx = new byte[n], sxi = new byte[n];
25       byte[] t = new byte[n << 1], s = new byte[3 * n];
26       t[0] = 2;                             // t = 2
27       BigArithmetic.Sqrt(x, x, n, t, n);    // x = sqrt(2)
28       BigArithmetic.Add(pi, t, x, n);       // pi = 2 + sqrt(2)
29       Array.Copy(pi, 1, pi, 0, n);
30       BigArithmetic.Sqrt(sx, sxi, n, x, n); // sx = sqrt(x)
31       Array.Copy(sx, y, n);                 // y = sqrt(x)
32       for (; ; )
33       {
34         BigArithmetic.Add(x, sx, sxi, n);       // x = sqrt(x) + 1/sqrt(x)
35         Array.Copy(x, 1, x, 0, n);
36         BigArithmetic.Divide(x, x, n, 2);       // x = x / 2
37         BigArithmetic.Sqrt(sx, sxi, n, x, n);   // sx = sqrt(x), sxi = 1/sqrt(x)
38         BigArithmetic.Multiply(t, y, n, sx, n); // t = y * sx
39         Array.Copy(t, 1, t, 0, n);
40         BigArithmetic.Add(t, t, sxi, n);        // t = y * sx + sxi
41         x[0]++;
42         y[0]++;
43         BigArithmetic.Inverse(s, n, y, n);      // s = 1 / (y + 1)
44         Array.Copy(t, 1, t, 0, n);
45         BigArithmetic.Multiply(y, t, n, s, n);  // y = t / (y + 1)
46         Array.Copy(y, 1, y, 0, n);
47         BigArithmetic.Multiply(t, x, n, s, n);  // t = (x + 1) / (y + 1)
48         int mm = t[1] - 1;                      // 若 t == 1 则收敛
49         int j = t[n] - mm;
50         if (j > 1 || j < -1)
51         {
52           for (j = 2; j < n; j++)
53           {
54             if (t[j] != mm)
55             {
56               Array.Copy(t, 1, t, 0, n);
57               BigArithmetic.Multiply(s, pi, n, t, n); // s = t * pi
58               Array.Copy(s, 1, pi, 0, n);             // pi = t * pi
59               break;
60             }
61           }
62           if (j < n) continue;
63         }
64         break;
65       }
66       return pi;
67     }
68   }
69 }

这个程序很简单，就是按照前面给出的公式进行计算而已。

下面是测试程序：

 1 using System;
 2 using System.Text;
 3 using System.Diagnostics;
 4 using Skyiv.Numeric;
 5 
 6 namespace Skyiv
 7 {
 8   sealed class Test
 9   {
10     static void Main(string[] args)
11     {
12       var stopWatch = Stopwatch.StartNew();
13       int digits = (args.Length > 0) ? int.Parse(args[0]) : 10000;
14       var bs = Pi.Compute(digits);
15       var sb = new StringBuilder();
16       sb.AppendFormat("Pi = {0}.{1}", bs[0], Environment.NewLine);
17       for (int i = 1; i <= (digits + 1) / 2; i++)
18       {
19         sb.Append(bs[i].ToString("D2"));
20         if (i % 25 == 0) sb.AppendLine();
21         else if (i % 5 == 0) sb.Append(' ');
22       }
23       stopWatch.Stop();
24       if ((digits + 1) / 2 % 25 != 0) sb.AppendLine();
25       sb.AppendFormat("DIGITS:{0:N0} ELAPSED:{1}", digits, stopWatch.Elapsed);
26       Console.Write(sb);
27     }
28   }
29 }

运行结果如下：

Pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
...
4610126483 6999892256 9596881592 0560010165 5256375678
DIGITS:10,000 ELAPSED:00:00:04.0998524

这个结果可以和我在2005年9月30日写的“计算圆周率的C#程序”这篇随笔中的结果相对照。