π是一个无数人追随的真正的神奇数字。我不是很清楚一个永远重复的无理数的迷人之处。在我看来,我乐于计算π,也就是计算π的值。因为π是一个无理数,它是无限的。这就意味着任何对π的计算都仅仅是个近似值。如果你计算100位,我可以计算101位并且更精确。迄今为止,有些人已经选拔出超级计算机来试图计算最精确的π。一些极值包括 计算π的5亿位。你甚至能从网上找到包含 π的一百亿位的文本文件(注意啦!下载这个文件可能得花一会儿时间,并且没法用你平时使用的记事本应用程序打开。)。对于我而言,如何用几行简单的python来计算π才是我的兴趣所在。
你总是可以 使用 math.pi 变量的 。它被 包含在 标准库中, 在你试图自己 计算它之前,你应该去使用它 。 事实上 , 我们将 用它来计算 精度 。作为 开始, 让我们看 一个 非常直截了当的 计算pi的 方法 。像往常一样,我将使用python 2.7,同样的想法和代码可能应用于不同的版本。我们将要使用的大部分算法来自 pi wikipedia page并加以实现。让我们看看下面的代码:
importsys
importmath
defmain(argv):
iflen(argv) !=1:
sys.exit(‘usage: calc_pi.py ‘)
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:])
这是个非常简单的脚本,你可以下载,运行,修改,和随意分享给别人。你能够看到类似下面的输出结果:
看到这些数字。不对! 我们输入的仅是 3.14,为什么我们得到了一些垃圾(junk)? 这是内存垃圾(memory junk)。 在坚果壳,python给你你想要的十进制数,再加上一点点额外的值。 只要精度小于垃圾号码开始时,它不会影响任何计算只要精度小于前面的垃圾号码(junk number)开始时。 您可以指定你想要多少位数的通过设置getcontext().prec 。我们试试。
很好。 现在让我们 试着用这个 来 看看我们是否能 与我们以前的 代码 有更好的 逼近 。 现在, 我通常 是反对 使用“ from library import * ” , 但在这种情况下, 它会 使代码 看起来更漂亮 。
importsys
importmath
fromdecimalimport*
defmain(argv):
iflen(argv) !=1:
sys.exit(‘usage: calc_pi.py ‘)
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:])
输出结果:
在代码中我们可以这样编写它:
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 ')
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 设置越大的值,你就会使计算越精确。让我们看看一些代码结果:
我们将只改变我们的变换公式,其余的代码将保持不变。点击这里下载python实现的贝拉公式。让我们看一看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
输出结果:
再一次,让我们看一下这个计算公式(假设我们有一个阶乘公式)。 点击这里可下载用 python 实现的 chudnovsky 公式。
下面是程序和输出结果:
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
所以我们有了什么结论?花哨的算法不会使机器浮点世界达到更高标准。我真的很期待能有一个比我们用求和公式时所能得到的更好的精度。我猜那是过分的要求。如果你真的需要用pi,就只需使用math.pi变量了。然而,作为乐趣和测试你的计算机真的能有多快,你总是可以尝试第一个计算出pi的百万位或者更多位是几。