圆周率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秒的时间。

 计算历史

日期	计算者	国籍	正确位数	详细纪录
前20世纪	未知	古巴比伦王国	1	π=3.125
前20世纪	未知	古印度	1	π=3.160493...
前12世纪	未知	中国	-	π=3
前6世纪中	圣经列王记上7章23节		-	π=3
前3世纪	阿基米德	古希腊	3	π=3.1418
公元前20年	维特鲁威	古罗马	1	π=3.125
公元前50年－公元前23年	刘歆	中国	1	π=3.1547
130年	张衡	中国	1	π=3.162277…
150年	未知	托勒密	3	π=3.141666…
250年	王蕃	中国	1	π=3.155555…
263年	刘徽	中国	5	π=3.14159
480年	祖冲之	中国	7	π=3.1415926
499年	阿耶波多	印度	3	π=3.1416
598年	婆罗摩笈多	印度	1	π=3.162277…
800年	花拉子米	乌兹别克	3	π=3.1416
12世纪	婆什迦罗第二	印度	4	π=3.14156
1220年	斐波那契	意大利	3	π=3.141818
1400年	Madhava		10	π=3.14159265359
1424年	Jamshid Masud Al Kashi		16
1573年	Valentinus Otho		6
1593年	弗朗索瓦·韦达	法国 [10]	9
1593年	Adriaan van Roomen		15
1596年	鲁道夫·范·科伊伦		20
1615年			32
1621年	威理博·司乃耳，范·科伊伦的学生		35
1665年	牛顿		16
1699年	Abraham Sharp		71
1700年	关孝和		10
1706年	John Machin		100
1706年	William Jones			引入希腊字母π
1719年	De Lagny		112	得出127位 前112位正确
1723年	建部贤弘		41
1730年	Kamata		25
1734年	莱昂哈德·欧拉			引入希腊字母π并肯定其普及性
1739年	松永良弼		50
1761年	约翰·海因里希·兰伯特			证明π是无理数
1775年	欧拉			指出π可能是超越数
1794年	Jurij Vega		136	得出140位小数 前136位正确
1794年	阿德里安-马里·勒让德			-
1841年	Rutherford		152	得出208位小数 前152位正确
1844年	Zacharias Dase及Strassnitzky		200
1847年	Thomas Clausen		248
1853年	Lehmann		261
1853年	William Rutherford		440
1855年	Richter		500
1874年	William Shanks		527	得出707位小数 前527位正确
1882年	Lindemann			证明π是超越数
1946年	D. F. Ferguson		620
1947年			710
1947年			808
1949年	J. W. Wrench爵士和L. R. Smith		2,037	首次使用计算机
1955年	J. W. Wrench爵士及L. R. Smith		3,089
1957年	G.E.Felton		7,480
1958年	Francois Genuys		10,000
1958年	G.E.Felton		10,020
1959年	Francois Genuys		16,167
1961年	IBM 7090 晶体管计算机		20,000
1961年	J. W. Wrench, Jr，及L. R. Smith		100,000
1966年			250,000
1967年			500,000
1974年			1,000,000
1981年	金田康正		2,000,000
1982年			4,000,000
1983年			8,000,000
1983年			16,000,000
1985年	Bill Gosper		17,000,000
1986年	David H. Bailey		29,000,000
1986年	金田康正		33,000,000
1986年			67,000,000
1987年			134,000,000
1988年			201,000,000
1989年	楚诺维斯基兄弟		480,000,000
1989年			535,000,000
1989年	金田康正		536,000,000
1989年	楚诺维斯基兄弟		1,011,000,000
1989年	金田康正		1,073,000,000
1992年			2,180,000,000
1994年	楚诺维斯基兄弟		4,044,000,000
1995年	金田康正和高桥大介		4,294,960,000
1995年			6,000,000,000
1996年	楚诺维斯基兄弟		8,000,000,000
1997年	金田康正和高桥大介		51,500,000,000
1999年			68,700,000,000
1999年			206,000,000,000
2002年	金田康正的队伍		1,241,100,000,000
2009年	高桥大介		2,576,980,370,000
2009年	法布里斯·贝拉	法国	2,699,999,990,000
2010年	近藤茂		5,000,000,000,000	[11]
2011年	IBM“蓝色基因” 超级电脑			π2的前60,000,000,000,000位二进制小数 [12]