Silicon Forest Math Series | RSA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Silicon Forest Math Series<br/>[Oregon Curriculum Network](http://4dsolutions.net/ocn/)\n",
    "\n",
    "## Introduction to Public Key Cryptography\n",
    "\n",
    "![Portland, Oregon](https://c6.staticflickr.com/6/5605/30351762741_ddbfbfdf71.jpg)\n",
    "\n",
    "Here in the Silicon Forest, we do not expect everyone to become a career computer programmer.  \n",
    "\n",
    "We do expect a lot of people will wish to program at some time in their career.  \n",
    "\n",
    "Coding skills give you the power to control machines and you might find appropriate and life-enhancing uses for this type of power.\n",
    "\n",
    "To help you get the flavor of coding, we leverage concepts we expect you're already getting through your math courses.  \n",
    "\n",
    "In moving from pre-computer math, to computer math, and back again, we develop important conceptual bridges.\n",
    "\n",
    "### Generating the Prime Numbers\n",
    "\n",
    "Lets look at a first such concept, that of a prime number.  \n",
    "\n",
    "The Fundamental Theorem of Arithmetic says every positive integer distills into a unique list of prime number factors.  Duplicates are allowed.\n",
    "\n",
    "But what are primes in the first place?  Numbers with no factors other than themselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31,\n",
      " 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,\n",
      " 79, 83, 89, 97, 101, 103, 107, 109,\n",
      " 113]\n"
     ]
    }
   ],
   "source": [
    "import pprint\n",
    "\n",
    "def primes():\n",
    "    \"\"\"generate successive prime numbers (trial by division)\"\"\"\n",
    "    candidate = 1\n",
    "    _primes_so_far = [2]     # first prime, only even prime\n",
    "    yield _primes_so_far[0]  # share it!\n",
    "    while True:\n",
    "        candidate += 2    # check odds only from now on\n",
    "        for prev in _primes_so_far:\n",
    "            if prev**2 > candidate:\n",
    "                yield candidate  # new prime!\n",
    "                _primes_so_far.append(candidate)\n",
    "                break\n",
    "            if not divmod(candidate, prev)[1]: # no remainder!\n",
    "                break                          # done looping\n",
    "                \n",
    "p = primes()  # generator function based iterator\n",
    "pp = pprint.PrettyPrinter(width=40, compact=True)\n",
    "pp.pprint([next(p) for _ in range(30)])  # next 30 primes please!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above algorithm is known as \"trial by division\".  \n",
    "\n",
    "Keep track of all primes discovered so far, and test divide them, in increasing order, into a candidate number, until: \n",
    "\n",
    "(A) either one of the primes goes evenly, in which case move on to the next odd \n",
    "\n",
    "or \n",
    "\n",
    "(B) until we know our candidate is a next prime, in which case yield it and append it to the growing list.\n",
    "\n",
    "If we get passed the 2nd root of the candidate, we conclude no larger factor will work, as we would have encountered it already as the smaller of the factor pair. \n",
    "\n",
    "Passing this 2nd root milestone triggers plan B.  Then we advance to the next candidate, ad infinitum.  \n",
    "\n",
    "Python pauses at each yield statement however, handing control back to the calling sequence, in this case a \"list comprehension\" containing a next() function for advancing to the next yield.\n",
    "\n",
    "![Brands of the PNW](https://c5.staticflickr.com/4/3006/3000202748_4a6a67e878.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coprimes, Totatives, and the Totient of a Number\n",
    "\n",
    "From here, we jump to the idea of numbers being coprime to one another.  A synonym for coprime is \"stranger.\"  Given two ordinary positive integers, they're strangers if they have no prime factors in common.  For that to be true, they'd have no shared factors at all (not counting 1).\n",
    "\n",
    "Guido van Rossum, the inventor of Python, gives us a pretty little implementation of what's known as Euclid's Method, an algorithm that's thousands of years old. It'll find the largest factor any two numbers have in common (gcd = \"greatest common divisor\").\n",
    "\n",
    "Here it is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n",
      "4\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "def gcd(a, b):\n",
    "    while b:\n",
    "        a, b = b, a % b\n",
    "    return a\n",
    "\n",
    "print(gcd(81, 18))\n",
    "print(gcd(12, 44))\n",
    "print(gcd(117, 17)) # strangers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does Euclid's Method work? That's a great question and one your teacher should be able to explain. First see if you might figure it out for yourself...  \n",
    "\n",
    "Here's one explanation:\n",
    "\n",
    "If a smaller number divides a larger one without remainder then we're done, and that will always happen when that smaller number is 1 if not before.  \n",
    "\n",
    "If there *is* a remainder, what then?  Lets work through an example.\n",
    "\n",
    "81 % 18 returns a remainder of 9 in the first cycle.  18 didn't go into 81 evenly but if another smaller number goes into both 9, the remainder, and 18, then we have our answer.  \n",
    "\n",
    "9 itself does the trick and we're done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "print(81 % 18) # 18 goes into \n",
    "print(18 % 9)  # so the new b becomes the answer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we had asked for gcd(18, 81) instead? 18 is the remainder (no 81s go into it) whereas b was 81, so the while loop simply flips the two numbers around to give the example above.\n",
    "\n",
    "The gcd function now gives us the means to compute *totients* and *totatives* of a number.  The totatives of N are the strangers less than N, whereas the totient is the number of such strangers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Totient of  100: 40\n",
      "Totient of 1000: 400\n"
     ]
    }
   ],
   "source": [
    "def totatives(N):\n",
    "    # list comprehension!\n",
    "    return [x for x in range(1,N) if gcd(x,N)==1]  # strangers only\n",
    "    \n",
    "def T(N):\n",
    "    \"\"\"\n",
    "    Returns the number of numbers between (1, N) that \n",
    "    have no factors in common with N:  called the \n",
    "    'totient of N' (sometimes phi is used in the docs)\n",
    "    \"\"\"\n",
    "    return len(totatives(N)) # how many strangers did we find?\n",
    "\n",
    "print(\"Totient of  100:\", T(100))\n",
    "print(\"Totient of 1000:\", T(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where to go next is in the direction of Euler's Theorem, a generalization of Fermat's Little Theorem.  The built-in pow(m, n, N) function will raise m to the n modulo N in an efficient manner. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Totient of 17: 16\n",
      " 1 [1]\n",
      " 8 [1, 2, 4, 8, 16, 15, 13, 9]\n",
      "16 [1, 3, 9, 10, 13, 5, 15, 11, 16, 14, 8, 7, 4, 12, 2, 6]\n",
      " 4 [1, 4, 16, 13]\n",
      "16 [1, 5, 8, 6, 13, 14, 2, 10, 16, 12, 9, 11, 4, 3, 15, 7]\n",
      "16 [1, 6, 2, 12, 4, 7, 8, 14, 16, 11, 15, 5, 13, 10, 9, 3]\n",
      "16 [1, 7, 15, 3, 4, 11, 9, 12, 16, 10, 2, 14, 13, 6, 8, 5]\n",
      " 8 [1, 8, 13, 2, 16, 9, 4, 15]\n",
      " 8 [1, 9, 13, 15, 16, 8, 4, 2]\n",
      "16 [1, 10, 15, 14, 4, 6, 9, 5, 16, 7, 2, 3, 13, 11, 8, 12]\n",
      "16 [1, 11, 2, 5, 4, 10, 8, 3, 16, 6, 15, 12, 13, 7, 9, 14]\n",
      "16 [1, 12, 8, 11, 13, 3, 2, 7, 16, 5, 9, 6, 4, 14, 15, 10]\n",
      " 4 [1, 13, 16, 4]\n",
      "16 [1, 14, 9, 7, 13, 12, 15, 6, 16, 3, 8, 10, 4, 5, 2, 11]\n",
      " 8 [1, 15, 4, 9, 16, 2, 13, 8]\n",
      " 2 [1, 16]\n"
     ]
    }
   ],
   "source": [
    "def powers(N):\n",
    "    totient = T(N)\n",
    "    print(\"Totient of {}:\".format(N), totient)\n",
    "    for t in totatives(N):\n",
    "        values = [pow(t, n, N) for n in range(totient + 1)]\n",
    "        cycle = values[:values.index(1, 1)]  # first 1 after initial 1\n",
    "        print(\"{:>2}\".format(len(cycle)), cycle)\n",
    "    \n",
    "powers(17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above we see repeating cycles of numbers, with the length of the cycles all dividing 16, the totient of the prime number 17. \n",
    "\n",
    "pow(14, 2, 17) is 9, pow(14, 3, 17) is 7, and so on, coming back around the 14 at pow(14, 17, 17) where 17 is 1 modulo 16.\n",
    "\n",
    "Numbers raised to any kth power modulo N, where k is 1 modulo the totient of N, end up staying the same number. For example, pow(m, (n \\* T(N)) + 1, N) == m for any n."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 1\n",
      "2 2\n",
      "3 3\n",
      "4 4\n",
      "5 5\n",
      "6 6\n",
      "7 7\n",
      "8 8\n",
      "9 9\n",
      "10 10\n",
      "11 11\n",
      "12 12\n",
      "13 13\n",
      "14 14\n",
      "15 15\n",
      "16 16\n"
     ]
    }
   ],
   "source": [
    "from random import randint\n",
    "\n",
    "def check(N):\n",
    "    totient = T(N)\n",
    "    for t in totatives(N):\n",
    "        n = randint(1, 10)\n",
    "        print(t, pow(t, (n * totient) + 1, N))\n",
    "\n",
    "check(17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In public key cryptography, RSA in particular, a gigantic composite N is formed from two primes p and q.  \n",
    "\n",
    "N's totient will then be (p - 1) \\* (q - 1). For example if N = 17 \\* 23 (both primes) then T(N) = 16 * 22."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 17\n",
    "q = 23\n",
    "T(p*q) == (p-1)*(q-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this totient, we'll be able to find pairs (e, d) such that (e * d) modulo T(N) == 1. \n",
    "\n",
    "We may find d, given e and T(N), by means of the Extended Euclidean Algorithm (xgcd below).\n",
    "\n",
    "Raising some numeric message m to the eth power modulo N will encrypt the message, giving c.  \n",
    "\n",
    "Raising the encrypted message c to the dth power will cycle it back around to its starting value, thereby decrypting it.\n",
    "\n",
    "c = pow(m, e, N)\n",
    "\n",
    "m = pow(c, d, N)\n",
    "\n",
    "where (e * d) % T(N) == 1.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "p = 37975227936943673922808872755445627854565536638199\n",
    "q = 40094690950920881030683735292761468389214899724061\n",
    "RSA_100 = p * q\n",
    "totient = (p - 1) * (q - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "# https://en.wikibooks.org/wiki/\n",
    "# Algorithm_Implementation/Mathematics/\n",
    "# Extended_Euclidean_algorithm\n",
    "\n",
    "def xgcd(b, n):\n",
    "    x0, x1, y0, y1 = 1, 0, 0, 1\n",
    "    while n != 0:\n",
    "        q, b, n = b // n, n, b % n\n",
    "        x0, x1 = x1, x0 - q * x1\n",
    "        y0, y1 = y1, y0 - q * y1\n",
    "    return  b, x0, y0\n",
    "\n",
    "# x = mulinv(b) mod n, (x * b) % n == 1\n",
    "def mulinv(b, n):\n",
    "    g, x, _ = xgcd(b, n)\n",
    "    if g == 1:\n",
    "        return x % n\n",
    "\n",
    "e = 3\n",
    "d = mulinv(e, totient)\n",
    "print((e*d) % totient)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22640069159164324404987585908\n"
     ]
    }
   ],
   "source": [
    "import binascii\n",
    "m = int(binascii.hexlify(b\"I'm a secret\"), 16)\n",
    "print(m)  # decimal encoding of byte string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11604682090980443295874366312042718286665199375418874910814662777671879383859486933312\n"
     ]
    }
   ],
   "source": [
    "c = pow(m, e, RSA_100) # raise to eth power\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22640069159164324404987585908\n"
     ]
    }
   ],
   "source": [
    "m = pow(c, d, RSA_100) # raise to dth power\n",
    "print(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "b\"I'm a secret\""
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "binascii.unhexlify(hex(m)[2:]) # m is back where we started."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What makes RSA hard to crack is that although N is public, d is not, and N's factors p and q have been thrown away. \n",
    "\n",
    "Because factoring N back into p, q is a super hard problem, if N is large enough, RSA remains a secure algorithm, and a favorite one, now that the patent has expired.\n",
    "\n",
    "The idea is when you want to encrypt a message for Alice, look up her public key N.  Only she has private key d, derived from her N at the time.  \n",
    "\n",
    "You may sign your message by raising it to your own secret key power, modulo your own public N, and she'll know for sure the message is from you, given she'll have your public key to decrypt it, once her own secret d has been applied.\n",
    "\n",
    "![Go By Train](https://c7.staticflickr.com/4/3343/3228050854_13115a5faf.jpg)\n",
    "\n",
    "### For Further Reading\n",
    "\n",
    "Linked Jupyter Notebooks:\n",
    "\n",
    "- [Generators and Coroutines](http://nbviewer.jupyter.org/github/4dsolutions/Python5/blob/master/Generators_and_Coroutines.ipynb) -- introduces Quadrays\n",
    "- [Pi Day Fun](http://nbviewer.jupyter.org/github/4dsolutions/Python5/blob/master/Pi%20Day%20Fun.ipynb) -- Ramanujan!\n",
    "- [Composition of Functions](http://nbviewer.jupyter.org/github/4dsolutions/Python5/blob/master/ComposingFunctions.ipynb) -- advanced Python\n",
    "- [STEM Mathematics](http://nbviewer.jupyter.org/github/4dsolutions/Python5/blob/master/STEM%20Mathematics.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}