比gmpy2内置log(5)更快的agm写法 - l1t1/note GitHub Wiki

需要用到内置const_pi, exp()函数。

import gmpy2
from gmpy2 import mpfr, exp, log, const_pi
import time

def accurate_ln_agm(x, p=1000):
    bp = int(p * 3.322)
    ctx = gmpy2.get_context()
    ctx.precision = int(bp) + 200  # 增加缓冲
    
    x = mpfr(x)
    s = gmpy2.ceil(p/2 * 2.31)
    
    def agm(a, b, tol=None):
        if tol is None:
            tol = mpfr(10)**(-bp-20)  # 更严格的容差
        for _ in range(100):  # 5百万位是43次
            a, b = (a + b)/2, gmpy2.sqrt(a * b)
            if abs(a - b) < tol:
                print(_, float(abs(a - b)));break
        return a #(a + b)/2
    
    bigx = exp(s)*x

    agm_val = agm(1, 4/bigx)
    result = const_pi() / (2 * agm_val) - s
    
    return mpfr(result, bp)
t=time.time()
ln5 = accurate_ln_agm(5, 1.6*10**6)
print("ln_agm", time.time()-t)
print(float(ln5),len(str(ln5)))
t=time.time()
exact = gmpy2.log(5)
print("gmpy2.log(5)", time.time()-t)
print(float(exact),len(str(exact)))
print(f"绝对误差: {abs(ln5 - exact):.3e}") 

输出

ln_agm 4.072713375091553
1.6094379124341003 1600036
gmpy2.log(5) 4.5735554695129395
1.6094379124341003 1600036
绝对误差: 1.994e-1600035