王德貴
圓周率π和e一樣,都是常用的數學常數,且都是無理數,也稱為超越數。在數學和物理課程中,經常用到π和e兩個常數,也經常會想,它們是怎么算出來的?怎么被發現的呢?今天我們來學習用Python求圓周率π的近似值。
其實圓周率求解方法很多,歸納起來有大致四類方法,故稱“求解圓周率四法”。我們選擇其中兩種來學習。
使用隨機數(或更常見的偽隨機數)來解決很多計算問題的方法,稱為蒙特卡洛法,它是統計模擬方法。當所求解問題是某種隨機事件出現的概率,或者是某個隨機變量的期望值時,通過某種“實驗”的方法,以這種事件出現的頻率估計這一隨機事件的概率,或者得到這個隨機變量的某些數字特征,并將其作為問題的解。主要步驟:構造或描述概率過程;實現從已知概率分布抽樣;建立各種估計量。
正方形內部有一個相切的圓,它們的面積之比是π/4(圖1)。

現在,在這個正方形內部,隨機產生n個點,計算它們與中心點的距離,并且判斷是否落在圓的內部。若這些點均勻分布,則圓周率 PI=4 * m/n, 其中m表示落到圓內投點數,n表示總的投點數。
顯然投點數目越多,數值越精確,誤差越小。于是我們采用循環的方法產生n個(0,1.0)區間上的隨機數,即第一象限內的點,然后判斷在單位圓內的個數m,求得圓周率近似值。
從前面的分析, 得出程序如下。這里用到了隨機函數模塊(random),利用random()生成一個(0,1.0)之間的隨機浮點數(圖2)。

單位圓標準方程為x+y=1:,條件判斷里,“x*x+y*y<=1”表示點(x,y)在圓上或圓內。投入1億個點,某次運行結果為3.1515122。
利用隨機數求圓周率,每次求得的近似值不一定相同。
所謂“割圓術”,是用圓內接正多邊形的面積去無限逼近圓面積并以此求取圓周率的方法。
魏晉時期數學家劉徽在《九章算術注》中為“割圓術”計算圓周率建立了嚴密的理論和完善的算法。約480年,南北朝時期的祖沖之在此基礎上算出圓周率在3.1415926和3.1415927之間,相當于精確到小數第7位,他是第一位將圓周率值計算到小數第7位的科學家。
用割圓術來求解圓周率。可以設置初始正多邊形為正三角形、正方形、正六邊形等等,我們用正六邊形,因為此時圓半徑和邊長相等,計算上簡便一些(圖3)。

如圖3所示,設圓半徑為1,這樣圓面積就是s=πr=π,計算出面積,即求出圓周率。基本思路就是無限的成倍分割,這樣新的正多邊形面積,就是在原正多邊形面積基礎上,每條邊與圓弧之間多了一個等腰三角形,對于正六邊形到正十二邊形的分割中,就是多了6個的?△ABC?面積。下面進行計算。

這樣依次計算下去,正多邊形的面積就會越來越接近圓面積,因此,正多邊形邊數越多,計算越精確。
這里用到了math模塊(數學基本運算庫),應用之前需要導入,有兩種導入方法:(1)import math;(2)from math import* 。兩種導入方法的區別,是Python等級考試四級內容,大家可以注意一下,本文不做贅述。sqrt()是開平方運算(圖4)。

運行程序,輸出結果如圖(圖5):

誤差已經很小了。當然如果邊數小于6,那就按六邊形算的了,所以取值要盡量的大,才會更接近圓周率的真值。
Python默認有效數字為17位,那么如何得到更精確的位數呢?
我們可以利用decimal模塊求解18-100位精度的有效數字,默認值為28位。decimal意思為十進制,提供了十進制浮點運算支持,主要是用來處理要求特別精確的小數。decimal所表示的數是完全精確的,而float浮點數不能使用decimal,因為float浮點數本身就不精確。
getcontext().prec = 50,即設置精度為50位(圖6)。結果如圖7。


從輸出可以看到第一行是28位有效數字,這是默認值。第三行是50位有效數字。第二和第四行,與設置精度無關,因為它們是浮點數轉換過來的,不準確。要想精度提高,必須將數值轉換為字符再輸出。
下面是割圓術基礎上,設定有效數字位數,來求近似值的程序(圖8)。

運行結果(圖9):

如果超過100位呢?那么后面就都用0補充了。更多位也不會有區別了(圖10)。

想要更精確的結果,必須設定小數點后位數,可以先放大多少倍,最后計算結束時再舍去多少位,以保證精度,這個過程稍難,程序如圖11。

運行結果可以看到,可以精確到小數點1000位(圖12)。

大家可以自己驗證,程序可以求出任意位數的小數。你注意到了嗎?程序是把計算結果,轉換為字符串后處理輸出的,這是因為輸出整數時,會自動轉換為科學記數法,不會顯示多于17位的有效數字。
SymPy模塊可以進行符號計算,可以定義符號變量,進行代數運算,以及微分運算、積分運算等。利用evalf()函數傳遞數據。比如要輸出1000位有效數字,則程序為(圖13):

結果如圖14:

從前面的分析和測試中,我們看到求解圓周率的方法很多,精度也不盡相同,一般來說Python默認的17位有效數字已經足夠了,如果需要高精度的結果,需要使用相應的方法。
通過圓周率的計算,我們也更深刻了解了圓周率的發展史。電子計算機的出現,把計算精度提高到了驚人的數量級。現在圓周率最高位數已計算到31.4萬億位,準確地說是31415926535897位,是2019年工程師愛瑪在谷歌云平臺的幫助下完成的。比2016年的紀錄又增加數萬億位。
其實π精度的追求更多的是人類對極限的追求,畢竟十位小數的π就足以使地球周界準確到一英寸以內,三十位小數的π便能使整個可見宇宙的四周準確到連最強大的顯微鏡都不能分辨的一個量。
如果你對此感興趣,可以繼續深入研究,求出更多位數,刷新圓周率計算史上的紀錄,成為一名數學家!