超越牛顿法？一个利用高阶导数逼近 π 的迭代构想

众所周知，圆周率 π\pi 是数学里最经典的常数之一。虽然我们平时背的是 3.1415926⋯3.1415926\cdots，但真正要高精度计算 π\pi，并不是靠“量圆”或者“背小数”，而是靠各种收敛极快的公式和迭代算法。

这篇文章主要想整理一下几类常见的求 π\pi 方法，最后介绍一个我自己瞎想出来的高阶迭代公式，请大佬轻喷。

#1. 最朴素的想法：从几何出发

最早的求 π\pi 方法，基本都和圆有关。

比如阿基米德的方法，就是用内接正多边形和外切正多边形夹逼圆周长。边数越多，多边形越接近圆，于是就能把 π\pi 夹在一个越来越小的区间里。

这个方法思想很漂亮：

内接多边形周长<圆周长<外切多边形周长\text{内接多边形周长}<\text{圆周长}<\text{外切多边形周长}

但它有个明显缺点：收敛速度不算快。想要很多位小数，需要把多边形边数搞得非常大。

所以它适合理解 π\pi，但不适合现代高精度计算。

#2. 无穷级数法

后来人们发现，π\pi 可以写成很多无穷级数。

最经典的例子是莱布尼茨公式：

π4=1−13+15−17+⋯\frac{\pi}{4}=1-\frac13+\frac15-\frac17+\cdots

也就是：

π=4∑k=0∞(−1)k2k+1\pi=4\sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}

这个公式非常简洁，而且很好证明，但问题也非常严重：收敛太慢。

比如要算很多位 π\pi，这个公式基本不现实。它的美感远大于实用价值。

#3. Machin 类公式

为了提高速度，可以利用反正切函数：

arctan⁡x=x−x33+x55−x77+⋯\arctan x=x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\cdots

当 xx 比较小时，这个级数收敛会快很多。

著名的 Machin 公式是：

π4=4arctan⁡15−arctan⁡1239\frac{\pi}{4}=4\arctan\frac15-\arctan\frac1{239}

因为 15\frac15 和 1239\frac1{239} 都比较小，所以实际收敛速度比莱布尼茨公式强得多。

很长一段时间里，Machin 类公式都是计算 π\pi 的重要工具。

#4. Ramanujan 与 Chudnovsky 公式

再往后，就出现了变态级别的公式。

Ramanujan 给出过很多极快收敛的 π\pi 公式，其中一类大概长这样：

1π=229801∑k=0∞(4k)!(1103+26390k)(k!)43964k\frac1\pi=\frac{2\sqrt2}{9801} \sum_{k=0}^{\infty} \frac{(4k)!(1103+26390k)}{(k!)^4 396^{4k}}

这个公式每一项能带来很多位正确数字。

现代计算 π\pi 最有名的公式之一是 Chudnovsky 公式：

1π=12∑k=0∞(−1)k(6k)!(13591409+545140134k)(3k)!(k!)36403203k+3/2\frac1\pi= 12\sum_{k=0}^{\infty} \frac{(-1)^k(6k)!(13591409+545140134k)} {(3k)!(k!)^3 640320^{3k+3/2}}

它收敛极快，是高精度计算 π\pi 的经典算法之一。

不过这类公式虽然快，但推导背景比较深，涉及模形式、椭圆函数等内容，不是正常人很容易自然想到的东西。

#5. 从“求 π\pi”变成“求方程的根”

其实还有一种思路：

因为

sin⁡π=0\sin \pi=0

所以求 π\pi，可以看成求方程

sin⁡x=0\sin x=0

的正根。

当然，sin⁡x=0\sin x=0 有很多根：

0,π,2π,3π,⋯0,\pi,2\pi,3\pi,\cdots

如果初值选在 π\pi 附近，比如 x0=3x_0=3，那么目标就是收敛到 π\pi。

这时自然可以想到牛顿法。

对

f(x)=sin⁡xf(x)=\sin x

使用牛顿法：

xn+1=xn−f(xn)f′(xn)x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

因为

f′(x)=cos⁡xf'(x)=\cos x

所以：

xn+1=xn−sin⁡xncos⁡xnx_{n+1}=x_n-\frac{\sin x_n}{\cos x_n}

也就是：

xn+1=xn−tan⁡xnx_{n+1}=x_n-\tan x_n

这个迭代在靠近 π\pi 时是二阶收敛的。

所谓二阶收敛，大概意思是：如果当前误差是 ene_n，下一步误差大约会变成 en2e_n^2 量级。也就是说，一旦接近正确答案，正确位数会增长得很快。

但我就想：牛顿法是二阶收敛，那有没有办法构造更高阶的迭代？

答案是：有的兄弟，有的

#6. 一个高阶求 π\pi 的想法

我的想法是，不直接对 sin⁡x\sin x 用牛顿法，而是考虑：

g(x)=1sin⁡xg(x)=\frac1{\sin x}

因为当 x→πx\to \pi 时，sin⁡x→0\sin x\to 0，所以 g(x)g(x) 在 x=πx=\pi 附近会趋向无穷。

换句话说，π\pi 不是 g(x)g(x) 的零点，而是 g(x)g(x) 的一个极点。

对这种函数，可以尝试利用高阶导数的比值构造迭代。

我想到的公式是：

