从一个讲座中看见了一个有趣的问题：祖冲之是怎么计算圆周率的，割圆法真的足够吗？并由此引出了数值计算中常见的几个提高阶数的技巧，值得记录一下。写着写着就刹不住了，将维基百科上的很多相关内容都记录了下来。

\(\pi\)的定义

\(\pi\)作为最重要的数学常数之一，关于它的定义方式数不胜数。

最初等的定义即为周长与直径的比 \[ \pi = \frac{C}{d} \]

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

可以利用积分公式来定义 \[ \pi = \int_{-1}^1 \frac{dx}{\sqrt{1-x^2}} \]

甚至还可以通过偏微分方程的特征值来定义，很多模型问题的特征值都是 \(n \pi\) 的形式。

\(\pi\)还有很多基于连分数的定义/展开式，由于不便于敲公式，这里略去。

\(\pi\)的计算历史

\(\pi\)的计算历史非常悠久，主要分为如下几个阶段：

  • 古代：自然产生的割圆法
  • 近代：
    • 16世纪开始的无穷级数法（包括梅钦类公式），基于一些含\(\pi\)的无穷级数公式
    • 统计模拟算法，即蒙特卡洛模拟，只是一类理论上可以计算的方法，实践中并不适合用来计算
  • 现代（计算机时代）：
    • 迭代算法，早期提出了迭代法，收敛速度比无穷级数快很多，有效位数在每一次迭代之后会增加数倍，但是占用内存较大
    • 快速收敛级数，以拉马努金的那些公式为代表，收敛速度和迭代算法相当，并且内存和复杂度更低，但是背后的代数背景很复杂
    • 阀门算法（1995年提出），与前面的所有算法不同，阀门算法每计算出一位数字，该数字就会像流过阀门的水一样不会再出现在后续的计算过程中，甚至有的算法可以计算十六进制下\(\pi\)的指定位数片段，称为位数萃取算法。

注意，近现代的算法都会很快超过浮点数的有效位数，因此在算法中对于大数的运算和存储需要额外的设计。

割圆法及加速技巧

割圆法的简要历史如下：

  • 第一个有记录的割圆法记录，来自古希腊的阿基米德，他得到的结果为\(3.1408 < \pi < 3.1429\)
  • 中国古代的割圆法源自曹魏时期的刘徽，他使用正3072边形得到了近似值\(\pi \approx 3.1416\)
  • 祖冲之使用割圆法和其它未知技巧，得到的结果为\(3.1415926 < \pi < 3.1415927\)，并提出约率\(\frac{22}7\)和密率\(\frac{142}{45}\)（如果不使用其它技巧，单纯依靠割圆法，需要使用正12288边形，这个量级的精确测量并不太可能）

割圆法的原理很简单，直接构造内接和外接正n边形即可得到 \[ 2 n \sin \frac{\pi}n < 2\pi < 2n \tan \frac{\pi}n \] 这里的正多边形边长数，通常从\(n=6\)开始，逐渐加倍到12, 24, 48, 96，192等。

TODO: 割圆法没有上面这么简答，实际上刘徽和祖冲之都用了一些加速技巧，有空再来整理一下。

定义\(\pi\)的上下界逼近数列 \[ A_n = n \sin \frac{\pi}n , B_n = n \tan \frac{\pi}n \]

基于泰勒展开可以得到 \[ \begin{aligned} A_n &= n \sin \frac{\pi}n = \pi - \frac{\pi^3}{3!}\left(\frac{1}n\right)^2 + \frac{\pi^5}{5!}\left(\frac{1}n\right)^4 - \frac{\pi^7}{7!}\left(\frac{1}n\right)^6 + \cdots \\ B_n &= n \tan \frac{\pi}n = \pi + \frac{\pi^3}{3}\left(\frac{1}n\right)^2 + \frac{2\pi^5}{15}\left(\frac{1}n\right)^4 + \frac{17\pi^7}{315}\left(\frac{1}n\right)^6 + \cdots \end{aligned} \] 因此它们都是\(\pi\)在\(O(n^{-2})\)量级的近似。

下面的计算结果是通过\(A_n\)和\(B_n\)直接获得的上下界，可以看出至少需要正12288边形才能得到3.1415926-3.1415927的估计。 

