cs1001.py , Tel Aviv University, Spring 2020

Recitation 8

We continued with another recurssion question. 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.

Takeaways:

  1. The probabilistic function is_prime, that uses Fermat's primality test, can be used to detect primes quickly and efficiently, but has a (very small) probability of error. Its time complexity is $O(n^3)$, where $n$ is the number of bits of the input.
  2. The DH protocol relies on two main principles: the following equality $(g^{a}\mod p)^b \mod p = g^{ab} \mod p $ and the (believed) hardness of the discrete log problem (given $g,p$, and $x = g^{a} \mod p$ finding $a$ is hard). Make sure you understand the details of the protocol.

Code for printing several outputs in one cell (not part of the recitation):

In [1]:
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

Recursion: another example

$N$-Queens

The $N$ queens problem is to determine how many possibilities are there to legally place $N$ queens on an $N$-by-$N$ chess board. Legally means no queen threatens another queen.

Solution: We build the solution incrementally, column by column. We maintain a partial solution (implemented as a list), which is initially empty.

The function $legal(partial, i)$ is given:

  • a partial solution placing the first $len(partial)$ queens on the leftmost columns
  • a positive integer $0\leq i<N$

It returns True if it is legal to place a queen in the next column on row $i$, given the partial placement $partial$

In [10]:
def legal(partial, i):
    ''' Can we place a queen in the next column in row i ? '''
    
    # any queens in the same row to the left ?
    left = [j for j in partial if j==i] 
    # diagonal up - left
    diag_up = [j for j in partial if j - partial.index(j) == i - len(partial)]
    # diagonal down - left
    diag_down = [j for j in partial if j + partial.index(j) == i + len(partial)]
    
    res = ( left == diag_up == diag_down == [])
    
    return res

The solution:

In [11]:
def queens (n, show = True):
    ''' how many ways to place n queens on an nXn board ? '''
    partial = [] # list representing partial placement of queens
    return queens_rec (n, partial, show)

def queens_rec (n, partial, show):
    ''' Given a list representing partial placement of queens ,
    can we legally extend it ? '''
    if len(partial) == n: # all n queens are placed legally
        return 1
    else:
        cnt = 0
        for i in range(n):
            # try to place a queen in row i of the next column
            if legal(partial, i):
                cnt += queens_rec(n, partial + [i], show)
        return cnt
In [12]:
print(queens(0))
print(queens(1))
print(queens(2))
print(queens(3))
print(queens(10))
1
1
0
0
724

Some intuition: assume we can find a solution for placing $k<N$ queens, how do we expand the solution to $k+1$? queens_rec returns the number of possible legal placements for $N$ queens, where $k$ are already placed at the leftmost columns and there are $N-k$ queens left to place. The recursive idea: Legally place queen number $(k+1)$ and recursively solve the problem, when there is one less queen to place.

Note that the complexity is $O(N!)$ (greater than $O(2^N)$)

Primes and Diffie-Hellman

Fermat's little theorem

Fermat's little theorem states that if $p$ is prime and $1 < a < p$, then $a^{p-1} (\textrm{mod}\ p) \equiv 1$

Compositeness witnesses

A witness is a piece of evidence that can be produced in order to prove a claim. In our case, the problem we are tackling is deciding whether a given number $m$ is prime or composite. We now describe three types of witnesses for compositeness:

  • A clear witness for compositeness can be a (not necessarily prime) factor of $m$. That is, a number $1 < b < m$ such that $b ~|~ m$ ($b$ divides $m$). Since $b$ is a non-trivial factor of $m$, $m$ is clearly composite.
  • Another, less trivial witness of compositeness is a GCD witness, that is a number $1 < b < m$ such that $\mathrm{gcd}(m,b) > 1$. Think, why is $b$ a witness of compositeness? Let $\mathrm{gcd}(m,b) = r$, then by definition $r ~|~ m$. As $r \leq b < m$ we get that $r$ is a factor of $m$ and thus $m$ is composite.
  • Fermat's primality test defines another set of witness of compositeness - a number $1 < a < m$ is a Fermat witness if $a^{m-1} (\textrm{mod}\ m) \not\equiv 1$

Why go through all of this work? Our tactic would be to randomly draw a number $b \in {1,\ldots, m - 1}$ and hope that $b$ is some witness of compositeness. Clearly, we'd like our witness pool to be as large as possible.

