感谢 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