大家好，这是我第一次写 Python，所以先挑个简单的题材来练习。
∞
使用的级数是 e = Σ (1/k!) = 1/0! + 1/1! + 1/2! + 1/3! + 1/4! + .....
k=0
先用 Stirling's approximation 估计大概需要计算前几项才能达到需要的精度
使用 Divide & Conquer 法，将级数对半切再对半切再对半切再对半切再对半切....
分开算成分数后再通分加起来。这样所要做的运算总量比起 for 循环从头加到尾
并没有变少，但是先算数字小的运算比较快，只有后期的巨大通分、最后一个巨大除法
以及转成十进制会特别慢。
大数运算的部份我使用了 gmpy2 module 的 mpz (整数运算) 及 mpfr (实数浮点运算)
我知道 Python 内建整数的大数运算，我也有试过，但速度还是比不上 mpz
至于 float，Python 似乎只有 53 bits 的 double float 精度，
所以要算到小数点以下十亿位的除法只能用 mpfr 了。
为了知道程式是否还有在跑，进度多少了，我做了一个进度条。
然而大部份的时间是花在进度 100% 之后的大除法和转成十进制，
所以大家看到 100% 就不动了请不要惊慌，它在 100% 之后还要跑很久
由于是第一次接触 Python，应该有很多地方写得不太正确，或是效率可以再改进，
请各位板友不吝指教。
☆ ☆ ☆
以下是程式码，
如果不加命令列参数的话，默认是计算十万位，所花时间不到一秒。
而在我的电脑上跑十亿位需要 52 分钟，最多消耗 10GB RAM
命令列参数像下面这样：
$ python3 e.py 1000000000
计算结果会存在纯文字档 e-py.txt 里面，十亿位的话档案大小约 1.2GB
如果是 UNIX 系的作业系统有 time 指令可以测速：
$ time python3 e.py 1000000000
为了方便大家剪贴下来实测，我也把程式码放在 ideone.com
但 ideone.com 不提供 gmpy2 module 因此无法直接在网站上执行。
# URL: https://ideone.com/TQluT9
#
# e.py - Calculate Eular's number e
#
import sys
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
old_p = 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():
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("e-py.txt", mode="w")
fd.write(" e = ")
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):
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)
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 = s(0, n_terms)
pf = mpfr(p);
qf = mpfr(q);
ef = gmpy2.div(p, q)
estr, exp, prec = mpfr.digits(ef)
estr = estr[0:d]
write_string(estr);
#
# main program
#
argc = len(sys.argv)
if (argc >= 2):
digits = int(sys.argv[1])
else:
digits = 100000
calc_e(digits)
# End of e.py