Shun Li's Notes

圆周率的数值计算

2.2k words
数值分析

圆周率的定义

作为最重要的数学常数之一，关于它的定义方式数不胜数，下面简要列举几种定义方式：

  1. 最初等的定义即为周长与直径的比

  2. 可以利用三角函数定义： 是正弦函数的最小正零点。 注意到三角函数的定义并不一定需要依赖于几何，它们可以独立依赖于微分方程的解来定义，从而避免循环定义。

  3. 可以利用积分公式来定义

  4. 还可以通过偏微分方程的特征值来定义，很多模型问题的特征值都是 的形式，这和第二种定义本质上是一样的。

  5. 还有很多基于连分数的定义/展开式，这里略去。

圆周率的计算历史

的计算历史非常悠久，大致可以分为如下几个阶段：

  1. 古代：自然产生的割圆法

  2. 近代：

    1. 16世纪开始的无穷级数法（包括梅钦类公式），基于一些含 的无穷级数公式

    2. 统计模拟算法，即蒙特卡洛模拟，只是一类理论上可以计算的方法，实践中并不适合用来计算

  3. 现代（计算机时代）：

    1. 迭代算法，早期提出了迭代法，收敛速度比无穷级数快很多，有效位数在每一次迭代之后会增加数倍，但是占用内存较大

    2. 快速收敛级数，以拉马努金的那些公式为代表，收敛速度和迭代算法相当，并且内存和复杂度更低，但是背后的代数背景很复杂

    3. 阀门算法（1995年提出），与前面的所有算法不同，阀门算法每计算出一位数字，该数字就会像流过阀门的水一样不会再出现在后续的计算过程中，甚至有的算法可以计算十六进制下 的指定位数片段，被称为位数萃取算法。

Remark

近现代的数值算法获得的有效精度会很快超过双精度/四精度浮点数的有效位数，因此在算法实现时对于 的存储和运算都需要进行额外的设计。

割圆法及加速技巧

割圆法的简要历史如下：

  1. 第一个有记录的割圆法记录，来自古希腊的阿基米德，他得到的结果为

  2. 中国古代的割圆法源自曹魏时期的刘徽，他使用正3072边形得到了近似值

  3. 祖冲之使用割圆法和其它未知的加速技巧，得到的结果为 ，并提出约率 和密率 （如果不使用其它技巧，单纯依靠割圆法，需要使用正12288边形，这个量级的精确测量并不太可能）

割圆法的原理很简单，直接构造内接和外接正n边形即可得到

这里的正多边形边长数，通常从 开始，逐渐加倍到12, 24, 48, 96，192等。

Remark

割圆法实际没有上面这么简单，根据文献分析，刘徽和祖冲之在实践中应该都使用了一些加速技巧。

定义 的上下界逼近数列

基于泰勒展开可以得到

因此它们都是 在 量级的近似。 通过 和 直接获得的上下界，至少需要正12288边形才能得到 的估计。

基于这两个直接测量可得的数列，有几种常见的提高阶数的技巧可以用来提升 的计算精度。这几种技巧可以用于加速割圆法的收敛，假设祖冲之使用了类似的某种提高收敛速度的技巧，那么获取同样的精度就不再需要正12288边形的计算，可以大约减少一个数量级。

  1. 加权平均：由于我们已知上下界数据，基于上界和下界某种加权平均，就可以得到 的 量级的近似

    这里的加权系数 可以通过精确公式推出，但是在实践中表达式可能是未知的，我们也可以通过已知计算数据近似获得，例如

    其中 是趋于 的数列，使用上述近似等式可以解得加权系数的近似值

    用近似的加权系数代入计算并不影响阶数

  2. 预估校正：在展开式中用已知数据近似替代精确值 ，可以得到 的 量级的近似

  3. 外推：外推是一种通用的加速技巧，利用展开式自身性质来消去自身的低阶项。 下面只使用 做外推，对于 也是类似的。 例如一阶外推可以得到 的 量级的近似

    继续二阶外推，可以得到 的 量级的近似

无穷级数法

在微积分之前的时代，数学家就已经发现了几个关于 的无穷乘积，例如

还有名为沃利斯乘积的表达式

在微积分建立之后，更多关于 的公式随之出现，牛顿自己就利用 计算了 的15位小数。例如对反正切进行泰勒展开

取 就可以得到

这个公式被称为格雷果里－莱布尼茨公式，形式非常简单，但是取 对应的收敛速度比较慢。

还有一些比格雷果里－莱布尼茨公式更快收敛的无穷级数，例如

有一类基于反正切但是形式比较独特的梅钦类公式（Machin-like formula），例如最经典的一个公式为

注意这个式子是精确成立的，在实际计算时通过对右侧的 进行泰勒展开计算，由于泰勒展开的位置更靠近 ，收敛速度更快。 梅钦类公式的推导过程基于三角函数的和差化积公式

假设满足 ，并且引入记号

那么就可以得到

反复使用这个式子就可以推导得到梅钦类公式。具体过程如下，首先准备

因此

最终得到

拉马努金（Ramanujan）提出了很多涉及 的神奇的级数公式，用到了模方程等代数背景，这些公式都可以被用来计算 ，并且收敛速度远远快于经典的无穷级数。 拉马努金至少一次性提出了14个相同类型的公式，最著名的一个公式如下

下面的楚德诺夫斯基（Chudnovsky）公式相当于对拉马努金公式的改进，每计算一项就能得到 的约14位有效数字，因而被用于突破圆周率数位的计算。

Remark

这两个公式的证明显然超出了本文的讨论范围。

迭代算法

有很多形式类似的迭代算法，这里记录一个典型的高斯-勒让德算法：

  1. 初始化

  2. 迭代计算

  3. 的近似值可以取为

Remark

高斯-勒让德算法并没有明确的迭代终止判定，选取一个足够大的迭代次数即可。高斯-勒让德算法的证明需要用到含有椭圆积分的恒等式，超出了本文的讨论范围。

关于这套算法，有一个更一般性的例子。

Example

给定两个正数 ，构造数列 和 ： ， ，递推关系为：

那么，这两个数列会收敛到同一个极限，称为 和 的算术几何平均值，记作 。

Proof

由归纳法容易证明： 递减， 递增，并且满足

又因为

因此，两个数列都收敛，并且收敛到同一个极限。

Note

可以用椭圆积分表示 的解析表达式

阀门算法

阀门算法是一类非常新颖的算法，在提出时曾震惊学界，因为它不再需要保存整个 ，而是像流水一样计算 的每一位数，并且后续不再需要这些数据，甚至可以直接计算指定片段而不需要先前的近似值，（位数萃取）只不过需要在十六进制意义下。

贝利—波尔温—普劳夫公式（BBP公式） 是阀门算法的典型

基于这个公式可以容易地设计出一种萃取算法来计算十六进制下的指定片段，因此这个公式主要被用来检查其它算法得到的新纪录的结果是否正确，将其它算法的结果转换到16进制并与之比较。

虽然阀门算法的思路非常新颖，但是 BBP 公式本身的证明并不困难：首先，对于 有下式成立

取几个不同的 就可以组合出 BBP 公式的右侧，即

BBP 公式得证。

数值分析