大家好，这是我第一次写 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/cn.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