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

前几天，我发表了一篇随笔：“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#程序”这篇随笔中的结果相对照。