计算机最原始的用途就是进行人类无法完成的复杂运算，算 PI 就是这样的运算之一。虽然算 PI 本身没有多大的实际意义，但是对于计算机爱好者来说作为一种编程的挑战，还是很有意思的。算 PI 看似简单，其实它还牵涉到一些有用的数学知识。

第一类算法：arctan 的级数展开

• PI/4 = 4 arctan(1/5) - arctan(1/239) (1)
• arctan(x) = x - x³/3 + x⁵/5 - x⁷/7 + .... (2)

很容易想到，要得到超高精度的 PI 值，实数在计算机中必须以数组的形式进行存取，数组的大小跟所需的有效位数成正比。在这个算法中，PI 的有效位数 n 随 (2) 的求和项数线性增加。而为计算 (2) 中的每一项，需要进行超高精度实数除以小整数(5², 239², 2k+1)的循环，循环所需次数也跟 n 成正比。所以，这个算法总的时间复杂度为 O(n²)。

这个算法的优点是简单，而且只需要进行整数运算。下面给出我写的算 PI 程序。在程序中，我采用了一些提高速度的措施：超高精度实数以数组的形式进行存取，数组元素的类型为 64 位整数(long long)，每个元素储存 12 个十进制位；对 x^k (x = 1/5, 1/239) 的头部和尾部的 0 的数量进行估计，只对非 0 的部分进行计算。

• pi.cpp C++ 源程序，在 Linux 下以 g++ pi.cpp -o pi -O2 编译
• pi.s 在 g++ 生成的汇编程序的基础上进行修改，速度更快，在 Linux 下以 g++ pi.s -o pi 编译

另外，还有许多跟 (1) 类似的式子，但不常用。例如：

• PI/4 = arctan(1/2) + arctan(1/3)
• PI/4 = 8 arctan(1/10) - arctan(1/239) - 4 arctan(1/515)

第二类算法：与 1/PI 有关的级数

• 1/PI = (sqrt(8) / 9801) sum_k=0~inf { [(4k)! (1103 + 26390k)] / [(k!)⁴ 396^4k] } (Ramanujan)
• 1/PI = (sqrt(10005) / 4270934400) sum_k=0~inf { [(6k)! (13591409 + 545140134k)] / [(k!)³ (3k)! (-640320)^3k] } (Chudnovsky)

以上两个级数(还有其它类似形式的级数，但不常用)比起 arctan 的泰勒级数要复杂得多。虽然仍然是线性收敛，总的时间复杂度也仍然是 O(n²)，但它们的收敛速度相当快， (Ramanujan) 每项可以增加 8 位有效数字， (Chudnovsky) 每项可以增加 14 位。

在这个算法中，除了要进行超高精度实数(数组形式)和小整数的运算外，还有一次超高精度实数的开方和倒数的运算，这需要用到 FFT(快速傅立叶变换)，在下文叙述。

第三类算法：算术几何平均值和迭代法

算术几何平均值(Arithmetic-Geometric Mean, AGM) M(a, b) 定义如下：

• a₀ = a, b₀ = b
  a_k = (a_k-1 + b_k-1) / 2, b_k = sqrt(a_k-1 b_k-1)
  M(a, b) = lim_k->inf a_k = lim_k->inf b_k

然后，由椭圆积分的一系列理论(抱歉，过程我不懂)可以推导出如下公式：

• a₀ = 1, b₀ = 1 / sqrt(2)
  1/PI = { 1 - sum_k=0~inf [2^k (a_k² - b_k²)] } / 2M(a₀, b₀)² (AGM)

根据这条公式可以制定适当的迭代算法。在迭代过程中，有效位数随迭代次数按 2 的指数增加，即每迭代一次有效位数乘 2。算法中的超高精度实数的乘、除、开方等运算需要使用 FFT，在下文叙述。综合考虑 FFT 的时间复杂度，整个算法的时间复杂度约为 O(n log(n)²)。

除了 (AGM) 以外，还有其它的迭代序列，它们具有同样的时间复杂度。例如下面的这个序列将按 4 的指数收敛到 1/PI：

• y₀ = sqrt(2) - 1, a₀ = 6 - 4 sqrt(2)
  y_k = [1 - sqrt(sqrt(1 - y_k-1⁴))] / [1 + sqrt(sqrt(1 - y_k-1⁴))], a_k = (1 + y_k)⁴ a_k-1 - 2^2k+1 y_k (1 + y_k + y_k²)
  1/PI = lim_k->inf a_k (Borwein)

FFT

如上所述，第二和第三类算法不可避免地要涉及超高精度实数(数组形式存取的多位数)的乘、除、开方等运算。多位数乘法如果按照常规方法来计算，逐位相乘然后相加，其时间复杂度将达到 O(n²)。使用 FFT 可大大减少计算量。

