# Lehmer's gcd algorithm;  revised version

DIGIT_BITS = 30
BASE = 1 << DIGIT_BITS

def nbits(n, correction = {
        '0': 4, '1': 3, '2': 2, '3': 2,
        '4': 1, '5': 1, '6': 1, '7': 1,
        '8': 0, '9': 0, 'a': 0, 'b': 0,
        'c': 0, 'd': 0, 'e': 0, 'f': 0}):
    """Number of bits in binary representation of the positive integer n,
    or 0 if n == 0.
    """
    if n < 0:
        raise ValueError("The argument to _nbits should be nonnegative.")
    hex_n = "%x" % n
    return 4*len(hex_n) - correction[hex_n[0]]

def lehmer_gcd(a, b):
    """Greatest common divisor of nonnegative integers a and b."""

    # initial reductions
    if a < b:
        a, b = b, a

    while b >= BASE:
        size = nbits(a) - DIGIT_BITS
        x, y = int(a >> size), int(b >> size)
        # single-precision arithmetic from here...
        A, B, C, D = 1, 0, 0, 1
        while True:
            if y+C == 0 or y+D == 0:
                break
            q = (x+A)//(y+C)
            if q != (x+B)//(y+D):
                break
            A, B, x, C, D, y = C, D, y, A-q*C, B-q*D, x-q*y
        # ...until here

        if B:
            a, b = A*a + B*b, C*a + D*b
        else:
            a, b = b, a % b

    a, b = int(b), int(a%b)
    # final single-precision gcd computation
    while b:
        a, b = b, a%b
    return a