[心得] 計算 e 到小數點下十億位
大家好，這是我第一次寫 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
--
桃樂絲: 可是, 如果你沒有頭腦, 為什麼會說話?
稻草人: ㄝ, 我也不知... 但是有些人沒有頭腦也能說超～多話呢。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 111.250.31.177 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Python/M.1614067863.A.B21.html
推
02/23 16:58,
3年前
, 1F
02/23 16:58, 1F
→
02/23 18:30,
3年前
, 2F
02/23 18:30, 2F
推
02/23 19:07,
3年前
, 3F
02/23 19:07, 3F
推
02/23 21:32,
3年前
, 4F
02/23 21:32, 4F
→
02/23 21:32,
3年前
, 5F
02/23 21:32, 5F
→
02/23 21:32,
3年前
, 6F
02/23 21:32, 6F
我是在 Debian Linux 10.7.0 (64-bit) 裡面跑的，Python 3 和 gmpy2 都是
Debian 附的 package
我改了幾行，讓程式顯示 gmpy2、GMP、MPFR library 版本，
以及檢查 precision 是否超過最大 precision
至於那位不願具名友人說的 mpfr.digits() 不存在的問題，我就真的沒頭緒了
https://ideone.com/aJQW4R
※ 編輯: Schottky (111.250.31.177 臺灣), 02/23/2021 22:23:12
推
02/23 22:27,
3年前
, 7F
02/23 22:27, 7F
→
02/23 22:33,
3年前
, 8F
02/23 22:33, 8F