Part 6: Ciphers and Key exchange

In this notebook, we introduce cryptography -- how to communicate securely over insecure channels. We begin with a study of two basic ciphers, the Caesar cipher and its fancier variant, the Vigenère cipher. The Vigenère cipher uses a key to turn plaintext (i.e., the message) into ciphertext (the coded message), and uses the same key to turn the ciphertext back into plaintext. Therefore, two parties can communicate securely if they -- and only they -- possess the key.

If the security of communication rests on possession of a common key, then we're left with a new problem: how do the two parties agree on a common key, especially if they are far apart and communicating over an insecure channel?

A clever solution to this problem was published in 1976 by Whitfield Diffie and Martin Hellman, and so it's called Diffie-Hellman key exchange. It takes advantage of modular arithmetic: the existence of a primitive root (modulo a prime) and the difficulty of solving the discrete logarithm problem.

This part complements Chapter 6 of An Illustrated Theory of Numbers.

Table of Contents

Ciphers

A cipher is a way of transforming a message, called the plaintext into a different form, the ciphertext, which conceals the meaning to all but the intended recipient(s). A cipher is a code, and can take many forms. A substitution cipher might simply change every letter to a different letter in the alphabet. This is the idea behind "Cryptoquip" puzzles. These are not too hard for people to solve, and are easy for computers to solve, using frequency analysis (understanding how often different letters or letter-combinations occur).

ASCII code and the Caesar cipher

Even though substitution ciphers are easy to break, they are a good starting point. To implement substitution ciphers in Python, we need to study the string type in a bit more detail. To declare a string variable, just put your string in quotes. You can use any letters, numbers, spaces, and many symbols inside a string. You can enclose your string by single quotes, like 'Hello' or double-quotes, like "Hello". This flexibility is convenient, if you want to use quotes within your string. For example, the string Prince's favorite prime is 1999 should be described in Python with double-quotes "Prince's favorite prime is 1999" so that the apostrophe doesn't confuse things.

Strings are indexed, and their letters can be retrieved as if the string were a list of letters. Python experts will note that strings are immutable while lists are mutable objects, but we aren't going to worry about that here.

In [1]:
W = "Hello"
print W
for j in range(len(W)):  # len(W) is the length of the string W.
    print W[j]  # Access the jth character of the string.
Hello
H
e
l
l
o

Each "letter" of a string again belongs to the string type. A string of length one is called a character.

In [3]:
print type(W)
print type(W[0])  # W[0] is a character.
<type 'str'>
<type 'str'>

Since computers store data in binary, the designers of early computers (1960s) created a code called ASCII (American Standard Code for Information Interchange) to associate to each character a number between 0 and 127. Every number between 0 and 127 is represented in binary by 7 bits (between 0000000 and 1111111), and so each character is stored with 7 bits of memory. Later, ASCII was extended with another 128 characters, so that codes between 0 and 255 were used, requiring 8 bits. 8 bits of memory is called a byte. One byte of memory suffices to store one (extended ASCII) character.

You might notice that there are 256 ASCII codes available, but there are fewer than 256 characters available on your keyboard, even once you include symbols like # and ;. Some of these "extra" codes are for accented letters, and others are relics of old computers. For example, ASCII code 7 (0000111) stands for the "bell", and older readers might remember making the Apple II computer beep by pressing Control-G on the keyboard ("G" is the 7th letter). You can look up a full ASCII table if you're curious.

Nowadays, the global community of computer users requires far more than 256 "letters" -- there are many alphabets around the world! So instead of ASCII, we can access over 100 thousand unicode characters. Scroll through a unicode table to see what is possible. In Python version 3.x, all strings are considered in Unicode, but in Python 2.7 (which we use), it's a bit trickier to work with Unicode.

Here we stay within ASCII codes, since they will suffice for basic English messages. Python has built-in commands chr and ord for converting from code-number (0--255) to character and back again.

In [4]:
chr(65)
Out[4]:
'A'
In [5]:
ord('A')
Out[5]:
65

The following code will produce a table of the ASCII characters with codes between 32 and 126. This is a good range which includes all the most common English characters and symbols on a U.S. keyboard. Note that ASCII code 32 corresponds to an empty space (an important character for long messages!)

In [6]:
for a in range(32,127):
    c = chr(a)
    print "ASCII %d is %s"%(a, c)
ASCII 32 is  
ASCII 33 is !
ASCII 34 is "
ASCII 35 is #
ASCII 36 is $
ASCII 37 is %
ASCII 38 is &
ASCII 39 is '
ASCII 40 is (
ASCII 41 is )
ASCII 42 is *
ASCII 43 is +
ASCII 44 is ,
ASCII 45 is -
ASCII 46 is .
ASCII 47 is /
ASCII 48 is 0
ASCII 49 is 1
ASCII 50 is 2
ASCII 51 is 3
ASCII 52 is 4
ASCII 53 is 5
ASCII 54 is 6
ASCII 55 is 7
ASCII 56 is 8
ASCII 57 is 9
ASCII 58 is :
ASCII 59 is ;
ASCII 60 is <
ASCII 61 is =
ASCII 62 is >
ASCII 63 is ?
ASCII 64 is @
ASCII 65 is A
ASCII 66 is B
ASCII 67 is C
ASCII 68 is D
ASCII 69 is E
ASCII 70 is F
ASCII 71 is G
ASCII 72 is H
ASCII 73 is I
ASCII 74 is J
ASCII 75 is K
ASCII 76 is L
ASCII 77 is M
ASCII 78 is N
ASCII 79 is O
ASCII 80 is P
ASCII 81 is Q
ASCII 82 is R
ASCII 83 is S
ASCII 84 is T
ASCII 85 is U
ASCII 86 is V
ASCII 87 is W
ASCII 88 is X
ASCII 89 is Y
ASCII 90 is Z
ASCII 91 is [
ASCII 92 is \
ASCII 93 is ]
ASCII 94 is ^
ASCII 95 is _
ASCII 96 is `
ASCII 97 is a
ASCII 98 is b
ASCII 99 is c
ASCII 100 is d
ASCII 101 is e
ASCII 102 is f
ASCII 103 is g
ASCII 104 is h
ASCII 105 is i
ASCII 106 is j
ASCII 107 is k
ASCII 108 is l
ASCII 109 is m
ASCII 110 is n
ASCII 111 is o
ASCII 112 is p
ASCII 113 is q
ASCII 114 is r
ASCII 115 is s
ASCII 116 is t
ASCII 117 is u
ASCII 118 is v
ASCII 119 is w
ASCII 120 is x
ASCII 121 is y
ASCII 122 is z
ASCII 123 is {
ASCII 124 is |
ASCII 125 is }
ASCII 126 is ~

Since we only work with the ASCII range between 32 and 126, it will be useful to "cycle" other numbers into this range. For example, we will interpret 127 as 32, 128 as 33, etc., when we convert out-of-range numbers into characters.

The following function forces a number into a given range, using the mod operator. It's a common trick, to make lists loop around cyclically.

In [7]:
def inrange(n,range_min, range_max):
    '''
    The input number n can be any integer.
    The output number will be between range_min and range_max (inclusive)
    If the input number is already within range, it will not change.
    '''
    range_len = range_max - range_min + 1
    a = n % range_len
    if a < range_min:
        a = a + range_len
    return a
In [8]:
inrange(13,1,10)
Out[8]:
3
In [9]:
inrange(17,5,50)
Out[9]:
17

Now we can implement a substitution cipher by converting characters to their ASCII codes, shuffling the codes, and converting back. One of the simplest substitution ciphers is called a Caesar cipher, in which each character is shifted -- by a fixed amount -- down the list. For example, a Caesar cipher of shift 3 would send 'A' to 'D' and 'B' to 'E', etc.. Near the end of the list, characters are shifted back to the beginning -- the list is considered cyclicly, using our inrange function.

Here is an implementation of the Caesar cipher, using the ASCII range between 32 and 126. We begin with a function to shift a single character.

In [10]:
def Caesar_shift(c, shift):
    '''
    Shifts the character c by shift units
    within the ASCII table between 32 and 126. 
    The shift parameter can be any integer!
    '''
    ascii = ord(c)
    a = ascii + shift # Now we have a number between 32+shift and 126+shift.
    a = inrange(a,32,126) # Put the number back in range.
    return chr(a)

Let's see the effect of the Caesar cipher on our ASCII table.

In [12]:
for a in range(32,127):
    c = chr(a)
    print "ASCII %d is %s, which shifts to %s"%(a, c, Caesar_shift(c,5)) # Shift by 5.
ASCII 32 is  , which shifts to %
ASCII 33 is !, which shifts to &
ASCII 34 is ", which shifts to '
ASCII 35 is #, which shifts to (
ASCII 36 is $, which shifts to )
ASCII 37 is %, which shifts to *
ASCII 38 is &, which shifts to +
ASCII 39 is ', which shifts to ,
ASCII 40 is (, which shifts to -
ASCII 41 is ), which shifts to .
ASCII 42 is *, which shifts to /
ASCII 43 is +, which shifts to 0
ASCII 44 is ,, which shifts to 1
ASCII 45 is -, which shifts to 2
ASCII 46 is ., which shifts to 3
ASCII 47 is /, which shifts to 4
ASCII 48 is 0, which shifts to 5
ASCII 49 is 1, which shifts to 6
ASCII 50 is 2, which shifts to 7
ASCII 51 is 3, which shifts to 8
ASCII 52 is 4, which shifts to 9
ASCII 53 is 5, which shifts to :
ASCII 54 is 6, which shifts to ;
ASCII 55 is 7, which shifts to <
ASCII 56 is 8, which shifts to =
ASCII 57 is 9, which shifts to >
ASCII 58 is :, which shifts to ?
ASCII 59 is ;, which shifts to @
ASCII 60 is <, which shifts to A
ASCII 61 is =, which shifts to B
ASCII 62 is >, which shifts to C
ASCII 63 is ?, which shifts to D
ASCII 64 is @, which shifts to E
ASCII 65 is A, which shifts to F
ASCII 66 is B, which shifts to G
ASCII 67 is C, which shifts to H
ASCII 68 is D, which shifts to I
ASCII 69 is E, which shifts to J
ASCII 70 is F, which shifts to K
ASCII 71 is G, which shifts to L
ASCII 72 is H, which shifts to M
ASCII 73 is I, which shifts to N
ASCII 74 is J, which shifts to O
ASCII 75 is K, which shifts to P
ASCII 76 is L, which shifts to Q
ASCII 77 is M, which shifts to R
ASCII 78 is N, which shifts to S
ASCII 79 is O, which shifts to T
ASCII 80 is P, which shifts to U
ASCII 81 is Q, which shifts to V
ASCII 82 is R, which shifts to W
ASCII 83 is S, which shifts to X
ASCII 84 is T, which shifts to Y
ASCII 85 is U, which shifts to Z
ASCII 86 is V, which shifts to [
ASCII 87 is W, which shifts to \
ASCII 88 is X, which shifts to ]
ASCII 89 is Y, which shifts to ^
ASCII 90 is Z, which shifts to _
ASCII 91 is [, which shifts to `
ASCII 92 is \, which shifts to a
ASCII 93 is ], which shifts to b
ASCII 94 is ^, which shifts to c
ASCII 95 is _, which shifts to d
ASCII 96 is `, which shifts to e
ASCII 97 is a, which shifts to f
ASCII 98 is b, which shifts to g
ASCII 99 is c, which shifts to h
ASCII 100 is d, which shifts to i
ASCII 101 is e, which shifts to j
ASCII 102 is f, which shifts to k
ASCII 103 is g, which shifts to l
ASCII 104 is h, which shifts to m
ASCII 105 is i, which shifts to n
ASCII 106 is j, which shifts to o
ASCII 107 is k, which shifts to p
ASCII 108 is l, which shifts to q
ASCII 109 is m, which shifts to r
ASCII 110 is n, which shifts to s
ASCII 111 is o, which shifts to t
ASCII 112 is p, which shifts to u
ASCII 113 is q, which shifts to v
ASCII 114 is r, which shifts to w
ASCII 115 is s, which shifts to x
ASCII 116 is t, which shifts to y
ASCII 117 is u, which shifts to z
ASCII 118 is v, which shifts to {
ASCII 119 is w, which shifts to |
ASCII 120 is x, which shifts to }
ASCII 121 is y, which shifts to ~
ASCII 122 is z, which shifts to  
ASCII 123 is {, which shifts to !
ASCII 124 is |, which shifts to "
ASCII 125 is }, which shifts to #
ASCII 126 is ~, which shifts to $

Now we can use the Caesar cipher to encrypt strings.

In [13]:
def Caesar_cipher(plaintext, shift):
    ciphertext = ''
    for c in plaintext:  # Iterate through the characters of a string.
        ciphertext = ciphertext + Caesar_shift(c,shift)  
    return ciphertext
In [14]:
print Caesar_cipher('Hello!  Can you read this?', 5)  # Shift forward 5 units in ASCII.
Mjqqt&%%Hfs%~tz%wjfi%ymnxD

As designed, the Caesar cipher turns plaintext into ciphertext by using a shift of the ASCII table. To decipher the ciphertext, one can just use the Caesar cipher again, with the negative shift.

In [15]:
print Caesar_cipher('Mjqqt&%%Hfs%~tz%wjfi%ymnxD', -5)  # Shift back 5 units in ASCII.
Hello!  Can you read this?

The Vigenère cipher

The Caesar cipher is pretty easy to break, by a brute force attack (shift by all possible values) or a frequency analysis (compare the frequency of characters in a message to the frequency of characters in typical English messages, to make a guess).

The Vigenère cipher is a variant of the Caesar cipher which uses an ecryption key to vary the shift-parameter throughout the encryption process. For example, to encrypt the message "This is very secret" using the key "Key", you line up the characters of the message above repeated copies of the key.

T h i s i s v e r y s e c r e t
K e y K e y K e y K e y K e y K e y K

Then, you turn everything into ASCII (or your preferred numerical system), and use the bottom row to shift the top row.

ASCII message 84 104 105 115 32 105 115 32 118 101 114 121 32 115 101 99 114 101 116
Shift 75 101 121 75 101 121 75 101 121 75 101 121 75 101 121 75 101 121 75
ASCII shifted 159 205 226 190 133 226 190 133 239 176 215 242 107 216 222 174 215 222 191
ASCII shifted in range 64 110 36 95 38 36 95 38 49 81 120 52 107 121 32 79 120 32 96

Finally, the shifted ASCII codes are converted back into characters for transmission. In this case, the codes 64,110,36,95, etc., are converted to the ciphertext "@n$_&$_&1Qx4ky Ox \`"

The Vigenère cipher is much harder to crack than the Caesar cipher, if you don't have the key. Indeed, the varying shifts make frequency analysis more difficult. The Vigenère cipher is weak by today's standards (see Wikipedia for a description of 19th century attacks), but illustrates the basic actors in a symmetric key cryptosystem: the plaintext, ciphertext, and a single key. Today, symmetric key cryptosystems like AES and 3DES are used all the time for secure communication.

Below, we implement the Vigenère cipher.

In [16]:
def Vigenere_cipher(plaintext, key):
    ciphertext = '' # Start with an empty string
    for j in range(len(plaintext)): 
        c = plaintext[j] # the jth letter of the plaintext
        key_index = j % len(key) # Cycle through letters of the key.
        shift = ord(key[key_index]) # How much we shift c by.
        ciphertext = ciphertext + Caesar_shift(c,shift) # Add new letter to ciphertext
    return ciphertext
In [17]:
print Vigenere_cipher('This is very secret', 'Key') # 'Key' is probably a bad key!!
@n$_&$_&1Qx4ky Ox `

The Vigenère cipher is called a symmetric cryptosystem, because the same key that is used to encrypt the plaintext can be used to decrypt the ciphertext. All we do is subtract the shift at each stage.

In [18]:
def Vigenere_decipher(ciphertext, key):
    plaintext = '' # Start with an empty string
    for j in range(len(ciphertext)): 
        c = ciphertext[j] # the jth letter of the ciphertext
        key_index = j % len(key) # Cycle through letters of the key.
        shift = - ord(key[key_index]) # Note the negative sign to decipher!
        plaintext = plaintext + Caesar_shift(c,shift) # Add new letter to plaintext
    return plaintext
In [19]:
Vigenere_decipher('@n$_&$_&1Qx4ky Ox `', 'Key')
Out[19]:
'This is very secret'
In [20]:
# Try a few cipher/deciphers yourself to get used to the Vigenere system.

The Vigenère cipher becomes an effective way for two parties to communicate securely, as long as they share a secret key. In the 19th century, this often meant that the parties would require an initial in-person meeting to agree upon a key, or a well-guarded messenger would carry the key from one party to the other.

Today, as we wish to communicate securely over long distances on a regular basis, the process of agreeing on a key is more difficult. It seems like a chicken-and-egg problem, where we need a shared secret to communicate securely, but we can't share a secret without communicating securely in the first place!