Now, if $m$ is a composite number, let $\mathrm{FACT}_m, \mathrm{GCD}_m, \mathrm{FERM}_m$ be the set of prime factors, GCD witnesses and Fermat witnesses for $m$'s compositeness. It is not hard to show that $$\mathrm{FACT}_m \subseteq \mathrm{GCD}_m \subseteq \mathrm{FERM}_m$$

But the real strength of Fermat's primality test comes from the fact that if $m$ is composite, then apart from very rare cases (where $n$ is a Carmichael number) it holds that $|\mathrm{FERM}_m| \geq m/2$. That is - a random number is a Fermat witness w.p. at least $1/2$.

A side note (for reference only) - Carmichael numbers are exactly the composite numbers $m$ where $\mathrm{GCD}_m = \mathrm{FERM}_m$

Every factor of a composite number is a Fermat's witness

Let $m$ be a composite number and write $m = ab$ for some $a,b>1$. We claim that $a$ is a Fermat witness. To see this, assume towards contradiction that $a^{m-1} (\textrm{mod}\ m) \equiv 1$, i.e. $a^{m-1} = c\cdot m + 1= c \cdot a \cdot b + 1$ for some $c \geq 1$.

Rearrange the above to get $a(a^{m-2} - c\cdot b) = 1$. However, $a > 1$ and $(a^{m-2} - c\cdot b) \in \mathbb{Z}$, a contradiction.

Primality test using Fermat's witness

We can use Fermat's little theorem in order test whether a given number $m$ is prime. That is, we can test whether we find a Fermat's witness $a\in\mathrm{FERM}_m$ for compositeness. 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 Fermat's witness.

In [9]:
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)$

In [14]:
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

Runtime analysis:

  • The main loop runs over $b$, dividing $b \to b/2$ at each iteration, so it runs $O(n)$ times.
  • In each iteration we do:
    • One operation of $b%2$ in $O(1)$ time
    • One operation of $b/2$ in $O(1)$ time
    • At most two multiplication and two modulu operations
    • Multiplication of two $n$ bit numbers runs in time $O(n^2)$
    • Modulu can be implemented by addition, division and multiplication: $a \textrm{ mod } b = a - (a // b) b$ and division runs in time $O(n^2)$ same as multiplication
    • Finally, the modulu operation keeps all numbers at most $n$ bits, thus the running time does not increase with each iteration
  • In total - $O(n^3)$

The probability of error:

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).

Failing over Carmichael numbers

In the rare case where we get a Carmichael number as input we have no guarantee on the performance of the test. In the following example we test whether the Carmichael number $N$ (which is composite) is prime or not using our test.

The number $N$ comes courtesy of Wikipedia

In [15]:
# p is a large prime:
p = 29674495668685510550154174642905332730771991799853043350995075531276838753171770199594238596428121188033664754218345562493168782883

# N is a Carmichael number (and is obviously composite):
N = p*(313*(p-1)+1)*(353*(p-1) + 1)

is_prime(N)
Out[15]:
True

Testing the prime number theorem: For a large n, a number of n bits is prime with a prob. of O(1/n)

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.

In [16]:
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
In [17]:
prob_prime(2, 10**4)
prob_prime(3, 10**4)
Out[17]:
1.0
Out[17]:
0.5011
In [18]:
prob_prime(100, 10**4)
Out[18]:
0.0133
In [19]:
prob_prime(200, 10**4)
Out[19]:
0.0063

Diffie Hellman from lecture

The protocol as code

In [8]:
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

Find a prime number

In [10]:
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:

In [21]:
import random
p = find_prime(10)
print(p)
g, a, b, x, y, key = DH_exchange(p)
g, a, b, x, y, key
971
Out[21]:
(226, 96, 658, 536, 69, 793)
In [22]:
print(pow(g, a, p))
print(pow(x, b, p))
536
793

Crack the Diffie Hellman code

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)$

In [24]:
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
In [28]:
g, a, b, x, y, key = DH_exchange(p)
print(a)
crack_DH(p, g, x)
617
Out[28]:
617

Trying to crack the protocol with a 100 bit prime

In [30]:
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)
888958161829799043985061475553
867441475552671949999208459311 553248002741623339226706792959 414752428718086349833513855252 225642142325534266491072467972 646307605600255401720101354566 203280842942198745051972342627
100000
200000
300000
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-30-6b3d13779c00> in <module>
      5 print(g, a, b, x, y, key)
      6 
----> 7 crack_DH(p, g, x)

<ipython-input-24-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)

In [31]:
a
Out[31]:
553248002741623339226706792959
In [32]:
a//100000/60/60/24/365
Out[32]:
1.75433790823701e+17
In [ ]: