Extended Euclidean algorithm

[Wiki]In arithmetic and computer programming, the extended Euclidean algorithm is an extension to the Euclidean algorithm, and computes, in addition to the greatest common divisor of integers a and b, also the coefficients of Bézout's identity, which are integers x and y such that ax+by=\gcd(a,b)

通常用于求modular multiplicative inverse，如17 x ≡ 1 (mod 43)，意思是存在s，t，使得43s + 17t = 1，求x

思路：c1, c2为p，a的系数，照着图中迭代

import math
import time


def gf_inv(a, p):
    r1 = p
    r2 = a
    c1 = 0
    c2 = 1

    while r1 > 1:
        q = r1 // r2
        r = r1 % r2
        c = c1 - c2 * q
        r1 = r2
        c1 = c2
        r2 = r
        c2 = c
    return c1 % p


def egcd(a, p):
    if a == 0:
        return (p, 0, 1)
    else:
        q = p // a
        r = p % a
        g, c1, c2 = egcd(r, a)
        return (g, c2 - q * c1, c1)


def mod_inv(a, p):
    g, x, y = egcd(a, p)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % p


def measure_time(func, *args):
    a, p = args
    start = time.time()
    if func.__name__ == 'pow':
        res = func(a, p - 2, p)
    else:
        res = func(a, p)
    end = time.time() - start
    print('{}: {}, runtime: {:<10f}'.format(func.__name__, res, end))


# 17 x ≡ 1 (mod 43)
# (3 ^ 6) x ≡ 1 (mod 31)
a_list = [17, 3 ** 6, math.ceil(math.sqrt(7919))]
p_list = [43, 31, 7919]
func = [gf_inv, mod_inv, pow]
for f in func:
    for i in range(len(a_list)):
        measure_time(f, a_list[i], p_list[i])

# runtime: iteration > recursion > built-in pow
# gf_inv: 38, runtime: 0.000011  
# gf_inv: 2, runtime: 0.000006  
# gf_inv: 4004, runtime: 0.000003  
# mod_inv: 38, runtime: 0.000010  
# mod_inv: 2, runtime: 0.000004  
# mod_inv: 4004, runtime: 0.000003  
# pow: 38, runtime: 0.000009  
# pow: 2, runtime: 0.000001  
# pow: 4004, runtime: 0.000002