Algorithm-ShorAlgorithm-Maths.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Mathematics Behind Shor's Algorithm\n",
    "\n",
    "Shor's algorithm is an algorithm that can factorize big numbers quickly. For an integer $N$, Shor's algorithm can find a factor of it in $O((\\log_{} N)^{3})$ which is an almost exponential speedup over classical algorithms. \n",
    "\n",
    "This algorithm brought a lot of interest into building quantum computers, because if an ideal quantum computer could be built, then the widely used public-key cryptography RSA scheme could be broken.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classical part\n",
    "\n",
    "Shor's algorithm is hybrid becasuse it mixes a classical part and a quantum part. In this notebook, we will discuss the mathematics behind the classical part and assume the quantum part as a black box. \n",
    "\n",
    "Firstly, we need to define what the algorithm does. It solves the factorization problem: given an intenger $N$, find two integers greater than one $P$ and $Q$ such that $PQ = N$, or state that $N$ is prime.\n",
    "\n",
    "Secondly, we need to define what the input of the algorithm is: it is $N$, an odd integer that is neither a prime nor the power of a prime. These assumptions are necessary for the algorithm to work, and when they are not respected factorizing is still easy.\n",
    "\n",
    "If $N$ is even, we can pick $P = 2$ and $Q = N/2$ so the even case is straightforward. \n",
    "\n",
    "If $N$ is prime, we can check with primality tests such as [Miller-Rabin's](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test) for primality in a reasonable time complexity so this case is also easier than solving general factorization.\n",
    "\n",
    "If $N$ is the power of a prime, we can check for every $2 \\leq k \\leq \\log_{3}N$ if $N^{1/k}$ is an integer. If it is, then $P = N^{1/k}$ and $Q = N^{(k-1)/k}$ is a solution and the case for powers of primes is solved.\n",
    "\n",
    "Below are an implementation of a naive primality checking and the check for an exact power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import gcd, log\n",
    "from random import randint\n",
    "\n",
    "def is_prime(N):\n",
    "    \"\"\"Returns if N is prime or not. Notice that this is not optimal,\n",
    "    there is faster primality testing algorithms e.g. Miller-Rabin\n",
    "    \"\"\"\n",
    "    \n",
    "    if N == 2:\n",
    "        return True  # only even prime\n",
    "    if N % 2 == 0 or N <= 1:\n",
    "        return False  # even numbers and 1 are not prime\n",
    "    \n",
    "    for i in range(3, N, 2):  # only try odd candidates\n",
    "        if i*i > N:\n",
    "            break  # we only need to check up to sqrt(N)\n",
    "        if N % i == 0:\n",
    "            return False  # found a factor\n",
    "    \n",
    "    return True\n",
    "\n",
    "def find_power_k(N):\n",
    "    \"\"\"Returns the smallest k > 1 such that N**(1/k) is an integer,\n",
    "    or 1 if there is no such k.\n",
    "    \"\"\"\n",
    "    \n",
    "    upper_bound = int(log(N)/log(3))\n",
    "    \n",
    "    for k in range(2, upper_bound + 1):\n",
    "        p = int(N**(1/k))\n",
    "        if p**k == N:\n",
    "            return k\n",
    "    \n",
    "    return 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum part\n",
    "\n",
    "Shor's algorithm speedup comes from the quantum subroutine for period finding. To find the smallest $r$ for $a^{r} \\equiv 1 \\pmod N$, the algorithm uses Phase Estimation algorithm and the Quantum Fourier Transform to make calculations faster than possible with a classical computer. We will assume that subroutine as a blackbox for now.\n",
    "\n",
    "Despite that, we can implement a slower classical version that gives the same result as that blackbox. We test our code with $a = 2$ and $N=15$, which has a solution $r = 4$ because $2^{4} = 16$ and $16 \\equiv 1 \\pmod {15} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def order_finding(a, N):\n",
    "    \"\"\"Returns the smallest r such that a^r = 1 mod N\n",
    "    Notice that this is a naive classic implementation and is\n",
    "    exponentially slower than the quantum version invented by Peter Shor.\n",
    "    \"\"\"\n",
    "    \n",
    "    i = 1\n",
    "    a_i = a % N\n",
    "    \n",
    "    while a_i != 1:\n",
    "        i += 1\n",
    "        a_i = (a_i * a) % N\n",
    "    \n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "order_finding(2, 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shor's algorithm\n",
    "\n",
    "Why focus on the smallest $r$ if it is easy to find an arbitrary one that works (e.g. $\\phi(N)$)? The fact is that the smallest $r$ has a nice property that is key in our analysis. First, let's rewrite our equation:\n",
    "\n",
    "$$\n",
    "a^{r} - 1 \\equiv 0 \\pmod N \\iff\n",
    "$$\n",
    "\n",
    "$$\n",
    "(a^{r/2} + 1)(a^{r/2} - 1) \\equiv 0 \\pmod N \\iff\n",
    "$$\n",
    "\n",
    "$$\n",
    "N \\mid (a^{r/2} + 1)(a^{r/2} - 1)\n",
    "$$\n",
    "\n",
    "Here, we assumed that $r$ is even, and that will be a requirement for our algorithm to work. Notice that $N \\nmid (a^{r/2} - 1)$, because that would imply $a^{r/2} \\equiv 1$ and violate our condition that $r$ was the smallest. Thus, the prime factors of $N$ are distributed among $a^{r/2} - 1$ and $a^{r/2} + 1$.\n",
    "\n",
    "There are two cases. In the first one, we're out of luck: if $N \\mid (a^{r/2} + 1)$, then it might be that all the factors are on $a^{r/2} + 1$ and we can conclude nothing about the prime factorization of $N$.\n",
    "\n",
    "In the second case, $N \\nmid (a^{r/2} + 1)$ and we can make a conclusion: because $N$ divides the product but not the numbers individually, at least one of $\\gcd(N, a^{r/2} + 1)$ or $\\gcd(N, a^{r/2} - 1)$ will be a non trivial factor of $N$. Our factorization is done.\n",
    "\n",
    "Thus, if finding $r$ is done quickly, we can try multiple values for $a$ until we found one that yields a factor. It can be shown that the probability that $a$ works is at least $1/2$, so with few attempts for $a$ we will find a factor. Shor's algorithm is concluded\n",
    "\n",
    "\n",
    "Given the math concepts behind Shor's algorithm, we can write the pseudo-code for the algorithm. Notice that in this step we are not worrying about the implementation of other parts of the algorithm: we assume that the classical parts of $\\gcd$, primality-testing, checking if the k-th root is an integer are implemented and available for use. We also assume that the quantum subroutine for period finding is available for use. This yields the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shor_algorithm(N):\n",
    "    \"\"\"Returns a pair of integers (P, Q) such that PQ = N for integer N\"\"\"\n",
    "    \n",
    "    if N % 2 == 0:  # even case\n",
    "        return (N//2, 2)\n",
    "    \n",
    "    if is_prime(N):  # prime case\n",
    "        return (N, 1)  # N is primes, factors cannot be found\n",
    "    \n",
    "    if find_power_k(N) > 1:  # prime power case\n",
    "        P = int(N**(1/find_power_k(N)))  # we find a k such that N**(1/k) is an integer\n",
    "        Q = N//P\n",
    "        return (P, Q)\n",
    "    \n",
    "    # Now we can assume that the criteria for Shor's algorithm is met\n",
    "    \n",
    "    while True:\n",
    "        # Non-deterministic, we will try multiple values for a\n",
    "        a = randint(2, N-1)  # pick random a\n",
    "        \n",
    "        if gcd(a, N) != 1:  # Lucky case: a and N are not coprime!\n",
    "            P = gcd(a, N)  # gcd yields a non-trivial factor\n",
    "            Q = N//P\n",
    "            return (P, Q)\n",
    "        \n",
    "        r = order_finding(a, N)  # quantum subroutine of the code\n",
    "        \n",
    "        if r % 2 == 0:\n",
    "            continue\n",
    "\n",
    "        P = gcd(a**(r//2) - 1, N)\n",
    "        if P != 1:\n",
    "            Q = N//P  # gcd yielded non trivial factor\n",
    "            return (P, Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing\n",
    "\n",
    "We can check if our code works! For example, let's take $N = {10013}$. That number has ${3}$ prime factors: ${17}$, ${19}$ and ${31}$. If we apply Shor's algorithm until we reach primality, we can find all of them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shor's algorithm found 10013 = 19 x 527 which is correct!\n"
     ]
    }
   ],
   "source": [
    "N = 10013\n",
    "P, Q = shor_algorithm(N)\n",
    "print(\n",
    "    \"Shor's algorithm found {} = {} x {} which is {}!\".format(N, P, Q, [\"incorrect\", \"correct\"][P*Q==N])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shor's algorithm found 527 = 17 x 31 which is correct again!\n"
     ]
    }
   ],
   "source": [
    "S, T = shor_algorithm(Q)\n",
    "print(\n",
    "    \"Shor's algorithm found {} = {} x {} which is {} again!\".format(Q, S, T, [\"incorrect\", \"correct\"][S*T==Q])\n",
    ")"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}