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)