sympy如何計算pi？

Question

sympy計算pi的數值背景是什麼？

我知道SymPy在後臺使用mpmath，這使得使用任意精度算法來執行計算成爲可能。這樣，一些特殊的常量，如e，pi，oo就被當作符號來處理，並且可以用任意的精度進行評估。

但是，Sympy如何確定任意小數位數？ Sympy如何在數字上做到這一點？

Answer 1

mpmath使用Chudnovsky公式計算pi（https://en.wikipedia.org/wiki/Chudnovsky_algorithm）。

Pi近似爲一個無限序列，其項減少爲（1/151931373056000）^ n，所以每個項增加大約14.18位數。這使得易於選擇多個術語N，從而達到期望的精度。

實際計算與整數運算來完成：即，爲了PREC位精度，* 2 ^（PREC）被計算pi的近似值。

下面是代碼，從mpmath/libmp/libelefun.py萃取（https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/libmp/libelefun.py）：

# Constants in Chudnovsky's series 
CHUD_A = MPZ(13591409) 
CHUD_B = MPZ(545140134) 
CHUD_C = MPZ(640320) 
CHUD_D = MPZ(12) 

def bs_chudnovsky(a, b, level, verbose): 
    """ 
    Computes the sum from a to b of the series in the Chudnovsky 
    formula. Returns g, p, q where p/q is the sum as an exact 
    fraction and g is a temporary value used to save work 
    for recursive calls. 
    """ 
    if b-a == 1: 
     g = MPZ((6*b-5)*(2*b-1)*(6*b-1)) 
     p = b**3 * CHUD_C**3 // 24 
     q = (-1)**b * g * (CHUD_A+CHUD_B*b) 
    else: 
     if verbose and level < 4: 
      print(" binary splitting", a, b) 
     mid = (a+b)//2 
     g1, p1, q1 = bs_chudnovsky(a, mid, level+1, verbose) 
     g2, p2, q2 = bs_chudnovsky(mid, b, level+1, verbose) 
     p = p1*p2 
     g = g1*g2 
     q = q1*p2 + q2*g1 
    return g, p, q 

@constant_memo 
def pi_fixed(prec, verbose=False, verbose_base=None): 
    """ 
    Compute floor(pi * 2**prec) as a big integer. 
    This is done using Chudnovsky's series (see comments in 
    libelefun.py for details). 
    """ 
    # The Chudnovsky series gives 14.18 digits per term 
    N = int(prec/3.3219280948/14.181647462 + 2) 
    if verbose: 
     print("binary splitting with N =", N) 
    g, p, q = bs_chudnovsky(0, N, 0, verbose) 
    sqrtC = isqrt_fast(CHUD_C<<(2*prec)) 
    v = p*CHUD_C*sqrtC//((q+CHUD_A*p)*CHUD_D) 
    return v 

這僅僅是普通的Python代碼，不同之處在於它依賴於一個額外的功能isqrt_fast()其計算平方根一個大整數。 MPZ是使用的大整數類型：gmpy.fmpz（如果可用），否則Python的內置long類型。裝飾器緩存計算值（在數值計算中通常需要重複使用pi，所以存儲它是有意義的）。

你可以看到它做一個基數轉換如下計算圓周率：

>>> pi_fixed(53) * 10**16 // 2**53 
mpz(31415926535897932) 

的關鍵技巧，使Chudnovsky式快速被稱爲二元分裂。無窮級數中的項滿足小系數的遞推關係（遞推方程可以在bs_chudnovsky函數的b-a == 1情況下看到）。不是按順序計算條件，而是將總和分成兩半;遞歸評估兩半，並將結果合併。最後，一個具有兩個大整數p和q，使得該系列的第一Ñ項之和恰好等於p/q。請注意，二進制拆分過程中不存在舍入誤差，這是該算法的一個非常好的功能;當計算平方根並在最後進行分割時，唯一的舍入發生。

大多數快速計算pi到高精度的程序都使用非常類似的策略，儘管有一些複雜的技巧可以進一步加速進程。

（注意：我是代碼的作者。）