伪·从零开始学算法 - 2.3 求圆周率

本文为圆周率日特辑~

简介

圆周率π是一个圆的周长与直径之比。即:π = C / d 。它是一个无限不循环小数,即无理数。

据说远古时期的人们根据经验以3作为圆周率的值(“径一周三”),这个值至今也能够用于一些对精度要求不高的场合(如手工做木桶)。但是,对于天体运行和面积测量等方面,3这个值未免太粗略。人们对于它的计算绝不仅限于此。

古埃及、古巴比伦、古印度、古代中国的遗迹中都能够找到关于圆周率的较精确的值,它们都能够达到3.1这个值。实际上,绝大多数情况下,那时的人还是把它当作有理数来计算的。总之,那时的计算可能仍然属于观测值加以拟合。

割圆术

割圆术是第一个有纪录、严谨计算π数值的算法。

一般来说,它是通过计算圆内接正多边形和/或圆外切正多边形的周长的方式来计算圆周率。随着正多边形的边数增长,其形状越接近圆,因此可以据它的周长或面积计算圆周率。理论上来说,它可以计算任意精度的圆周率的值。

目前所知最早这么做的是约公元前250年的阿基米德。阿基米德的算法是在计算圆的外切正六边形及内接正六边形的边长,以此计算π的上限及下限,之后再将六边形变成十二边形,继续计算边长……,一直计算到正96边形为止。由此,他计算出π的值在223/71和22/7(3.140845...和3.142857...)之间。

此后的托勒密等人也使用割圆术来计算圆周率,但他们是单向迫近,不如阿基米德的双向迫近精确。

三国时期的中国,刘徽也创立了割圆术。与之前不同的是,他的理论是数学史上最严谨完备简洁的割圆术:使用极限论来论证割圆术的正确性;双向迫近。

割之弥细,所失弥少。割之又割,以至于不可割,则与圆周合体而无所失矣。

与阿基米德相似,刘徽也是从正六边形开始计算的。不同的是,刘徽以面积来计算。他自己通过分割圆为192边形,计算出圆周率在3.141024与3.142704之间,取其近似值157/50。

刘徽割圆术原理

在下图中可以发现,如果计原多边形(黄色部分)面积为S1,新的多边形(黄色+绿色部分)面积为S2,圆面积为S,那么新的多边形与原多边形的面积之差(绿色部分)为S2 - S1,矩形ABCD的面积和为2(S2 - S1)。则有:

S2 < S < S2 + (S2 - S1)

S2 < S < S2 + 2(S2 - S1)

刘徽圆周率不等式示意图

后来刘徽发明一种快捷算法,通过新的多边形与原多边形的面积之差近似成等比级数来计算,可以只用96边形得到和1536边形同等的精确度,从而得令他自己满意的3.1416。

南北朝的祖冲之,运用刘徽的算法,继续分割到12288边形,求得24576边形的面积,得3.14159261864 < π < 3.141592706934,精确到了小数点后七位,保持了世界最准确圆周率达900年之久。此后他又算出来两个比较简单的近似值:约率(22/7 > π)和密率(355/113 > π)。

至于割圆术的算法,简要解释如下:

割圆术

设上图中圆半径为1,弦心距为hn;正n边形的边长为xn,面积为Sn。由勾股定理,得

hn = sqrt[1 - (xn / 2)2]

x2n = sqrt[(xn)2 + (1 = hn)2]

(n ≥ 6)

纯文字形式有点乱,可以看这个

得x6 = 1。

而且我们可以得到,正2n边形面积等于正n边形面积加n个等腰三角形的面积,即

S2n = Sn + nxn(1 - hn) / 2 (n ≥ 6)

随着n的增大,S2n不断趋近于π。

于是我们可以据此编制割圆术的算法:

割圆术(1)

据此求在正n边形下,圆周率的下限和上限:

割圆术(2)

割圆术在16~17世纪之前一直都是精确计算圆周率的解决方法。奥地利天文学家克里斯托夫·格林伯格在1630年用1040边形计算到第38位小数,至今这仍是利用多边形算法可以达到最准确的结果(当然,是人工的情况下,我在测试算法的时候,仅用53.6 ms就能计算到小数点后100位)。

无穷级数

16世纪及17世纪时,π的计算开始改用无穷级数的计算方式。但是,公式并不是唯一的。1735年,欧拉解决了巴塞尔问题,得到了关于π的一个非常精妙的公式(可能很多人都知道):

算法如下:

欧拉的公式推导的算法

i越大,结果越精确。

此外,还有莱布尼茨公式、尼拉卡莎级数等一系列级数。但是许多级数的收敛速度很慢,莱布尼茨公式要到500000项之后,才会精确到π的第五位小数。


莱布尼茨公式
尼拉卡莎级数

印度数学家斯里尼瓦瑟·拉马努金在1914年发表了许多与π相关的公式,这些公式十分新颖,极为优雅而又颇具数学深度,收敛速度也非常快。
以下为一例:

第一位使用拉马努金公式计算π并取得进展的是比尔·高斯珀,他在1985年算得了小数点后一千七百万位。

拉马努金公式开创了现代数值近似算法的先河。,此后波尔文兄弟和楚德诺夫斯基兄弟进一步发展了这类算法。后者于1987年提出了楚德诺夫斯基公式,如下所示:

此公式每计算一项就能得到π的约14位数值,因而被用于突破圆周率的数位的计算。利用该公式,有人已经计算到小数点后第一万亿位。

但是很显然,这些公式已经变得非常复杂,难以记忆(虽然没有多少人需要记这个的)。

迭代算法

二十世纪中期计算机技术的发展、革新再次引发了计算π位数的热潮,除了之前的无穷级数,利用迭代算法计算圆周率的方法也被提出。

迭代算法因为收敛速度比无穷级数快很多,在1980年代以后广为使用。无穷级数随着项次的增加,一般来说正确的位数也会增加几位,但迭代算法每多一次计算,正确的位数会呈几何级数增长。

有一个比较有名的算法是高斯-勒让德算法,于1975年被理查德·布伦特和尤金·萨拉明独立发现。日本筑波大学于2009年8月17日宣布利用此算法计算出π小数点后2576980370000(2.58万亿)位数字。知名的电脑性能测试程序Super PI也使用此算法。

该算法的步骤如下:

第一步:设置初始值:

第二步:反复执行以下步骤直到an与bn之间的误差到达所需精度:

则π的近似值为:

流程图如下:

高斯-勒让德算法

它以迅速收敛著称,算法每执行一步正确位数就会加倍,只需25次迭代即可产生π的4500万位正确数字。不过,它的缺点是内存密集。

蒙特卡罗方法

前面介绍的是靠计算得出的方法,但这个却与概率与统计有关。

蒙特卡罗方法(Monte Carlo method)是以概率统计理论为指导的一类非常重要的数值计算方法,通过进行大量重复试验计算事件发生的频率,按照大数定律,可以求得π的近似值。

布丰投针问题就是其中一个应用的例子。当一枚长度为l的针随机地往一个画满间距为t(l ≤ t)的平行线的平面上抛掷n次, 如果针与平行直线相交了m次,那么当n充分大时就可根据以下公式算出π的近似值:π≈(2nl)/(mt)。

布丰投针问题

另一个例子是随机地往内切四分之一圆的正方形内抛掷大量的点,落在四分之一圆内的点的数量与抛掷点的总量的比值会近似等于π/4。

蒙特卡罗方法

我们以第二个例子为例,如果要编制算法,我们可以使用解析几何的方法,把“随机地往正方形内抛掷大量的点”转化为“定义一个点,该点的横坐标和纵坐标取[0,1]范围的随机数(两坐标不一定相同),重复若干次”,把“落在四分之一圆内”这个条件转化为“点是否在x2+y2=1内”。于是得下面的流程图:

蒙特卡罗方法流程图

它的缺点也很明显:随机、不精确。

结语

还有很多算法,在此我就不一一列举了。

最后我放上一张图,以展示π的计算的发展。

当数学家发现新的算法、电脑变得普及时,π的已知小数位急剧增加。注意垂直坐标使用了对数坐标

一般而言,π值并不需要过于精确便能够满足大部分的数学运算的需求。按照约尔格·阿恩特(Jörg Arndt)及克里斯托夫·黑内尔(Christoph Haenel)的计算,39个数位已足够运算绝大多数的宇宙学的计算需求,因为这个精确度已能够将可观测宇宙圆周的精确度准确至一个原子大小。但是大家仍然对π的计算乐此不疲,这有很多原因,但计算π的本身让我想起了一个故事:

有人问英国登山家乔治·马洛里为何想要攀登珠穆朗玛峰。他回答说:“因为山就在那儿。”

参考资料

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,324评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,303评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 150,192评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,555评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,569评论 5 365
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,566评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,927评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,583评论 0 257
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,827评论 1 297
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,590评论 2 320
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,669评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,365评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,941评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,928评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,159评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,880评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,399评论 2 342

推荐阅读更多精彩内容

  • 我在写1.1节的时候本来是要写这个的,但是突然就忘了……就作为一节来写吧。 顺便说一下,1946年的今天,世界上第...
    阿啊阿吖丁阅读 3,120评论 1 0
  • 勾股定理 圓 圓形的概念的形成,是人類認知歷史上的一大里程碑。 圓周率 定义1 定义2 以圆形半径为边长作一正方形...
    光剑书架上的书阅读 1,616评论 0 4
  • 14年的夏天,我去了一趟火车站,记忆里儿时的旅行都是从这里出发,仿佛只有到了这里,才能与风为伴,抵达远方。呆滞了一...
    Daisy沐阅读 123评论 0 0
  • 不知从什么时候起, 不再虚荣自负, 不再好大喜功, 喜欢微笑着去面对所有的得失。 曾经以为非我莫属的东西,原来并不...
    54谭小姐阅读 302评论 0 1
  • 随着人囗的迅速增长,工农业的迅速发展,我们生态环境的污染也随之而来。工厂排放的污水,废气废物所产生的污染,人们乱砍...
    王姚文阅读 166评论 0 0