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)