1# Chinese Remainder Theorem
  2def crt(a_i, m_i):
  3    if len(a_i)!=len(m_i):
  4        raise Exception("error, a_i and m_i must be of the same length")
  5
  6    M=1
  7    for i in range(len(m_i)):
  8        M = M * m_i[i]
  9
 10    x = 0
 11    for i in range(len(a_i)):
 12        M_i = M/m_i[i]
 13        y_i = Integer(mod(M_i^-1, m_i[i]))
 14        x = x + a_i[i] * M_i * y_i
 15    return mod(x, M)
 16
 17# gcd, using Binary Euclidean algorithm
 18def gcd(a, b):
 19    g=1
 20    # random_elementove powers of two from the gcd
 21    while mod(a, 2)==0 and mod(b, 2)==0:
 22        a=a/2
 23        b=b/2
 24        g=2*g
 25    # at least one of a and b is now odd
 26    while a!=0:
 27        while mod(a, 2)==0:
 28            a=a/2
 29        while mod(b, 2)==0:
 30            b=b/2
 31        # now both a and b are odd
 32        if a>=b:
 33            a = (a-b)/2
 34        else:
 35            b = (b-a)/2
 36
 37    return g*b
 38
 39def gcd_recursive(a, b):
 40    if mod(a, b)==0:
 41        return b
 42    return gcd_recursive(b, mod(a,b))
 43
 44
 45# Extended Euclidean algorithm
 46# Inputs: a, b
 47# Outputs: r, x, y, such that r = gcd(a, b) = x*a + y*b
 48def egcd(a, b):
 49    s=0
 50    s_=1
 51    t=1
 52    t_=0
 53    r=b
 54    r_=a
 55    while r!=0:
 56        q = r_ // r
 57        (r_,r) = (r,r_ - q*r)
 58        (s_,s) = (s,s_ - q*s)
 59        (t_,t) = (t,t_ - q*t)
 60    
 61    d = r_
 62    x = s_
 63    y = t_
 64    
 65    return d, x, y
 66
 67def egcd_recursive(a, b):
 68    if b==0:
 69        return a, 1, 0
 70
 71    g, x, y = egcd_recursive(b, a%b)
 72    return g, y, x - (a//b) * y
 73
 74# Inverse modulo N, using the Extended Euclidean algorithm
 75def inv_mod(a, N):
 76    g, x, y = egcd(a, N)
 77    if g != 1:
 78        raise Exception("inv_mod err, g!=1")
 79    return mod(x, N)
 80
 81# Tests
 82
 83#####
 84# Chinese Remainder Theorem tests
 85a_i = [5, 3, 10]
 86m_i = [7, 11, 13]
 87assert crt(a_i, m_i) == 894
 88
 89a_i = [3, 8]
 90m_i = [13, 17]
 91assert crt(a_i, m_i) == 42
 92
 93
 94#####
 95# gcd, using Binary Euclidean algorithm tests
 96assert gcd(21, 12) == 3
 97assert gcd(1_426_668_559_730, 810_653_094_756) == 1_417_082
 98
 99assert gcd_recursive(21, 12) == 3
100
101#####
102# Extended Euclidean algorithm tests
103assert egcd(7, 19) == (1, -8, 3)
104assert egcd_recursive(7, 19) == (1, -8, 3)
105
106#####
107# Inverse modulo N tests
108assert inv_mod(7, 19) == 11