Remarkably, this secret-sharing problem can be solved with some modular arithmetic tricks. This is the subject of the next section.

Exercises

  1. A Caesar cipher was used to encode a message, with the resulting ciphertext: 'j!\'1r$v1"$v&&+1t}v(v$2'. Use a loop (brute force attack) to figure out the original message.

  2. Imagine that you encrypt a long message (e.g., 1000 words of standard English) with a Vigenère cipher. How might you detect the length of the key, if it is short (e.g. 3 or 4 characters)?

  3. Consider running a plaintext message through a Vigenère cipher with a 3-character key, and then running the ciphertext through a Vigenère cipher with a 4-character key. Explain how this is equivalent to running the original message through a single cipher with a 12-character key.

Key exchange

Now we study Diffie-Hellman key exchange, a remarkable way for two parties to share a secret without ever needing to directly communicate the secret with each other. Their method is based on properties of modular exponentiation and the existence of a primitive root modulo prime numbers.

Primitive roots and Sophie Germain primes

If $p$ is a prime number, and $GCD(a,p) = 1$, then recall Fermat's Little Theorem: $$a^{p-1} \equiv 1 \text{ mod } p.$$ It may be the case that $a^\ell \equiv 1$ mod $p$ for some smaller (positive) value of $\ell$ however. The smallest such positive value of $\ell$ is called the order (multiplicative order, to be precise) of $a$ modulo $p$, and it is always a divisor of $p-1$.

The following code determines the order of a number, mod $p$, with a brute force approach.

In [21]:
def mult_order(a,p):
    '''
    Determines the (multiplicative) order of an integer
    a, modulo p.  Here p is prime, and GCD(a,p) = 1.
    If bad inputs are used, this might lead to a 
    never-ending loop!
    '''
    current_number = a % p
    current_exponent = 1
    while current_number != 1:
        current_number = (current_number * a)%p
        current_exponent = current_exponent + 1
    return current_exponent
    
In [24]:
for j in range(1,37):
    print "The multiplicative order of %d modulo 37 is %d"%(j,mult_order(j,37))  
    # These orders should all be divisors of 36.
