圆周率的近似计算方法综述

序言人们很早就知道圆的周长与直径之比是一个常数，数学家们把这一比率用希腊字母兀来表示，称之为圆周率。圆周率兀是科技领域中最直观和最主要的常数，它是一个极其驰名的数。在日常生活中人们经常与兀接触，并且从有文字记载开始，圆周率就引进了外行人和学者们的兴趣，古今中外许多科学家在兀值计算上献出了自己的智慧和劳动，甚至奉献了自己的一生。因此，准确计算圆周率的值，不仅直接涉及到兀值计算时的需要，而且通过圆周率的数值计算促进了数学的发展。兀值的计算伴随着人类的进步而发展，作为一个非常重要的常数，它最早是解决有关圆的计算问题，所以，求出它的尽量准确的近似值，就是一个极其迫切的问题了。早在二千多年前，古希腊著名数学家阿基米德第一个用科学方法度量圆的周长，得出圆周长与直径之比（圆周率）为3.14；我国杰出数学家刘徽（公元前3世纪）提出震惊中外的“割圆术”求出圆周率的近似值为3.1416；南北朝伟大科学家祖冲之又进一步将圆周率计算在介于3.1415926与3.1615927之间的8位可靠数字。直至1882年德国数学家林德曼证明了兀不仅是一个无理数，而且是一个超越数，给几千年来对兀的认识历史划上了一个句号……在一般工程应用中，对兀值的精度只要求十几位，但是在某些特殊场合需要高精度的圆周率兀值。在信息技术发展迅速的今天，尤其是电脑的发明以来，人们对兀的计算位数大大增加，如今，借助大型计算机对兀有效的计算位数已达小数点后的27000亿位；同时兀的计算也已成为验证超大型计算机计算效率和工作可靠性的一种有效手段。尽管目前数学家已经将兀值计算出小数点后27000亿位，但是，人们对兀的研究还没有完，始终都在追求计算出更为准确的兀值，兀值里仍有许多未解的谜团。现在，圆周率的准确程度在一定程度上反映了一个地区和时代的数学水平，因此，兀的值还要继续计算下去。本文通过利用割圆术、韦达公式、级数加速法、拉马努金公式、迭代法等近似计算方法的介绍和计算实验，来综合表述圆周率兀的计算方法。圆周率的起源及早期发展1.1圆周率简介圆周率是代表圆周长和直径的比例的一个常数（约等于3.1415926）。在日常生活中，通常都用3.14来代表圆周率去进行计算。早期的圆周率没有确定的字母表示，直至1600年，英国威廉•奥托兰特首先使用兀表示圆周率，1737年欧拉在其著作中使用兀。后来被数学家广泛接受，一直沿用至今。圆周率不仅是一个无理数，而且还是一个超越数。早在1767年，兰伯特就证明了兀是一个无理数；1794年，勒让德证明了兀也是无理数；1882年，林德曼证明了兀是超越数。早期是通过实验对兀值进行估算的，第一个用科学方法寻求圆周率数值的人是阿基米德，他在《圆的度量》（公元前3世纪）中用圆内接和外切正多边形的周长确定圆周长的上下界，从正六边形开始，逐次加倍计算到正96边形，得到（3+（10/71））＜兀＜（3+（1/7）），开创了圆周率计算的几何方法（亦称古典方法，或阿基米德方法），得出精确到小数点后两位的兀值。中国数学家刘徽在注释《九章算术》（263年）时只用圆内接正多边形就求得兀的近似值，也得出精确到两位小数的兀值，他的方法被后人称为割圆术。他用割圆术一直算到圆内接正192边形，得出兀R打0（约为3.16），这被称为“徽率”。南北朝时代著名数学家祖冲之进一步得出精确到小数点后7位的兀值（约5世纪下半叶），给出不足近似值3.1415926和过剩近似值3.1415927，还得到两个近似分数值，密率355/113和约率22/7。他的辉煌成就比欧洲至少早了1000年。1.2早期的圆周率数学中的圆，溯源到上古的时候，就引起了人类的探索。《墨经》书中说它是“一中同长也”（“一中”即一个中心或中点。“一中同长”就是到一个心的点的距离都相等，是对圆的定义）。成语说：“不以规矩，不成方圆。”等到人们知道了比例的概念之后，人们自然关顾圆周的长度与圆的直径之间一定的比例-2-常数。尽管圆有大有小，但对一个圆来说，其周长/与直径d之间的比例常数就是圆周率兀。历史上曾采用过圆周率的多种近似值，早期大都是通过实验而得到的结果。古代东方常粗略地用3作为兀的值。我们可以在《旧约•历代志下》第四章（4：2）看到：“他又造一铜海，样式是圆的，径10肘，高5肘，围30肘。”这说明，当时的希伯来人近似以3作为圆周长与直径之比。这相当于拿圆的内接正六边形的周长近似圆的周长。我国第一部《周髀算经》中，就记载有圆〃周三径一〃这一结论。在我国，木工师傅有两句从古流传下来的口诀：叫做：〃周三径一，方五斜七〃，意思是说，直径为1的圆，周长大约是3,边长为5的正方形，对角线之长约为7。这正反映了早期人们对圆周率兀和克这两个无理数的粗略估计。东汉时期官方还明文规定圆周率取3为计算面积的标准。在历史上，兀从粗略的近似3开始，有不少数学家都对圆周率作出过研究。魏晋时，刘徽曾用使正多边形的边数逐渐增加去逼近圆周的方法（即“割圆术”），求得兀的近似值3.1416。汉朝时，张衡得出n的平方除以16等于5/8,即兀等于10的开方（约为3.162）。虽然这个值不太准确，但它简单易理解，所以也在亚洲风行了一阵。王蕃（229-267）发现了另一个圆周率值，这就是3.156,但没有人知道他是如何求出来的。公元5世纪，祖冲之和他的儿子以正24576边形，求出圆周率约为355/113,和真正的值相比，误差小于八亿分之一。这个纪录在一千年后才给打破。约在公元530年，印度数学大师阿耶波多利用384边形的周长，算出圆周率约为3.1415926＜兀＜3.1415927。欧洲斐波那契算出圆周率约为3.1418。圆周率的近似计算历程2.1圆周率的早期计算2.1.1实验时期通过实验对兀值进行估算，这是计算兀的的第一阶段。这种对兀值的估算基本上都是以观察或实验为根据，是基于对一个圆的周长和直径的实际测量而得出的。早期的人们还使用了其它的粗糙方法。如古埃及、古希腊人曾用谷粒摆在圆形上，以数粒数与方形对比的方法取得数值。或用匀重木板锯成圆形和方形以秤量对比取值……由此，得到圆周率的稍好些的值。如古埃及人应用了约四千年的4(8/9)2=3.1605。在印度，公元前六世纪，曾取兀=血=3.162。在我国东、西汉之交，新朝王莽令刘歆制造量的容器一一律嘉量斛。刘歆在制造标准容器的过程中就需要用到圆周率的值。为此，他大约也是通过做实验，得到一些关于圆周率的并不划一的近似值。现在根据铭文推算，其计算值分别取为3.1547,3.1992,3.1498,3.2031比径一周三的古率已有所进步。人类的这种探索的结果，当主要估计圆田面积时，对生产没有太大影响，但以此来制造器皿或其它计算就不合适了。蒲丰在《或然性算术实验》一书中，提出了用实验方法计算兀。这个实验方法的操作很简单：找一根粗细均匀，长度为d的细针，并在一张白纸上画上一组间距为l的平行线(方便起见，常取l二d/2)，然后一次又一次地将小针任意投掷在白纸上。这样反复地投多次，数数针与任意平行线相交的次数，于是就可以得到兀的近似值。因为蒲丰本人证明了针与任意平行线相交的概率为p=2l/兀d。利用这一公式，可以用概率方法得到圆周率的近似值。在一次实验中，他选取l二d/2，然后投针2212次，其中针与平行线相交704次，这样求得圆周率的近似值为2212/704=3.142。当实验中投的次数相当多时，就可以得到兀的更精确的值。1850年，沃尔夫在投掷5000多次后，得到兀的近似值为3.1596。目前宣称用这种方法得到最好结果的是意大利人拉兹瑞尼。在1901年，他重复这项实验，作了3408次投针，求得兀的近似值为3.1415929。2.1.2几何法时期一一割圆法凭直观推测或实物度量，来计算兀值的实验方法所得到的结果是相当粗略的。因此，古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。阿基米德真正使圆周率计算建立在科学的基础上。他是科学地研究这一常数的第一个人，是他首先提出了一种能够借助数学过程而不是通过测量的、能够把兀的值精确到任意精度的方法。由此，开创了圆周率计算的第二阶段。阿基米德求圆周率的更精确近似值的方法，体现在他的一篇论文《圆的度量》之中。在这一书中，阿基米德第一次用上、下界来确定兀的近似值，他用-4-几何方法证明了“圆周长与圆直径之比小于3+(1/7)而大于3+(10/71)”，他还提供了误差的估计。重要的是，这种方法从理论上而言，能够求得圆周率的更准确的值。到公元150年左右，希腊天文学家托勒密得出兀=3.1416，取得了自阿基米德以来的巨大进步。图2.1割圆术在我国，数学家刘徽率在《九章算术》方田章“圆田术”注中，提出割圆术作为计算圆的周长、面积以及圆周率的基础。他从圆内接正六边形出发，将边数逐渐加倍，如图2.1，设圆面积为5o半径为,，圆内接正〃边形边长为匕，周长为Ln,面积为5。将边数加倍后，得到圆内接正2n边形，其边长、周长、面积分别记为七、匕、以。当ln已知，用勾股定理求出七。即如图所示得TOC\o"1-5"\h\z\o"CurrentDocument"22:],222\o"CurrentDocument"l2n=AE7AC+CE2=\：(2仃-(_仃]2(览)求得了内接正n边形的周长Ln，即可求得正2n边形的面积：._1_lr1_S=n(一AB•OE)=n•—二一L•r2n222n(2.2)刘徽割圆术还注意到，如果在内接n边形的每条边上作一高为CE的矩形，就可以证明\o"CurrentDocument"S2n<S0<S2n+('2〃一"(2.3)由此，他从正六边形一值计算到192边形，得出兀^3.14，通常称为“徽率”。南北朝时期的祖冲之计算出了圆周率数值的上下限：3.1415926〈兀<3.1415927，由于史料上没有关于祖冲之推算圆周率“正数”方-5-法的记载，一般认为这个正数的获得是沿用了刘徽的割圆术。如按刘徽割圆术从正六边形出发连续算到正24576边形时，恰好得到这一结果。用Mathmatic计算圆内接6144边形的结果如下：C：iDanmeatsand割圆术冲]=1；a[n_];=Sqrt[2-Sqrt[4-a[n-1]A2]]i碧t[n]:=*2Ao*a[n]•(*多边形周长s[n]:=3t2An*a[n];n=Input[11inputthemdy旨[^x2n)tauribernPrint正七砧七吨破的边七是：七注何]PrMM咽峰正'6再度*蜩七是：，,W[b[n]]]Print["SrtS正七62七哽外嬲是;%N[a[n]]]正乱44迫观鲜是：“-)2十)2十i2十小十］2十］2十(2十J?十由十龙U内隹三614顾邮七是:5.283：513内措正61440fiBKft：3.心9C：iDanmeatsand2.2圆周率的经典计算公式2.2.1基本计算1数值积分法j14/Jax=兀⑴由定积分01+x2计算出该积分的数值，即可得到兀的近似值。x=a+_ai(i=0,1,A,n)⑵将区间［a，b］n等分，则分点in，计算定积分S=jbf(x)dxa。利用定积分的几何意义，可以将小曲边梯形的面积近似地用矩形、梯形来代替，就有了梯形公式、矩形公式:①矩形公式左矩形公式S牝…习f(X)nii=0右矩形公式S牝土£f(x)nii=1中矩形公式S牝…云f(Xi+Xi+1)n2i=0梯形公式SMjp£f(X)+f(X0)+f(Xn)n|_.1，2由Mathmatic编程取n为1000计算(见附录1(A)),由计算可知，中矩形公式取得的结果最接近圆周率值。将小曲边梯形的曲边用三次抛物线来代替，就有了辛普森公式。辛普森(Simpson)公式Sm^―^f(x)+f(x)+2习f(x)+疙f(二)]

6n0n.2Li=1i=0」由Mathmatic编程取n为1000计算(见附录I(B))，由附录I(A)、(B)得，辛普森公式计算更为精确。2蒙特卡罗算法(MonteCarlo)在直角坐标系中，以O(0,0),A(1,0),8(0,1),C(1,1)为定点作一正方形，面积设其中m个点落入扇形区域内。则为V1。以原点为圆心,半径为1在该正方形内作扇形，面积S2=I。在该正方形内随机投入n个点，4m兀M，n由Mathmatic计算分别取n为1000,1001,1002得结果计算如下:

设其中m个点落入扇形区域内。则由Mathmatic计算朋蒙特卡罗算法一成ornirx蒙特卡罗算灌:、DomingtE；aitd夹七七&卒1红豆2如8.般牛-6£1匚81：3其乎，宪11_Inp-Litz[''Plca.seinput:thenvLml>«x:11;"]7HI=O；For[i_=11,i>0zi--rm.=m+If[Random[]A2-i-:Fta_iidom[]*2色1rLO]]FNT4*m./n,30]n.=1000;3.L92c<-jaoaaooooooooooooooooc'ooix=Iiiput-[11P^leaLBeinputtinenu2rJben.:11]m=0;For[1=afi>0z1——,m=m+If[Random[]姑2+Random口八2W1,LO]];M[4*m./nf30]n.==1001;3.27912037912067912387312067911-Inptit[R,Plea.seinputtliem_u-irL>ezt?n:F,];in=0;For[□_=i>05.--fm.=m+Ifz[Random[J2-i-:Rja-iidom[J|2<.1rLO]]/N[4*m./n,30]n.=1002;3.213-5728542914171656686£2£74£5125%-<mi』|.2.2.2级数法】、莱布尼兹级数（1674年发现）:七霜1844年，数学家达什在无计算机的情况下用此公式计算出了兀的前200位小数。误差估计式为：rn数。误差估计式为：rn:1-1+1-1+…+(-1)n-1I3572n-1)1<2n-1根据莱布尼茨级数公式及误差公式计算结果如下图所示：莱布尼茨缱甄nb□叵区|匚：迎口fummtwmdSe11i理旦，红豆£008.研阡翎1匚813C98姑桌面茉布尼我nPrint[MnMf''jt”,■'误差“]Table[{n,N[4Sum[(-1)"j/(2j+1),{j,0,n]].10],N[Pi-4Sun[(-1)Aj/(2j+1),(j.E!!}]]}「{n「100.1000.200}]HTaBleFonnii二误差IDO3.1S1493401-0.009900753003.1445L4904-0.003322255003.143588660-C'.0C'L9960L7003-143C'J918G-0,00L4^6539003.142702531-0.0C'L1098B125%▲<2.欧拉的两个级数（1748年发现）k=1k=0(2k+1)2S3欧拉为甄nb□回1™.2C：\DocumentsandSettingsH2OO8.m-6ElC8L3C98A\^欧拉骚号吝k=0Print["n，1,11JT'\''设差叮V■Table[{nrN[Sqrt[6Sum[1/七{'L,n}]],10],N[Pi-N[Sqrt[6Sum[1/「2,fj,1,□}]]]]},E，100,1000,200}]TalDleForiri=n7T误差IDO3.1320765320.00951612-30C3.1384132450.003173415003.1396B4L230.001908537003.14022314f0.0013G3519003.1405320310.00106062I)2Print[F,n”，"■7T”,”说差■叮VTableF{nrNrSartFSSluhTI/f2i+11A2..0,n}11r10,N[Pi-N[Sqrt[8Sum[1/[2j-1)f{j，OrE]]]]},(n.100.1000.200}]//TableFormnyr误差LDC3.13B4335140.00315314■2002.14O5249S90.001057S950C3.140S5724C0.00063541370C3.141L3dMl0.0J045411290C3.1412333450.000353305L2E%▲<莱布尼兹级数和欧拉的这两个级数的收敛速度较慢。3.基于arctanx的级数由泰勒级数arctanx=工('"*22*1当x—1时有—=工(1"即为莱布尼兹TOC\o"1-5"\h\z2k+142k+1k=0k=0级数。当1x1的值越接近于0,级数收敛的速度越快。令x=tana=5，1仁2tana2x5a=arctan—tan2a===—\o"CurrentDocument"5，1-tan2a1-x212，因此，tan4a=2tan2a1-tan22a121201—"11193=4a-『，就非常接近于012014V—1Vtan4a-11191tanp==一七-=——1+tan4a1+120239H9由此，英国天文学教授JohnMachin得出Machin公式兀=16a一4P=16arctan上一arctan—5239=16^空.上/空.-±-k02k+152k+1上*+12392k+1他利用这个公式计算到了100位的圆周率。由此原理，可以得到高斯公式兀=48arcta^—+32arctan—20arctan—!—TOC\o"1-5"\h\z1857239斯托梅尔公式11,…1兀=24arctan—+8arctan+4arctan—857239类似公式兀—=5cot-1(6)-cot-1兀(99、(452761)、2543/-cot-1(1393)—=6cot-1(8)-cot-1(452761)、2543/-cot-1(1393)兀—=8cot-1(10)-2cot-1兀—=12cot-1(18)+3cot-1(70)+5cot-1(99)+8cot-1(307)4兀■—=5cot-1(7)+2cot-1—=8cot-1(10)-cot-1(239)-4cot-1(515)4-=12cot-1(18)+8cot-1(99)+3cot-1(239)+8cot-1(307)42.2.3迭代法1、1593年，韦达给出2<2v'2+巨\：'2+\：'2+克

