π是一個(gè)無(wú)數(shù)人追隨的真正的神奇數(shù)字。我不是很清楚一個(gè)永遠(yuǎn)重復(fù)的無(wú)理數(shù)的迷人之處。在我看來(lái)，我樂(lè)于計(jì)算π，也就是計(jì)算π的值。因?yàn)棣惺且粋€(gè)無(wú)理數(shù)，它是無(wú)限的。這就意味著任何對(duì)π的計(jì)算都僅僅是個(gè)近似值。如果你計(jì)算100位，我可以計(jì)算101位并且更精確。迄今為止，有些人已經(jīng)選拔出超級(jí)計(jì)算機(jī)來(lái)試圖計(jì)算最精確的π。一些極值包括 計(jì)算π的5億位。你甚至能從網(wǎng)上找到包含 π的一百億位的文本文件（注意啦！下載這個(gè)文件可能得花一會(huì)兒時(shí)間，并且沒(méi)法用你平時(shí)使用的記事本應(yīng)用程序打開。）。對(duì)于我而言，如何用幾行簡(jiǎn)單的Python來(lái)計(jì)算π才是我的興趣所在。
你總是可以 使用 math.pi 變量的 。它被 包含在 標(biāo)準(zhǔn)庫(kù)中， 在你試圖自己 計(jì)算它之前，你應(yīng)該去使用它 。 事實(shí)上 ， 我們將 用它來(lái)計(jì)算 精度 。作為 開始， 讓我們看 一個(gè) 非常直截了當(dāng)?shù)?計(jì)算Pi的 方法 。像往常一樣，我將使用Python 2.7，同樣的想法和代碼可能應(yīng)用于不同的版本。我們將要使用的大部分算法來(lái)自Pi WikiPedia page并加以實(shí)現(xiàn)。讓我們看看下面的代碼：
 

importsys
importmath
 
defmain(argv):
 
  iflen(argv) !=1:
    sys.exit('Usage: calc_pi.py <n>')
 
  print'\nComputing Pi v.01\n'
   
  a=1.0
  b=1.0/math.sqrt(2)
  t=1.0/4.0
  p=1.0
     
  foriinrange(int(sys.argv[1])):
    at=(a+b)/2
    bt=math.sqrt(a*b)
    tt=t-p*(a-at)**2
    pt=2*p
     
    a=at;b=bt;t=tt;p=pt
     
  my_pi=(a+b)**2/(4*t)
  accuracy=100*(math.pi-my_pi)/my_pi
     
  print"Pi is approximately: "+str(my_pi)
  print"Accuracy with math.pi: "+str(accuracy)
   
if__name__=="__main__":
  main(sys.argv[1:])

這是個(gè)非常簡(jiǎn)單的腳本，你可以下載，運(yùn)行，修改，和隨意分享給別人。你能夠看到類似下面的輸出結(jié)果： 

 你會(huì)發(fā)現(xiàn)，盡管 n 大于4 ,我們逼近 Pi 精度卻沒(méi)有多大的提升。 我們可以猜到即使 n的值更大，同樣的事情（pi的逼近精度沒(méi)有提升）依舊會(huì)發(fā)生。幸運(yùn)的是，有不止一種方法來(lái)揭開這個(gè)謎。使用 Python Decimal （十進(jìn)制）庫(kù)，我們可以就可以得到更高精度的值來(lái)逼近Pi。讓我們來(lái)看看庫(kù)函數(shù)是如何使用的。這個(gè)簡(jiǎn)化的版本，可以得到多于11位的數(shù)字 通常情況小Python 浮點(diǎn)數(shù)給出的精度。下面是Python Decimal 庫(kù)中的一個(gè)例子 ：

wpid-python_decimal_example-2013-05-28-12-54.png

看到這些數(shù)字。不對(duì)！ 我們輸入的僅是 3.14，為什么我們得到了一些垃圾(junk)？ 這是內(nèi)存垃圾(memory junk)。 簡(jiǎn)單點(diǎn)說(shuō)，Python給你你想要的十進(jìn)制數(shù)，再加上一點(diǎn)點(diǎn)額外的值。 只要精度小于垃圾數(shù)，它不會(huì)影響任何計(jì)算。通過(guò)設(shè)置getcontext().prec 你可以的到你想要的位數(shù) 。我們?cè)囋嚒?

看到這些數(shù)字。不對(duì)！ 我們輸入的僅是 3.14，為什么我們得到了一些垃圾(junk)？ 這是內(nèi)存垃圾(memory junk)。 簡(jiǎn)單點(diǎn)說(shuō)，Python給你你想要的十進(jìn)制數(shù)，再加上一點(diǎn)點(diǎn)額外的值。 只要精度小于垃圾數(shù)，它不會(huì)影響任何計(jì)算。通過(guò)設(shè)置getcontext().prec 你可以的到你想要的位數(shù) 。我們?cè)囋嚒?

很好。 現(xiàn)在讓我們 試著用這個(gè) 來(lái) 看看我們是否能 與我們以前的 代碼 有更好的 逼近 。 現(xiàn)在， 我通常 是反對(duì) 使用“ from library import * ” ， 但在這種情況下， 它會(huì) 使代碼 看起來(lái)更漂亮 。
 