The multiplicative order of 1 modulo 37 is 1
The multiplicative order of 2 modulo 37 is 36
The multiplicative order of 3 modulo 37 is 18
The multiplicative order of 4 modulo 37 is 18
The multiplicative order of 5 modulo 37 is 36
The multiplicative order of 6 modulo 37 is 4
The multiplicative order of 7 modulo 37 is 9
The multiplicative order of 8 modulo 37 is 12
The multiplicative order of 9 modulo 37 is 9
The multiplicative order of 10 modulo 37 is 3
The multiplicative order of 11 modulo 37 is 6
The multiplicative order of 12 modulo 37 is 9
The multiplicative order of 13 modulo 37 is 36
The multiplicative order of 14 modulo 37 is 12
The multiplicative order of 15 modulo 37 is 36
The multiplicative order of 16 modulo 37 is 9
The multiplicative order of 17 modulo 37 is 36
The multiplicative order of 18 modulo 37 is 36
The multiplicative order of 19 modulo 37 is 36
The multiplicative order of 20 modulo 37 is 36
The multiplicative order of 21 modulo 37 is 18
The multiplicative order of 22 modulo 37 is 36
The multiplicative order of 23 modulo 37 is 12
The multiplicative order of 24 modulo 37 is 36
The multiplicative order of 25 modulo 37 is 18
The multiplicative order of 26 modulo 37 is 3
The multiplicative order of 27 modulo 37 is 6
The multiplicative order of 28 modulo 37 is 18
The multiplicative order of 29 modulo 37 is 12
The multiplicative order of 30 modulo 37 is 18
The multiplicative order of 31 modulo 37 is 4
The multiplicative order of 32 modulo 37 is 36
The multiplicative order of 33 modulo 37 is 9
The multiplicative order of 34 modulo 37 is 9
The multiplicative order of 35 modulo 37 is 36
The multiplicative order of 36 modulo 37 is 2

A theorem of Gauss states that, if $p$ is prime, there exists an integer $b$ whose order is precisely $p-1$ (as big as possible!). Such an integer is called a primitive root modulo $p$. For example, the previous computation found 12 primitive roots modulo $37$: they are 2,5,13,15,17,18,19,20,22,24,32,35. To see these illustrated (mod 37), check out this poster (yes, that is blatant self-promotion!)

For everything that follows, suppose that $p$ is a prime number. Not only do primitive roots exist mod $p$, but they are pretty common. In fact, the number of primitive roots mod $p$ equals $\phi(p-1)$, where $\phi$ denotes Euler's totient. On average, $\phi(n)$ is about $6 / \pi^2$ times $n$ (for positive integers $n$). While numbers of the form $p-1$ are not "average", one still expects that $\phi(p-1)$ is a not-very-small fraction of $p-1$. You should not have to look very far if you want to find a primitive root.

The more difficult part, in practice, is determining whether a number $b$ is or is not a primitive root modulo $p$. When $p$ is very large (like hundreds or thousands of digits), $p-1$ is also very large. It is certainly not practical to cycle all the powers (from $1$ to $p-1$) of $b$ to determine whether $b$ is a primitive root!

The better approach, sometimes, is to use the fact that the multiplicative order of $b$ must be a divisor of $p-1$. If one can find all the divisors of $p-1$, then one can just check whether $b^d \equiv 1$ mod $p$ for each divisor $d$. This makes the problem of determining whether $b$ is a primitive root just about as hard as the problem of factoring $p-1$. This is a hard problem, in general!

But, for the application we're interested in, we will want to have a large prime number $p$ and a primitive root mod $p$. The easiest way to do this is to use a Sophie Germain prime $q$. A Sophie Germain prime is a prime number $q$ such that $2q + 1$ is also prime. When $q$ is a Sophie Germain prime, the resulting prime $p = 2q + 1$ is called a safe prime.

Observe that when $p$ is a safe prime, the prime decomposition of $p-1$ is $$p-1 = 2 \cdot q.$$ That's it. So the possible multiplicative orders of an element $b$, mod $p$, are the divisors of $2q$, which are $$1, 2, q, \text{ or } 2q.$$ In order to check whether $b$ is a primitive root, modulo a safe prime $p = 2q + 1$, we must check just three things: is $b \equiv 1$, is $b^2 \equiv 1$, or is $b^q \equiv 1$, mod $p$? If the answer to these three questions is NO, then $b$ is a primitive root mod $p$.

In [26]:
def is_primroot_safe(b,p):
    '''
    Checks whether b is a primitive root modulo p,
    when p is a safe prime.  If p is not safe,
    the results will not be good!
    '''
    q = (p-1) / 2   # q is the Sophie Germain prime
    if b%p == 1:  # Is the multiplicative order 1?
        return False
    if (b*b)%p == 1:  # Is the multiplicative order 2?
        return False
    if pow(b,q,p) == 1:  # Is the multiplicative order q?
        return False
    return True  # If not, then b is a primitive root mod p.

This would not be very useful if we couldn't find Sophie Germain primes. Fortunately, they are not so rare. The first few are 2, 3, 5, 11, 23, 29, 41, 53, 83, 89. It is expected, but unproven that there are infinitely many Sophie Germain primes. In practice, they occur fairly often. If we consider numbers of magnitude $N$, about $1 / \log(N)$ of them are prime. Among such primes, we expect about $1.3 / \log(N)$ to be Sophie Germain primes. In this way, we can expect to stumble upon Sophie Germain primes if we search for a bit (and if $\log(N)^2$ is not too large).

The code below tests whether a number $p$ is a Sophie Germain prime. We construct it by simply testing whether $p$ and $2p+1$ are both prime. We use the Miller-Rabin test (the code from the previous Python notebook) in order to test whether each is prime.

In [27]:
from random import randint # randint chooses random integers.

