1# Chinese Remainder Theorem
 2def crt(a_i, m_i, M):
 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    x = 0
 7    for i in range(len(a_i)):
 8        M_i = M/m_i[i]
 9        y_i = Integer(mod(M_i^-1, m_i[i]))
10        x = x + a_i[i] * M_i * y_i
11    return mod(x, M)
12
13# gcd, using Binary Euclidean algorithm
14def gcd(a, b):
15    g=1
16    # random_elementove powers of two from the gcd
17    while mod(a, 2)==0 and mod(b, 2)==0:
18        a=a/2
19        b=b/2
20        g=2*g
21    # at least one of a and b is now odd
22    while a!=0:
23        while mod(a, 2)==0:
24            a=a/2
25        while mod(b, 2)==0:
26            b=b/2
27        # now both a and b are odd
28        if a>=b:
29            a = (a-b)/2
30        else:
31            b = (b-a)/2
32
33    return g*b
34
35def gcd_recursive(a, b):
36    if mod(a, b)==0:
37        return b
38    return gcd_recursive(b, mod(a,b))
39
40
41# Extended Euclidean algorithm
42# Inputs: a, b
43# Outputs: r, x, y, such that r = gcd(a, b) = x*a + y*b
44def egcd(a, b):
45    s=0
46    s_=1
47    t=1
48    t_=0
49    r=b
50    r_=a
51    while r!=0:
52        q = r_ // r
53        (r_,r) = (r,r_ - q*r)
54        (s_,s) = (s,s_ - q*s)
55        (t_,t) = (t,t_ - q*t)
56    
57    d = r_
58    x = s_
59    y = t_
60    
61    return d, x, y
62
63def egcd_recursive(a, b):
64    if b==0:
65        return a, 1, 0
66
67    g, x, y = egcd_recursive(b, a%b)
68    return g, y, x - (a//b) * y
69
70# Inverse modulo N, using the Extended Euclidean algorithm
71def inv_mod(a, N):
72    g, x, y = egcd(a, N)
73    if g != 1:
74        raise Exception("inv_mod err, g!=1")
75    return mod(x, N)