圆周率PI 几千年的进展

原来一直没弄明白的圆周率，突然想给小孩子讲圆的周长怎么计算的，结果PI都不会，计算啥！原来PI计算了几千年 才有的结果。

1、蒙特卡洛法：非常直观易懂，但效率极低

蒙特卡洛是欧洲著名赌城，编程中采用的蒙特卡洛法就是指用随机数的方法。

设想一个正方形，边长a=2，则其内接圆的半径r=1，则圆和正方形的面积比是：

Sy/Sz=πr2/a2=π/4

因此,

π=4Sy/Sz

我们用蒙特卡洛法来获得面积比Sy/Sz：让计算机产生很多的随机点(x,y)，其中x、y都是[0,1]之间的随机数，可以用判断 (x2+y2) 是否小于1的方法来确定是否在圆内，则圆内的点数代表圆面积，总点数代表总面积，两者相除就是面积比Sy/Sz。

2、割圆迭代法：也很直观易懂，效率还不错

现在我们用刘徽和祖冲之的割圆术试试。

原理是这样的：设想一个圆，半径r=1，内接一个正n边形，边长为a，那么当这个多边形的边数n越大，它的周长就越接近圆的周长。所以我们就可以用多边形的周长代替圆的周长并应用圆的周长公式得到：

π=na/2

圆的内接正n边形边长a可以迭代计算：

a[n]=sqrt(2−2sqrt(1−(a[n/2]/2)2))

我们从正六边形开始，它的边长a=r=1，然后计算正12边形、正24边形、正48边形...逐渐逼近圆。

内接正n边形           π计算结果               
6                   3.000000000000      
12                  3.105828541230      
24                  3.132628613281      
48                  3.139350203047      
96                  3.141031950891      
192                 3.141452472285      
384                 3.141557607912      
768                 3.141583892149      
1536                3.141590463237      
3072                3.141592106043      
6144                3.141592516588      
12288               3.141592618641      
24576               3.141592645321      
49152               3.141592645321      
98304               3.141592645321  

3、梅钦级数法：不直观但也易懂，效率略高

大Boss梅钦（John Machin）用他的神奇公式亲手把π算到了小数点后100位，今天，我们借助Python编程，可以远远超过这个精度。

梅钦公式：

π=16arctan(1/5)−4arctan(1/239)

mathematica

Clear["Global`*"]
expi0 = N[Pi, 50];
ex1 = 16 ArcTan [1/5] - 4 ArcTan [1/239]
N[ex1, 100]
arctan[x_] := Normal[Series[ArcTan [x], {x, 0, 30}]]
ex2 = ArcTan[1/5];
Print["原函数：", N[ex2, 30]]
ex3 = arctan[a]
a = 1/5;
Print["级数化简后：", N[ex3, 30]]
ex5 = ex2 - ex3;
N[ex5, 30]
exout = 16 arctan [x] - 4 arctan [y]
x = 1/5;
y = 1/239;
expi0
expi1 = N[exout, 50]
expi0 - expi1

4. 拉马努金 也算是神公式  计算也比较快

计算12次，只能到8位的精度 

5. Chudnovsky在拉马努金公式 下再优化  

 只用计算6次就能到99位精度

ex1 = Sum[((-1)^
      k*(6 k)!*(13591409 + 
       545140134 k))/((3 k)!*(k!)^3*640320^(3 k + 3/2)), {k, 0, 
   n}];   (* Chudnovsky *)
chud[n_] := 
 N[1/(12*Sum[((-1)^
          k*(6 k)!*(13591409 + 
           545140134 k))/((3 k)!*(k!)^3*640320^(3 k + 3/2)), {k, 0, 
       n}]), 100];
a2 = chud[6]

6.目前最优化方法  高斯 - 勒让德算法Gauss\[Dash]Legendre algorithm

 

 只用计算6次 就能精确到小数后171位

Clear["Global`*"]
precision = 10000;
Print["==>正确的圆周率："]
pi1 = N[Pi, precision]
a0 = 1;    (*不能指定小数位 不然后面不能指定精度*)  
b0 = 1/Sqrt[2];
t0 = 1/4;
p0 = 1;  (* 高斯-勒让德算法Gauss\[Dash]Legendre algorithm *)
 
Do[
 a1 = (a0 + b0)/2;
 b1 = Sqrt[a0 *b0];
 t1 = t0 - p0 (a0 - a1)^2;
 p1 = 2 p0;
 a0 = a1;
 b0 = b1;
 t0 = t1;
 p0 = p1;
 rst = (a1 + b1)^2/(4 t1);
 Print["==>", n, "  ", N[rst - pi1, 30]], {n, 1, 11, 1}]
Print["==>计算的结果：(精度)", precision, "位  ", N[rst, precision]]

==>6  -2.3085807149343902668213207344*10^-171

计算11次 结果到了5583位精度了。 普通电脑2秒的时间。

 计算历史