设有复数数组 a[k] 和 b[k] (k=0~n-1)，正向和反向的离散傅立叶变换(DFT)定义如下： (i = sqrt(-1))

• b = FFT_forward(a) : b[k] = sum_j=0~n-1 ( a[j] e^-i*j*2PI*k/n ) (3)
• b = FFT_backward(a) : b[k] = (1/n) sum_j=0~n-1 ( a[j] e^i*j*2PI*k/n ) (4)

(3) 和 (4) 中的 (1/n) 可以放在任何一个式子中，也可以拆成 (1/sqrt(n)) 同时放在两个式子中，目的是保证正向和反向傅立叶变换以后不会相差一个因子。

当 n 的所有素因子均为小整数，尤其是当 n 为 2 的整数次幂的时候，使用适当的算法经过仔细的协调，可以避免多余的计算，使离散傅立叶变换 (3) 和 (4) 减少至 O(n log(n)) 的时间复杂度，即所谓的快速傅立叶变换(FFT)。具体的细节请查阅相关书籍。下面给出我写的一段 FFT 程序，仅供参考。另外也有已经开发的 FFT 函数库，例如 FFTW ，可以直接使用。

• fft.cpp FFT 的 C++ 源程序

利用 FFT，要计算 n₁ 位和 n₂ 位的两个多位数乘法，可以这样进行：开辟两个长度为 n(n>=n₁+n₂，取 2^m 最佳) 的复数数组，将两个多位数从低位到高位分别填入，高位补 0。对两个数组分别进行正向傅立叶变换。将得到的两个变换后的数组的对应项相乘，然后进行反向傅立叶变换，最后得到一个结果数组。由于傅立叶变换是在复数域中进行的，因此还要对结果数组进行取整和进位，才能得到最终的乘积。

值得留意的是傅立叶变换的精度问题。我们知道，在计算机中实数用单精度数或双精度数表示，它们会存在一定的误差。在计算多位数乘法时，n 往往是一个很大的数字，傅立叶变换过程中需要对数组的每一项进行求和，如何保证精度带来的误差不会因为求和而超出允许的范围？我的观点是必须使用双精度实数，而且由于统计特性，精度带来的误差在求和过程中不会很大，一般不会影响计算的正确性。如果需要保证计算的正确性，我想到两种检查方法。第一种是取模验算。例如，如果乘数和被乘数对 17 的模分别是 8 和 6，那么积对 17 的模就应该是 14。第二种是检查运算结果中浮点数偏离整数的最大值。如果偏差只有比如 10^-3 量级，我们可以认为这个尺度的乘法运算很安全；如果偏差达到 0.5，说明运算已经出错了；如果偏差达到 0.1 量级，那也比较危险，也许换个别的乘数和被乘数就溢出了。

多位数的倒数和开方可以通过牛顿迭代求根法转化为乘法运算。例如，要计算 x = 1/a ，根据牛顿迭代法令 f(x) = 1/x - a ，可以得到以下迭代序列：

• x₀ ~= 1/a
  x_k = x_k-1 - f(x_k-1)

  /f'(x_k-1) = 2x_k-1 - ax_k-1² (5)

要计算 x = sqrt(a) ，可以先计算 x = 1 / sqrt(a) ，令 f(x) = 1/x² - a ，可以得到以下迭代序列：

• x₀ ~= 1 / sqrt(a)
  x_k = x_k-1 - f(x_k-1)/f'(x_k-1) = (3/2)x_k-1 - (1/2)ax_k-1³ (6)

(5) 和 (6) 均以 2 的指数收敛到所求结果。还存在其它更复杂一些的迭代序列，它们以更高的指数收敛，在此不提。不过需要提醒的是，跟 (AGM) 不同，这里 (5) 和 (6) 中的 x₀ 只是 1/a 和 1 / sqrt(a) 的约值，在前几次的迭代中不必进行满 n 位数的乘法运算，因而可以减少计算量。

示例程序

作为 AGM 和 FFT 算 PI 的完整过程演示，这里是我新写的算 PI 程序。很遗憾，我的程序比网上可以找到的其它算 PI 程序慢大约 100 倍，而且消耗更多的内存。:-( 目前还不清楚它的瓶颈所在。不管怎么说，它总算是我的第一个 AGM 和 FFT 算 PI 程序，祝贺！:-)

• faint-pi-1.0.3.tar.gz C 源程序，在 Linux 下以 gcc -std=c99 *.c -o pi -lm -O3 编译
• pi.txt.gz 32M 位 PI

综述

以上介绍了三类算 PI 的算法。第一类算法的速度最慢，基本上已经过时。后两类算法的速度目前相差不大，最为常用。迭代法虽然在时间复杂度上有理论上的优势，但实现起来较为复杂，实际上也并不见得比 1/PI 级数法快。

来源：http://elephant.net.cn/articles/pi.php