Python快速计算高精度圆周率
################################################################
# Computing pi by Binary Splitting Algorithm with GMP library. #
################################################################
# 经测试 计算三千两百万为PI所需时间如下
# CPU 计算时间 文件流写入时间
# Intel(R) Core(TM) i7-4720HQ CPU @ 3.40GHz 363.33s 30.08s
# Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz 343.81s 27.82s
# Intel(R) Xeon(R) Gold 6148 CPU @ 2.40GHz 278.65s 8.54s
import gmpy2
import numpy as np
import time
# 在gmpy2中，mpz高精度整形， mpfr为高精度浮点
# 首先封装PQT类，定义P、Q、T为mpz类型数据
class PQT:
def __init__(self):
self.P = gmpy2.mpz(0)
self.Q = gmpy2.mpz(0)
self.T = gmpy2.mpz(0)
# 定义Chudnovsky类
class Chudnovsky:
# 初始化变量
def __init__(self, digits):
self.DIGITS = digits # 要计算的位数
self.A = gmpy2.mpz(13591409) # 常数A，mpz
self.B = gmpy2.mpz(545140134) # 常数B，mpz
self.C = gmpy2.mpz(640320) # 常数C，mpz
self.D = gmpy2.mpz(426880) # 常数D，mpz
self.E = gmpy2.mpz(10005) # 常数E，mpz
self.aa = gmpy2.mpz(53360) # 常数a，用来计算DIGITS_PER_TERM
self.aa *= 53360 # mpz类型与int类型加减乘后依然为mpz类型
self.aa *= 53360
self.aa = gmpy2.log(self.aa) # 常数aa = log(53360^3)
# 常数DIGITS_PER_TERM，用来计算递归的次数
self.DIGITS_PER_TERM = self.aa / gmpy2.log(10)
# 常数C3_24(C^3/24)，mpz。C为mpz但是与int做除法后，会自动变为mpfr类型，所以需要强制类型转换为mpz
self.C3_24 = gmpy2.mpz(self.C * self.C * self.C / 24)
# N,mpz。递归计算PQT时只用计算到int(DIGITS/DIGITS_PER_TERM)为止即可
self.N = int(self.DIGITS / self.DIGITS_PER_TERM)
# 最后输出pi的精度。在gmpy2中，精度的单位是bit，所以十进制输出每一位需要乘log2(10)
self.PREC = int(self.DIGITS * np.log2(10))
# 构建一棵二叉树，递归遍历每个节点，复杂度为O(n)。
# 目的是得到P(0,N)、Q(0,N)、T(0,N)，用于直接求得pi
# P、Q、T定义见PDF
# p_k = (2*k - 1)(6*k - 1)(6*k - 1)
# q_k = (k^3) * (C^3) / 24
# a_k = (-1)^k * (A + B*k)
def compPQT(self, n1, n2):
pqt = PQT() # PQT对象
if n1 + 1 == n2: # 递归的出口
pqt.P = 2 * n2 - 1 # P(n2-1, n2) = a_n2
pqt.P *= 6 * n2 - 1
pqt.P *= 6 * n2 - 5
pqt.Q = self.C3_24 * n2 * n2 * n2 # Q(n2-1, n2) = q_n2
pqt.T = (self.A + self.B * n2) * pqt.P # T(n2-1, n2) = a_n2 * p_n2
if (n2 & 1) == 1: # (-1)^n2
pqt.T = - pqt.T
else:
m = int((n1 + n2) / 2) ############
pqt1 = self.compPQT(n1, m) ## 二叉树 ##
pqt2 = self.compPQT(m, n2) ############
pqt.P = pqt1.P * pqt2.P # P(n1, n2) = P(n1, m) * P(m, n2)
pqt.Q = pqt1.Q * pqt2.Q # Q(n1, n2) = Q(n1, m) * Q(m, n2)
pqt.T = pqt1.T * pqt2.Q + pqt1.P * pqt2.T # T(n1, n2) = T(n1, m) * Q(m, n2) + P(n1, m) * T(m, n2)
return pqt
#####################################################
# _____ Q(0, N) #
# pi = 426880 · √10005 · ——————————————————————— #
# A·Q(0, N) + T(0, N) #
#####################################################
def compPi(self):
print("**** PI Computation ( " + str(self.DIGITS) + " digits )")
t0 = time.time()
pqt = self.compPQT(0, self.N) # 计算(0, N)的PQT
gmpy2.get_context().precision = self.PREC # 设置mpfr的精度为PREC
pi = self.D * gmpy2.sqrt(gmpy2.mpfr(self.E)) * pqt.Q # 计算pi的分子
pi /= (self.A * pqt.Q + pqt.T) # 除以分母得到pi
t1 = time.time()
print('TIME (COMPUTE):' + str(t1-t0))
file = open('', 'w') # 将pi写入文件保存
file.write(str(pi))
t2 = time.time()
print('TIME (WRITE):' + str(t2-t1))
print('LAST 100 DIGITS:')
return str(pi)[-100:]
digit = int(input('Input digits:'))
main = Chudnovsky(digit)
print(main.compPi())