三个计算ln(2)的python程序 - l1t1/note GitHub Wiki

1.利用AGM算法

这是在博客园：计算自然对数的快速算法基础上，让DeepSeek译成利用gmpy2库的Python, 并且将原文预计算常量π、ln(10), 改写为用DeepSeek给出的Brent-Salamin方法（Fast Multiple-Precision Evaluation of Elementary Functions RICHARD P. BRENT）AGM计算π函数、前面介绍过的二进制分裂+阶乘逆计算e函数，最终计算ln(2)的算法出自（T. Sasaki and Y. Kanada, Practically fast multiple-precision evaluation of log(x), J. Inf. Process. 5 (1982), 247–250.），从而实现不依赖高精度库内置常量的程序。因为是拼凑的，所以格式有点乱。

import gmpy2
from gmpy2 import mpz, mpfr, get_context
import time 
import math
def agm_pi(prec: int) -> mpfr:
    """完全修正的AGM计算π实现"""
    #ctx = gmpy2.get_context()
    #ctx.precision = int(1.1 * prec)  # 额外10%精度缓冲
    
    # 初始化
    a = mpfr(1)
    b = 1 / gmpy2.sqrt(mpfr(2))
    t = mpfr(0.25)
    p = mpfr(1)
    
    # 动态计算阈值
    threshold = mpfr(2)**(-int(prec * 0.9))  # 关键修正
    
    for _ in range(1000000):#int(log2(prec)) + 5):
        a_new = (a + b) / 2
        t -= p * (a - a_new)**2
        a, b, p = a_new,  gmpy2.sqrt(a*b), 2*p
        
        # 动态收敛判断
        if abs(a - b) < threshold:
            break
    
    pi_approx = (a + b)**2 / (4 * t)
    return pi_approx


def calculate_e_terms(l, r):
    if l == r:
        return mpz(1), mpz(l)
    mid = (l + r) // 2
    p1, q1 = calculate_e_terms(l, mid)
    p2, q2 = calculate_e_terms(mid + 1, r)
    return p1 * q2 + p2, q1 * q2

def calculate_e_fast(n, target_digits):
    """优化版：纯整数运算生成小数部分"""
    # 计算所需二进制精度
    ln_n_fact = float(n * math.log(n) - n + 0.5*math.log(2*math.pi*n))
    prec_bits = int((ln_n_fact + math.log(n)) / math.log(2)) + 128
    get_context().precision = prec_bits

    p, q = calculate_e_terms(2, n)
    power = mpz(10)**target_digits
    decimal_part = (p * power) // q  # 关键优化点
    
    return mpfr("2."+str(decimal_part))

def find_required_terms(digits):
    """计算所需最小n（使用斯特林公式）"""
    ln10 = math.log(10)
    n = 1
    while True:
        ln_n_fact = n * math.log(n) - n + 0.5*math.log(2*math.pi*n)
        if (ln_n_fact - n * ln10) > digits * ln10:
            return n
        n += 1
# 初始化计算环境
def init_ln(precision=1000):
    """初始化高精度对数计算环境"""
    # 设置计算精度（比特位数 = 十进制位数 * log2(10) + 安全余量）
    get_context().precision = int(precision * 3.33)
    prec= int(precision * 3.33)
    # 定义全局常量
    global pi2, ln10, eps2, eps1, e
    pi2 = agm_pi(prec)*2#gmpy2.const_pi() * 2
    ln10 = gmpy2.log(10)
    #e=gmpy2.exp(1)
    t=time.time()
    #e = (1 + 1/mpfr(2)**(prec//2)) ** mpfr(2)**(prec//2)
    e=calculate_e_fast(find_required_terms(precision),precision )
    get_context().precision = int(precision * 3.33)

    #error=e-gmpy2.exp(1)
    #print(f"  误差: {error:.3e}")
    print(float(e), time.time()-t, "s")
    #exit()
    eps2 = mpfr(10) ** (-precision)
    eps1 = eps2 * eps2

def log(x):
    """计算自然对数ln(x)"""
    if x <= 0:
        raise ValueError("Must be positive")
    if x == 1:
        return mpfr(0)
    
    x = mpfr(x)
    k = 0

    
    # 范围调整（将x规范到(0.01, 0.1]区间）
    while x > 0.1:
        x /= e
        k += 1
    while x <= 0.01:
        x *= e
        k -= 1
    return k - negative_log(x)#k * ln10 - negative_log(x)

def negative_log(q):
    """辅助函数：计算(0.01, 0.1]区间内的对数"""
    q = mpfr(q)
    r = q
    s = q
    n = q
    q2 = q * q
    q1 = q2 * q
    
    p = True
    while (n * q1) > eps1:
        n *= q1
        s += n
        q1 *= q2
        r += n if (p := not p) else -n
    
    u = 1 - 2 * r
    v = 1 + 2 * s
    t = u / v
    
    a = mpfr(1)
    b = gmpy2.sqrt(1 - t * t * t * t)
    
    while (a - b) > eps2:
        t = (a + b) / 2
        b = gmpy2.sqrt(a * b)
        a = t
    
    return pi2 / (a + b) / v / v

def log10(x):
    """计算常用对数log10(x)"""
    return log(x) / ln10