def Miller_Rabin(p, base):
    '''
    Tests whether p is prime, using the given base.
    The result False implies that p is definitely not prime.
    The result True implies that p **might** be prime.
    It is not a perfect test!
    '''
    result = 1
    exponent = p-1
    modulus = p
    bitstring = bin(exponent)[2:] # Chop off the '0b' part of the binary expansion of exponent
    for bit in bitstring: # Iterates through the "letters" of the string.  Here the letters are '0' or '1'.
        sq_result = result*result % modulus  # We need to compute this in any case.
        if sq_result == 1:
            if (result != 1) and (result != exponent):  # Note that exponent is congruent to -1, mod p.
                return False  # a ROO violation occurred, so p is not prime
        if bit == '0':
            result = sq_result 
        if bit == '1':
            result = (sq_result * base) % modulus
    if result != 1:
        return False  # a FLT violation occurred, so p is not prime.
    
    return True  # If we made it this far, no violation occurred and p might be prime.

def is_prime(p, witnesses=50):  # witnesses is a parameter with a default value.
    '''
    Tests whether a positive integer p is prime.
    For p < 2^64, the test is deterministic, using known good witnesses.
    Good witnesses come from a table at Wikipedia's article on the Miller-Rabin test,
    based on research by Pomerance, Selfridge and Wagstaff, Jaeschke, Jiang and Deng.
    For larger p, a number (by default, 50) of witnesses are chosen at random.
    '''
    if (p%2 == 0): # Might as well take care of even numbers at the outset!
        if p == 2:
            return True
        else:
            return False 
    
    if p > 2**64:  # We use the probabilistic test for large p.
        trial = 0
        while trial < witnesses:
            trial = trial + 1
            witness = randint(2,p-2) # A good range for possible witnesses
            if Miller_Rabin(p,witness) == False:
                return False
        return True
    
    else:  # We use a determinisic test for p <= 2**64.
        verdict = Miller_Rabin(p,2)
        if p < 2047:
            return verdict # The witness 2 suffices.
        verdict = verdict and Miller_Rabin(p,3)
        if p < 1373653:
            return verdict # The witnesses 2 and 3 suffice.
        verdict = verdict and Miller_Rabin(p,5)
        if p < 25326001:
            return verdict # The witnesses 2,3,5 suffice.
        verdict = verdict and Miller_Rabin(p,7)
        if p < 3215031751:
            return verdict # The witnesses 2,3,5,7 suffice.
        verdict = verdict and Miller_Rabin(p,11)
        if p < 2152302898747:
            return verdict # The witnesses 2,3,5,7,11 suffice.
        verdict = verdict and Miller_Rabin(p,13)
        if p < 3474749660383:
            return verdict # The witnesses 2,3,5,7,11,13 suffice.
        verdict = verdict and Miller_Rabin(p,17)
        if p < 341550071728321:
            return verdict # The witnesses 2,3,5,7,11,17 suffice.
        verdict = verdict and Miller_Rabin(p,19) and Miller_Rabin(p,23)
        if p < 3825123056546413051:
            return verdict # The witnesses 2,3,5,7,11,17,19,23 suffice.
        verdict = verdict and Miller_Rabin(p,29) and Miller_Rabin(p,31) and Miller_Rabin(p,37)
        return verdict # The witnesses 2,3,5,7,11,17,19,23,29,31,37 suffice for testing up to 2^64. 
    
In [28]:
def is_SGprime(p):
    '''
    Tests whether p is a Sophie Germain prime
    '''
    if is_prime(p):  # A bit faster to check whether p is prime first.
        if is_prime(2*p + 1):  # and *then* check whether 2p+1 is prime.
            return True

Let's test this out by finding the Sophie Germain primes up to 100, and their associated safe primes.

In [29]:
for j in range(1,100):
    if is_SGprime(j):
        print j, 2*j+1
2 5
3 7
5 11
11 23
23 47
29 59
41 83
53 107
83 167
89 179

Next, we find the first 100-digit Sophie Germain prime! This might take a minute!

In [30]:
test_number = 10**99 # Start looking at the first 100-digit number, which is 10^99.
while not is_SGprime(test_number):
    test_number = test_number + 1
print test_number
1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000044239

In the seconds or minutes your computer was running, it checked the primality of almost 90 thousand numbers, each with 100 digits. Not bad!

The Diffie-Hellman protocol

When we study protocols for secure communication, we must keep track of the communicating parties (often called Alice and Bob), and who has knowledge of what information. We assume at all times that the "wire" between Alice and Bob is tapped -- anything they say to each other is actively monitored, and is therefore public knowledge. We also assume that what happens on Alice's private computer is private to Alice, and what happens on Bob's private computer is private to Bob. Of course, these last two assumptions are big assumptions -- they point towards the danger of computer viruses which infect computers and can violate such privacy!

The goal of the Diffie-Hellman protocol is -- at the end of the process -- for Alice and Bob to share a secret without ever having communicated the secret with each other. The process involves a series of modular arithmetic calculations performed on each of Alice and Bob's computers.

The process begins when Alice or Bob creates and publicizes a large prime number p and a primitive root g modulo p. It is best, for efficiency and security, to choose a safe prime p. Alice and Bob can create their own safe prime, or choose one from a public list online, e.g., from the RFC 3526 memo. Nowadays, it's common to take p with 2048 bits, i.e., a prime which is between $2^{2046}$ and $2^{2047}$ (a number with 617 decimal digits!).

For the purposes of this introduction, we use a smaller safe prime, with about 256 bits. We use the SystemRandom functionality of the random package to create a good random prime. It is not so much of an issue here, but in general one must be very careful in cryptography that one's "random" numbers are really "random"! The SystemRandom function uses chaotic properties of your computer's innards in order to initialize a random number generator, and is considered cryptographically secure.

In [34]:
from random import SystemRandom # Import the necessary package.

r = SystemRandom().getrandbits(256)
print "The random integer is ",r
print "with binary expansion",bin(r)  #  r is an integer constructed from 256 random bits.
print "with bit-length ",len(bin(r)) - 2  # In case you want to check.  Remember '0b' is at the beginning.
The random integer is  45778098513937897927757416377531402057253432819631942923254557322726928037628
with binary expansion 0b110010100110101011111000100110110111111001010011100110010011101001000001011010000011100101011010001000010100000100010011110001111111101011101001010011101000101100101111110010010110001111010101111100001010111000110111101110000000001111001000011011011111100
with bit-length  255
In [40]:
def getrandSGprime(bitlength):
    '''
    Creates a random Sophie Germain prime p with about 
    bitlength bits.
    '''
    while True:
        p = SystemRandom().getrandbits(bitlength) # Choose a really random number.
        if is_SGprime(p):
            return p    

The function above searches and searches among random numbers until it finds a Sophie Germain prime. The (possibly endless!) search is performed with a while True: loop that may look strange. The idea is to stay in the loop until such a prime is found. Then the return p command returns the found prime as output and halts the loop. One must be careful with while True loops, since they are structured to run forever -- if there's not a loop-breaking command like return or break inside the loop, your computer will be spinning for a long time.

In [41]:
q = getrandSGprime(256) # A random ~256 bit Sophie Germain prime
p = 2*q + 1 # And its associated safe prime

print "p is ",p  # Just to see what we're working with.
print "q is ",q
p is  122707501782650049149641770078281232537400641066909203053473747068579132177267
q is  61353750891325024574820885039140616268700320533454601526736873534289566088633

Next we find a primitive root, modulo the safe prime p.

In [42]:
def findprimroot_safe(p):
    '''
    Finds a primitive root, 
    modulo a safe prime p.
    '''
    b = 2  # Start trying with 2.
    while True:  # We just keep on looking.
        if is_primroot_safe(b,p):
            return b
        b = b + 1 # Try the next base.  Shouldn't take too long to find one!
In [43]:
g = findprimroot_safe(p)
print g
2

The pair of numbers $(g, p)$, the primitive root and the safe prime, chosen by either Alice or Bob, is now made public. They can post their $g$ and $p$ on a public website or shout it in the streets. It doesn't matter. They are just tools for their secret-creation algorithm below.

Alice and Bob's private secrets

Next, Alice and Bob invent private secret numbers $a$ and $b$. They do not tell anyone these numbers. Not each other. Not their family. Nobody. They don't write them on a chalkboard, or leave them on a thumbdrive that they lose. These are really secret.

But they don't use their phone numbers, or social security numbers. It's best for Alice and Bob to use a secure random number generator on their separate private computers to create $a$ and $b$. They are often 256 bit numbers in practice, so that's what we use below.

In [45]:
a = SystemRandom().getrandbits(256)  # Alice's secret number
b = SystemRandom().getrandbits(256)  # Bob's secret number

print "Only Alice should know that a = %d"%(a)
print "Only Bob should know that b = %d"%(b)

print "But everyone can know p = %d and g = %d"%(p,g)
Only Alice should know that a = 114480947817489678722698132041311177496755779684221373891734317238736713023199
Only Bob should know that b = 97717430435546061486743031500808224958316500057061257004787704865563088657970
But everyone can know p = 122707501782650049149641770078281232537400641066909203053473747068579132177267 and g = 2

Now Alice and Bob use their secrets to generate new numbers. Alice computes the number $$A = g^a \text{ mod } p,$$ and Bob computes the number $$B = g^b \text{ mod } p.$$

In [46]:
A = pow(g,a,p)  # This would be computed on Alice's computer.
B = pow(g,b,p)  # This would be computed on Bob's computer.

Now Alice and Bob do something that seems very strange at first. Alice sends Bob her new number $A$ and Bob sends Alice his new number $B$. Since they are far apart, and the channel is insecure, we can assume everyone in the world now knows $A$ and $B$.

In [47]:
print "Everyone knows A = %d and B = %d."%(A,B)
Everyone knows A = 75920269441495140152557761658482862313742531041516025298431576783015820029854 and B = 68601777633459389691993124698389618527684118371645684321716024809600796197217.

