gcd算法实现

[HJimmyK]

实现任意精度的gcd算法，对不同输入长度下的数都有良好的性能。

[HJimmyK]

Lehmer算法已被实现，在绝大多数测试中，Lehmer算法比basecase快4-5左右

[HJimmyK]

这是lehmer算法的python代码，实际中的实现基本一致。

def lehmer_gcd(a, b):
    """
    使用 Lehmer 算法计算两个非负大整数 a 和 b 的最大公约数。
    要求 a >= 0, b >= 0，并且不同时为 0。
    """
    if a == 0:
        return b
    if b == 0:
        return a

    # 确保 a >= b
    if a < b:
        a, b = b, a

    while b != 0:
        # 如果 a 和 b 都很小，直接用欧几里得算法收尾
        if a.bit_length() <= 64 or b.bit_length() <= 64:
            while b:
                a, b = b, a % b
            break

        # 取 a 和 b 的高位（例如前 64 位）
        k = a.bit_length() - 64
        if k < 0:
            k = 0
        x = a >> k
        y = b >> k

        # 初始化矩阵系数
        A, B, C, D = 1, 0, 0, 1

        while True:
            if y == 0:
                break
            q = x // y
            
            # 更新矩阵
            A, B, C, D = C, D, A - q * C, B - q * D
            x, y = y, x - q * y
            # 如果新 y 太小，停止
            if y < (1 << 32):  # 保证下一次 q 还可靠
                break

        if C == 0:
            # 退化为普通欧几里得一步
            a, b = b, a % b
        else:
            # 应用变换
            new_a = A * a + B * b
            new_b = C * a + D * b
            a, b = new_a, new_b

        # 保持 a >= b >= 0
        if a < b:
            a, b = b, a
        if b < 0:
            b = -b

    return a

import math
import random
# 测试
random.seed(12345)
for i in range(10000):
    c = random.getrandbits(608)
    a = random.getrandbits(1601)*c
    b = random.getrandbits(1607)*c
    r1 = lehmer_gcd(a, b)
    r2 = math.gcd(a, b)
    if(r1 != r2):
        print("error")