# 使用示例
if __name__ == "__main__":
    # 初始化（精度设为1000位）
    init_ln(100000)
    
    # 测试计算
    num = mpfr(2)#"2.718281828459045")
    result=log(num)
    print(str(result)[:100002])
    #print(f"ln({num}) = {len(str(result))}")
    #print(f"log10(100) = {log10(100)}")

    #exact = gmpy2.log(mpfr(2))
    #print(f"与内置函数差异: {abs(result - exact):.3e}")    
    # 验证1000位精度
    #high_prec_num = mpfr("3.141592653589793238462643383279")
    #print(f"ln(π) 前50位: {str(log(high_prec_num))[:50]}...")

2.利用Zuniga-II二进制分裂法

本程序和下一个程序均来自论文Fast Log() calculation for Arbitrary Precision number, 由DeepSeek一次性翻译完成，除了因为这个高精度库函数的精度采取的是二进制长度，做了*3.33的处理，其余一字未改。

import gmpy2
from gmpy2 import mpz, mpfr, get_context
import math
from functools import lru_cache

def compute_ln2_zuniga(precision):
    """基于Zuniga-II算法的ln(2)高精度计算"""
    # 1. 参数初始化
    ctx = get_context()
    kmax = int(math.ceil(precision * math.log(10) / math.log(3888))) + 10  # 安全增量
    workprec = int(precision + 2 * math.log(kmax + 1))  # 增强工作精度
    
    # 2. 设置GMP上下文
    ctx.precision = int(workprec ) + 256  # 二进制精度+额外缓冲

    # 3. 定义二进制分裂递归函数
    @lru_cache(maxsize=None)
    def binary_split(a, b):
        if a + 1 == b:
            b6 = 6 * b
            bb = mpz(b)
            
            # 计算三项式系数
            q = mpz(216) * mpz(b6 - 1) * mpz(b6 - 5)
            r = mpz(b) * mpz(2 * b - 1)
            p = mpz(1794) * bb - mpz(297)
            return p, q, r
        
        # 区间分裂
        mid = (a + b) // 2
        pL, qL, rL = binary_split(a, mid)
        pR, qR, rR = binary_split(mid, b)
        
        # 合并结果 (p = pL*qR + pR*rL)
        return pL*qR + pR*rL, qL*qR, rL*rR

    # 4. 执行主计算
    p, q, _ = binary_split(0, kmax)
    q *= mpz(2)  # 最终修正 q *= 2

    # 5. 高精度转换
    ctx.precision = int(precision * 1.5)  # 最终精度缓冲
    result = mpfr(p) / mpfr(q)
    
    # 精度锁定
    ctx.precision = int(precision * 3.33)
    return mpfr(result, precision)

# 验证测试
if __name__ == "__main__":
    precision = int(100000*3.33)  # 测试1000位精度
    ln2 = compute_ln2_zuniga(precision)
    exact = gmpy2.log(2)
    
    print(f"计算结果[前50位]: {str(ln2)[:50]}")
    print(f"标准值差异: {abs(ln2 - exact):.3e}")
    print(f"有效位数: {-math.log10(abs(float(ln2 - exact))):.1f}")

3.利用Xiao二进制分裂法

import gmpy2
from gmpy2 import mpz, mpfr
import math

def compute_ln2_xiao(precision):
    """基于Xiao的二进制分裂算法计算ln(2)的gmpy2实现"""
    # 1. 计算所需项数kmax和工作精度
    kmax = int(math.ceil(precision * math.log(10) / math.log(450000)))
    workprec = int(math.ceil(precision + 1 + 1 * math.log(kmax)))
    
    # 2. 初始化gmpy2变量
    ctx = gmpy2.get_context()
    ctx.precision = int(precision)  # 二进制位数 ≈ 十进制位数*log2(10)
    
    # 3. 定义递归的二进制分裂函数
    def binary_splitting(a, b):
        if a + 1 == b:
            b4 = 4 * b
            b10 = 10 * b
            bb = mpz(b)
            
            # 计算q = 1440*(10b-1)(10b-3)(10b-7)(10b-9)
            q = mpz(1440) * mpz((b10-1)*(b10-3)) * mpz((b10-7)*(b10-9))
            
            # 计算r = b(2b-1)(4b-1)(4b-3)
            r = mpz(b*(2*b-1)) * mpz((b4 - 1) * (b4 - 3))
            
            # 计算p = (3927264b³ - 4300512b² + 1209726b - 81891)
            p = ((mpz(3927264) * bb - mpz(4300512)) * bb + mpz(1209726)) * bb - mpz(81891)
            return p, q, r
        
        # 分裂区间[a,b] → [a,mid]和[mid,b]
        mid = (a + b) // 2
        p_left, q_left, r_left = binary_splitting(a, mid)
        p_right, q_right, r_right = binary_splitting(mid, b)
        
        # 合并结果: p = p_left*q_right + p_right*r_left
        p = p_left * q_right + p_right * r_left
        q = q_left * q_right
        r = r_left * r_right
        return p, q, r
    
    # 4. 执行主计算
    p, q, r = binary_splitting(0, kmax)
    q *= mpz(4)  # q *= 4
    
    # 5. 转换为高精度浮点数
    fp = mpfr(p, workprec)
    fq = mpfr(q, workprec)
    result = fp / fq
    
    # 调整最终精度
    ctx.precision = int(precision)
    return mpfr(result, precision)

# 测试代码
if __name__ == "__main__":
    precision = int(100000 *3.33) # 100位十进制精度
    ln2 = compute_ln2_xiao(precision)
    print(f"ln(2) [前{precision}位]: {len(str(ln2))}") #[:30]}")#precision+2]}")
    print(f"与gmpy2内置函数差异: {abs(ln2 - gmpy2.log(2)):.3e}")

三个程序计算10万位ln(2)均低于1秒。