王德貴

圓周率π和e一樣，都是常用的數(shù)學(xué)常數(shù)，且都是無(wú)理數(shù)，也稱(chēng)為超越數(shù)。在數(shù)學(xué)和物理課程中，經(jīng)常用到π和e兩個(gè)常數(shù)，也經(jīng)常會(huì)想，它們是怎么算出來(lái)的？怎么被發(fā)現(xiàn)的呢？今天我們來(lái)學(xué)習(xí)用Python求圓周率π的近似值。

其實(shí)圓周率求解方法很多，歸納起來(lái)有大致四類(lèi)方法，故稱(chēng)“求解圓周率四法”。我們選擇其中兩種來(lái)學(xué)習(xí)。

使用隨機(jī)數(shù)（或更常見(jiàn)的偽隨機(jī)數(shù)）來(lái)解決很多計(jì)算問(wèn)題的方法，稱(chēng)為蒙特卡洛法，它是統(tǒng)計(jì)模擬方法。當(dāng)所求解問(wèn)題是某種隨機(jī)事件出現(xiàn)的概率，或者是某個(gè)隨機(jī)變量的期望值時(shí)，通過(guò)某種“實(shí)驗(yàn)”的方法，以這種事件出現(xiàn)的頻率估計(jì)這一隨機(jī)事件的概率，或者得到這個(gè)隨機(jī)變量的某些數(shù)字特征，并將其作為問(wèn)題的解。主要步驟：構(gòu)造或描述概率過(guò)程;實(shí)現(xiàn)從已知概率分布抽樣;建立各種估計(jì)量。

正方形內(nèi)部有一個(gè)相切的圓，它們的面積之比是π/4（圖1）。

現(xiàn)在，在這個(gè)正方形內(nèi)部，隨機(jī)產(chǎn)生n個(gè)點(diǎn)，計(jì)算它們與中心點(diǎn)的距離，并且判斷是否落在圓的內(nèi)部。若這些點(diǎn)均勻分布，則圓周率 PI=4 * m/n， 其中m表示落到圓內(nèi)投點(diǎn)數(shù)，n表示總的投點(diǎn)數(shù)。

顯然投點(diǎn)數(shù)目越多，數(shù)值越精確，誤差越小。于是我們采用循環(huán)的方法產(chǎn)生n個(gè)（0，1.0）區(qū)間上的隨機(jī)數(shù)，即第一象限內(nèi)的點(diǎn)，然后判斷在單位圓內(nèi)的個(gè)數(shù)m，求得圓周率近似值。

從前面的分析， 得出程序如下。這里用到了隨機(jī)函數(shù)模塊（random），利用random（）生成一個(gè)（0，1.0）之間的隨機(jī)浮點(diǎn)數(shù)（圖2）。

單位圓標(biāo)準(zhǔn)方程為x+y=1：，條件判斷里，“x*x+y*y<=1”表示點(diǎn)（x，y）在圓上或圓內(nèi)。投入1億個(gè)點(diǎn)，某次運(yùn)行結(jié)果為3.1515122。

利用隨機(jī)數(shù)求圓周率，每次求得的近似值不一定相同。

所謂“割圓術(shù)”，是用圓內(nèi)接正多邊形的面積去無(wú)限逼近圓面積并以此求取圓周率的方法。

魏晉時(shí)期數(shù)學(xué)家劉徽在《九章算術(shù)注》中為“割圓術(shù)”計(jì)算圓周率建立了嚴(yán)密的理論和完善的算法。約480年，南北朝時(shí)期的祖沖之在此基礎(chǔ)上算出圓周率在3.1415926和3.1415927之間，相當(dāng)于精確到小數(shù)第7位，他是第一位將圓周率值計(jì)算到小數(shù)第7位的科學(xué)家。

用割圓術(shù)來(lái)求解圓周率?？梢栽O(shè)置初始正多邊形為正三角形、正方形、正六邊形等等，我們用正六邊形，因?yàn)榇藭r(shí)圓半徑和邊長(zhǎng)相等，計(jì)算上簡(jiǎn)便一些（圖3）。

如圖3所示，設(shè)圓半徑為1，這樣圓面積就是s=πr=π，計(jì)算出面積，即求出圓周率?；舅悸肪褪菬o(wú)限的成倍分割，這樣新的正多邊形面積，就是在原正多邊形面積基礎(chǔ)上，每條邊與圓弧之間多了一個(gè)等腰三角形，對(duì)于正六邊形到正十二邊形的分割中，就是多了6個(gè)的?△ABC?面積。下面進(jìn)行計(jì)算。

這樣依次計(jì)算下去，正多邊形的面積就會(huì)越來(lái)越接近圓面積，因此，正多邊形邊數(shù)越多，計(jì)算越精確。

