计算圆周率的Chudnovsky公式

mathematica

\[\frac{640320^{3/2}}{12\pi}=\frac{426880\sqrt{10005}}{\pi}=\sum^\infty_{k=0}\frac{(6k)!(545140134k+13591409)}{(3k)!(k!)^3\left(-640320\right)^{3k}}\]


https://en.wikipedia.org/wiki/Chudnovsky_algorithm
从维基百科偷来的

mathematica

还有用gauss-勒让德    AGM-method计算圆周率的，但是估计消耗的内存非常巨大，以致于用上面的办法了

mathematica

拉马努金圆周率公式：
\[
\frac{1}{\pi}=\frac{2\sqrt{2}}{99^2}\sum^\infty_{k=0}\frac{(4k)!}{(k!)^4}\frac{(26390k+1103)}{396^{4k}}
\]

mathematica

Super PI is a single threaded benchmark that calculates pi to a specific number of digits. It uses the Gauss-Legendre algorithm and is a Windows port of a program used by Yasumasa Kanada in 1995 to compute pi to 2³² digits
http://www.superpi.net/About/

mathematica

PiFast : the fastest windows program to compute pi
http://numbers.computation.free. ... Program/pifast.html

mathematica

In Part 3 we managed to calculate 1,000,000 decimal places of π with Machin's arctan formula. Our stated aim was 100,000,000 places which we are going to achieve now!  Fun with Maths and Python  This is a having fun with maths and python article. See the introduction for important information!  We have still got a long way to go though, and we'll have to improve both our algorithm (formula) for π and our implementation.  The current darling of the π world is the Chudnovsky algorithm which is similar to the arctan formula but it converges much quicker. It is also rather complicated. The formula itself is derived from one by Ramanjuan who's work was extraordinary in the extreme. It isn't trivial to prove, so I won't! Here is Chudnovsky's formula for π as it is usually stated:
https://www.craig-wood.com/nick/articles/pi-chudnovsky/

mathematica

https://arxiv.org/pdf/1809.00533.pdf
A detailed proof of the Chudnovsky formula with means of basic complex analysis

mathematica

http://www.numberworld.org/y-cru ... .html#pi_chudnovsky

似乎在这个地方

nyy

1. (*拉马努金计算圆周率*)
   
2. Clear["Global`*"];(*清除所有变量*)
   
3. aa=Sum[(4k)!*(1103+26390k)/((k!)^4*396^(4k)),{k,0,1000}]//Simplify
   
4. pi=9801/(2*Sqrt[2])/aa
   
5. N[pi-Pi,1000*8]
   

复制代码

计算结果
9.03899782943515430057174757283191068463940040255532974114*10^-7992
因此前面的7991项应该是对的，然后7991/(1000+1)=7.98301698302，也就是平均每项增加7.98个有效数字。
计算到1000，但是由于下标是从0开始的，所以是1001项。

nyy

nyy 发表于 2023-12-21 11:06
计算结果
9.03899782943515430057174757283191068463940040255532974114*10^-7992
因此前面的7991项应该 ...

1. (*Chudnovsky公式计算圆周率*)
   
2. Clear["Global`*"];(*清除所有变量*)
   
3. aa=Sum[(6k)!*(545140134k+13591409)/((3k)!*(k!)^3*(-640320)^(3k)),{k,0,1000}]//Simplify;
   
4. pi=426880*Sqrt[10005]/aa;
   
5. N[pi-Pi,1000*15+10]
   

复制代码

计算结果
-5.2991698741646948426044878234657120328437939087364571033028283320584\
5678250355866207553044563563292959758700266461440817956947591344489822\
7345254380289970993962021582852574098349058844609078300876164962780648\
5297186654543759362327170886404835309478522345959515603456366477085389\
4570169858117902511119500908775907035111977407572937540033696321321140\
4368756280389497603794896888152285180954901797338694204228811273853627\
8755825621023131343593465980435120018279395073716876861696636790478206\
2890643665518965870130201999599313626114234352329597948145084851595722\
0910496520132025334553425938764684093574176330847273635161072260885014\
9656446256673600604299943239862948121533750795201181299611132502681975\
5078209636485580404835050882736729175872875322728445031161384715786553\
9316037193667580285488783578233502014067745866495196510981445704642149\
249991007549275878145253*10^-14197

大约对了14196项，14196/1001=14.1818181818，平均每项增加14.18个有效数字。
这个只是统计意义上的。并不是严格的数学证明