{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RIVEST SHAMIR ADLEMAN (RSA)\n",
    "\n",
    "Let $p,q$ be two large primes of comparable size. Put $n = pq$. This computation takes about $O(\\log_2(n)^2)$ many flops if you use elementary school method of multiplying these numbers (it is actually $O(\\log_2(n)\\log\\log_2(n))$ many flops if you employ Fast Fourier Transform multiplication algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However trying to factor $n$ into $p$ and $q$ using the brute force method will take $O(\\max\\{p,q\\})$ modular divisions. which is on the order of $O(\\sqrt{n})$ if $p$ and $q$ are close enough to one another. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "48782919860664620234782399461351856881494516091785396502382566480860785545607256493096155801868780530190395411670"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "163298476128374618237491823746912837468127467125648723654817356418765*298734691328746192837461923874619287461983278"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 2, 3, 13]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def BruteForceFactor(n):\n",
    "    primelist = []\n",
    "    d = 2\n",
    "    while n!= 1:\n",
    "        if d*d>n:\n",
    "            primelist.append(n)\n",
    "            return primelist\n",
    "        elif n%d == 0:\n",
    "            primelist.append(d)\n",
    "            n = n//d\n",
    "        else:\n",
    "            d = d + 1\n",
    "    return primelist\n",
    "\n",
    "BruteForceFactor(156)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 3, 3607, 3803]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(123456789)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1234567891]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(1234567891)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 2, 2, 2, 2, 3, 5, 323339, 3347983, 2375923237887317]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(1234567891011121314151617181920)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Euler $\\phi$ Function\n",
    "\n",
    "We are interested in the order of the group $(\\mathbb{Z}/n\\mathbb{Z})^*$. That is, what is the number $$\\phi(n) := \\#\\{1\\leq a \\leq n: \\operatorname{gcd}(a,n) = 1\\}.$$\n",
    "\n",
    "If $n = p$ is a prime, we know that every integer less than $p$ is relatively prime to $p$, and hence $\\phi(p) = p-1$. \n",
    "\n",
    "If $n = p^k$ a power of a prime, then the only integers that are not relatively prime to $n$ are multiples of $p$, and there are $n/p = p^{k-1}$ many of them. Hence $\\phi(p^k) = p^k - p^{k-1} = p^k(1 - 1/p)$.\n",
    "\n",
    "Lastly by Chinese Remainder theorem we note that if $n = ab$ with $\\operatorname{gcd}(a,b) = 1$, then there is an isomorphism of groups\n",
    "$$\n",
    "    (\\mathbb{Z}/n\\mathbb{Z})^* \\cong (\\mathbb{Z}/a\\mathbb{Z})^* \\times (\\mathbb{Z}/b\\mathbb{Z})^*.\n",
    "$$\n",
    "There are $\\phi(n)$ elements on the left, and $\\phi(a)\\phi(b)$ elements on the right.\n",
    "\n",
    "Because of this we say that $\\phi$ is a multiplicative function and for a general $n$ of the form,\n",
    "$$\n",
    "n = p_1^{k_1} p_2^{k_2} \\cdots p_r^{k_r}\n",
    "$$\n",
    "we deduce that\n",
    "$$\n",
    " \\phi(n) = (p_1^{k_1} - p_1^{k_1-1}) \\cdots (p_r^{k_r} - p_r^{k_r-1}) = n \\prod_{p | n}  \\left(1 - \\frac 1p\\right)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def EulerPhi(n, factors = None):\n",
    "    if factors == None:\n",
    "        primelist = BruteForceFactor(n)\n",
    "    else:\n",
    "        primelist= factors\n",
    "    result = 1\n",
    "    while len(primelist)>1:\n",
    "        ult = primelist[-1]\n",
    "        penult = primelist[-2]\n",
    "        if ult == penult:\n",
    "            result = result*ult\n",
    "            primelist.pop()\n",
    "        else:\n",
    "            result = result*(ult -1)\n",
    "            primelist.pop()\n",
    "    return result*(primelist[-1]-1)\n",
    "\n",
    "EulerPhi(101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EulerPhi(25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1234567890"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EulerPhi(1234567891)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EulerPhi(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fastexp(x,n,m):\n",
    "    z = 1\n",
    "    while n!=0:\n",
    "        if n%2==1:\n",
    "            z = (z*x)%m\n",
    "        x = (x*x)%m\n",
    "        n = n//2\n",
    "    return z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponentiation for hiding information\n",
    "\n",
    "Let $p = 1234567891$ be a prime, and choose a number $m$ as the message to be sent.\n",
    "\n",
    "\n",
    "Computing $m^{123876} \\pmod{p}$ gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1010044275"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 1234567891\n",
    "m = 58008 #a message for the other side\n",
    "c = fastexp(m,1234567, p)\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1f9a79f7710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "ks = range(2,640)\n",
    "y = [fastexp(m,k,p) for k in ks]\n",
    "s = [5 for i in ks]\n",
    "\n",
    "plt.scatter(ks,y,s)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cracking $m$ from $c$ with what is above.\n",
    "\n",
    "Actually it is not very hard to solve $m$ given $c$ using the fact that $p$ is prime, and we have $$|(\\mathbb{Z}/p\\mathbb{Z})^*| = p-1$$ for a prime number $p$. Hence we have that $$ a^{p-1} \\equiv 1 \\pmod{p}.$$ This is  because the order of any element $a$ divides the order of the group (that is $p-1$). So $a^p \\equiv a \\pmod{p}$ for all integers $a$, or in fact $a^{e} \\equiv a \\pmod p$ for any $e \\equiv 1 \\pmod {p-1}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let us bring the things we have written again. \n",
    "\n",
    "def bezout(n,m):\n",
    "    (g,u,v) = (n,1,0)\n",
    "    (e,s,t) = (m,0,1)\n",
    "    r = g%e\n",
    "    while r != 0:\n",
    "        q = g//e\n",
    "        (eprime, sprime,tprime) = (e,s,t)\n",
    "        (e,s,t) = (r,u - q*s,v  - t*q)\n",
    "        (g,u,v) = (eprime,sprime,tprime)\n",
    "        r = g%e\n",
    "    return (e,s,t)\n",
    "\n",
    "def gcd(n,m):\n",
    "    (g,u,v) = bezout(n,m)\n",
    "    return g\n",
    "\n",
    "def inverse(x,n):\n",
    "    (g,u,v) = bezout(x,n)\n",
    "    if g!=1:\n",
    "        return None\n",
    "    else:\n",
    "        return (u%n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18033013"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e = 1234567\n",
    "d = inverse(e, p-1)\n",
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(e*d)%1234567890"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9798"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 10403\n",
    "fastexp(fastexp(fastexp(fastexp(fastexp(fastexp(7580,5,10403),6,n),7,n),8,n),9,n),10,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "58008"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fastexp(c,d,p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "101"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(9797,10403)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2^ 0  = 1 mod 7\n",
      "2^ 1  = 2 mod 7\n",
      "2^ 2  = 4 mod 7\n",
      "2^ 3  = 1 mod 7\n",
      "2^ 4  = 2 mod 7\n",
      "2^ 5  = 4 mod 7\n",
      "2^ 6  = 1 mod 7\n",
      "2^ 7  = 2 mod 7\n",
      "2^ 8  = 4 mod 7\n",
      "2^ 9  = 1 mod 7\n",
      "2^ 10  = 2 mod 7\n",
      "2^ 11  = 4 mod 7\n",
      "2^ 12  = 1 mod 7\n",
      "2^ 13  = 2 mod 7\n",
      "2^ 14  = 4 mod 7\n"
     ]
    }
   ],
   "source": [
    "for k in range(15):\n",
    "    print(\"2^\", k,\" =\", fastexp(2,k,7), 'mod 7')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got $m$ back!!! Why was that, well here is the reason\n",
    "$$ c^{d} \\equiv (m^e)^d \\equiv m^{ed} \\pmod p$$ and $ed \\equiv 1 \\pmod {p-1}$ and therefore $$c^{d} \\equiv m^{ed} \\equiv m \\pmod {p}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What if the modulus $n$ is not prime. \n",
    "\n",
    "Then we can do calculation in $G = (\\mathbb{Z}/n\\mathbb{Z})^*$ however we don't necessarily know $|G| = \\phi(n)$. If we do, then the above method also solves the equation for $m$:\n",
    "$$ m^e = \\equiv c \\pmod{n}.$$\n",
    "\n",
    "\n",
    "Knowing $\\phi(n)$ is easy if you know the factorization of $n$, which means that we are looking for numbers that are hard to factor. The smallest such numbers are of the form\n",
    "$$ n = pq$$\n",
    "for primes $p$ and $q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[12345678910111]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q = 12345678910111\n",
    "BruteForceFactor(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15241578775018915845901"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = p*q\n",
    "n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dont even try BruteForceFactor(n), it takes too long!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import randint as rdint\n",
    "\n",
    "def RSAkeys(n,factors):\n",
    "    e = rdint(n//4,n//2)\n",
    "    while not all([gcd(e, p-1)==1 for p in factors]):\n",
    "        e = rdint(n//4,n//2)\n",
    "    order = EulerPhi(n,factors)\n",
    "    d = inverse(e,order)\n",
    "    return (e,d)\n",
    "\n",
    "def RSAencrypt(m, n, e):\n",
    "    c = fastexp(m,e,n)\n",
    "    return c\n",
    "\n",
    "def RSAdecrypt(c,n, d):\n",
    "    m = fastexp(c,d,n)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3937139939893578395651, 15184922774872748407451)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(e,d) = RSAkeys(n,[p,q])\n",
    "(e,d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15241578762672002367900"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(p-1)*(q-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(e*d)%((p-1)*(q-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7399584266783735548498"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = RSAencrypt(58008, n, e)\n",
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "58008"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RSAdecrypt(c,n,d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The RSA protocole:\n",
    "\n",
    "Alice chooses two primes $p,q$ and computes $n = pq$. Choses an $e$ such that $\\operatorname{gcd}(e,(p-1)(q-1))$. Computes $d$ such that $ed \\equiv 1 \\pmod{\\phi(n)}$. \n",
    "\n",
    "Alice publishes $e, n$ and keeps $d$ to herself.\n",
    "\n",
    "Bob has a message $m$, and computes $c \\equiv m^e \\pmod{n}$ at the privacy of his own office and sends Alice the cipher $c$. \n",
    "\n",
    "Alice (since only she has knowledge of the secret key $e$) computes $c^e \\pmod{n}$ which is equivalent to $m$. \n",
    "\n",
    "## Factorization is as hard as finding $\\phi(n)$.\n",
    "\n",
    "Note that if $n = pq$ then $\\phi(n) = (p-1)(q-1) = pq - (p + q) + 1= n - (p + q) + 1$. Thus if we know $n$ and $\\phi(n)$ then we can easily calculate $p + q = n- \\phi(n) +1$, and knowing both $pq$ and $p+q$ we solve the quadratic equation $$X^2 + (p + q) X + pq = 0$$ in order to get $p$ and $q$. \n",
    "\n",
    "## The solution is unique\n",
    "\n",
    "We are given $c \\equiv m^{e} \\pmod{n}$ and for $\\operatorname{gcd}(e, (p-1)(q-1))=1$. So we are trying to find the $e$-th root of $c$ modulo $n$. We said that with $d$ given as $ed \\equiv 1 \\pmod \\phi(n)$ then $x \\equiv c^d\\pmod{n}$ is a solution to $x^e \\equiv c \\pmod{n}$. It is also the only solution:\n",
    "\n",
    "Assume $x$ were a solution, and that $de = 1 +  k\\phi(n) $ for some integer $k$, \n",
    "$$ \\begin{align*} x &\\equiv x^{de - k\\phi(n) }\\pmod{n} \\\\ &\\equiv (x^e)^d (x^{\\phi(n)})^{-k} \\pmod{n}\\\\\n",
    "&\\equiv c^d 1^{-k} \\equiv c^d \\pmod{n}\\end{align*} $$\n",
    "\n",
    "<i>(Note that if the exponent $e$ is not relatively prime to the order of the group, it is not necessarily true that $x^e \\equiv c \\pmod{n}$ is solvable for generic $c$. For example not every integer modulo $p$ is a square.)</i>\n",
    "\n",
    "\n",
    "## Is this really as hard as factorization\n",
    "\n",
    "It is generally assumed (especially by public) that the difficulty of breaking RSA is equivalent to factorization. Technically what one needs is to solve the equation $x^e \\equiv c \\pmod{n}$, which may or may not be strictly easier than the problem of factorizing $n$. The expert opinion is a bit more complicated:\n",
    "\n",
    "- D. Boneh and R. Venkatesan. <i><b> Breaking RSA may not be equivalent to factoring</b></i>. In Advances in Cryptology—EUROCRYPT ’98 (Espoo), volume 1403 of Lecture Notes in Comput. Sci., pages 59–71. Springer, Berlin, 1998.\n",
    "- Brown, D.R.L. J Cryptol <i><b>Breaking RSA may be as difficult as factoring.</b></i> Journal of Cryptology January 2016, Volume 29, Issue 1, pp 220–241.\n",
    "-     Divesh Aggarwal, Ueli Maurer, <i> <b>Breaking RSA Generically Is Equivalent to Factoring </b></i> Joux A. (eds) Advances in Cryptology - EUROCRYPT 2009. EUROCRYPT 2009. Lecture Notes in Computer Science, vol 5479. Springer, Berlin, Heidelberg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FERMAT's Test for Primality\n",
    "\n",
    "If $p$ is a prime number then $(\\mathbb{Z}/p\\mathbb{Z})^*$ has $p-1$ integers and therefore the order of any nonzero element $a \\pmod{p}$ divides $p-1$. In particular,\n",
    "$$ a^{p-1} \\equiv 1 \\pmod{p}.$$\n",
    "Equivalently\n",
    "$$ a^{p}\\equiv  a \\pmod{p}$$ for all integers $a$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "211367123"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 329487987\n",
    "fastexp(2,n, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 19, 5780491]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $2^n \\not\\equiv 2 \\pmod{n}$ we know that $n$ is not a prime. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 12801\n",
    "fastexp(2,n,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 17, 251]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9693"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fastexp(3,n,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen even though $n = 12801$ is not a prime, $2^n \\equiv 2 \\pmod{n}$. Yet the base $3$ comes to save the day, and is a witness that $n$ is composite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def CompositeWitness(a,n):\n",
    "    return fastexp(a,n,n)!=(a%n)\n",
    "\n",
    "CompositeWitness(2,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CompositeWitness(3,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see of all the composite numbers, how many of them are caught by various Fermat witnesses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4, 6, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 22, 24, 25, 26, 27, 28, 30, 32]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def CompositeList(N):\n",
    "    return list(filter(lambda n: len(BruteForceFactor(n))>1, range(2,N)))\n",
    "\n",
    "comp = CompositeList(10000)\n",
    "len(comp)\n",
    "comp[:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[341,\n",
       " 561,\n",
       " 645,\n",
       " 1105,\n",
       " 1387,\n",
       " 1729,\n",
       " 1905,\n",
       " 2047,\n",
       " 2465,\n",
       " 2701,\n",
       " 2821,\n",
       " 3277,\n",
       " 4033,\n",
       " 4369,\n",
       " 4371,\n",
       " 4681,\n",
       " 5461,\n",
       " 6601,\n",
       " 7957,\n",
       " 8321,\n",
       " 8481,\n",
       " 8911]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(filter(lambda n: not CompositeWitness(2,n),comp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[561, 1105, 1729, 2465, 2701, 2821, 6601, 8911]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(filter(lambda n: not (CompositeWitness(2,n) or CompositeWitness(3,n)),comp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[561, 1105, 1729, 2465, 2821, 6601, 8911]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def CarmichaelNumbers(N):\n",
    "    #Returns all the composite numbers n less than N such that\n",
    "    #All a < n are non-witnesses for the compositeness of n\n",
    "    comp = CompositeList(N)\n",
    "    a = 2\n",
    "    while (a < comp[-1]) and len(comp)>0:\n",
    "        comp = list(filter(lambda n: not CompositeWitness(a,n), comp))\n",
    "        a = a + 1\n",
    "    return list(comp)\n",
    "\n",
    "CarmichaelNumbers(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Carmichael Numbers are the wrenches in this scheme to determine primes.\n",
    "\n",
    "Even though $561 = 3\\times 11\\times 17$, it is true that $a^n \\equiv a \\pmod{n}$ for all integers $a$. So we have a sufficient method for testing primality.\n",
    "\n",
    "<b>Theorem.</b> All Charmichael numbers are of the form $n = p_1 p_2 \\cdots p_k$ with $k>3$ such that the prime divisors $p_i$ are all distinct and $(p_i -1) | (n-1)$ for all $i$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 11, 17]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(561)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5412\n"
     ]
    }
   ],
   "source": [
    "X = range(10**8,10**8 + 100000)\n",
    "X = filter(lambda n: not(CompositeWitness(2,n) or CompositeWitness(3,n)),X)\n",
    "X = list(X)\n",
    "print(len(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100017223 is not a prime, but pretends to be\n"
     ]
    }
   ],
   "source": [
    "for p in X:\n",
    "    if len(BruteForceFactor(p))>1:\n",
    "        print(str(p) + ' is not a prime, but pretends to be')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Miller Rabin Test Comes to the Rescue\n",
    "\n",
    "Let us test whether $n$ is a prime using the following algorithm. First factorize as $$n -1 = 2^k q$$ with $q$ an odd integer. Let $L_n$ be the set of nonzero integers $a$ less than $n$ such that either $$a^{2^kq} \\equiv 1 \\pmod{n}$$ or in the set $$a^q, a^{2q}, a^{2^2q}, \\ldots, a^{2^{k-1}q}$$ of successive squares they are all congruent to $1$ or one element is congruent to $-1$ modulo $n$. \n",
    "\n",
    "If $n$ is prime then $L_n = (\\mathbb{Z}/N\\mathbb{Z})^*$. \n",
    "\n",
    "If $n$ is not prime then $|L_n| \\leq (n-1)/4$. This means that for $n$ composite $\\%75$ of the bases $a$ are Witnesses for the compositeness of $n$ using the Miller Rabin primality test. <i> This means that a probabilistic argument can be made.</i>\n",
    "\n",
    "<i>Partial Proof</i>: If $n=p$ is a prime then since $a^{n-1} \\equiv 1 \\pmod {n}$ then eitber all integers in the sequence are equivalent to $1$ or in this list of succesive squares one integer $b \\not \\equiv 1 \\pmod{p}$ and $b^2 \\equiv 1 \\pmod{p}$. This implies that $p | b^2 -1$ hence $p | (b-1)$ or $p | b + 1$. The first case would mean $b \\equiv 1 \\pmod{p}$, and hence we have $b \\equiv -1 \\pmod{p}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def powersOfTwo(n):\n",
    "    \"\"\"\n",
    "    Returns (k,q) such that n = 2^k q\n",
    "    \"\"\"\n",
    "    q = n\n",
    "    k = 0\n",
    "    while q%2==0:\n",
    "        k = k + 1\n",
    "        q = q//2\n",
    "    return (k,q)\n",
    "\n",
    "def MillerRabinWitness(a,n):\n",
    "    \"\"\"\n",
    "    Tests whether a is a witness for the compositeness of n\n",
    "    \"\"\"\n",
    "    g = gcd(n,a)\n",
    "    if g > 1 and g < n:\n",
    "        return True\n",
    "    (k,q) = powersOfTwo(n-1)\n",
    "    b = fastexp(a,q,n)\n",
    "    if b==1:\n",
    "        return False\n",
    "    for i in range(k):\n",
    "        if b==n-1:\n",
    "            return False\n",
    "        b = (b*b)%n\n",
    "    return True\n",
    "\n",
    "from math import log as log\n",
    "from random import randint as randint\n",
    "\n",
    "def MillerRabinPrimalityTest(n, precision = None):\n",
    "    if precision == None:\n",
    "        precision = min(1/n,1/1000)\n",
    "    if n==2:\n",
    "        return True\n",
    "    if n%2==0:\n",
    "        return False\n",
    "    A = int(log(precision,3/4))+1\n",
    "    for trial in range(A):\n",
    "        a = randint(2,n-1)\n",
    "        if MillerRabinWitness(a,n):\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "MillerRabinPrimalityTest(105,0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1000000000000000000000000000057, 1000000000000000000000000000099]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def FindPrimes(N,h):\n",
    "    #Returns the list of (Miller-Rabin)--primes in the interval [N, N + h)\n",
    "    X = range(N, N + h)\n",
    "    X = list(filter(MillerRabinPrimalityTest, X))\n",
    "    return X\n",
    "\n",
    "FindPrimes(1000000000000000000000000000000,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deterministic Miller-Rabin Test,\n",
    "\n",
    "It can be shown that if one assumes the Riemann hypothesis, then if $n$ is not a prime, then an integer $a$ will be a witness for the compositeness of $a$ such that $$a \\leq 2\\ln(n).$$ This gives us a deterministic test for primality (dependent on Generalized Riemann Hypothesis). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pollard $p-1$ factorization algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(101, 103)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Pollardp(n, B):\n",
    "    a = 2\n",
    "    k = max(1,B//1000) #this is a optimization technique, we don't check the gcd everytime\n",
    "    for j in range(2,B):\n",
    "        a = fastexp(a,j,n)\n",
    "        if j%k ==0:\n",
    "            g = gcd(a-1,n)\n",
    "            if g >1 and g<n:\n",
    "                return (g,n//g)\n",
    "    return None\n",
    "\n",
    "Pollardp(101*103,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p =  223456789121 q =  223456790033\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "49932936808059655630993"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = 10**7\n",
    "\n",
    "def BSmooth(n,B):\n",
    "    #Checks whether n is B smooth or not\n",
    "    for prime in range(2,B):\n",
    "        while n%prime==0:\n",
    "            n=n//prime\n",
    "    if n ==1:\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "\n",
    "def BSmoothFinder(alist, B):\n",
    "    for n in alist:\n",
    "        if BSmooth(n-1,B):\n",
    "            return n\n",
    "    return None\n",
    "    \n",
    "primes = FindPrimes(223456789101,1000)\n",
    "q = primes.pop()\n",
    "p = BSmoothFinder(primes,B)\n",
    "print('p = ',p, 'q = ', q)\n",
    "n = p*q\n",
    "n "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(223456789121, 223456790033)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Pollardp(n,B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 2, 2, 2, 2, 2, 2, 5, 31, 71, 158633]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(p-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 2, 2, 2, 7, 643, 3102877]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(q-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Factorization Via Difference of Squares\n",
    "\n",
    "If we find integers such that $a^2 \\equiv b^2 \\pmod n$ with $n = pq$ then we can use that since $$pq | (a-b) (a + b) $$\n",
    "we may have that $p | (a + b)$ and $b | (a - b) $.\n",
    "\n",
    "How to find such numbers. \n",
    "\n",
    "Let $p$ and $q$ be as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 914387"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1869 ^2 = [2, 2, 2, 2, 3, 5, 5, 5, 5, 5, 5] mod n\n",
      "1909 ^2 = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 5, 11] mod n\n",
      "3387 ^2 = [3, 5, 5, 5, 11, 11, 11] mod n\n",
      "4586 ^2 = [3, 3, 5, 11] mod n\n",
      "5023 ^2 = [2, 2, 2, 2, 2, 2, 2, 5, 7, 11, 11] mod n\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt as sqrt\n",
    "a = int(sqrt(n)) + 1\n",
    "sample =0\n",
    "while sample <5:\n",
    "    x = (a*a)%n\n",
    "    if BSmooth(x,13):\n",
    "        print(a,\"^2 =\", BruteForceFactor(x), \"mod n\")\n",
    "        sample = sample +1\n",
    "    a = a + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9835 164255\n"
     ]
    }
   ],
   "source": [
    "a = (1869*1909*3387)%n\n",
    "b = (2**9*3*5**5*11**2)%n\n",
    "print(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1103"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(a - b,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "829"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(a + b,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1103 * 829 ==n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 9873417"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3144 ^2 = [3, 7, 7, 7, 11] mod n\n",
      "3177 ^2 = [2, 2, 2, 3, 7, 7, 11, 17] mod n\n",
      "3487 ^2 = [2, 2, 2, 7, 7, 7, 7, 7, 17] mod n\n",
      "3555 ^2 = [2, 2, 2, 2, 2, 2, 3, 7, 11, 11, 17] mod n\n",
      "4578 ^2 = [2, 3, 5, 5, 5, 5, 17, 19] mod n\n",
      "4806 ^2 = [2, 3, 7, 13, 17, 19, 19] mod n\n",
      "4833 ^2 = [3, 5, 7, 7, 17, 17, 17] mod n\n",
      "5499 ^2 = [2, 3, 3, 5, 5, 5, 5, 5, 11] mod n\n",
      "5655 ^2 = [2, 3, 3, 3, 11, 11, 19, 19] mod n\n",
      "6288 ^2 = [2, 2, 3, 7, 7, 7, 11] mod n\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt as sqrt\n",
    "a = int(sqrt(n)) + 1\n",
    "sample =0\n",
    "while sample <10:\n",
    "    x = (a*a)%n\n",
    "    if BSmooth(x,20):\n",
    "        print(a,\"^2 =\", BruteForceFactor(x), \"mod n\")\n",
    "        sample = sample +1\n",
    "    a = a + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6315897 3814236\n"
     ]
    }
   ],
   "source": [
    "a = (3144*3177*3487)%n\n",
    "b = (2**3*3**7**5*11*17)%n\n",
    "print(a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(a + b, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3, 23, 143093]"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BruteForceFactor(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "407 ^2 = [3, 3, 3, 3, 5, 7] mod n\n",
      "412 ^2 = [2, 3, 3, 5, 7, 11] mod n\n",
      "435 ^2 = [7, 7, 7, 7, 11] mod n\n",
      "484 ^2 = [2, 3, 3, 3, 3, 3, 3, 7, 7] mod n\n",
      "638 ^2 = [3, 3] mod n\n",
      "755 ^2 = [2, 2, 2, 2, 11] mod n\n",
      "757 ^2 = [2, 2, 2, 2, 2, 2, 2, 5, 5] mod n\n",
      "763 ^2 = [2, 2, 2, 2, 2, 5, 7, 11] mod n\n",
      "777 ^2 = [2, 2, 2, 5, 7, 11, 11] mod n\n",
      "814 ^2 = [2, 2, 3, 3, 3, 3, 5, 7] mod n\n"
     ]
    }
   ],
   "source": [
    "n = 81407\n",
    "a = int(sqrt(n)) + 1\n",
    "sample =0\n",
    "while sample <10:\n",
    "    x = (a*a)%n\n",
    "    if BSmooth(x,12):\n",
    "        print(a,\"^2 =\", BruteForceFactor(x), \"mod n\")\n",
    "        sample = sample +1\n",
    "    a = a + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8788"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(503**2)%n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70085"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(313*484)%n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "34398"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2*3**3*7**2*13)%n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 78005 \n",
    "b = 34398"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "81407"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "112403"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(a - b,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = (2*3**3*5*7**4*11*13)%n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = (313*407*412*435)%n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "61444 14835\n"
     ]
    }
   ],
   "source": [
    "print(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "641"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gcd(a + b,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}