1
2
3
4
5
6
7
8
9
10
11
12
[    6] 3.000000000000 < pi < 3.464101615138
[   12] 3.105828541230 < pi < 3.215390309173
[   24] 3.132628613281 < pi < 3.159659942098
[   48] 3.139350203047 < pi < 3.146086215131
[   96] 3.141031950891 < pi < 3.142714599645
[  192] 3.141452472285 < pi < 3.141873049980
[  384] 3.141557607912 < pi < 3.141662747057
[  768] 3.141583892148 < pi < 3.141610176605
[ 1536] 3.141590463228 < pi < 3.141597034322
[ 3072] 3.141592105999 < pi < 3.141593748771
[ 6144] 3.141592516692 < pi < 3.141592927385
[12288] 3.141592619365 < pi < 3.141592722039

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

加权平均

由于我们已知上下界数据，基于上界和下界某种加权平均，就可以得到\(\pi\)的\(O(n^{-4})\)量级的近似 \[ C_n = \frac23 A_n + \frac13 B_n = \pi + \frac{\pi^5}{20} \left(\frac{1}n\right)^4 + \frac{\pi^7}{56} \left(\frac{1}n\right)^6 + \cdots \] 这里的加权系数\(\alpha = \frac23\)既可以通过精确公式推出，实践中通常表达式是未知的，我们也可以通过已知计算数据近似获得，例如 \[ \begin{aligned} \alpha =& \alpha_n A_n + (1-\alpha_n) B_n = \alpha_{2n} A_{2n} + (1-\alpha_{2n}) B_{2n} \approx \alpha_{n} A_{2n} + (1-\alpha_{n}) B_{2n} \end{aligned} \] 其中\(\{\alpha_n\}\)是趋于\(\alpha\)的数列，使用上述近似等式可以解得加权系数的近似值 \[ \alpha_n^* = \frac{B_{2n} - B_n}{(B_{2n} - B_n) - (A_{2n} - A_n)} = \frac23 + \frac{\pi^2}8 \left(\frac{1}n\right)^2 + \frac{\pi^4}{128} \left(\frac{1}n\right)^4 + \cdots \] 用近似的加权系数代入计算并不影响阶数 \[ C_n^* = \alpha_n^* A_n + (1-\alpha_n^*) B_n = \pi - \frac{\pi^5}{80} \left(\frac{1}n\right)^4 - \frac{3\pi^7}{1792}\left(\frac{1}n\right)^6 + \cdots \]

下面的计算结果是使用近似加权系数得到的 

1
2
3
4
5
6
7
8
9
10
11
12
[    6] pi ~ 3.138532233545
[   12] pi ~ 3.141406484423
[   24] pi ~ 3.141581097518
[   48] pi ~ 3.141591932576
[   96] pi ~ 3.141592608546
[  192] pi ~ 3.141592650775
[  384] pi ~ 3.141592653414
[  768] pi ~ 3.141592653579
[ 1536] pi ~ 3.141592653589
[ 3072] pi ~ 3.141592653590
[ 6144] pi ~ 3.141592653590
[12288] pi ~ 3.141592653590

预估校正

在展开式中用已知数据近似替代精确值\(\pi\)，可以得到\(\pi\)的\(O(n^{-4})\)量级的近似 \[ \begin{aligned} A_n^c &= A_n + \frac{A_n^3}{3!}\left(\frac{1}n\right)^2 = \pi -\frac{3\pi^5}{40}\left(\frac{1}n\right)^4 +\frac{\pi^7}{56}\left(\frac{1}n\right)^6 +\cdots \\ B_n^c &= B_n - \frac{B_n^3}{3}\left(\frac{1}n\right)^2 = \pi -\frac{\pi^5}{5}\left(\frac{1}n\right)^4 -\frac{4\pi^7}{21}\left(\frac{1}n\right)^6 + \cdots \end{aligned} \]

计算结果如下，每一行依次为两个数列的校正 

1
2
3
4
5
6
7
8
9
10
11
12
[    6] pi ~ 3.125000000000 | 3.079201435678
[   12] pi ~ 3.140503718291 | 3.138438763306
[   24] pi ~ 3.141523757576 | 3.141405133041
[   48] pi ~ 3.141588334395 | 3.141581076807
[   96] pi ~ 3.141592383434 | 3.141591932254
[  192] pi ~ 3.141592636702 | 3.141592608541
[  384] pi ~ 3.141592652534 | 3.141592650775
[  768] pi ~ 3.141592653524 | 3.141592653414
[ 1536] pi ~ 3.141592653586 | 3.141592653579
[ 3072] pi ~ 3.141592653590 | 3.141592653589
[ 6144] pi ~ 3.141592653590 | 3.141592653590
[12288] pi ~ 3.141592653590 | 3.141592653590

