Shor's Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Course Project - Shor's Algorithm\n",
    "\n",
    "CS251 Quantum Information Science, 2018 @ ACM Honors Class, SJTU\n",
    "<script src=\"https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML\" type=\"text/javascript\"></script>\n",
    "\n",
    "## Lecturer\n",
    "\n",
    "Prof. **Runyao Duan**\n",
    "\n",
    "\n",
    "## Group Member\n",
    "\n",
    "- **Zetian Jiang** (516030910577)\n",
    "- **Zhanghao Wu** (516030910593)\n",
    "- **Zhou Fan** (516030910574)\n",
    "\n",
    "\n",
    "## Abstract\n",
    "\n",
    "We implemented Shor's Algorithm in this project. The quantum part of the algorithm, i.e. the phase estimation subroutine, is implemented in Q# programming language, and we encapsulated this part to provide a flexible callable interface. In this way one can easily call our quantum subroutine in programming languages other than Q# or C#, for instance, Python. With this cross-language feature, we were able to implement the classical part of Shor's algorithm in Python, which is presented later in this Jupyter notebook. In addition, we introduced parallel computation in the implementation of Shors' algorithm, running multiple order finding subroutine simultaneously to accelerate the process. Finally, we wrote a test script to factorize random numbers from 1 to 100, and the implementation passed the test.\n",
    "\n",
    "\n",
    "## Adavanced Feature\n",
    "\n",
    "Our implementation of Shor's algorithm has the following two advanced features described in the course task documentation:\n",
    "\n",
    "- **Cross-language interoperability**: one can call our Q# subroutine in languages other than C#, such as Python\n",
    "- **Parallel computation**: we introduced parallel computation in the implementation of Shor's algorithm\n",
    "\n",
    "## Requirement\n",
    "\n",
    "- conda\n",
    "- .Net core sdk >= 2.0\n",
    "\n",
    "## How to run\n",
    "\n",
    "1. To install the dependencies, run command `conda env create -f ./environment.yml`\n",
    "2. Activate the python environment, run command `activate qsharp-samples`\n",
    "3. Run command `cd ./myShorLib/; dotnet build`, which will build the dynamic library of the Q#\n",
    "4. Run command `cd ../Classical` to go into the classical part of the codes.\n",
    "5. You can now run `python shor.py` to run the Shor's Algorithm and the program will wait for your input of the integer to be factorized. Or you can run `python test.py` to run the Shor's Algorithm to factorize integers from 1 to 100, and output the result to ./test.out\n",
    "\n",
    "\n",
    "## Directory Structure \n",
    "\n",
    "### Classical Part\n",
    "\n",
    "[gcd.py](./gcd.py): This code calculate the greatest common divider of two given number a and b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd(a, b):\n",
    "    if b > a:\n",
    "        a, b = b, a\n",
    "    while b > 0:\n",
    "        a = a % b\n",
    "        a, b = b, a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[fastPow.py](./Classical/fastPow.py): This code calculate the answer of $a^b\\pmod N$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fastPow(a, b, N):\n",
    "    ans = 1\n",
    "    while (b > 0):\n",
    "        if b % 2:\n",
    "            ans = ans * a % N\n",
    "        b = b // 2\n",
    "        a = a * a\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[miller_robin.py](./Classical/miller_robin): This code test that whether a given number $N$ is a prime number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def miller_robin(n):\n",
    "    if n == 2 or n == 3:\n",
    "        return True\n",
    "    if n % 2 == 0:\n",
    "        return False\n",
    "\n",
    "    t = 10\n",
    "    q = 0\n",
    "    m = n - 1\n",
    "    while m % 2 == 0:\n",
    "        q += 1\n",
    "        m /= 2\n",
    "    for _ in range(t):\n",
    "        a = random.randint(2, n-2)\n",
    "        x = fastPow.fastPow(a, m, n)\n",
    "        if x == 1:\n",
    "            continue\n",
    "        j = 0\n",
    "        while j < q and x != n-1:\n",
    "            x = (x * x) % n\n",
    "            j += 1\n",
    "        if j >= q:\n",
    "            return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[power.py](./Classical/power.py): This code test whether a given number $N$ is in the form of $a^b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power(N):\n",
    "    def isPower(l, r, s, N):\n",
    "        if (l > r):\n",
    "            return -1\n",
    "        mid = (l + r) / 2\n",
    "        ans = fastPow.fastPowBool(mid, s, N)\n",
    "        if (ans == N):\n",
    "            return mid\n",
    "        elif (ans < N):\n",
    "            return isPower(mid+1, r, s, N)\n",
    "        else:\n",
    "            return isPower(l, mid-1, s, N)\n",
    "\n",
    "    s = int(math.floor(math.log(N, 2))) + 1\n",
    "    r = int(math.floor(math.sqrt(N))) + 1\n",
    "    for i in range(2, s):\n",
    "        ans = isPower(2, r, i, N)\n",
    "        if ans != -1:\n",
    "            return ans\n",
    "    return -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[shor.py](./Classical/shor.py): This code is the exactly we combine the tools together and factorize a given number $N$.\n",
    "* *shor* is the interface that can be called to find factorization of $N$\n",
    "* *order_finding* is the function that can be called to find possible order for the given number pair $(x, n)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class myShor():\n",
    "    def __init__(self, precision, thread_num):\n",
    "        self.Precision = precision\n",
    "        self.Thread_Num = thread_num\n",
    "\n",
    "    def order_finding(self, x, n):\n",
    "        qsim = qsharp.QuantumSimulator()\n",
    "\n",
    "        def phase_estimation(x, n):\n",
    "            tmp = int(qsim.run(PhaseEstimation, x, n,\n",
    "                               self.Precision).result().clr_object)\n",
    "            theta = float(tmp) / (2**self.Precision)\n",
    "            return theta\n",
    "\n",
    "        for _ in range(2):\n",
    "            theta = phase_estimation(x, n)\n",
    "            if theta == 0:\n",
    "                print(\"========\\nOrder Finding for: x={}, N={}\\nGet theta: {}\\nNo r estimated\\n\".format(\n",
    "                    x, n, theta))\n",
    "                continue\n",
    "\n",
    "            r = fraction.denominator(theta, n)\n",
    "            print(\"========\\nOrder Finding for: x={}, N={}\\nGet theta: {}\\nEstimate r: {}\\n\".format(\n",
    "                x, n, theta, r))\n",
    "            for i in r:\n",
    "                m = fastPow.fastPow(x, i, n)\n",
    "                if m == 1:\n",
    "                    return i\n",
    "        return -1\n",
    "\n",
    "    def shor(self, n):\n",
    "        if miller_robin.miller_robin(n):\n",
    "            return (1, n)\n",
    "        else:\n",
    "            tmp = power.power(n)\n",
    "            if tmp != -1:\n",
    "                return (tmp, n // tmp)\n",
    "            else:\n",
    "                if (n % 2 == 0):\n",
    "                    return (2, n // 2)\n",
    "                while True:\n",
    "                    # Parrel computing for some random x\n",
    "                    xlist = random.sample(range(3, n - 1), self.Thread_Num)\n",
    "                    g = [gcd.gcd(x, n) for x in xlist]\n",
    "                    for idx, g in enumerate(g):\n",
    "                        if (g != 1):\n",
    "                            # ======= For debug ===========\n",
    "                            # while gcd.gcd(xlist[idx], n) != 1:\n",
    "                            #     newx = random.randint(3, n - 1)\n",
    "                            #     xlist[idx] = newx\n",
    "                            # ======= In Real Quantum Computer =========\n",
    "                            return (g, n // g)\n",
    "\n",
    "                    print(\"======== Order Finding Started ========\")\n",
    "                    threadPool = ThreadPool(processes=self.Thread_Num)\n",
    "                    results = []\n",
    "                    for x in xlist:\n",
    "                        results.append(threadPool.apply_async(\n",
    "                            self.order_finding, args=(x, n)))\n",
    "                    threadPool.close()\n",
    "                    threadPool.join()\n",
    "                    results = [r.get() for r in results]\n",
    "\n",
    "                    for r in results:\n",
    "                        if r == -1:\n",
    "                            continue\n",
    "                        if (r % 2 == 0):\n",
    "                            s = fastPow.fastPow(x, r // 2, n)\n",
    "                            if (s != 1 and s != n-1):\n",
    "                                g1 = gcd.gcd(s+1, n)\n",
    "                                g2 = gcd.gcd(s-1, n)\n",
    "                                if (g1 != 1):\n",
    "                                    return (g1, n // g1)\n",
    "                                elif (g2 != 1):\n",
    "                                    return (g2, n // g2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[fraction.py](./Classical/fraction.py): This code estimate the denominator of possible fraction representing the given decimal $x$. The fra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def denominator(x, N):\n",
    "    seq = []\n",
    "    if (x < 1):\n",
    "        x = 1 / x\n",
    "    k = 10\n",
    "    for i in range(10):\n",
    "        seq.append(math.floor(x))\n",
    "        if (abs(x - math.floor(x)) < 1e-9):\n",
    "            k = i + 1\n",
    "            break\n",
    "        x = 1 / (x - math.floor(x))\n",
    "    ans = []\n",
    "    ans2 = []\n",
    "    for i in range(k):\n",
    "        up = 1\n",
    "        down = seq[i]\n",
    "        for j in range(i):\n",
    "            t = seq[i-1-j]\n",
    "            d = down\n",
    "            down = down * t + up\n",
    "            up = d\n",
    "        ans.append(int(down))\n",
    "        ans2.append(int(up))\n",
    "    return ans\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[test.py](./Classical/test.py): This code test our implementaion of shor's algorithm for numbers from 1 to 100. And our implementation handles the task very well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alg = myShor(6, 3)\n",
    "with open('test.out', 'w', encoding='utf8') as outfile:\n",
    "    for i in range(1, 100):\n",
    "        print(alg.shor(i), file=outfile)\n",
    "        outfile.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantum Part\n",
    "[Operation.qs](./myShorLib/Operation.qs): This implement the quantum part of the algorithm\n",
    "* *PhaseEstimation* is the interface for quantum phase estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "operation PhaseEstimation(\n",
    "    x : Int, \n",
    "    N : Int, \n",
    "    precision : Int) : Int\n",
    "{\n",
    "    body\n",
    "    {\n",
    "        let eigenLength = BitSize(N);\n",
    "        mutable res = 0;\n",
    "\n",
    "        using (eigenvector = Qubit[eigenLength]){\n",
    "            X(eigenvector[0]);\n",
    "\n",
    "            using (target = Qubit[precision]){\n",
    "\n",
    "                PhaseEstimationImpl(x, N, target, eigenvector);\n",
    "\n",
    "                for (i in 0..(precision - 1)) {\n",
    "                    set res = res * 2;\n",
    "                    if (M(target[i]) == One) {\n",
    "                        set res = res + 1;\n",
    "                    }\n",
    "                }\n",
    "                ResetAll(target);\n",
    "                ResetAll(eigenvector);\n",
    "            }\n",
    "        }\n",
    "\n",
    "        return res;\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* *PhaseEstimationImpl* implements the details for phase estimation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "operation PhaseEstimationImpl (\n",
    "    x : Int,\n",
    "    N : Int,\n",
    "    target : Qubit[],\n",
    "    eigenvector : Qubit[]) : ()\n",
    "{\n",
    "    body {\n",
    "        let targetLength = Length(target);\n",
    "\n",
    "        for (idx in 0..(targetLength - 1)) {\n",
    "            H(target[idx]);\n",
    "        }\n",
    "\n",
    "        mutable power = 1;\n",
    "        for (idx in 0..(targetLength - 1)) {\n",
    "            (Controlled ConstructU) ([target[targetLength - 1 -idx]], (x, N, power, eigenvector)); \n",
    "            set power = power * 2;\n",
    "        }\n",
    "\n",
    "        (InverseFT)(target);\n",
    "        // (Adjoint QFT)(BigEndian(target));\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* $U$ is constructed by the code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "operation ConstructU(\n",
    "    x : Int,\n",
    "    modulus : Int, \n",
    "    power : Int, \n",
    "    target : Qubit[]) : ()\n",
    "{\n",
    "    body {\n",
    "        ModularMultiplyByConstantLE(\n",
    "            PowerMod(x, power, modulus), \n",
    "            modulus,\n",
    "            LittleEndian(target)\n",
    "        ); \n",
    "    }\n",
    "    adjoint auto\n",
    "\n",
    "    controlled auto\n",
    "    adjoint controlled auto \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* *Inverse fourier transformation* is implemented as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "operation InverseFT(qs : Qubit[]) : () \n",
    "{\n",
    "    body {\n",
    "        let qLength = Length(qs);\n",
    "        for (i in 0..(qLength - 1)) {\n",
    "            for (j in 0..(i-1)) {\n",
    "                (Controlled R1Frac) ([qs[j]], (1, i - j, qs[i]));\n",
    "            }\n",
    "            H(qs[i]);\n",
    "        }\n",
    "\n",
    "        for (i in 0..qLength/2-1) {\n",
    "            SWAP(qs[i],qs[qLength-1-i]);\n",
    "        }\n",
    "\n",
    "    }\n",
    "    adjoint auto\n",
    "    controlled auto\n",
    "    controlled adjoint auto\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Process of Shor's Algorithm\n",
    "\n",
    "The input is a number $N$ to be factorized.\n",
    "1. Check if $N$ is a prime number, this can be done with some classical algorithm such as [Miller-Robin](./Classical/miller-robin.py) primality test. If $N$ is prime number, the factorization is no longer needed.\n",
    "2. Check if $N$ is a power of a prime number, that is, if there exists a prime number $a$ such that $N = a^k$. If so, then $a$ is a factor of $N$, then we have the factorization of $N$. And the checking process can be found in [power](./Classical/power.py)\n",
    "3. Now the number $N$ can be represented as $$N = p_1^{k_1} p_2^{k_2} ... p_n^{k_n}.$$ Perform the following operations until we find a factorization of $N$:\n",
    "   1. Generate a random number $x < N$.\n",
    "   2. Check if $x$ is coprime to $N$, using the [gcd](./Classical/gcd.py), if not, then the greatest common divisor of $x$ and $N$ is a factor of $N$. Then we have a factorization of $N$.\n",
    "   3. Now we know that $x$ and $N$ are coprime to each other. There exists an order $r$ of $x$ such that $$x^r \\equiv 1 \\pmod N.$$ To get the order $r$, the [quantum phase estimation](./myShorLib/Operation.qs) comes in, the details of which will be discussed in next part of this report.\n",
    "   4. The quantum phase estimation returns a number $\\theta = \\frac{s}{r}$ in decimal with certain precision. Then a [continued fraction decomposition](./Classical/fraction.py) is needed to estimate the exact value of $r$. \n",
    "   5. The steps to estimate $r$ from $\\theta$ are as follows. Present $\\theta$ as $$\\frac{1}{a_0 + \\frac{1}{a_1 + \\dots \\frac{1}{a_n + b}}}.$$ Regardless of b, we use $$\\frac{1}{a_0 + \\frac{1}{a_1 + \\dots \\frac{1}{a_n}}} = \\frac{c_n}{d_n}$$ as its estimation. Check each $d_n \\leq N$ whether $d_n$ is the order of x.\n",
    "   6. Unfortunately, because of the limited precision of the phase returned by the quantum algorithm, the randomness of the program and the limitation of the continued fraction decomposition for reducible fractions, we may not get a proper $r$. \n",
    "   7. If the order $r$ is odd, try to run the phase estimation again or try another $x$. Else we have $$(x^{\\frac{r}{2}} + 1)(x^{\\frac{r}{2}} - 1) = x^r - 1 \\equiv 0 \\pmod N.$$ If $x^{\\frac{r}{2}} + 1 \\pmod N$ and $x^{\\frac{r}{2}} - 1 \\pmod N$ are both not zero, then we have a factorization of $N$. If not, try to run the phase estimation again or try another $x$. And here is the exact place where parallel computing plays a role.\n",
    "\n",
    "## Quantum Phase Estimation\n",
    "\n",
    "1. Input two register. The first one contains $t$ qubits in state $\\left|{0}\\right\\rangle$,where $t$ is the precision of the estimation. The second one is eigenvector of U (actually it can be any vector that is in the space spanned by the eignevectors of U, and in the code it is an initialized vector with $\\left|{1}\\right\\rangle$ at the first and $\\left|{0}\\right\\rangle$ at other position). Apply the circuit shown in Figure 1 to get each bit of $\\phi$.\n",
    "2.  The circuit begins by applying a Hadamard transform to the first register, followed by application of controlled-U operations on the second register, with U raised to successive powers of two. \n",
    "![](./others/pe2.png)\n",
    "*Figure 1: The first stage of the phase estimation procedure*\n",
    "3. Apply inverse Fourier transform circuit on first register while doing nothing on second one. The Fourier transform circuit is shown in Figure 2 and the normalization factors and swap gates which should be applied at the circuit are not shown. \n",
    "The circuit begins by applying a Hadamard transform , followed by application of controlled-$R_i$ operations where operation $R_k$'s matrix representation is $\\begin{pmatrix}1 & 0 \\\\0 & e^{2 \\pi i / 2^k}\\end{pmatrix}$. Finally it ends with a swap gate to swap $\\left|{\\phi_i}\\right\\rangle$ and $\\left|{\\phi_{n-i+1}}\\right\\rangle$.\n",
    "    ![](./others/qft.png)\n",
    "*Figure 2: Quantum Fourier Transform circuit*\n",
    "4. The final state of first register is $\\left|{\\theta}\\right\\rangle$ where $\\theta$ is a t bit estimation of a element in $\\{\\frac{s}{r} | $r is the order$, s=0 \\dots r-1\\}$. Measure the first register and get a decimal $\\theta$.\n",
    "\n",
    "## Cross-language Interoperability\n",
    "\n",
    "To encapsulate the quantum subroutine implemented in Q#, we separate the Q# operations as a standalone .NET standard library project. And we compiled this project to generate the library assemblies. Then one can use operations in that library in languages other than C# by adding a reference to the assemblies file. For instance, to use our Q# library in Python, one could use the Python.NET package to load the assemblies, and one should also add references to the assemblies from the Q# development kit released by Microsoft as dependencies.\n",
    "\n",
    "The following example shows how to load our assemblies in Python. One should first have Python.NET package installed in Python environment in order to use the clr module. One can use the command `pip install pythonnet` to install Python.NET package. It is also needed to add the path to the assemblies file in the system path such that the `clr` module could find the assemblies. Another requirement is Python library `qsharp`, which provides a QuantumSimulator class which wraps the Microsoft.Quantum.Simulation.Simulators.QuantumSimulator .NET class provided with the Quantum Development Kit. This wrapper provides a few useful convienence features that we will use later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('./myShorLib/bin/Release/netstandard2.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qsharp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import clr\n",
    "clr.AddReference(\"Microsoft.Quantum.Canon\")\n",
    "clr.AddReference('myShorLib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here `myShorLib` is the quantum library we implemented in Q#. Having added reference to the assemblies, we can now use operations in the library. To call the operation, we could just import it just as a normal Python package."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parallel computation\n",
    "\n",
    "In Shor's algorithm, we use the quantum phase estimation subroutine to find the order $r$ of a random number $x < N$, i.e. a positive number $r$ such that\n",
    "$$x^r \\equiv 1 \\pmod N$$\n",
    "We introduced parallel computation in this process, that is, performing the order finding subroutine for multiple random numbers simultaneously with multi-threading technique. This could accelerate the process of Shor's algorithm when $N$ is large. With the help of the cross-language interoperability of our implementation of the phase estimation subroutine, we can easily implement the parallel part in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threadPool = ThreadPool(processes=self.Thread_Num)\n",
    "results = []\n",
    "for x in xlist:\n",
    "    results.append(threadPool.apply_async(\n",
    "        self.order_finding, args=(x, n)))\n",
    "threadPool.close()\n",
    "threadPool.join()\n",
    "results = [r.get() for r in results]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, when we calls the code to factorize 21. The code will calculate several $x$ at the same time, and get the result in a short time. The process is shown in the figure below.\n",
    "![](./others/parallel.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:qsharp-samples]",
   "language": "python",
   "name": "conda-env-qsharp-samples-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}