圆周率的计算史源远流长,从中国古代的祖冲之到现代科学计算器,将圆周率\(\pi\)精确到小数点后的位数是人类一直在追求的科学事业,发展出来的计算公式也是层出不穷:
下面用while循环写出Machain公式计算圆周率的代码:
def arctg(x,N): #(x,N=5)则代表N的缺省值是5,即调用函数时不给值就用5
n=1
m=0
while n<=N: #控制计算N次,可以提高N以提高计算精度
m+=(-1)**(n-1)*x**(2*n-1)/(2*n-1)
n+=1
return m
pi=4*(4*arctg(1/5,5)-arctg(1/239,5))
N=5
print('当N=%d时,圆周率的近似值为%.20f'%(N,pi))
但是当N大于11时会发现输出的pi值没有发生变化了
PS C:\python> python '.\ln(x)''.py'
当N=11时,圆周率的近似值为3.14159268240439937259
PS C:\python> python '.\ln(x)''.py'
当N=11时,圆周率的近似值为3.14159268240439937259
出现这种情况的原因时计算机只能存储整数,所以实数都是约数,这样浮点运算会有误差
所以我们可以测试一下python对长整数的支持:
x=12345679101112131415161718192021222324252627282930
print(x)
PS C:\python> python '.\ln(x)''.py'
12345679101112131415161718192021222324252627282930
说明python可以准确地输出长整数,因此我们考虑将浮点数转化为整数来保证计算的精度,例如我们计算1/239到小数点后30位:
precision=30 #精确到小数点后的位数
largen=10**(precision+5) #一般多延5位保证取舍不要影响精度
x=(1*largen)//239 #用扩大后的数整除239
m=x//(10**5) #去掉多延的5位
print('当精度为小数点后%d位时,输出的结果为0.00%d'%(precision,m))
我们先对Machain公式进行适当的变换:
这样做的目的是让两个数列相同项数的数先相减,再求和,用除以公比来找到下一项的方式设计程序
接下来就按照上面的例子,对相应的代码进行优化:
precision=100 #精确到小数点后的位数
largen=10**(precision+5) #一般多延5位保证取舍不要影响精度
N=precision*2 #将循环次数设置为两倍的精确度
term_1=4*largen//5 #用整除得到第一组数值
term_2=-largen//239
m=term_1+term_2 #计算得到第一个括号的值
for i in range(N-1): #不要忘了i是从0开始取值的
n=i+2
term_1=term_1//25 #得到下一个括号里的值,注意全都要用整除
term_2=term_2//239**2
m+=(term_1+term_2)*(-1)**(n-1)/(2*n-1)
pi=4*m//10**5 #去掉多延的5位
pi=pi%10**(precision) #除以10的100次方取余数是为了将输出写成3.几的形式
print('当精度为%d时,圆周率计算结果为3.%d'%(precision,pi))
这样我们就得到了我们想要的结果
PS C:\python> python '.\ln(x)''.py'
当精度为100时,圆周率计算结果为3.1415926535897923187905000975459403616237497995798551546679370505449715292938495673272074805945827328
在上个随笔提出精确度的疑问过后,这次对其中一个小部分(如何用整数实现浮点数的精确度问题)进行了探索和记录。刚好今晚考完高数就来赶博客了,比起高数还是目前的大计爱我啊QWQ
继续加油吧!