We reviewed some properties of prime numbers and used them for primality testing. We reviewed the Diffie-Hellman protocol for finding a shared secret key and also tried to crack it. Then, we discussed OOP and special methods that allow us to define operators in Python.
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
Fermat's little theorem: if $p$ is prime and $1 < a < p$, then $a^{p-1} (\textrm{mod}\ p) \equiv 1$
On the other hand: If $p$ is composite, write $p = ab$ for some $a,b>1$. Now, assume toward contradiction that $a^{ab-1} (\textrm{mod}\ ab) \equiv 1$. It follows that $ ab ~|~a^{ab-1} -1$ and therefore $a~|~a^{ab-1} -1$. However, as clearly $a ~|~a^{ab-1}$ we must have $a~|~1$, a contradiction.
We can use Fermat's little theorem in order test whether a given number is prime. Note that if the number has $n$ bits than testing all possible $a$-s will require $O(2^n)$ iterations (a lot!).
Instead, we will try 100 random $a$-s in the range and see if one works as a witness.
import random
def is_prime(m, show_witness=False):
""" probabilistic test for m's compositeness """''
for i in range(0, 100):
a = random.randint(1, m - 1) # a is a random integer in [1..m-1]
if pow(a, m - 1, m) != 1:
if show_witness: # caller wishes to see a witness
print(m, "is composite", "\n", a, "is a witness, i=", i + 1)
return False
return True
For $a,b,c$ of at most $n$ bits each, time complexity of modpower is $O(n^3)$
def modpower(a, b, c):
""" computes a**b modulo c, using iterated squaring """
result = 1
while b > 0: # while b is nonzero
if b % 2 == 1: # b is odd
result = (result * a) % c
a = (a*a) % c
b = b // 2
return result
First, notice that if the function says that an imput number $m$ is not prime, then it is true. The function can make a mistake only is the case where a number $m$ is not prime, and is excidentally categorized by the function as prime. This can happen if all $100$ $a$'s that the function tried were not witnesses.
A quick computation shows that if $m$ is not a Charmichael number then at least $\frac{1}{2}$ of all possible $a$s are witnesses, so in almost all cases the probability for error is $(\frac{1}{2})^{100}$ (this is extremely low).
We decide on the size of the sample (to avoid testing all possible $2^{n-1}$ numbers of $n$ bits) and test whether each number we sample is prime. Then we divide the number of primes with the size of the sample.
def prob_prime(n, sample):
cnt = 0
for i in range(sample):
m = random.randint(2**(n-1), 2**n - 1)
cnt += is_prime(m)
return cnt / sample
prob_prime(2, 10**4)
prob_prime(3, 10**4)
1.0
0.5057
prob_prime(100, 10**4)
0.0148
prob_prime(200, 10**4)
0.0081
Diffie Hellman from lecture
def DH_exchange(p):
""" generates a shared DH key """
g = random.randint(1, p - 1)
a = random.randint(1, p - 1)# Alice's secret
b = random.randint(1, p - 1)# Bob's secret
x = pow(g, a, p)
y = pow(g, b, p)
key_A = pow(y, a, p)
key_B = pow(x, b, p)
#the next line is different from lecture
return g, a, b, x, y, key_A #key_A=key_B
def find_prime(n):
""" find random n-bit long prime """
while(True):
candidate = random.randrange(2**(n-1), 2**n)
if is_prime(candidate):
return candidate
Demostration:
import random
p = find_prime(10)
print(p)
g, a, b, x, y, key = DH_exchange(p)
g, a, b, x, y, key
557
(247, 88, 353, 367, 2, 206)
print(pow(g, a, p))
print(pow(x, b, p))
367 206
There is no known way to find $a$ efficiently, so we try the naive one: iterating over all $a$-s and cheking whether the equation $g^a \mod p = x$ holds for them.
If we found $a'$ that satisfies the condition but is not the original $a$, does it matter?
The time complexity of crack_DH is $O(2^nn^3)$
def crack_DH(p, g, x):
''' find secret "a" that satisfies g**a%p == x
Not feasible for large p '''
for a in range(1, p - 1):
if a % 100000 == 0:
print(a) #just to estimate running time
if pow(g, a, p) == x:
return a
return None #should never get here
#g, a, b, x, y, key = DH_exchange(p)
print(a)
crack_DH(p, g, x)
88
88
The algorithm crack_DH can return a different private key ($a$) than the one chosen by Alice, i.e. - $crack\_DH(p,g,x) = a' \neq a$, however, in this case we have, by definition of the cracking algorithm: $g^{a'} \textrm{ mod } p = g^{a} \textrm{ mod } p$, thus:
$$y^{a'} \textrm{ mod } p = \left(g^{b}\right)^{a'} \textrm{ mod } p = \left(g^{a'}\right)^{b} \textrm{ mod } p = \left(g^{a} \textrm{ mod } p \right)^{b} \textrm{ mod } p = g^{ab} \textrm{ mod } p$$I.e. - we can still compute the shared secret!
import random
p = find_prime(100)
print(p)
g, a, b, x, y, key = DH_exchange(p)
print(g, a, b, x, y, key)
crack_DH(p, g, x)
972762959355461077758092346931 711085486200126306605788797672 301337617630740085784318464474 953829050690741734973058379801 673125365307331009336129574105 220291338338430064459432090221 203856719177576397001814724913 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-78-6b3d13779c00> in <module>() 5 print(g, a, b, x, y, key) 6 ----> 7 crack_DH(p, g, x) <ipython-input-76-9a63fe3af2e8> in crack_DH(p, g, x) 3 Not feasible for large p ''' 4 for a in range(1, p - 1): ----> 5 if a % 100000 == 0: 6 print(a) #just to estimate running time 7 if pow(g, a, p) == x: KeyboardInterrupt:
Analyzing the nubmer of years it will take to crack the protocol if $a$ is found at the end (assuming iterating over 100000 $a$s takes a second)
a
301337617630740085784318464474
a//100000/60/60/24/365
9.55535317195396e+16
A link to special method names in python: http://getpython3.com/diveintopython3/special-method-names.html
We represent a polynomial as a list of coefficients, where the $i$'th element in the list is the coefficient of $x^i$
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
f = Polynomial([1, 1])
g = Polynomial([2, 3, 4])
f
g
(1*x^0) + (1*x^1)
(2*x^0) + (3*x^1) + (4*x^2)
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
f = Polynomial([1, 1])
f.degree()
f.evaluate(5)
1
6
p = Polynomial([1,0,1])
q = Polynomial([1,0,1])
p == q
False
Without it, Python compares memory adderesses
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
f = Polynomial([1, 1])
f2 = Polynomial([1, 1])
g = Polynomial([2, 3, 4])
f == g
f == f2
f.__eq__(f2)
False
True
True
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
def __add__(self, other):
assert isinstance(other, (Polynomial, int))
if isinstance(other, int):
terms = [c for c in self.coeffs]
terms[0] += other
return Polynomial(terms)
#else, other is a polynomial
size = min(self.degree(), other.degree()) + 1
terms = [self.coeffs[i] + other.coeffs[i] for i in range(size)]
terms += self.coeffs[size : self.degree() + 1]
terms += other.coeffs[size : other.degree() + 1]
#remove leading zeros
last = len(terms) - 1
while last > 0 and terms[last] == 0:
last -= 1
return Polynomial(terms[:last+1])
f = Polynomial([1, 1])
g = Polynomial([2, 3, 4])
f+g
f+100
100+f
(3*x^0) + (4*x^1) + (4*x^2)
(101*x^0) + (1*x^1)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-85-886b2663ae8b> in <module>() 46 f+g 47 f+100 ---> 48 100+f TypeError: unsupported operand type(s) for +: 'int' and 'Polynomial'
Since there is no support for addition with a Polynomial in the class int, python needs a definition of right addition for int in class Polynomial
When executing the command $f + 3$, python will first search for an __add__ method in class int that can support a second operand of type polynomial (self = int, other = polynomial). If such method does not exist, python will then search for an __radd__ method in class polynomial that can support a second operand of type int (self= polynomial, other = int). If such method does not exist, an exception is thrown.
Since + is commutative, then we set __radd__ to be identical to __add__
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
def __add__(self, other):
assert isinstance(other, (Polynomial, int))
if isinstance(other, int):
terms = [c for c in self.coeffs]
terms[0] += other
return Polynomial(terms)
#else, other is a polynomial
size = min(self.degree(), other.degree()) + 1
terms = [self.coeffs[i] + other.coeffs[i] for i in range(size)]
terms += self.coeffs[size : self.degree() + 1]
terms += other.coeffs[size : other.degree() + 1]
#remove leading zeros
last = len(terms) - 1
while last > 0 and terms[last] == 0:
last -= 1
return Polynomial(terms[:last+1])
__radd__ = __add__ #to allow int+Polynomial
f = Polynomial([1, 1])
1 + f
(2*x^0) + (1*x^1)
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
def __add__(self, other):
assert isinstance(other, (Polynomial, int))
if isinstance(other, int):
terms = [c for c in self.coeffs]
terms[0] += other
return Polynomial(terms)
#else, other is a polynomial
size = min(self.degree(), other.degree()) + 1
terms = [self.coeffs[i] + other.coeffs[i] for i in range(size)]
terms += self.coeffs[size : self.degree() + 1]
terms += other.coeffs[size : other.degree() + 1]
#remove leading zeros
last = len(terms) - 1
while last > 0 and terms[last] == 0:
last -= 1
return Polynomial(terms[:last+1])
__radd__ = __add__ #to allow int+Polynomial
def __neg__(self):
return Polynomial([-co for co in self.coeffs])
f = Polynomial([1, 1])
-f
(-1*x^0) + (-1*x^1)
Using the fact we already have __neg__ and __add__ we define $f-g = f + (-g)$
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
def __add__(self, other):
assert isinstance(other, (Polynomial, int))
if isinstance(other, int):
terms = [c for c in self.coeffs]
terms[0] += other
return Polynomial(terms)
#else, other is a polynomial
size = min(self.degree(), other.degree()) + 1
terms = [self.coeffs[i] + other.coeffs[i] for i in range(size)]
terms += self.coeffs[size : self.degree() + 1]
terms += other.coeffs[size : other.degree() + 1]
#remove leading zeros
last = len(terms) - 1
while last > 0 and terms[last] == 0:
last -= 1
return Polynomial(terms[:last+1])
__radd__ = __add__ #to allow int+Polynomial
def __neg__(self):
return Polynomial([-co for co in self.coeffs])
def __sub__(self, other):
assert isinstance(other, (int, Polynomial))
return (self + (-other))
f = Polynomial([1, 1])
g = Polynomial([2, 3, 4])
g - f
g - 1
(1*x^0) + (2*x^1) + (4*x^2)
(1*x^0) + (3*x^1) + (4*x^2)
g
(2*x^0) + (3*x^1) + (4*x^2)
Using the fact we already have __neg__ and __sub__ we define $f-g = -(g-f)$
since oprator - is not commutative, then we cannot set __rsub__ to be identical to __sub__
Note that the following implementation of __rsub__ is promlematic, because it will keep calling itself (infinitely)
def __rsub__(self, other):
return other-self #PROBLEMATIC!
class Polynomial():
def __init__(self, coeffs):
self.coeffs = coeffs
def __repr__(self):
terms = [" + (" + str(self.coeffs[k]) + "*x^" + \
str(k) + ")" \
for k in range(len(self.coeffs)) \
if self.coeffs[k] != 0]
terms = "".join(terms)
if terms == "":
return "0"
else:
return terms[3:] #discard leftmost '+'
def degree(self):
return len(self.coeffs) - 1
def evaluate(self, x):
result = 0
for i, c in enumerate(self.coeffs):
result += c * pow(x, i)
return result
def __eq__(self, other):
assert isinstance(other, Polynomial)
return self.coeffs == other.coeffs
def __add__(self, other):
assert isinstance(other, (Polynomial, int))
if isinstance(other, int):
terms = [c for c in self.coeffs]
terms[0] += other
return Polynomial(terms)
#else, other is a polynomial
size = min(self.degree(), other.degree()) + 1
terms = [self.coeffs[i] + other.coeffs[i] for i in range(size)]
terms += self.coeffs[size : self.degree() + 1]
terms += other.coeffs[size : other.degree() + 1]
#remove leading zeros
last = len(terms) - 1
while last > 0 and terms[last] == 0:
last -= 1
return Polynomial(terms[:last+1])
__radd__ = __add__ #to allow int+Polynomial
def __neg__(self):
return Polynomial([-co for co in self.coeffs])
def __sub__(self, other):
assert isinstance(other, (int, Polynomial))
return (self + (-other))
#__rsub__ = __sub__ #not what we want...
def __rsub__(self, other):
# other - self
return -(self - other) #why not just other-self ?
f = Polynomial([1, 1])
g = Polynomial([2, 3, 4])
1 - f
1 - g
(-1*x^1)
(-1*x^0) + (-3*x^1) + (-4*x^2)