偽·從零開始學算法 - 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。由勾股定理，得

得x6 = 1。

而且我們可以得到，正2n邊形面積等於正n邊形面積加n個等腰三角形的面積，即

隨著n的增大，S(2n)不斷趨近於π。

於是我們可以據此編制割圓術的算法：

割圓術(1)

據此求在正n邊形下，圓周率的下限和上限：

割圓術(2)

割圓術在16~17世紀之前一直都是精確計算圓周率的解決方法。奧地利天文學家克裡斯託夫·格林伯格在1630年用10^40邊形計算到第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]範圍的隨機數（兩坐標不一定相同），重複若干次」，把「落在四分之一圓內」這個條件轉化為「點是否在x^2+y^2=1內」。於是得下面的流程圖：

蒙特卡羅方法流程圖

它的缺點也很明顯：隨機、不精確。

結語

還有很多算法，在此我就不一一列舉了。

最後我放上一張圖，以展示π的計算的發展。

（當數學家發現新的算法、電腦變得普及時，π的已知小數位急劇增加。注意垂直坐標使用了對數坐標）

一般而言，π值並不需要過於精確便能夠滿足大部分的數學運算的需求。按照約爾格·阿恩特（Jrg Arndt）及克裡斯託夫·黑內爾（Christoph Haenel）的計算，39個數位已足夠運算絕大多數的宇宙學的計算需求，因為這個精確度已能夠將可觀測宇宙圓周的精確度準確至一個原子大小。但是大家仍然對π的計算樂此不疲，這有很多原因，但計算π的本身讓我想起了一個故事：

有人問英國登山家喬治·馬洛裡為何想要攀登珠穆朗瑪峰。他回答說：「因為山就在那兒。」

參考資料

圓周率 - 維基百科，自由的百科全書圓周率計算年表 - 維基百科，自由的百科全書割圓術 (劉徽) - 維基百科，自由的百科全書普通高中課程標準實驗教科書（必修） 數學（A版） 必修3. 人民教育出版社高斯-勒讓德算法 - 維基百科，自由的百科全書