作者Schottky (顺风相送)
看板Python
标题[心得] 计算π到小数点下十亿位
时间Tue Feb 23 18:28:22 2021
上一篇的程式稍微改变一下可以用来计算圆周率π (根本改超多)
使用 Chudnovsky 公式
https://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
R(a,b)*P(a,b)
定义 S(a,b) = -------------
Q(a,b)
(6b)! (6a)!
R(a,b) = ----- ÷ ----- (第 a+1 项到第 b 项提出来的乘数)
(3b)! (3a)!
b!
Q(a,b) = (----)^3*(640320)^(3*(b-a)) (第 a+1 项到第 b 项的分母)
a!
P(a,b) = 第 a+1 项到第 b 项的 (13591409+545140134*k) 分子部份通分後的和
接着就和上一篇差不多，把级数对半切再对半切再对半切再对半切......
最後再一层一层通分合并起来即可。
☆ ☆ ☆
其实我还没跑完十亿位的计算，估计要一两个小时 (被打)
用我的电脑计算一亿位大约六分半钟，消耗 1.2GB RAM
(注：後来跑完了，耗时 1 小时 28 分钟，RAM 用量最大值 12GB)
以下是程式码，没给命令列参数的话预设一样是十万位数。
如果要计算十亿位，命令列参数像这样指定：
$ python3 pi.py 1000000000
注意要用 python3，假如出现语法错误 SyntaxError: invalid syntax 可能是使用到
python 2，可以用 python --version 确认版本
UNIX 系的作业系统有 time 指令可以测速
$ time python3 pi.py 1000000000
程式码一样放在 ideone.com 方便大家复制贴上，
但因为网站不提供 gmpy2 module 而无法执行。
#
https://ideone.com/jAiHdM
#
# pi.py - Calculate Pi
#
import sys
import math
import gmpy2
from gmpy2 import mpfr
from gmpy2 import mpz
#
# Global Variables
#
count = 0
total = 0
old_p = 0
#
# Show Progress
#
def progress():
global count, old_p, total
p = int(math.floor(1000*count/total+0.5))
if (p > old_p):
old_p = p
g = int(math.floor(72.5*count/total+0.5))
for c in range(72):
if (c<g):
print("H", sep="", end="")
else:
print("-", sep="", end="")
print(" ", p/10, "%\r", sep="", end="", flush=True)
if (count == total):
print("\n", sep="", end="")
#
# Write digit string
#
def write_string(digit_string):
fd = open("pi-py.txt", mode="w")
fd.write(" pi = ")
for c in range(len(digit_string)):
if ((c != 1) and (c % 50 == 1)):
fd.write("\t")
fd.write(digit_string[c])
if (c == 0):
fd.write(".")
elif ((c % 1000) == 0):
fd.write(" << ")
fd.write(str(c))
fd.write("\r\n")
elif ((c % 500) == 0):
fd.write(" <\r\n")
elif ((c % 50) == 0):
fd.write("\r\n")
elif ((c % 5) == 0):
fd.write(" ")
# Final new-line
if ((c%50) != 0):
fd.write("\r\n")
fd.close()
#
# Recursive funcion.
#
def s(a, b, max):
global count
m = math.ceil((a + b) / 2)
if (b - a == 1):
if (a == 0):
r = 120 # 6!
q = mpz(640320**3)
p = gmpy2.sub( gmpy2.mul(q, 13591409),
gmpy2.mul(r, 13591409+545140134) )
else:
r = mpz(8 * (a*6+1) * (a*6+3) * (a*6+5))
q = mpz((b*640320)**3)
if ((b%2) == 0):
p = gmpy2.mul(mpz(13591409 + b*545140134), r)
else:
p = gmpy2.mul(mpz(-13591409 - b*545140134), r)
else:
p1, q1, r1 = s(a, m, max)
p2, q2, r2 = s(m, b, max)
# Merge
p = gmpy2.add( gmpy2.mul(p1, q2), gmpy2.mul(p2, r1) )
q = gmpy2.mul(q1, q2)
if (b != max):
r = gmpy2.mul(r1, r2)
else:
r = 0
count += 1
progress()
return p, q, r;
#
# Calculate e
#
def calc_pi(digits):
global total
d = digits+1
n_terms = math.ceil(d*math.log(10)/(3*math.log(53360)))
precision = math.ceil(d * math.log2(10)) + 4;
print("d = ", d, ", n = ", n_terms, ", precision = ", precision, sep="")
gmpy2.get_context().precision = precision
gmpy2.get_context().emax = 1000000000000;
print("Real precision =", gmpy2.get_context().precision)
total = n_terms * 2 - 1 # Max progress
p, q, r = s(0, n_terms, n_terms)
q = gmpy2.mul(q, 426880)
print("Recursion done. Processing division.")
ef = gmpy2.div(q, p)
ef = gmpy2.mul(ef, gmpy2.sqrt(10005))
print("Division done. Converting to digits.")
estr, exp, prec = mpfr.digits(ef)
estr = estr[0:d]
print("Writing to file.")
write_string(estr);
#
# main program
#
argc = len(sys.argv)
if (argc >= 2):
digits = int(sys.argv[1])
else:
digits = 100000
calc_pi(digits)
# End of pi.py
--
桃乐丝: 可是, 如果你没有头脑, 为什麽会说话?
稻草人: ㄝ, 我也不知... 但是有些人没有头脑也能说超～多话呢。
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 111.250.31.177 (台湾)
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Python/M.1614076108.A.9EC.html
※ 编辑: Schottky (111.250.31.177 台湾), 02/23/2021 19:25:34