用其他高阶迭代法达到半秒计算根号二小数点后百万位 - l1t1/note GitHub Wiki

在调试泰勒法计算根号二时发现一个有意思的现象：用10进制整数幂做BASE，不如用2进制整数幂，即使要加上最后返回时x*10进制整数幂//2进制整数幂,也是利大于弊，这在x86上更明显。

在阅读论文时发现还有很不错的牛顿迭代法的优化版本，DeepSeek提供了一种二次收敛法，我自己找到了一种halley法，它们与泰勒法相比，共同的特点是：迭代次数非常少，从几千、万余次（取决于初始佩尔值）降低到几十次，相关初始值要求也不高，随便给个0.7都行，但每次迭代计算量很大，耗时很长。采用2进制整数幂做BASE后，两种方法分别缩短了三分之二和一半时间，AMD机器上达到0.54和0.32s。现在单进程速度已经超出了预期，也超过了煞费苦心调整并行块大小的泰勒法，数学才是王道，在它面前，一切程序逻辑优化都是微不足道的，当然是在GMP库提供的高效高精度计算库的基础上，只不过没有用库函数gmpy2.root(mpfr(2), 2)本身。

DeepSeek的贡献在于，把原来的浮点乘除改写成基于BASE的整数乘除。但它提供的二次收敛法，迭代结果是负二分之根号二，我调换初始值的分子分母后就变成正的二分之根号二。

直接上最后优化的代码，值得一提的是，DeepSeek还对二进制乘除做了移位优化，AMD上的实验证明并没有提高效率，估计GMP库自动转化了。移位操作还特别容易引起优先级问题，要加很多层括号，可读性差，最后没有采纳（在ARM CPU上这种改写还是值得的）。

def halley_integer_correct2(digits):
    #bs=int(digits/0.301)
    BASE = mpz(2)**int(digits/0.301)#mpz(10)**digits
    #if BASE >mpz(10)**(digits+10):print("BASE >mpz(10)**digits",len(str(BASE)), len(str(mpz(10)**digits)))
    # 初始值使用更高精度的佩尔解避免首次迭代失效
    m, n = mpz(9369319), mpz(6625109)
    x = (m * BASE) // n  # 初始值≈1.414213562
    #x=x*(1+2*(2-x*x)/(3*x*x+2))
    for _ in range(int(mpz(digits).bit_length() // 2 + 2)):  # 收敛更快，减少迭代次数
        oldx=x 
        x_sq = (x * x) // BASE
        numerator = 2 * (2 * BASE - x_sq)  # 2*(2 - x^2)
        denominator = 3 * x_sq + 2 * BASE  # 3x^2 + 2
        correction = (BASE + (numerator * BASE // denominator) ) 
        x = (x * correction) // BASE

    return (x*mpz(10)**digits)//BASE 
    
def sqrt2_by_quadratic_opt(digits):
    t=time.time()
    # 初始化佩尔方程解 (m,n)
    m = mpz(9369319)
    n = mpz(6625109)
    
    # 设置计算精度 (BASE = 10^digits)
    bin_shift = int(digits / 0.301)+1   # 安全余量
    BASE = mpz(2)**bin_shift
    DEC_BASE = mpz(10)**digits
    x = (n * BASE) // m  # 初始近似值
    #print(x)
    # 二次收敛迭代
    for _ in range(int(mpz(digits).bit_length() + 1)):
        #x_sq = (x * x) // BASE
        #x = (x * (3 * BASE - 2 * x_sq)) // (2 * BASE)
        x = (x * (BASE*3 - 2*x**2//BASE)) // (2*BASE)  
    print("二次收敛opt", time.time()-t, "iter", int(mpz(digits).bit_length() + 1))
    return (x* DEC_BASE) >> (bin_shift-1)#x*mpz(2)