为丘德诺夫斯基算法提速

请问如下的代码还可以优化,提高速度吗


```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的算法,可以分享一下吗,为圆周率日做准备

提速肯定可以,只是多少的问题。从代码来看,各种的代码都被封装了,这个要优化恐怕只有从被封装的核心代码开始了。

数学我不懂啊。。。。如果计算机内存足够,用高斯-勒让德算法如何?