外推

外推是一种通用的技巧，利用展开式自身性质来消去自身的低阶项。

下面只使用\(A_n\)做外推，对于\(B_n\)也是类似的。 例如一阶外推可以得到\(\pi\)的\(O(n^{-4})\)量级的近似 \[ E_n = \frac43 A_{2n} - \frac13 A_n = \pi - \frac{\pi^5}{480} \left(\frac{1}n\right)^4 + \frac{\pi^7}{16128} \left(\frac{1}n\right)^6 + \cdots \] 继续二阶外推，可以得到\(\pi\)的\(O(n^{-6})\)量级的近似 \[ \widetilde{E}_n = \frac{16}{15} E_{2n} - \frac{1}{15}E_n = \pi - \frac{\pi^7}{322560} \left(\frac{1}n\right)^6 + \cdots \]

计算结果如下，每一行依次为一阶外推和二阶外推的结果 

1
2
3
4
5
6
7
8
9
10
11
12
[    6] pi ~ 3.141104721640 | 3.141592453898
[   12] pi ~ 3.141561970632 | 3.141592650458
[   24] pi ~ 3.141590732969 | 3.141592653541
[   48] pi ~ 3.141592533505 | 3.141592653589
[   96] pi ~ 3.141592646084 | 3.141592653590
[  192] pi ~ 3.141592653121 | 3.141592653590
[  384] pi ~ 3.141592653560 | 3.141592653590
[  768] pi ~ 3.141592653588 | 3.141592653590
[ 1536] pi ~ 3.141592653590 | 3.141592653590
[ 3072] pi ~ 3.141592653590 | 3.141592653590
[ 6144] pi ~ 3.141592653590 | 3.141592653590
[12288] pi ~ 3.141592653590 | 3.141592653590

无穷级数法

在微积分发明之前，数学家发现了几个关于\(\pi\)的无穷乘积，例如 \[ \frac{2}{\pi} = \frac{\sqrt{2}}2 \cdot \frac{\sqrt{2+\sqrt{2}}}2 \cdot \frac{\sqrt{2+\sqrt{2+\sqrt{2}}}}2 \cdots \] 还有名为沃利斯乘积的表达式 \[ \frac{\pi}2 = \frac21 \cdot \frac23 \cdot \frac43 \cdot \frac45 \cdot \frac65 \cdot \frac67 \cdot \frac87 \cdot \frac89 \cdots \]

格雷果里-莱布尼茨公式

在微积分建立之后，许多关于\(\pi\)的公式随之出现，例如牛顿自己就基于\(\arcsin\)计算了\(\pi\)的15位小数。可以通过对反正切进行泰勒展开 \[ \arctan z = z - \frac{z^3}3 + \frac{z^5}5 - \frac{z^7}7 + \cdots \] 据此可以计算\(\pi\)，例如取\(z=1\)就可以得到 \[ \frac{\pi}4 = 1 - \frac13 + \frac15 - \frac17 + \cdots \] 这个公式被称为格雷果里－莱布尼茨公式，形式非常简单，但是取\(z=1\)的收敛速度非常慢。

还有一些比格雷果里－莱布尼茨公式更快收敛的无穷级数，例如 \[ \pi = 3 + \frac{4}{2 \times 3 \times 4} - \frac{4}{4 \times 5 \times 6} + \frac{4}{6 \times 7 \times 8} + \cdots \] 还有 \[ \frac{\pi}2 = 1 + \frac13 + \frac13 \frac25 + \frac13 \frac25 \frac37 + \cdots \]

梅钦类公式

基于反正切，有一类形式比较独特的梅钦类公式（Machin-like formula），例如最经典的一个公式 \[ \frac{\pi}4 = 4\arctan \frac15 - \arctan \frac{1}{239} \] 注意这个式子是精确成立的，在实际计算时对右侧泰勒展开，由于泰勒展开的位置更靠近0，因此收敛更快。

