OI

π和e的无穷级数计算

Posted by wyj on June 13, 2023

曾经我一直以为算$\pi$的那些无穷级数是很蠢的,因为就算收敛速度足够快(如线性收敛),每算级数的一项到$n$位就需要至少$O(n)$的时间,还要算$O(n)$项才能得出正确的前$n$位,这样就至少$O(n^2)$了;需要使用高斯-勒让德迭代这样的二阶收敛的方法才能有$\tilde{O}(n)$的复杂度,从而让超高位数的计算变为可能。现在我才知道Binary splitting这一神奇方法的存在,它能用$\tilde{O}(n)$的时间算一个无穷级数的前$n$项和(当然,不是什么级数都可以,但相当广泛的一类级数——超几何级数都能被这样计算)。

e的无穷级数

这里说的就是最最普通的级数$e=\sum_{n=0}^{\infty}\dfrac{1}{n!}$。在数学里我们关心的一般是级数的收敛性,但是计算时需要关心的是有限和与$e$之间的误差。众所周知,

\[\begin{aligned} e-\sum_{k=0}^{n}\dfrac{1}{k!}&=\sum_{k=n+1}^{\infty}\dfrac{1}{k!}\\ &\le \sum_{k=n+1}^{\infty}\dfrac{1}{(n+1)!(n+1)^{k-n-1}}\\ &=\dfrac{1}{n\cdot n!} \end{aligned}\]

并且可以用Stirling公式估算$\log(n!)=O(n\log n)$,因此算$n$项就有$O(n\log n)$位的精度,收敛速度是超线性的,但是并没有比线性快很多,硬算级数是肯定算不了的。

我先是被wiki上的介绍误导了,以为Binary splitting的分治就是计算出$\sum_{k=a}^{b}\dfrac{1}{k!}$的分母和分子(它乱写什么$S(a,b)=\sum_{n=a}^{b}\dfrac{p_n}{q_n}=\dfrac{P(a,b)}{Q(a,b)}$),然后做有理数加法。结果我一想,这不是更慢了吗,我肯定要把每个$S(a,a)$都算出来,那不是等于还要把级数的每一项都算一遍?总不可能真这么简单,否则就不会为这个算法专门起一个名字了。实际上对于$e$来讲,分治要算的是$S(a,b)=\dfrac{P(a,b)}{Q(a,b)}=(a-1)!\sum_{k=a}^{b}\dfrac{1}{k!}$。其中$S(a,a)=\dfrac{1}{a}$,是可以$O(1)$得到的有理数。用$S(a,m)$和$S(m+1,b)$算出$S(a,b)$是非常简单的,因为归纳易证$Q(a,b)=\prod_{k=a}^{b}k$,有

\[(a-1)!\sum_{k=m+1}^{b}\dfrac{1}{k!}=\dfrac{(a-1)!}{m!}\dfrac{P(m+1,b)}{Q(m+1,b)}=\dfrac{P(m+1,b)}{Q(a,m)Q(m+1,b)}\]

然后做有理数加法就行了。

复杂度很好分析,为了算$n$位,需要$m=O(\dfrac{n}{\log n})$项。而分治过程中的数大小不超过$m!$,位数是$O(n)$级别的,因此乘法是$O(n\log n)$,加上分治就是$O(n\log^2 n)$。这个方法算$e$的前556万位(用级数的前$10^6$项)的速度和使用高斯-勒让德迭代算$\pi$的前$2^{20}$位差不多快,都是1.5s左右(同样是我用C++的高精度计算库gmp实现的,这里是代码

π的无穷级数

这里我用的是最最古老的公式之一——梅钦(Machin)公式:

\[\dfrac{\pi}{4}=4\arctan \dfrac{1}{5}-\arctan\dfrac{1}{239}\]

证明:注意到$(239-i)(5+i)^4=114244+114244i$,两边取辐角,证毕。

这个公式的年代非常久远,可表现也并不逊色。我们用泰勒展开计算$\arctan \dfrac{1}{z}$,其中$z\in\mathbb{N},\ z>1$:

\[\arctan \dfrac{1}{z}=\sum_{k=0}^{\infty}\dfrac{(-1)^k}{(2k+1)z^{2k+1}}\]

这是典型的振荡型无穷级数$S_n=\sum_{k\le n}a_k$,且满足$S_{2n}<S_{\infty},S_{2n+1}>S_{\infty}$.于是$\vert S_{n}-S_{\infty}\vert \le \vert a_{n+1}\vert =\dfrac{1}{(2n+3)z^{2n+3}}$,而$z>1$,因此是线性收敛的。(实际上这些推导都没有说另外一面,即这个估计不会比真实情况差太多;但这也是很好说理的,不过有点烦,我懒得了)

可以用类似的操作来分治:$S(a,b)=\dfrac{P(a,b)}{Q(a,b)}=\sum_{k=a}^{b}\dfrac{(-1)^k}{(2k+1)z^{2(k-a)+1}}$,边界情况为$S(a,a)=\dfrac{(-1)^a}{(2a+1)z}$,在$z=O(1)$时是可以$O(1)$算出的。为了能合并,同时还要算出$z^{2(b-a+1)}$。具体的计算过程就不再列出了,这里是代码;我更想强调的是这个过程与进制转换的相似性:这相当于把一个$z$进制的小数转成了10进制,区别不过在于小数的每一位被改成了有理数!而进制转换的分治算法是我很久之前就知道的(可能是还在搞OI的时候?记不清了),与这里的分治其实完全相同,这样来看我这么多年一直不知道无穷级数的快速高精度求法是很奇怪的。

由于是线性收敛的,需要$O(n)$项来算前$n$位。分母是$O(n)!!z^{2O(n)}$级别的,因此位数是$O(n\log n)$;乘法复杂度是$O(n\log^2 n)$,因此总复杂度是$O(n\log^3 n)$。可以看到这个复杂度比高斯-勒让德迭代多了一个$\log$,因此算一百万位的用时是高斯-勒让德迭代的3~4倍。应该和普通的分治FFT一样,可以当区间太小时转为使用暴力来优化常数,但我现在不关心这个了。

梅钦公式是相当古老的,现在有许多无穷级数的表现都比它好得多,但也都只是线性收敛而已,因此只是常数优化,复杂度还是$O(n\log^3 n)$,似乎仍不如高斯-勒让德迭代。然而这类计算无穷级数的算法作为分治算法(且两个子问题互不依赖),相当适合并行化,因此更适用于打破计算$\pi$位数的世界记录。从2009年到现在,世界纪录都是用Chudnovsky algorithm计算的,是一个类似的无穷级数。

关于更通用的Binary splitting方法和其余一些常数的计算,参见多次达成计算$\pi$位数世界记录的程序y-cruncher的介绍:Binary Splitting Recursion Library