In [8]:
# Extended Euclidean algorithm (builtin)
# xgcd(a,b) returns a triple, where the first number is gcd, and the
# next two numbers are coefficients c,d such that ac + bd =1  
xgcd(3,5)
Out[8]:
(1, 2, -1)
In [9]:
3*2 + 5*(-1)
Out[9]:
1
In [10]:
# modular exponentiation is important for public key cryptography
(2^5) % 19
Out[10]:
13
In [11]:
(17^100000) % 19 # this works, but if 100000 is replaced by a huge number, an error will occur
Out[11]:
17
In [56]:
def exp_mod_slow(a,e,b):
    """Return (a^e) % b, reducing at every stage."""
    prod = 1
    for i in range(e):
        prod = (prod*a) % b  # reduce at every stage, rather than just at end
    return prod

# runtime (for fixed a,b) is about e, which exponential in number of digits of e
In [44]:
exp_mod_slow(17,10000,19)
Out[44]:
17
In [13]:
# In-class Exercise: compute a^(2^n) mod b
In [14]:
def exp_mod_powtwo(a,n,b):
    """Returns a^(2^n) % b, using the repeated squaring method."""
    prod = a
    for i in range(n):
        prod = (prod)^2 % b
    return prod

#runtime (for fixed a,b) is about n (exp_mod_slow would take time 2^n)
In [21]:
exp_mod_powtwo(3,1000,23)  # fast even to compute 3^(2^1000) mod 23
Out[21]:
3
In [54]:
bin(235)  # builtin function outputs binary expansion; ignore first two characters of output
Out[54]:
'0b11101011'
In [55]:
# kth bit of binary expansion of 235
k=1
bin(235)[-1-k]
Out[55]:
'1'
In [49]:
length = len(bin(235))-2
bin(235)[-1-(length-1)]  # gets the highest bit in binary expansion
Out[49]:
'1'
In [45]:
def exp_mod(a,e,b):
    """Return a^e % b, computed using repeated squaring."""
    
    bin_rep=bin(e) # first find binary representation of e 
    prod = 1    
    cur_pow_two = a # this will run through a^(2^0), a^(2^1), a^(2^2),...
    for i in range(len(bin_rep)-2): 
        if bin_rep[-1-i]=='1':
            prod = prod * cur_pow_two % b # multiply by cur_pow_two if the corresponding bit is 1
        cur_pow_two = (cur_pow_two)^2 % b  
    return prod

# runtime (for fixed a,b) is approx number of bits of e
In [46]:
exp_mod(2,10,18)
Out[46]:
16
In [47]:
exp_mod(2,10000000,19) # this is very fast
Out[47]:
17
In [48]:
exp_mod_slow(2,10000000,19) # takes several seconds
Out[48]:
17
In [49]:
exp_mod(3, 3^1000 + 17^1000, 23) # very fast; this would produce an error with exp_mod_slow
Out[49]:
9
In [50]:
# Fermat's little theorem: a^(p-1) mod p = 1
# if p prime, and a is not a multiple of p
In [51]:
# Example
exp_mod(3,16,17)
Out[51]:
1
In [52]:
# Euler's theorem: a^(phi(n)) = 1 mod n, 
# n is any number, a relatively prime to n 
# What is phi(n)? (For n prime, phi(n)=n-1)
# phi(n) = # of positive integer less than n 
# that are relatively prime to n
In [53]:
# Example: phi(10) = 4
exp_mod(3,4,10)
Out[53]:
1