這里用到了math模塊（數(shù)學(xué)基本運(yùn)算庫(kù)），應(yīng)用之前需要導(dǎo)入，有兩種導(dǎo)入方法：（1）import math;（2）from math import* 。兩種導(dǎo)入方法的區(qū)別，是Python等級(jí)考試四級(jí)內(nèi)容，大家可以注意一下，本文不做贅述。sqrt（）是開(kāi)平方運(yùn)算（圖4）。

運(yùn)行程序，輸出結(jié)果如圖（圖5）：

誤差已經(jīng)很小了。當(dāng)然如果邊數(shù)小于6，那就按六邊形算的了，所以取值要盡量的大，才會(huì)更接近圓周率的真值。

Python默認(rèn)有效數(shù)字為17位，那么如何得到更精確的位數(shù)呢？

我們可以利用decimal模塊求解18-100位精度的有效數(shù)字，默認(rèn)值為28位。decimal意思為十進(jìn)制，提供了十進(jìn)制浮點(diǎn)運(yùn)算支持，主要是用來(lái)處理要求特別精確的小數(shù)。decimal所表示的數(shù)是完全精確的，而float浮點(diǎn)數(shù)不能使用decimal，因?yàn)閒loat浮點(diǎn)數(shù)本身就不精確。

getcontext（）.prec = 50，即設(shè)置精度為50位（圖6）。結(jié)果如圖7。

從輸出可以看到第一行是28位有效數(shù)字，這是默認(rèn)值。第三行是50位有效數(shù)字。第二和第四行，與設(shè)置精度無(wú)關(guān)，因?yàn)樗鼈兪歉↑c(diǎn)數(shù)轉(zhuǎn)換過(guò)來(lái)的，不準(zhǔn)確。要想精度提高，必須將數(shù)值轉(zhuǎn)換為字符再輸出。

下面是割圓術(shù)基礎(chǔ)上，設(shè)定有效數(shù)字位數(shù)，來(lái)求近似值的程序（圖8）。

運(yùn)行結(jié)果（圖9）：

如果超過(guò)100位呢？那么后面就都用0補(bǔ)充了。更多位也不會(huì)有區(qū)別了（圖10）。

想要更精確的結(jié)果，必須設(shè)定小數(shù)點(diǎn)后位數(shù)，可以先放大多少倍，最后計(jì)算結(jié)束時(shí)再舍去多少位，以保證精度，這個(gè)過(guò)程稍難，程序如圖11。

運(yùn)行結(jié)果可以看到，可以精確到小數(shù)點(diǎn)1000位（圖12）。

大家可以自己驗(yàn)證，程序可以求出任意位數(shù)的小數(shù)。你注意到了嗎？程序是把計(jì)算結(jié)果，轉(zhuǎn)換為字符串后處理輸出的，這是因?yàn)檩敵稣麛?shù)時(shí)，會(huì)自動(dòng)轉(zhuǎn)換為科學(xué)記數(shù)法，不會(huì)顯示多于17位的有效數(shù)字。

SymPy模塊可以進(jìn)行符號(hào)計(jì)算，可以定義符號(hào)變量，進(jìn)行代數(shù)運(yùn)算，以及微分運(yùn)算、積分運(yùn)算等。利用evalf（）函數(shù)傳遞數(shù)據(jù)。比如要輸出1000位有效數(shù)字，則程序?yàn)椋▓D13）：

結(jié)果如圖14：

從前面的分析和測(cè)試中，我們看到求解圓周率的方法很多，精度也不盡相同，一般來(lái)說(shuō)Python默認(rèn)的17位有效數(shù)字已經(jīng)足夠了，如果需要高精度的結(jié)果，需要使用相應(yīng)的方法。

通過(guò)圓周率的計(jì)算，我們也更深刻了解了圓周率的發(fā)展史。電子計(jì)算機(jī)的出現(xiàn)，把計(jì)算精度提高到了驚人的數(shù)量級(jí)?，F(xiàn)在圓周率最高位數(shù)已計(jì)算到31.4萬(wàn)億位，準(zhǔn)確地說(shuō)是31415926535897位，是2019年工程師愛(ài)瑪在谷歌云平臺(tái)的幫助下完成的。比2016年的紀(jì)錄又增加數(shù)萬(wàn)億位。

其實(shí)π精度的追求更多的是人類(lèi)對(duì)極限的追求，畢竟十位小數(shù)的π就足以使地球周界準(zhǔn)確到一英寸以?xún)?nèi)，三十位小數(shù)的π便能使整個(gè)可見(jiàn)宇宙的四周準(zhǔn)確到連最強(qiáng)大的顯微鏡都不能分辨的一個(gè)量。

如果你對(duì)此感興趣，可以繼續(xù)深入研究，求出更多位數(shù)，刷新圓周率計(jì)算史上的紀(jì)錄，成為一名數(shù)學(xué)家！