作者Schottky (顺风相送)
看板Python
标题[心得] 计算ｅ到小数点下十亿位 ─ 超进化版
时间Fri Feb 26 03:08:49 2021
感谢 raagi 与快乐的朋友们指导，我把程式里不少地方都修改过了，分号也都删了，
执行速度快了将近三倍，接近 C 语言做同样运算的速度了。
过程中发现：
1. 要注意 gmpy2 module 的版本，我用 2.1.0a4 才有 mpfr.digits() 可以用，
但 2.0.8 却没有 mpfr.digits()
因此增加了 gmpy2, GMP, MPFR 的版本显示，方便诊断问题。
2. Windows 的 Python 不知为何 MPFR precision 偏低，低於十亿位的要求，
这次也特别在执行运算前检查这个问题并警告。
程式码一样顺便放在
https://ideone.com/3B6wdb 方便大家复制贴上
#!/usr/bin/env python3
#
# e-mpz.py - Calculate Eular's number e
#
import sys
import time
import math
import gmpy2
from gmpy2 import mpfr
from gmpy2 import mpz
#
# Constants used in Stirling's approximation
#
E = float(2.718281828459045235360287)
pi = float(3.141592653589793238462643)
C = math.log10(2*pi) / 2
#
# Global Variables
#
count = 0
total = 0
grad = 0
step = 0
#
# Stirling's approximation
#
def logfactorial(n):
return (C + math.log10(n)/2 + n*(math.log10(n)-math.log10(E)))
#
# Estimate how many terms in the serie sould be calculated.
#
def terms(digits):
upper = 2
lower = 1
while (logfactorial(upper)<digits):
upper <<= 1
else:
lower = upper/2
while ((upper-lower) > 1):
n = (upper+lower)/2
if (logfactorial(n) > digits):
upper = n
else:
lower = n
return n
#
# Show Progress
#
def progress_init(max):
global count, total, grad, step
total = max
count = 0
step = int(total / 1000)
grad = int(step / 2)
def progress():
global count, total, grad, step
if (count > grad):
grad += step
g = int(math.floor(72.5*count/total+0.5))
p = int(math.floor(1000.5*count/total+0.5))
msg = "H" * g + "-" * (72-g) + " " + str(p/10) + "%\r"
if (grad > total):
msg += "\n"
print(msg, sep="", end="", flush=True)
#
# Write digit string
#
def write_string(digit_string):
fd = open("e-py.txt", mode="w")
fd.write(" e = ")
fd.write(digit_string[0])
fd.write(".")
for c in range(1, len(digit_string)-1, 50):
if (c != 1):
fd.write("\t")
fd.write(digit_string[c:c+50])
if ((c % 1000) == 951):
fd.write(" << ")
fd.write(str(c+49))
fd.write("\r\n")
elif ((c % 500) == 451):
fd.write(" <\r\n")
else:
fd.write("\r\n")
# Final new-line
fd.write("\r\n")
fd.close()
#
# Recursive funcion.
#
def s(a, b):
global count
m = math.ceil((a + b) / 2)
if (a == b):
q = mpz(1)
if (a == 0):
p = mpz(1)
else:
p = mpz(0)
elif (b - a == 1):
if (a == 0):
p = mpz(2)
q = mpz(1)
else:
p = mpz(1)
q = mpz(b)
else:
p1, q1 = s(a, m)
p2, q2 = s(m, b)
# Merge
p = gmpy2.add(gmpy2.mul(p1, q2), p2)
q = gmpy2.mul(q1, q2)
count += 1
progress()
return p, q
#
# Calculate e
#
def calc_e(digits):
global total
d = digits+1
n_terms = int(terms(d))
precision = math.ceil(d * math.log2(10)) + 4
print("d = ", d, ", n = ", n_terms, ", precision = ", precision)
print("gmpy2 version:", gmpy2.version())
print("MP version:", gmpy2.mp_version())
print("MPFR version:", gmpy2.mpfr_version())
max_precision = gmpy2.get_max_precision()
print("max_precision =", max_precision)
max_emax = gmpy2.get_emax_max()
print("max_emax =", max_emax)
if (max_precision < precision):
print("Error! Max precision is too small! Program terminated.")
return
gmpy2.get_context().precision = precision
gmpy2.get_context().emax = max_emax
print("Real precision = ", gmpy2.get_context().precision)
progress_init(n_terms * 2 - 1) # Initialize progress bar
start_time = time.monotonic_ns()
p, q = s(0, n_terms)
end_time = time.monotonic_ns()
elapsed = (end_time - start_time) / 1000000000
print("Recursion:", elapsed, "seconds.")
start_time = time.monotonic_ns()
pf = mpfr(p)
qf = mpfr(q)
ef = gmpy2.div(pf, qf)
end_time = time.monotonic_ns()
elapsed = (end_time - start_time) / 1000000000
print("Grand division:", elapsed, "seconds.")
start_time = time.monotonic_ns()
estr, exp, prec = mpfr.digits(ef)
estr = estr[0:d]
end_time = time.monotonic_ns()
elapsed = (end_time - start_time) / 1000000000
print("Convert to decimal digits:", elapsed, "seconds.")
start_time = time.monotonic_ns()
write_string(estr)
end_time = time.monotonic_ns()
elapsed = (end_time - start_time) / 1000000000
print("Write file:", elapsed, "seconds.")
#
# main program
#
if __name__ == '__main__':
argc = len(sys.argv)
if (argc >= 2):
digits = int(sys.argv[1])
else:
digits = 100000
calc_e(digits)
# End of e-mpz.py
--
桃乐丝: 可是, 如果你没有头脑, 为什麽会说话?
稻草人: ㄝ, 我也不知... 但是有些人没有头脑也能说超～多话呢。
--
※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 118.165.64.143 (台湾)
※ 文章网址: https://webptt.com/cn.aspx?n=bbs/Python/M.1614280134.A.B3D.html
1F：推 TuCH: pi的精度会影响e的精度吗 02/26 08:18
π和ｅ是独立的两支程式，不会互相影响哟！
但是 gmpy2 module (或者说 MPFR library) 本身是有精度 (precision) 上限
只不过在我的运作环境 (Debian Linux) 中，这个上限非常高，
足足有 4,611,686,018,427,387,903 bits，在撞到上限前应该只需要担心 RAM 不够
不要说 RAM 了，连硬碟都没见过这麽大的硬碟 (四百万TB)
至於 Windows 环境，我自己并没有真的在 Windows 中测试过，说不定是我误会了
不知道各位能不能帮我试试看，现在这个版本在 Windows 能不能跑到十亿位？
2F：→ OrzOGC: 好强,我只会写粪code...QQ 02/26 20:25
感谢称赞 :)
能达到要求就是好 code (拇指)
※ 编辑: Schottky (111.250.54.51 台湾), 02/28/2021 06:08:43