Now Alice, on her private computer, computes $B^a$ mod $p$. She can do that because everyone knows $B$ and $p$, and she knows $a$ too.

Similarly, Bob, on his private computer, computes $A^b$ mod $p$. He can do that because everyone knows $A$ and $p$, and he knows $b$ too.

Alice and Bob do not share the results of their computations!

In [48]:
print pow(B,a,p)  # This is what Alice computes.
85348116833196709219191884338872059190642395205540062913372343974475490625996
In [49]:
print pow(A,b,p)  # This is what Bob computes.
85348116833196709219191884338872059190642395205540062913372343974475490625996

Woah! What happened? In terms of exponents, it's elementary. For $$B^a = (g^{b})^a = g^{ba} = g^{ab} = (g^a)^b = A^b.$$ So these two computations yield the same result (mod $p$, the whole way through).

In the end, we find that Alice and Bob share a secret. We call this secret number $S$. $$S = B^a = A^b.$$

In [50]:
S = pow(B,a,p)  # Or we could have used pow(A,b,p)
print S
85348116833196709219191884338872059190642395205540062913372343974475490625996

This common secret $S$ can be used as a key for Alice and Bob to communicate hereafter. For example, they might use $S$ (converted to a string, if needed) as the key for a Vigenère cipher, and chat with each other knowing that only they have the secret key to encrypt and decrypt their messages.

In [59]:
#  We use the single-quotes for a long string, that occupies multiple lines.
#  The backslash at the end of the line tells Python to ignore the newline character.
#  Imagine that Alice has a secret message she wants to send to Bob.  
#  She writes the plaintext on her computer. 

plaintext = '''Did you hear that the American Mathematical Society has an annual textbook sale? \
            It's 40 percent off for members and 25 percent off for everyone else.'''
In [60]:
# Now Alice uses the secret S (as a string) to encrypt.  
ciphertext = Vigenere_cipher(plaintext, str(S))
print ciphertext
# Alice sends the following ciphertext to Bob, over an insecure channel.
|[email protected]>G=7:GW 2JB7G2K>65FU&@9?:EJV;;ES5FR5CBO9ATL;MN3FBAYM8E;rTXQQVXSSQYVWP#GXMQmaXI:F7>[email protected]<@[email protected]@[email protected]?5VdnQD9J69CHYG;;W<DLP<I;LSFG;S:EE7d
In [61]:
# When Bob receives the ciphertext, he decodes it with the secret S again.
print Vigenere_decipher(ciphertext, str(S))
Did you hear that the American Mathematical Society has an annual textbook sale?             It's 40 percent off for members and 25 percent off for everyone else.

To have confidence in this protocol, one needs to be convinced that their secret is truly secret! The public has a lot of information: they know

  1. the prime $p$,
  2. the primitive root $g$,
  3. the number $A = g^a$ (mod $p$), and
  4. the number $B = g^b$ (mod $p$).

If the public could figure out either $a$ or $b$, they could figure out the secret (by raising $A^b$ or $B^a$ like Alice and Bob did).

This is called the discrete logarithm problem. If we know the value of $g^a$ mod $p$, can we figure out the possible value(s) of $a$? If this were ordinary arithmetic, we would say that $a = \log_g(A)$. But this is modular arithmetic, and there's no easy way to figure out such logarithms. The values of $g^a$ tend to bounce all over the place, mod $p$, especially since we chose $a$ to be pretty large (256 bits!).

The security of the Diffie-Hellman protocol, i.e., the security of Alice and Bob's shared secret, depends on the difficulty of the discrete logarithm problem. When $p$ is a large (e.g. 2048 bits) safe prime, and $a$ and $b$ are suitably large (roughly 256 bits), there seems to be no way to solve the discrete logarithm problem mod $p$ in any reasonable amount of time. Someday we might have quantum computers to quickly solve discrete logarithm problems, and the Diffie-Hellman protocol will not be secure. But for now, Alice and Bob's secret key seems safe.

Exercises

  1. How many Sophie Germain primes are there between 1 and 1000000? What proportion of primes in this range are Sophie Germain primes?

  2. It is expected that there are infinitely primes $p$ such that 2 is a primitive root mod $p$. Study the density of such primes. For example, among the primes up to 1000, how often is $2$ a primitive root? Does this density seem to change?

  3. Adapt Diffie-Hellman to work with a group of three parties who wish to share a common secret. Hint: the common secret will have the form $g^{abc}$, and other exponents like $g^a$, $g^b$, $g^c$, $g^{ab}$, $g^{bc}$, $g^{ac}$ will be public information.

  4. Sadly, Alice and Bob have agreed to use the primitive root $g = 3$ and the prime $p = 65537$. Listening into their conversation, you intercept the following: $A = 40360$ and $B = 21002$ and the ciphertext is $;6HWD;P5LVJ99W+EH9JVx=I:V7ESpGC^. If you know that they use a protocol with a Vigenère cipher, with key equal the string associated to their secret $S$, what is the plaintext message? Hint: you should be able to solve the discrete logarithm problem with a brute force attack.