优享资讯 | 世界上收敛最快的计算 π 的公式是什么？

知乎热榜 ( ) • 2024-04-03 21:48

世界上收敛最快的计算 π 的公式是什么？

hqztrue的回答

反对高票 @酱紫君的回答。其中对其他答主算法的批判是不对的，充满了ZZ（无意冒犯，只是觉得++YY==ZZ很逗）。关注这位答主挺久了，但这次的前半篇回答私以为有失水准。

为避免误导群众，特开本回答。接下来对酱紫罗列的理由进行逐条反驳。

“但凡自己编程试下，就会发现不可能算到一万位”：

编程试了，一万位很轻松，见后文。

2. “首先你这堆根号套根号, 你怎么控制精度? 记住每开一次根号, 精度减半”：

根据数值分析小常识，开根并不是一个可怕的事情，反倒性质良好。所谓的“精度减半”听起来“减半”有点吓人，其实只是有效数字减少1位，这么想就不慌了。上文公式中只包含两层嵌套根号，精度完全可以接受。

3. “然后就是一堆阶乘, 你自己算算展开 100 项是多少位的数字除多少位的数字”：

想算一万位 $\pi$ 的话，这里只需要求出和式大约前两百项，对应阶乘只需要算到几千位，非常小了。数量级上如果我们想算 $n$ 位 $\pi$ ，那只需要算 $\tilde{O}(n)$ 长度的阶乘。

4. “这堆东西全部乘起来你再加起来还能保留有效精度就有鬼了....”：

之前说了开根完全不是问题。根据误差分析小常识，乘除、正数加法也不是问题。真正比较危险的其实是对正负项的求和，不过公式里的所有数值用有效数字 $\tilde{O}(n)$ 位的浮点数计算就好，没有致命的精度问题，花足够的时间就可以计算到任意精度。

(p.s. 我相信世界上没有鬼)

5. “不止这个答案, 其他所有用级数一项一项加起来, 里面还含有阶乘的答案, 全部都是错的”：

可以说它们收敛速度（指实际计算速度，这里可能有一些定义的分歧）比较慢所以不好，但拿精度作为理由批判它们就不对了。

6. “这里所有出现的数字都是整数, 这个才是可编程计算的无精度损失版本”：

高精度浮点数计算完全不是可怕的事情，合理分析控制精度就不会出问题，不要陷入误区。强求整数运算是不必要的，并不会让式子变得更优美，也不是是否可以编程计算的判断准则。

神奇的是根据对酱紫过往高质量回答的印象，我不认为我这里说的数值分析小常识有哪个是他不知道的，所以只能理解为偶发性的失误？

接下来是写码验证环节。不过要算出正确的 $\pi$ ，首先得把公式给弄对。调查了下发现网上的这个公式十有八九都抄错了，例如原答主 @NGC13009 的式子里写错了好几个地方（多出了两个额外的系数、 $n$ 应该乘在 $B$ 上面而不是 $A$ 上面、 $C$ 的指数处理不对，疑似原因是与另一组参数混淆），这个网页把 $\sqrt{-C^3}$ 的项乘除弄反了，并漏掉了一个 $\sqrt{5}$ 。好不容易才在这里找到了正确的式子，这也太坑了吧：

$\frac{1}{\pi}=\sum_{n=0}^{+\infty} \frac{(6 n) !}{(3 n) !(n !)^3} \cdot \frac{{\color{red}A}+{\color{green}B}n}{{\color{blue}C}^{3n}\cdot \sqrt{{-{\color{blue}C}}^3}},\\$ 其中常数 $A,B,C$ 参见上文原答主的截图。

最后用sagemath可以很简单地写出代码，基本上对着公式照抄就行，还是很方便的，注意把精度设置对就好。分析一下那个公式的计算复杂度：开根、乘除、加减都是近线性复杂度的，需要求前 $O(n)$ 项和，所以总复杂度是 $\tilde{O}(n^2)$ ，不算快但也没有号称的那么慢。

测试下发现仅需0.03s就能算到一万位，花6s可以算到十万位。虽然和更快的算法相比仍然是意料之中的龟速（毕竟是个接近平方复杂度的算法），但至少证明了算到一万位还是很轻松的。

代码如下：

import time

t1=time.time()
L=10000
N=log(10)/log(2)*L+100
R=RealField(N)
sqrt5=R(5).sqrt()

A=(63365028312971999585426220+28337702140800842046825600*sqrt5
    +384*sqrt5*(10891728551171178200467436212395209160385656017
    +4870929086578810225077338534541688721351255040*sqrt5).sqrt())
B=(7849910453496627210289749000+3510586678260932028965606400*sqrt5
    +2515968*R(3110).sqrt()*(6260208323789001636993322654444020882161
    +2799650273060444296577206890718825190235*sqrt5).sqrt())
C=(-214772995063512240-96049403338648032*sqrt5-1296*sqrt5
    *(10985234579463550323713318473+4912746253692362754607395912*sqrt5).sqrt())

inv_pi=0
fac=cur=1
C3=C^3
C0=(-C3).sqrt()

T=L//50+5
for n in range(T):
    v=cur*(A+B*n)/(fac^3*C0)
    inv_pi+=v
    C0*=C3
    fac*=n+1
    cur*=8*(n*6+1)*(n*6+3)*(n*6+5)

pi=1/inv_pi
print('time=',time.time()-t1)
print(str(pi)[:L+1])

p.s. 1.5k赞的原回答公式有误，虽然评论区有人指出但是似乎没人提到正确的式子；1.3k赞的酱紫回答存在严重的知识错误但似乎没太多人发现，看起来大家还是看着一堆数学公式字很多好像感觉很厉害就点个赞了，没有仔细验证。上古知乎遗风去哪了... 我字多你们赞不赞我？

以上。