#!/usr/bin/env python
# coding: utf-8

# # 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](http://illustratedtheoryofnumbers.com).

# ## Table of Contents
# 
# - [Ciphers](#cipher)
# - [Key exchange](#keyexchange)
# 

# <a id='cipher'></a>

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


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


# Since computers store data in binary, the designers of early computers (1960s) created a code called [**ASCII**](https://en.wikipedia.org/wiki/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](http://www.asciitable.com/) 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](https://unicode-table.com/en/) 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)


# In[5]:


ord('A')


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


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


# In[9]:


inrange(17,5,50)


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


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


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


# ### 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](https://en.wikipedia.org/wiki/Vigen%C3%A8re_cipher#Cryptanalysis) 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](https://en.wikipedia.org/wiki/Advanced_Encryption_Standard) and [3DES](https://en.wikipedia.org/wiki/Triple_DES) 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!!


# 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')


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

# <a id='keyexchange'></a>

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


# 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](https://fineartamerica.com/featured/epicycles-modulo-37-martin-weissman.html) (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


# 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


# 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](https://tools.ietf.org/html/rfc3526).  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.


# 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


# 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


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


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


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


# In[49]:


print pow(A,b,p)  # This is what Bob computes.


# 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


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


# In[61]:


# When Bob receives the ciphertext, he decodes it with the secret S again.
print Vigenere_decipher(ciphertext, str(S))


# 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](https://en.wikipedia.org/wiki/Artin%27s_conjecture_on_primitive_roots) 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.
# 
# 

# In[ ]: