请问如下的代码还可以优化,提高速度吗
```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('pi.txt', '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())
这里附上源代码他人的链接 https://blog.csdn.net/ws995339251/article/details/85680710?ops_request_misc=&request_id=&biz_id=102&utm_term=Python%20%E9%AB%98%E7%B2%BE%E5%BA%A6%E8%AE%A1%E7%AE%97%E5%9C%86%E5%91%A8%E7%8E%87&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduweb~default-0-85680710.142%5Ev73%5Epc_new_rank,201%5Ev4%5Eadd_ask,239%5Ev2%5Einsert_chatgpt
最后,这是丘德诺夫斯基算法,如果大家有想法改进,欢迎提出宝贵的建议!
如果大家有更好的方法,可以不修改的基础上重新设计
或者各位谁知道y-cruncher的算法,可以分享一下吗,为圆周率日做准备
提速肯定可以,只是多少的问题。从代码来看,各种的代码都被封装了,这个要优化恐怕只有从被封装的核心代码开始了。
数学我不懂啊。。。。如果计算机内存足够,用高斯-勒让德算法如何?