Python implementation of (extended) Lehmer's GCD algorithm

lgcd.py

import math, random

def lgcd(x, y):
    if x < y:
        return lgcd(y, x)

    shift = max(x.bit_length() // 64, y.bit_length() // 64)
    xbar = x >> (shift * 64)
    ybar = y >> (shift * 64)
    while y > 2**64:
        print("oprand", x, y)
        a, b, c, d = 1, 0, 0, 1
        while ybar + c != 0 and ybar + d != 0:
            q = (xbar + a) // (ybar + c)
            p = (xbar + b) // (ybar + d)
            if q != p:
                break

            a, c = c, a - q*c
            b, d = d, b - q*d
            xbar, ybar = ybar, xbar - q*ybar

        if b == 0:
            x, y = y, x % y
        else:
            print("reduced", a, b, c, d)
            x, y = a*x + b*y, c*x + d*y

    return math.gcd(x, y)

def xgcd(x, y):
    xneg, yneg = -1 if x < 0 else 1, -1 if y < 0 else 1
    x, y = abs(x), abs(y)

    # it's maintained that r = s * x + t * y, last_r = last_s * x + last_t * y
    last_r, r = x, y
    last_s, s = 1, 0
    last_t, t = 0, 1

    while r > 0:
        q = last_r // r
        last_r, r = r, last_r - q*r
        last_s, s = s, last_s - q*s
        last_t, t = t, last_t - q*t

    return last_r, last_s * xneg, last_t * yneg

def lxgcd(x, y):
    if x < y:
        g, cy, cx = xgcd(y, x)
        return g, cx, cy

    shift = max(x.bit_length() // 64, y.bit_length() // 64)
    xbar = x >> (shift * 64)
    ybar = y >> (shift * 64)

    last_s, s = 1, 0
    last_t, t = 0, 1

    while y > 2**64:
        print("oprand", x, y)
        a, b, c, d = 1, 0, 0, 1
        while ybar + c != 0 and ybar + d != 0:
            q = (xbar + a) // (ybar + c)
            p = (xbar + b) // (ybar + d)
            if q != p:
                break

            a, c = c, a - q*c
            b, d = d, b - q*d
            xbar, ybar = ybar, xbar - q*ybar

        if b == 0:
            q = x // y
            x, y = y, x % y
            last_s, s = s, last_s - q*s
            last_t, t = t, last_t - q*t
        else:
            print("reduced", a, b, c, d)
            x, y = a*x + b*y, c*x + d*y
            last_s, s = a*last_s + b*s, c*last_s + d*s
            last_t, t = a*last_t + b*t, c*last_t + d*t

    # x = last_s * X + last_t * Y
    # y = s * X + t * Y
    # notice that here x, y could be negative
    g, cx, cy = xgcd(x, y)

    # g = cx * x + cy * y = (cx * last_s + cy * s) * X + (cx * last_t + cy * t) * X
    return g, cx * last_s + cy * s, cx * last_t + cy * t

# some tests
x = random.randrange(0, 2**190)
y = random.randrange(0, 2**190)
g = lgcd(x, y)
assert(g == math.gcd(x, y))
print("gcd =", g)

assert(xgcd(x, y)[0] == math.gcd(x, y))

g, cx, cy = lxgcd(x, y)
assert(g == math.gcd(x, y))
assert(cx*x + cy*y == g)
print("gcd =", g)