大家好，這是我第一次寫 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://webptt.com/m.aspx?n=bbs/Python/M.1614067863.A.B21.html 我是在 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