——=・・这一不寻常的公式是兀的最早分析表达式。其推导过程是:sint=2cos!sinl=4cos—cosLsin—244=8cos—cosLcos^sin—2488=A对于任意的N，总有sint=2nsinHcos-^-2nn=12n，sint2N.——sint2nn=1lNt—11cos——2n令nT3时，有^^£=rfcos-^-

tn=1取t=—，可得2—=cos—cos8cos16A=H兀cos——2n+1丸cos—=8丸

cos—=16兀cos—=TOC\o"1-5"\h\z42cos—+丸cos—=8丸

cos—=164_22+\：2:cos—+1卜2+e+1■822由归纳法得:—cos—2n+1«2+l''2+A+"'2+“2（n重根号）由公式(1)和(2)可得韦达公式28兀=cos—兀n=1=m2n+1n=1■2+\；‘2+A+(2+巨2<2\;2+t2i/2+%；2+22+、2+\2+^2人甚至在今天，这个公式的优美也会令我们赞叹不已。它表明仅仅借助数字2,通过一系列的加、乘、除和开平方就可算出兀值。2迭代公式⑴1989年，Borwein发现下列收敛于—的迭代公式:兀=4；1-y4n-11-zn1+zn=6-4巨y=、•‘2_10zna0a=(1+y)4a-22n=4；1-y4n-11-zn1+zn=6-4巨迭代误差有估计式1兀———an<164ne-2兀4ny=-—ny=-—nyn1(1+dn仞亏+印2)2a0=1/2a=y2a一5”-1+侦七1(yn12-2yn1+5)⑵1996年，Bailey发现另一个收敛于-1的迭代公式:兀y=5(j5-2),c=(2-—)2,d=(三)0(nyn-1)e=d(7+c)2-3d3+7+c)25迭代误差有估计式兀-一<165ne-标an2.3PC机计算1Ramanujan公式9801兀=2,，2若(4n)!(1103+26390n)44n(n!)4994nn=0该公式是由印度数学家拉马努金(1887—1920)给出的，是计算兀的一个极其有效的公式。1985年，数学家比尔高•斯帕伊利用该公式在计算机上计算出兀的1750万位小数。该级数收敛速度非常快，级数每增加一项，大约可以提高8位小数的精度。2改进的计算公式Chudnovsky112寸(-1)n(6n)!13591409+545140134n—=乙兀6403202n=°(〃!)3(3〃)！6403203n该级数每增加一项，大约可以提高14位小数的精度。1999年9月，日本东京大学教授金田康正和其助手用时37小时21分，计算出了兀的2061.5843亿位小数，检验用时46小时7分钟。3二分法计算算法原理：f(x)二sin(x)在乂二兀点取的零点，又有函数在3与4之间只有一个零点。因此，可对区间［3,4］进行二分，逐渐达到兀。a(1)=3;b(1)=4;k=1;iff(a(k))*f((a(k)/2+b(k))/2)<0a(k+1)=a(k);b(k+1)=(a(k)+b(k))/2;k=k+1;elsea(k+1)=(a(k)+b(k))/2;b(k+1)=b(k);k=k+1;当k=n时结束。x=(a(n)+b(n))/2.算法程序：Functionxx=pi_1(n)a(1)=3;b(1)=4;k=1;fori=1;ifsin(a(k))*sin((a(k)+b(k))/2)<0a(k+1)=a(k);b(k+1)=(a(k)+b(k))/2;k=k+1;elsea(k+1)=(a(k)+b(k))/2;b(k+1)=b(k);k=k+1;endendxx=(a(k)+b(k))/2;执行程序，选择迭代次数为20，可以得到兀=3.14159250259399.计算兀的方法还有很多本文仅列具几种常见并且典型的计算方法。最后不得不提一句的是，为什么对于兀的位数的竞争会持续不断呢？这也应该是相当一部分人的问题吧，可能是因为分数诞生的缘故吧！圆周率计算的最新记录及意义1946年，世界第一台计算机ENIAC制造成功，标志着人类历史迈入了电脑时代。电脑的出现使计算有了根本性的变革。1949年，ENIAC根据梅钦公式计算到2035（一说是2037）位小数，包括准备和整理时间在内仅用了70小时。随着计算机的发展一日千里，圆周率兀值的纪录也就被频频打破。1973年，有人就把圆周率算到了小数点后100万位，并将结果印成一本二百页厚的书，可谓世界上最枯燥无味的书了。1989年兀值突破10亿大关，1995年10月超过64亿位。1999年9月30日，《文摘报》报道，日本东京大学教授金田康正已求到2061.5843亿位的小数值。如果将这些数字打印在A4大小的复印纸上，令每页印2万位数字，那么，这些纸摞起来将高达五六百米。2002年12月6日，金田康正利用一台超级计算机，计算出圆周率小数点后一兆二千四百一十一亿位数，改写了他本人两年前创造的纪录。据悉，金田教授还与日立制作所的员工合作，利用目前计算能力居世界第二十六位的超级计算机，使用新的计算方法，耗时四百多个小时，才计算出新的数位，比他一九九九年九月计算出的小数点后二千六百一十一位提高了六倍。圆周率小数点后第一兆位数是二，第一兆二千四百一十一亿位数为五。如果一秒钟读一位数，大约四万年后才能读完。2009年8月日本筑波大学科学家利用超级计算机耗时29小时计算出圆周率兀的小数点后25770亿位，取得了吉尼斯世界纪录。但该纪录在保持了不到一年的时间里，就被法国巴黎市的软件工程师法布里斯•贝拉德打破。贝拉德使用家用台式电脑，运用比日本科学家有效20倍的改良后的查德诺夫斯基方程算法，耗时131天时间计算出了圆周率小数点后27000亿位这个最新精确值，如果以每秒钟一个数字的速度朗读，这个最新纪录至少要花费49000年时间才能朗读完。[10]圆周率兀对于计算各种数量，例如体积、面积、周长以及任何与圆、圆柱、圆锥、球有重要的作用。但随着数学的不断发展，兀的应用不再局限于求圆的面积和周长，椭圆、萁舌线、旋轮线等面积公式中也都出现了兀值。此外，一-15-些函数的定义，积分的计算，指数的构成等都要用到兀。把兀的值计算到小数点后那么远的位置，并不仅仅是为了创造什么计算的纪录，而是有着更多的理由。一是为了确保兀的“正态”统计信息。如果在一个实数的小数展开中，所有数字以相同频率出现，我们就说它是“简单正态的”；如果相同长度的数字串以相同频率出现，我们就说它是“正态的”。我们还不知道兀（甚至其）是正态的或简单正态的。从1949年ENIAC开始，关于兀的神奇计算就为着解决这个统计特征。从兀的那些展开看，似乎说明它也许是正态的。Shanks在1873年计算的带错的707位又似乎说明兀不是简单正态的。兀是正态的还是非正态的，这个问题当然不可能由计算机来解决。我们面临着一个典型的问题，它需要更高的数学才智，不能单靠计算。这类问题的存在，至少应该能部分医治“计算机万能”病，它在今天似乎颇为流行。不论普通大众还是学数学的同学，似乎越来越感觉任何数学问题都可以通过足够聪明的电子机器来解决。