[心得] 計算ｅ到小數點下十億位 ─ 超進化版
感謝 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://www.ptt.cc/bbs/Python/M.1614280134.A.B3D.html
推
02/26 08:18,
3年前
, 1F
02/26 08:18, 1F
π和ｅ是獨立的兩支程式，不會互相影響喲！
但是 gmpy2 module (或者說 MPFR library) 本身是有精度 (precision) 上限
只不過在我的運作環境 (Debian Linux) 中，這個上限非常高，
足足有 4,611,686,018,427,387,903 bits，在撞到上限前應該只需要擔心 RAM 不夠
不要說 RAM 了，連硬碟都沒見過這麼大的硬碟 (四百萬TB)
至於 Windows 環境，我自己並沒有真的在 Windows 中測試過，說不定是我誤會了
不知道各位能不能幫我試試看，現在這個版本在 Windows 能不能跑到十億位？
→
02/26 20:25,
3年前
, 2F
02/26 20:25, 2F
感謝稱讚 :)
能達到要求就是好 code (拇指)
※ 編輯: Schottky (111.250.54.51 臺灣), 02/28/2021 06:08:43