We have reviewed some properties of prime numbers and used them for primality testing. We revied the Diffie-Hellman protocol for finding a shared secret key and also tried to crack it. At the second part of the recitation, 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$
Equivalently: if $m$ is not a prime then there exists $1 < a < m$ such that $a^{m-1} (\textrm{mod}\ m) \not\equiv 1$. Such a number $a$ is called a witness to the fact that $m$ is not prime.
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. According to the Miller-Rabin theorem $\frac{3}{4}$ of all possible $a$s are witnesses, so the probability for error is $(\frac{1}{4})^{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**5)
prob_prime(200, 10**4)
0.0075
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)
p
g,a,b,x,y,key = DH_exchange(p)
g,a,b,x,y,key
883
(196, 96, 702, 747, 46, 386)
pow(g, a, p)
pow(x, b, p)
747
386
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
import random
p = find_prime(100)
# p
g,a,b,x,y,key = DH_exchange(p)
g,a,b,x,y,key
# crack_DH(p,g,x)
(585908982447390826849189842691, 541120442351353498570334017705, 1930761781644675344995819292, 70659930174800164907189704268, 779461075719798735415516947073, 1045758344420625426709522664848)
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)
541120442351353498570334017705//100000//60//60//24//365
171588166651240962
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
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
(3*x^0) + (4*x^1) + (4*x^2)
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)
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):
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)