使用快速傅里叶变换计算圆周率

前几天，我发表了一篇随笔：“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   {
     //= pp.780-781, mppi, 20.6 任意精度的运算
     /// <summary>
     /// 计算圆周率到小数点后 digits 位数字。
     /// 计算结果保存在返回的字节数组中，每个字节存放两个十进制数字，字节数组的第一个元素是“03”。
     /// 字节数组的长度可能大于 (digits + 1) / 2 + 1，只保证小数点后前 digits 个十进制数字是准确的。
     /// </summary>
     /// <param name="digits">小数点后的十进制数字个数</param>
     /// <returns>存放圆周率的字节数组</returns>
     public static byte[] Compute(int digits)
     {
       if (digits < 0) throw new ArgumentOutOfRangeException("digits", "can't less than zero");
       int n = Math.Max(5, (digits + 1) / 2 + 2);
       byte[] pi = new byte[n + 1];
       byte[] x = new byte[n + 1], y = new byte[n << 1];
       byte[] sx = new byte[n], sxi = new byte[n];
       byte[] t = new byte[n << 1], s = new byte[3 * n];
       t[0] = 2;                             // t = 2
       BigArithmetic.Sqrt(x, x, n, t, n);    // x = sqrt(2)
       BigArithmetic.Add(pi, t, x, n);       // pi = 2 + sqrt(2)
       Array.Copy(pi, 1, pi, 0, n);
       BigArithmetic.Sqrt(sx, sxi, n, x, n); // sx = sqrt(x)
       Array.Copy(sx, y, n);                 // y = sqrt(x)
       for (; ; )
       {
         BigArithmetic.Add(x, sx, sxi, n);       // x = sqrt(x) + 1/sqrt(x)
         Array.Copy(x, 1, x, 0, n);
         BigArithmetic.Divide(x, x, n, 2);       // x = x / 2
         BigArithmetic.Sqrt(sx, sxi, n, x, n);   // sx = sqrt(x), sxi = 1/sqrt(x)
         BigArithmetic.Multiply(t, y, n, sx, n); // t = y * sx
         Array.Copy(t, 1, t, 0, n);
         BigArithmetic.Add(t, t, sxi, n);        // t = y * sx + sxi
         x[0]++;
         y[0]++;
         BigArithmetic.Inverse(s, n, y, n);      // s = 1 / (y + 1)
         Array.Copy(t, 1, t, 0, n);
         BigArithmetic.Multiply(y, t, n, s, n);  // y = t / (y + 1)
         Array.Copy(y, 1, y, 0, n);
         BigArithmetic.Multiply(t, x, n, s, n);  // t = (x + 1) / (y + 1)
         int mm = t[1] - 1;                      // 若 t == 1 则收敛
         int j = t[n] - mm;
         if (j > 1 || j < -1)
         {
           for (j = 2; j < n; j++)
           {
             if (t[j] != mm)
             {
               Array.Copy(t, 1, t, 0, n);
               BigArithmetic.Multiply(s, pi, n, t, n); // s = t * pi
               Array.Copy(s, 1, pi, 0, n);             // pi = t * pi
               break;
             }
           }
           if (j < n) continue;
         }
         break;
       }
       return pi;
     }
   }
 }

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

下面是测试程序：

 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   {
     static void Main(string[] args)
     {
       var stopWatch = Stopwatch.StartNew();
       int digits = (args.Length > 0) ? int.Parse(args[0]) : 10000;
       var bs = Pi.Compute(digits);
       var sb = new StringBuilder();
       sb.AppendFormat("Pi = {0}.{1}", bs[0], Environment.NewLine);
       for (int i = 1; i <= (digits + 1) / 2; i++)
       {
         sb.Append(bs[i].ToString("D2"));
         if (i % 25 == 0) sb.AppendLine();
         else if (i % 5 == 0) sb.Append(' ');
       }
       stopWatch.Stop();
       if ((digits + 1) / 2 % 25 != 0) sb.AppendLine();
       sb.AppendFormat("DIGITS:{0:N0} ELAPSED:{1}", digits, stopWatch.Elapsed);
       Console.Write(sb);
     }
   }
 }

运行结果如下：

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#程序”这篇随笔中的结果相对照。