梅钦类公式的推导过程基于三角函数的和差化积公式 \[ \begin{aligned} \sin(\alpha+\beta) &= \sin \alpha \cos \beta + \cos \alpha \sin \beta\\ \cos(\alpha+\beta) &= \cos \alpha \cos \beta - \sin \alpha \sin \beta\\ \end{aligned} \] 假设满足\(|\alpha + \beta|\le \frac{\pi}2\)，并且引入记号 \[ \begin{aligned} a_1 = \sin(\alpha),\,\,b_1 = \cos(\alpha)\\ a_2 = \sin(\beta),\,\,b_2 = \cos(\beta) \end{aligned} \] 那么就可以得到 \[ \begin{aligned} & \tan(\arctan \frac{a_1}{b_1} + \arctan \frac{a_2}{b_2}) = \tan(\alpha + \beta) = \frac{a_1 b_2 + a_2 b_1}{b_1 b_2 - a_1 a_2} \\ \Rightarrow & \arctan \frac{a_1}{b_1} + \arctan \frac{a_2}{b_2} = \arctan \left( \frac{a_1 b_2 + a_2 b_1}{b_1 b_2 - a_1 a_2} \right) \end{aligned} \] 反复使用这个式子就可以推导得到梅钦类公式，具体过程如下，首先准备 \[ \begin{aligned} &\arctan \frac15 + \arctan \frac15 = \arctan \frac{5+5}{25-1} = \arctan \frac{5}{12} \\ &\arctan \frac{5}{12} + \arctan \frac{5}{12} = \arctan \frac{60+60}{144-25} = \arctan \frac{120}{119} \end{aligned} \] 因此 \[ \begin{aligned} 4 \arctan \frac15 - \frac{\pi}4 =& \arctan \frac{120}{119} + \arctan \frac{-1}{1} = \arctan \frac{120-119}{119+120} = \arctan \frac{1}{239} \end{aligned} \] 最终得到 \[ \frac{\pi}4 = 4\arctan \frac15 - \arctan \frac{1}{239} \]

快速收敛级数

拉马努金（Ramanujan）提出了很多涉及\(\pi\)的神奇的级数公式，用到了模方程等代数背景，这些公式都可以被用来计算\(\pi\)，并且收敛速度远远快于经典的无穷级数。

拉马努金至少一次性提出了14个相同类型的公式，最著名的一个公式如下 \[ \frac{1}{\pi} = \frac{2\sqrt 2}{9801} \sum_{k=0}^\infty \frac{(4k)! (1103+26390k)}{(k!)^4 396^{4k}} \]

下面的楚德诺夫斯基（Chudnovsky）公式相当于拉马努金公式的改进，每计算一项就能得到\(\pi\)的约14位有效数字，因而被用于突破圆周率数位的计算。 \[ \frac{1}{\pi} = 12 \sum_{k=0}^\infty \frac{(-1)^k (6k)! (13591409+545140134k)}{(3k)!(k!)^3 640320^{3k+3/2}} \]

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

迭代算法

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

  1. 初始化：\(a_0 = 1\)，\(b_0 = \frac{1}{\sqrt{2}}\)，\(t_0 = \frac14\)，\(p_0=1\)
  2. 迭代计算：\(a_{n+1} = \frac12 (a_n + b_n)\)，\(b_{n+1} = \sqrt{a_n b_n}\)，\(t_{n+1} = t_n - p_n(a_n-a_{n+1})^2\)，\(p_{n+1} = 2p_n\)
  3. \(\pi\)的近似值为：\(\pi \approx (a_n + b_n)^2/(4t_n)\)

注：

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

阀门算法

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

贝利－波尔温－普劳夫公式（BBP公式）是其中的典型 \[ \pi = \sum_{k=0}^\infty \frac{1}{16^k} \left( \frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6} \right) \] 基于这个公式可以容易地设计出一种萃取算法来计算指定片段，因此这个公式主要被用来检查其它算法得到的新纪录的结果是否正确，将其它算法的结果转换到16进制并与之比较。

虽然算法的思路非常新颖，但是公式本身的证明并不困难：首先，对于\(n < 8\)有下式成立 \[ \int_0^{1/\sqrt{2}} \frac{x^{n-1}}{1-x^8}\,dx = \int_0^{1/\sqrt{2}} \sum_{k=0}^\infty x^{n-1+8k}\,dx = \frac{1}{2^{n/2}} \sum_{k=0}^\infty \frac{1}{16^k} \frac{1}{8k+n} \] 因此BBP公式右边变为 \[ \sum_{k=0}^\infty \frac{1}{16^k} \left( \frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6} \right) = \int_0^{1/\sqrt{2}} \frac{4 \sqrt{2} - 8 x^3 - 4\sqrt{2} x^4 - 8 x^5}{1-x^8}\,dx = \pi \] 得证。