importsys
importmath
fromdecimalimport*
 
defmain(argv):
 
  iflen(argv) !=1:
    sys.exit('Usage: calc_pi.py <n>')
 
  print'\nComputing Pi v.01\n'
   
  a=Decimal(1.0)
  b=Decimal(1.0/math.sqrt(2))
  t=Decimal(1.0)/Decimal(4.0)
  p=Decimal(1.0)
     
  foriinrange(int(sys.argv[1])):
    at=Decimal((a+b)/2)
    bt=Decimal(math.sqrt(a*b))
    tt=Decimal(t-p*(a-at)**2)
    pt=Decimal(2*p)
     
    a=at;b=bt;t=tt;p=pt
     
  my_pi=(a+b)**2/(4*t)
  accuracy=100*(Decimal(math.pi)-my_pi)/my_pi
     
  print"Pi is approximately: "+str(my_pi)
  print"Accuracy with math.pi: "+str(accuracy)
   
if__name__=="__main__":
  main(sys.argv[1:])

 
輸出結(jié)果: 

 好了。我們更準(zhǔn)確了，但看起來(lái)似乎有一些舍入。從n = 100和n = 1000，我們有相同的精度?，F(xiàn)在怎么辦？好吧，現(xiàn)在我們來(lái)求助于公式。到目前為止，我們計(jì)算Pi的方式是通過(guò)對(duì)幾部分加在一起。我從DAN 的關(guān)于Calculating Pi 的文章中發(fā)現(xiàn)一些代碼。他建議我們用以下3個(gè)公式：

    Bailey–Borwein–Plouffe 公式
   Bellard的公式
    Chudnovsky 算法

讓我們從Bailey–Borwein–Plouffe 公式開始。它看起來(lái)是這個(gè)樣子： 

 在代碼中我們可以這樣編寫它：
 

import sys
import math
from decimal import *
 
def bbp(n):
  pi=Decimal(0)
  k=0
  while k < n:
    pi+=(Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))
    k+=1
  return pi
 
def main(argv):
 
    if len(argv) !=2:
    sys.exit('Usage: BaileyBorweinPlouffe.py <prec> <n>')
     
  getcontext().prec=(int(sys.argv[1]))
  my_pi=bbp(int(sys.argv[2]))
  accuracy=100*(Decimal(math.pi)-my_pi)/my_pi
 
  print"Pi is approximately "+str(my_pi)
  print"Accuracy with math.pi: "+str(accuracy)
   
if __name__=="__main__":
  main(sys.argv[1:])

 
拋開“ 包裝”的代碼，BBP（N）的功能是你真正想要的。你給它越大的N和給 getcontext().prec 設(shè)置越大的值，你就會(huì)使計(jì)算越精確。讓我們看看一些代碼結(jié)果：

這有許多數(shù)字位。你可以看出，我們并沒(méi)有比以前更準(zhǔn)確。所以我們需要前進(jìn)到下一個(gè)公式，貝拉公式，希望能獲得更好的精度。它看起來(lái)像這樣： 

 我們將只改變我們的變換公式，其余的代碼將保持不變。點(diǎn)擊這里下載Python實(shí)現(xiàn)的貝拉公式。讓我們看一看bellards(n)：
 

def bellard(n):
  pi=Decimal(0)
  k=0
  while k < n:
    pi+=(Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1)+Decimal(1)/(10*k+9)-Decimal(64)/(10*k+3)-Decimal(32)/(4*k+1)-Decimal(4)/(10*k+5)-Decimal(4)/(10*k+7)-Decimal(1)/(4*k+3))
    k+=1
  pi=pi*1/(2**6)
  return pi

   哦，不，我們得到的是同樣的精度。好吧，讓我們?cè)囋嚨谌齻€(gè)公式， Chudnovsky 算法，它看起來(lái)是這個(gè)樣子： 

   再一次，讓我們看一下這個(gè)計(jì)算公式（假設(shè)我們有一個(gè)階乘公式）。 點(diǎn)擊這里可下載用 python 實(shí)現(xiàn)的 Chudnovsky 公式。

下面是程序和輸出結(jié)果：
 

def chudnovsky(n):
  pi=Decimal(0)
  k=0
  while k < n:
    pi+=(Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))*(13591409+545140134*k)/(640320**(3*k)))
    k+=1
  pi=pi*Decimal(10005).sqrt()/4270934400
  pi=pi**(-1)
  return pi

    所以我們有了什么結(jié)論？花哨的算法不會(huì)使機(jī)器浮點(diǎn)世界達(dá)到更高標(biāo)準(zhǔn)。我真的很期待能有一個(gè)比我們用求和公式時(shí)所能得到的更好的精度。我猜那是過(guò)分的要求。如果你真的需要用PI，就只需使用math.pi變量了。然而，作為樂(lè)趣和測(cè)試你的計(jì)算機(jī)真的能有多快，你總是可以嘗試第一個(gè)計(jì)算出Pi的百萬(wàn)位或者更多位是幾。

最新評(píng)論