[心得] 計算π到小數點下十億位
上一篇的程式稍微改變一下可以用來計算圓周率π (根本改超多)
使用 Chudnovsky 公式
https://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
R(a,b)*P(a,b)
定義 S(a,b) = -------------
Q(a,b)
(6b)! (6a)!
R(a,b) = ----- ÷ ----- (第 a+1 項到第 b 項提出來的乘數)
(3b)! (3a)!
b!
Q(a,b) = (----)^3*(640320)^(3*(b-a)) (第 a+1 項到第 b 項的分母)
a!
P(a,b) = 第 a+1 項到第 b 項的 (13591409+545140134*k) 分子部份通分後的和
接著就和上一篇差不多，把級數對半切再對半切再對半切再對半切......
最後再一層一層通分合併起來即可。
☆ ☆ ☆
其實我還沒跑完十億位的計算，估計要一兩個小時 (被打)
用我的電腦計算一億位大約六分半鐘，消耗 1.2GB RAM
(註：後來跑完了，耗時 1 小時 28 分鐘，RAM 用量最大值 12GB)
以下是程式碼，沒給命令列參數的話預設一樣是十萬位數。
如果要計算十億位，命令列參數像這樣指定：
$ python3 pi.py 1000000000
注意要用 python3，假如出現語法錯誤 SyntaxError: invalid syntax 可能是使用到
python 2，可以用 python --version 確認版本
UNIX 系的作業系統有 time 指令可以測速
$ time python3 pi.py 1000000000
程式碼一樣放在 ideone.com 方便大家複製貼上，
但因為網站不提供 gmpy2 module 而無法執行。
# https://ideone.com/jAiHdM
#
# pi.py - Calculate Pi
#
import sys
import math
import gmpy2
from gmpy2 import mpfr
from gmpy2 import mpz
#
# Global Variables
#
count = 0
total = 0
old_p = 0
#
# 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("pi-py.txt", mode="w")
fd.write(" pi = ")
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, max):
global count
m = math.ceil((a + b) / 2)
if (b - a == 1):
if (a == 0):
r = 120 # 6!
q = mpz(640320**3)
p = gmpy2.sub( gmpy2.mul(q, 13591409),
gmpy2.mul(r, 13591409+545140134) )
else:
r = mpz(8 * (a*6+1) * (a*6+3) * (a*6+5))
q = mpz((b*640320)**3)
if ((b%2) == 0):
p = gmpy2.mul(mpz(13591409 + b*545140134), r)
else:
p = gmpy2.mul(mpz(-13591409 - b*545140134), r)
else:
p1, q1, r1 = s(a, m, max)
p2, q2, r2 = s(m, b, max)
# Merge
p = gmpy2.add( gmpy2.mul(p1, q2), gmpy2.mul(p2, r1) )
q = gmpy2.mul(q1, q2)
if (b != max):
r = gmpy2.mul(r1, r2)
else:
r = 0
count += 1
progress()
return p, q, r;
#
# Calculate e
#
def calc_pi(digits):
global total
d = digits+1
n_terms = math.ceil(d*math.log(10)/(3*math.log(53360)))
precision = math.ceil(d * math.log2(10)) + 4;
print("d = ", d, ", n = ", n_terms, ", precision = ", precision, sep="")
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, r = s(0, n_terms, n_terms)
q = gmpy2.mul(q, 426880)
print("Recursion done. Processing division.")
ef = gmpy2.div(q, p)
ef = gmpy2.mul(ef, gmpy2.sqrt(10005))
print("Division done. Converting to digits.")
estr, exp, prec = mpfr.digits(ef)
estr = estr[0:d]
print("Writing to file.")
write_string(estr);
#
# main program
#
argc = len(sys.argv)
if (argc >= 2):
digits = int(sys.argv[1])
else:
digits = 100000
calc_pi(digits)
# End of pi.py
--
桃樂絲: 可是, 如果你沒有頭腦, 為什麼會說話?
稻草人: ㄝ, 我也不知... 但是有些人沒有頭腦也能說超～多話呢。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 111.250.31.177 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Python/M.1614076108.A.9EC.html
※ 編輯: Schottky (111.250.31.177 臺灣), 02/23/2021 19:25:34