xn+1=xn+(p−1)dp−2dxp−2(1sin⁡x)∣x=xndp−1dxp−1(1sin⁡x)∣x=xnx_{n+1} = x_n+ (p-1) \frac{ \left.\dfrac{d^{p-2}}{dx^{p-2}}\left(\dfrac1{\sin x}\right)\right|_{x=x_n} }{ \left.\dfrac{d^{p-1}}{dx^{p-1}}\left(\dfrac1{\sin x}\right)\right|_{x=x_n} }

其中 p≥2p\ge 2。

如果写得紧凑一点，就是：

xn+1=xn+(p−1)g(p−2)(xn)g(p−1)(xn),g(x)=1sin⁡xx_{n+1} = x_n+ (p-1) \frac{ g^{(p-2)}(x_n) }{ g^{(p-1)}(x_n) }, \quad g(x)=\frac1{\sin x}

并希望有：

lim⁡n→∞xn=π\lim_{n\to\infty}x_n=\pi

#7. 为什么这个公式看起来有道理？

在 x=πx=\pi 附近，令误差：

e=x−πe=x-\pi

因为

sin⁡x=sin⁡(π+e)=−sin⁡e\sin x=\sin(\pi+e)=-\sin e

而当 ee 很小时：

sin⁡e≈e\sin e\approx e

所以：

sin⁡x≈−e\sin x\approx -e

于是：

1sin⁡x≈−1e\frac1{\sin x}\approx -\frac1e

也就是说，g(x)=1sin⁡xg(x)=\frac1{\sin x} 在 π\pi 附近的主要部分很像：

−1x−π-\frac1{x-\pi}

而对

h(x)=1x−πh(x)=\frac1{x-\pi}

求导，有：

h(k)(x)=(−1)kk!(x−π)−k−1h^{(k)}(x)=(-1)^k k!(x-\pi)^{-k-1}

于是高阶导数的比值大概会把 x−πx-\pi 提取出来。

这就说明，用

g(p−2)(x)g(p−1)(x)\frac{g^{(p-2)}(x)}{g^{(p-1)}(x)}

这种结构，确实有机会构造出一个把 xx 往 π\pi 推的高阶迭代。

直观地说，这个公式不是像牛顿法那样看函数的切线，而是利用更多阶导数的信息，试图一次性榨出更多误差项。

#8. 特殊情况：当 p=2p=2 时

当 p=2p=2，公式变成：

xn+1=xn+g(xn)g′(xn)x_{n+1} = x_n+ \frac{ g(x_n) }{ g'(x_n) }

其中：

g(x)=1sin⁡xg(x)=\frac1{\sin x}

计算导数：

g′(x)=−cos⁡xsin⁡2xg'(x)=-\frac{\cos x}{\sin^2 x}

所以：

g(x)g′(x)=1sin⁡x−cos⁡xsin⁡2x−sin⁡xcos⁡x−tan⁡x\frac{g(x)}{g'(x)} = \frac{\frac1{\sin x}}{-\frac{\cos x}{\sin^2 x}} -\frac{\sin x}{\cos x} -\tan x

因此：

xn+1=xn−tan⁡xnx_{n+1}=x_n-\tan x_n

刚好退化成普通牛顿法。

这点我感觉还挺有意思：也就是说，这个公式可以看作牛顿法的一种高阶推广。

#9. 这个公式的优点和问题

这个公式理论上的优点是：

如果 pp 越大，它利用的导数阶数越高，局部收敛阶也可能越高。也就是说，在已经离 π\pi 比较近的时候，它可能比普通牛顿法收敛得更猛。

但是它也有很现实的问题：

第一，高阶导数不好算。

1sin⁡x\frac1{\sin x} 的高阶导数会越来越复杂，实际计算时并不一定划算。

第二，高阶迭代不等于高效算法。

一个算法快不快，不只看迭代次数，还要看每一步的计算量。如果一步里要算特别复杂的高阶导数，那么总耗时可能反而更大。

第三，初值很重要。

这种迭代应该是局部收敛的，也就是初始值要离 π\pi 比较近。比如 x0=3x_0=3 这种就比较自然。如果初值乱选，可能跑到别的根，甚至发散。

所以我不敢说它一定能超过 Chudnovsky 这种顶级公式，但从迭代法角度看，它至少是一个挺有意思的高阶构造。

#10. 总结

求 π\pi 的方法有很多：

从阿基米德的几何夹逼，到莱布尼茨级数，再到 Machin 公式、Ramanujan 公式、Chudnovsky 公式，思路越来越抽象，收敛也越来越快。

而本文最后这个公式，是从另一个角度出发：

把求 π\pi 看成求 sin⁡x=0\sin x=0 的根，再把根转化成 1sin⁡x\frac1{\sin x} 的极点，然后利用高阶导数比值构造迭代。

公式如下：

xn+1=xn+(p−1)dp−2dxp−2(1sin⁡x)∣x=xndp−1dxp−1(1sin⁡x)∣x=xnx_{n+1} = x_n+ (p-1) \frac{ \left.\dfrac{d^{p-2}}{dx^{p-2}}\left(\dfrac1{\sin x}\right)\right|_{x=x_n} }{ \left.\dfrac{d^{p-1}}{dx^{p-1}}\left(\dfrac1{\sin x}\right)\right|_{x=x_n} }

并猜想在合适初值下：

lim⁡n→∞xn=π\lim_{n\to\infty}x_n=\pi

其中 pp 越大，理论上的局部收敛阶越高。