{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "628c747b-084c-4fb6-843e-7c7f4cbf07a1",
   "metadata": {},
   "source": [
    "# Factoring Integers With Shor's Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0627c4-95bd-47b9-9ab5-e8374d25508e",
   "metadata": {},
   "source": [
    "The most famous application of quantum computers is factoring integers using Shor's algorithm. This algorithm is particularly significant because an efficient factorization algorithm could potentially break modern asymmetric encryption schemes, such as RSA.\n",
    "\n",
    "For small integers, this quantum algorithm can be simulated on classical computers. The main challenge in classical implementation lies in the order-finding algorithm. We will first introduce the classical implementation of this algorithm as a preliminary step, and then proceed to explain the quantum order-finding algorithm.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "148568d8",
   "metadata": {},
   "source": [
    "First let's install some libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "486fb3d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install contfrac==1.0.0 -q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "88619b35-624e-4853-8d2e-6d9952390e5e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from math import gcd, log2, ceil\n",
    "import numpy as np\n",
    "import random\n",
    "import cudaq\n",
    "from cudaq import *\n",
    "import fractions\n",
    "import matplotlib.pyplot as plt\n",
    "import contfrac\n",
    "\n",
    "if cudaq.num_available_gpus() > 0 and cudaq.has_target(\"nvidia\"):\n",
    "    cudaq.set_target(\"nvidia\")\n",
    "else:\n",
    "    print(\"CUDA or GPU support is unavailable. Running with CPU simulator. Performance may be significantly reduced.\")\n",
    "    cudaq.set_target(\"qpp-cpu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8523c6f-d729-4d5e-8df9-1780c343a58e",
   "metadata": {},
   "source": [
    "## Shor's algorithm\n",
    "The factoring problem for an integer $N$ aims to find integers $a_1$, $a_2$ which are factors of $N$.  This problem is interesting when $N$ is not prime and the solutions $a_1$ and $a_2$ are non-trivial.  In other words, we will attempt to find integers $a_1$ and $a_2$ satisfying $N=a_1\\cdot a_2$ with $a_1\\neq 1$ and $a_2\\neq 1$.\n",
    "\n",
    "The algorithm consists of two ideas:\n",
    "\n",
    "* Reduce the problem of factoring the integer $N$ to the order-finding problem.\n",
    "* Solve the order-finding problem: \n",
    "Given integers $a$ and $N$ so that $a$ and $N$ share no common factors (i.e., the greatest common divisor of $a$ and $N$ is 1), find the smallest positive integer which satisfies $a^r \\equiv 1 \\mod N$.  This value $r$ is refered to as the *order of* $a \\mod N$.\n",
    "\n",
    "These two ideas are combined in the following steps in Shor's algorithm:\n",
    "\n",
    "0. Rule out the easy case that $N$ is even\n",
    "1. Select a random integer $a$ between $2$ and $N-1$\n",
    "2. Check to see if $a$ is a factor of $N$ (if so we're done!)\n",
    "3. Find the order of $a \\mod N$, (i.e., find $r$ so that $a^r\\equiv 1 \\mod N$)\n",
    "4. Check to see if $a^{\\frac{r}{2}}-1$ or $a^{\\frac{r}{2}}+1$ are integers and share common, non-trivial, divisors with $N$ (if so we're done!)\n",
    "5. If no factor is found, repeat steps 1 through 4 above.\n",
    "\n",
    "The function `shors_algorithm` below outlines these steps for both the classical and the quantum implementation of Shor's algorithm. For the purposes of demonstration, we will also control the initial random integer selected in step 1 so that we can investigate cases in which the selected integer produces common divisors of $N$ in step 4 and others in which step 4 produces no common factors.\n",
    "\n",
    "\n",
    "**Note about terminology:** Some literature refers to the \"period-finding problem\" in Shor's algorithm.  The order-finding problem, as we have described above, can be recast as finding the period of the function $f(x) = a^x\\mod N$, by noticing that the period of $f(x)$ is one more than the order of $a\\mod N$. The period finding problem is more general than the order-finding problem since it aims to find the period of any periodic function, not just modular exponentiation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7dd05729",
   "metadata": {},
   "outputs": [],
   "source": [
    "def shors_algorithm(N, initial, quantum):\n",
    "    \"\"\" Factor N using Shor's algorithm with initial starting value and choice of \n",
    "    using either classical or quantum approach for the period finding step\n",
    "    Parameters\n",
    "    ----------\n",
    "    N: int \n",
    "        Integer to be factored (assumed to be non-prime and >1)\n",
    "    initial: int \n",
    "        Initial choice of the random integer between 2 and N-1\n",
    "    quantum: boolean\n",
    "        True if we will call the quantum order-finding algorithm, and false if we call the classical one for step 3.   \n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    a1, a2: int\n",
    "        Non-trivial factors of N\n",
    "    \"\"\"\n",
    "\n",
    "    # 0. Check to see if N is even.\n",
    "    if N % 2 == 0:\n",
    "        divisor1 = 2\n",
    "        divisor2 = N // 2\n",
    "        print(\"Found factors:\", divisor1, divisor2)\n",
    "        return (divisor1, divisor2)\n",
    "\n",
    "    attempts = [initial]\n",
    "    while (True):\n",
    "        # 1. Select a random integer between 2 and N-1\n",
    "        if len(attempts) == 1:\n",
    "            a = initial\n",
    "        else:\n",
    "            a = random.choice(\n",
    "                [n for n in range(N - 1) if n not in attempts and n != 1])\n",
    "\n",
    "        # 2. See if the integer selected in step 1 happens to factor N\n",
    "        print(\"Trying a =\", a)\n",
    "        divisor1 = gcd(a, N)\n",
    "        if divisor1 != 1:\n",
    "            divisor2 = N // divisor1\n",
    "            print(\"Found factors of N={} by chance: {} and {}\".format(\n",
    "                N, divisor1, divisor2))\n",
    "            return (divisor1, divisor2)\n",
    "\n",
    "        # 3. Find the order of a mod N (i.e., r, where a^r = 1 (mod N))\n",
    "        if quantum == True:\n",
    "            r = find_order_quantum(a, N)\n",
    "        else:\n",
    "            r = find_order_classical(a, N)\n",
    "        print(\"The order of a = {} is {}\".format(a, r))\n",
    "\n",
    "        # 4. If the order of a is found and it is\n",
    "        # * even and\n",
    "        # * not a^(r/2) = -1 (mod N),\n",
    "        # then test a^(r/2)-1 and a^(r/2)+1 to see if they share factors with N.\n",
    "        # We also want to rule out the case of finding the trivial factors: 1 and N.\n",
    "        divisor1, divisor2 = test_order(a, r, N)\n",
    "        if (divisor1 != 0):  # test_order will return a 0 if no factor is found\n",
    "            print(\"Found factors of N = {}: {} and {}\".format(\n",
    "                N, divisor1, divisor2))\n",
    "            return divisor1, divisor2\n",
    "\n",
    "        # 5. Repeat\n",
    "        print(\"retrying...\")\n",
    "        attempts.append(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed682e81",
   "metadata": {},
   "source": [
    "Let's first explore the idea of reducing the factoring problem into the order-finding problem. This gets captured in step 4 described above.  In this step, we have already established that $a$ and $N$ share no factors other than $1$ (i.e., $\\gcd(a,N)=1$) and we have found $r$, the order of $a\\mod N.$  With this information we know that $a^r \\equiv 1 \\mod N.$  Rewritten in another way:\n",
    "$$\n",
    "a^r -1 \\equiv 0\\mod N\n",
    "\\tag{1}.\n",
    "$$ \n",
    " \n",
    " If $r$ is even, we can rewrite $a^r$ as $x^2$ where $x=a^\\frac{r}{2}$. Next, we can factor equation $(1)$ using the identity $x^2-1 = (x-1)(x+1)$:\n",
    " $$ \n",
    " (a^{\\frac{r}{2}} - 1)(a^{\\frac{r}{2}} + 1) \\equiv 0\\mod N\n",
    "\\tag{2}.\n",
    "$$  \n",
    "\n",
    "If, in addition, the equation\n",
    "$$\n",
    "a^{\\frac{r}{2}}  \\not\\equiv -1\\mod N\n",
    "\\tag{3}\n",
    "$$\n",
    "is satisfied, then  at least one of the terms $a^{\\frac{r}{2}} - 1$ or $a^{\\frac{r}{2}} + 1$ must share a common factor with $N$. [Peter Shor](https://arxiv.org/abs/quant-ph/9508027) demonstrated that there is greater than a $50\\%$ chance of randomly selecting a value for $a$ satisfying these properties.\n",
    "\n",
    "The code block below defines a function that tests whether $r$ is even and whether equation $(3)$ is satisfied and searches for a non-trivial factor of $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "id": "6c8eed12",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_order(a, r, N):\n",
    "    \"\"\"Checks whether or not a^(r/2)+1 or a^(r/2)-1 share a non-trivial common factor with N\n",
    "    Parameters\n",
    "    ----------\n",
    "    a: int\n",
    "    r: int\n",
    "    N: int\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    int, int factors of N, if found; 0,0 otherwise  \n",
    "    \"\"\"\n",
    "\n",
    "    if r != None and (r % 2 == 0) and a**r % N == 1:\n",
    "        if (a**(int(r / 2))) % N != -1:\n",
    "            possible_factors = [gcd(r - 1, N), gcd(r + 1, N)]\n",
    "            for test_factor in possible_factors:\n",
    "                if test_factor != 1 and test_factor != N:\n",
    "                    return test_factor, N // test_factor\n",
    "    # period did not produce a factor\n",
    "    else:\n",
    "        print('No non-trivial factor found')\n",
    "        return 0, 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a43af433-905c-4070-94b3-0f6be6f67f6b",
   "metadata": {},
   "source": [
    "### Solving the order-finding problem classically\n",
    "\n",
    "The key component of Shor's algorithm is an efficient quantum algorithm to find the order $r$ of $a \\mod N$. While a straightforward classical algorithm can solve this problem, it is notably inefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "id": "14f4b0d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_order_classical(a, N):\n",
    "    \"\"\"A naive classical method to find the order r of a (mod N).\n",
    "    Parameters\n",
    "    ----------\n",
    "    a: int\n",
    "        an integer in the interval [2,N-1]\n",
    "    N: int\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    r: int \n",
    "        Period of a^x (mod N)\n",
    "    \"\"\"\n",
    "    assert 1 < a and a < N\n",
    "    r = 1\n",
    "    y = a\n",
    "    while y != 1:\n",
    "        y = y * a % N\n",
    "        r += 1\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aae6b36d",
   "metadata": {},
   "source": [
    "Let's see how Shor's algorithm works on a few examples using the classical order-finding problem.  Notice how often our choice of value for the initial guess for $a$ produces factors of $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "3f75b314-9f0c-40c5-a7bc-58ae235fc183",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trying a = 42\n",
      "Found factors of N=123 by chance: 3 and 41\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(3, 41)"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_integer = 123  #edit this value to try out a few examples\n",
    "# Edit the value in the line below to try out different initial guesses for a.\n",
    "# What happens when you choose a = 42 for the integer 123?\n",
    "# What happens when you choose a = 100 for the integer 123?\n",
    "initial_value_to_start = 42  # edit this value; it should be less than my_integer\n",
    "\n",
    "shors_algorithm(my_integer, initial_value_to_start, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "275efced-b10e-4438-9b5e-f3534079304b",
   "metadata": {},
   "source": [
    "### Solving the order-finding problem with a quantum algorithm\n",
    "\n",
    "\n",
    "The Fourier transform is a classical computation that provides a more efficient algorithm than the one encoded in `find_order_classical` for identifying the period of $f(x) = a^x\\mod N$. The quantum version of the Fourier transform is the central idea of Shor's Algorithm. This efficient quantum solution derives the period from a measurement of $n = \\lceil log2(N) \\rceil$ qubits. \n",
    "\n",
    "The image below outlines the quantum kernel used to find the order of $a$. Notice the last step involves applying the Inverse Quantum Fourier Transform."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba9fbd7",
   "metadata": {},
   "source": [
    "![](images/shor_circuit.png)\n",
    "\n",
    "Figure. Circuit diagram for finding the phase of the modular multiplication gate $U|x\\rangle = |a^x\\mod N\\rangle$, which will be used to compute the order of $a$ modulo $N$.  The number of qubits in the control register determines the accuracy of estimating the phase of $U$. The size of the work register depends on $N$. The goal of this section is to code the  diagram as a kernel named `phase_kernel`. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7eff3b6",
   "metadata": {},
   "source": [
    "We will need to create a quantum kernel for the Inverse Quantum Fourier Transform.  Additionally we'll need to create a kernel for modular multiplication: $g(y) = ay \\mod N$, which can be repeatedly applied $x$-times to $y=1$ to carry our modular exponentation $f(x)=a^x\\mod N$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "567e104e",
   "metadata": {},
   "source": [
    "#### Inverse quantum Fourier transform\n",
    "In the code block below we define a kernel for the quantum Fourier transform and then use `cudaq.adjoint` to create a kernel for the inverse quantum Fourier transform used in the quantum order-finding algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "0e04a0a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define kernels for the quantum Fourier transform and the inverse quantum Fourier transform\n",
    "@cudaq.kernel\n",
    "def quantum_fourier_transform(qubits: cudaq.qview):\n",
    "    qubit_count = len(qubits)\n",
    "    # Apply Hadamard gates and controlled rotation gates.\n",
    "    for i in range(qubit_count):\n",
    "        h(qubits[i])\n",
    "        for j in range(i + 1, qubit_count):\n",
    "            angle = (2 * np.pi) / (2**(j - i + 1))\n",
    "            cr1(angle, [qubits[j]], qubits[i])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def inverse_qft(qubits: cudaq.qview):\n",
    "    cudaq.adjoint(quantum_fourier_transform, qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39db9e6a",
   "metadata": {},
   "source": [
    "#### Quantum kernels for modular exponentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b8b8884",
   "metadata": {},
   "source": [
    "\n",
    "While Shor's algorithm provably can factor numbers faster than known classical techniques, the resources required for implementing Shor's algorithm are hefty. A full-scale implementation to factor an $L$-bit number could require a quantum kernel with $5L+1$ qubits to achieve accuracy for the continued fractions part of the algoirthm and as many as $72L^3$ quantum gates for the modular exponentiaion [(Beckman, Chari, Devabhaktuni, & Preskill, 1996)](https://arxiv.org/pdf/quant-ph/9602016).  Both of these are well beyond the capabilities of modern quantum hardware for numbers of interest. Advancements to reduce the number of gates and qubits required focus on optimizing the kernel for $f(x) = a\\cdot x \\mod N$ for properties of the number to be factored. The difficulty is to create efficient quantum kernels (in terms of qubit and gate count) that compute $a\\cdot x \\mod{N}$.  \n",
    "\n",
    "For the purposes of this demonstration, we will consider only the order-finding problem for the values of $a$ to be either $4$ or $5$ with $N=21$. We'll be using the quantum circuits depicted in [this paper](https://arxiv.org/abs/2103.13855) and [this report](https://physlab.org/wp-content/uploads/2023/05/Shor_s_Algorithm_23100113_Fin.pdf), respectively. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a48d6fa1",
   "metadata": {},
   "source": [
    "##### The case N = 21 and a = 5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "c4b8e2e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cudaq.kernel\n",
    "def modular_mult_5_21(work: cudaq.qview):\n",
    "    \"\"\"\"Kernel for multiplying by 5 mod 21\n",
    "    based off of the circuit diagram in \n",
    "    https://physlab.org/wp-content/uploads/2023/05/Shor_s_Algorithm_23100113_Fin.pdf\n",
    "    Modifications were made to change the ordering of the qubits\"\"\"\n",
    "    x(work[0])\n",
    "    x(work[2])\n",
    "    x(work[4])\n",
    "\n",
    "    swap(work[0], work[4])\n",
    "    swap(work[0], work[2])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def modular_exp_5_21(exponent: cudaq.qview, work: cudaq.qview,\n",
    "                     control_size: int):\n",
    "    \"\"\" Controlled modular exponentiation kernel used in Shor's algorithm\n",
    "    |x> U^x |y> = |x> |(5^x)y mod 21>\n",
    "    \"\"\"\n",
    "    x(work[0])\n",
    "    for exp in range(control_size):\n",
    "        ctrl_qubit = exponent[exp]\n",
    "        for _ in range(2**(exp)):\n",
    "            cudaq.control(modular_mult_5_21, ctrl_qubit, work)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c644786",
   "metadata": {},
   "source": [
    "Verify in the code block below that the kernels defined above do in fact carry out multiplication and exponentiation with $a=5$ and $N=21$ by changing the `iterations` variable below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "id": "e93fcf39",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     ╭───╮╭───╮      \n",
      "q0 : ┤ x ├┤ x ├─╳──╳─\n",
      "     ╰───╯╰───╯ │  │ \n",
      "q1 : ───────────┼──┼─\n",
      "     ╭───╮      │  │ \n",
      "q2 : ┤ x ├──────┼──╳─\n",
      "     ╰───╯      │    \n",
      "q3 : ───────────┼────\n",
      "     ╭───╮      │    \n",
      "q4 : ┤ x ├──────╳────\n",
      "     ╰───╯           \n",
      "\n",
      "Measurement results from sampling: { 10100:200 }\n",
      "\n",
      "For x = 1, 5^x mod 21 = 5\n",
      "For x = 1, the computed result of the circuit is 5\n"
     ]
    }
   ],
   "source": [
    "# Demonstrate iterated application of 5y mod 21 where y = 1\n",
    "@cudaq.kernel\n",
    "def demonstrate_mod_exponentiation(iterations: int):\n",
    "    qubits = cudaq.qvector(5)\n",
    "    x(qubits[0]\n",
    "     )  # initalizes the qubits in the state for y = 1 which is |10000>\n",
    "    for _ in range(iterations):\n",
    "        modular_mult_5_21(qubits)\n",
    "\n",
    "\n",
    "shots = 200\n",
    "\n",
    "# The iterations variable determines the exponent in 5^x mod 21.\n",
    "# Change this value to verify that the demonstrate_mod_exponentiation\n",
    "# kernel carries out the desired calculation.\n",
    "iterations = 1\n",
    "\n",
    "print(cudaq.draw(demonstrate_mod_exponentiation, iterations))\n",
    "\n",
    "results = cudaq.sample(demonstrate_mod_exponentiation,\n",
    "                       iterations,\n",
    "                       shots_count=shots)\n",
    "\n",
    "print(\"Measurement results from sampling:\", results)\n",
    "\n",
    "# Reverse the order of the most probable measured bit string\n",
    "# and convert the binary string to an integer\n",
    "integer_result = int(results.most_probable()[::-1], 2)\n",
    "\n",
    "print(\"For x = {}, 5^x mod 21 = {}\".format(iterations, (5**iterations) % 21))\n",
    "print(\"For x = {}, the computed result of the circuit is {}\".format(\n",
    "    iterations, integer_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff07579c",
   "metadata": {},
   "source": [
    "##### The case N = 21 and a = 4:\n",
    "\n",
    "The example below is one where the work register has been optimized to use fewer gates and qubits than would have been necessary through a straightforward implementation of modular multiplication as seen in the previous case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "69963075",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cudaq.kernel\n",
    "def modular_exp_4_21(exponent: cudaq.qview, work: cudaq.qview):\n",
    "    \"\"\" Controlled modular exponentiation kernel used in Shor's algorithm\n",
    "     |x> U^x |y> = |x> |(4^x)y mod 21>\n",
    "     based off of the circuit diagram in https://arxiv.org/abs/2103.13855\n",
    "     Modifications were made to account for qubit ordering differences\"\"\"\n",
    "    swap(exponent[0], exponent[2])\n",
    "    # x = 1\n",
    "    x.ctrl(exponent[2], work[1])\n",
    "\n",
    "    # x = 2\n",
    "    x.ctrl(exponent[1], work[1])\n",
    "    x.ctrl(work[1], work[0])\n",
    "    x.ctrl([exponent[1], work[0]], work[1])\n",
    "    x.ctrl(work[1], work[0])\n",
    "\n",
    "    # x = 4\n",
    "    x(work[1])\n",
    "    x.ctrl([exponent[0], work[1]], work[0])\n",
    "    x(work[1])\n",
    "    x.ctrl(work[1], work[0])\n",
    "    x.ctrl([exponent[0], work[0]], work[1])\n",
    "    x.ctrl(work[1], work[0])\n",
    "    swap(exponent[0], exponent[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84a197e2",
   "metadata": {},
   "source": [
    "Now we are ready to define the `phase_kernel` to carry out the instructions in the circuit diagram drawn at the beginning of this section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "206bcac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "@cudaq.kernel\n",
    "def phase_kernel(control_register_size: int, work_register_size: int, a: int,\n",
    "                 N: int):\n",
    "    \"\"\" \n",
    "    Kernel to estimate the phase of the modular multiplication gate |x> U |y> = |x> |a*y mod 21> for a = 4 or 5\n",
    "    \"\"\"\n",
    "\n",
    "    qubits = cudaq.qvector(control_register_size + work_register_size)\n",
    "    control_register = qubits[0:control_register_size]\n",
    "    work_register = qubits[control_register_size:control_register_size +\n",
    "                           work_register_size]\n",
    "\n",
    "    h(control_register)\n",
    "\n",
    "    if a == 4 and N == 21:\n",
    "        modular_exp_4_21(control_register, work_register)\n",
    "    if a == 5 and N == 21:\n",
    "        modular_exp_5_21(control_register, work_register, control_register_size)\n",
    "\n",
    "    inverse_qft(control_register)\n",
    "\n",
    "    # Measure only the control_register and not the work_register\n",
    "    mz(control_register)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "id": "ca9b570a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     ╭───╮                                                                    »\n",
      "q0 : ┤ h ├──●────●────●───●──●────────────────────────────────────────────────»\n",
      "     ├───┤  │    │    │   │  │                                                »\n",
      "q1 : ┤ h ├──┼────┼────┼───┼──┼───●────●────●───●──●───●────●────●───●──●──────»\n",
      "     ├───┤  │    │    │   │  │   │    │    │   │  │   │    │    │   │  │      »\n",
      "q2 : ┤ h ├──┼────┼────┼───┼──┼───┼────┼────┼───┼──┼───┼────┼────┼───┼──┼───●──»\n",
      "     ├───┤╭─┴─╮  │    │   │  │ ╭─┴─╮  │    │   │  │ ╭─┴─╮  │    │   │  │ ╭─┴─╮»\n",
      "q3 : ┤ x ├┤ x ├──┼────┼───╳──╳─┤ x ├──┼────┼───╳──╳─┤ x ├──┼────┼───╳──╳─┤ x ├»\n",
      "     ╰───╯╰───╯  │    │   │  │ ╰───╯  │    │   │  │ ╰───╯  │    │   │  │ ╰───╯»\n",
      "q4 : ────────────┼────┼───┼──┼────────┼────┼───┼──┼────────┼────┼───┼──┼──────»\n",
      "               ╭─┴─╮  │   │  │      ╭─┴─╮  │   │  │      ╭─┴─╮  │   │  │      »\n",
      "q5 : ──────────┤ x ├──┼───┼──╳──────┤ x ├──┼───┼──╳──────┤ x ├──┼───┼──╳──────»\n",
      "               ╰───╯  │   │         ╰───╯  │   │         ╰───╯  │   │         »\n",
      "q6 : ─────────────────┼───┼────────────────┼───┼────────────────┼───┼─────────»\n",
      "                    ╭─┴─╮ │              ╭─┴─╮ │              ╭─┴─╮ │         »\n",
      "q7 : ───────────────┤ x ├─╳──────────────┤ x ├─╳──────────────┤ x ├─╳─────────»\n",
      "                    ╰───╯                ╰───╯                ╰───╯           »\n",
      "\n",
      "################################################################################\n",
      "\n",
      "                                                                            »\n",
      "────────────────────────────────────────────────────────────────────────────»\n",
      "                                                                            »\n",
      "────────────────────────────────────────────────────────────────────────────»\n",
      "                                                                            »\n",
      "──●────●───●──●───●────●────●───●──●───●────●────●───●──●───●────●────●───●─»\n",
      "  │    │   │  │ ╭─┴─╮  │    │   │  │ ╭─┴─╮  │    │   │  │ ╭─┴─╮  │    │   │ »\n",
      "──┼────┼───╳──╳─┤ x ├──┼────┼───╳──╳─┤ x ├──┼────┼───╳──╳─┤ x ├──┼────┼───╳─»\n",
      "  │    │   │  │ ╰───╯  │    │   │  │ ╰───╯  │    │   │  │ ╰───╯  │    │   │ »\n",
      "──┼────┼───┼──┼────────┼────┼───┼──┼────────┼────┼───┼──┼────────┼────┼───┼─»\n",
      "╭─┴─╮  │   │  │      ╭─┴─╮  │   │  │      ╭─┴─╮  │   │  │      ╭─┴─╮  │   │ »\n",
      "┤ x ├──┼───┼──╳──────┤ x ├──┼───┼──╳──────┤ x ├──┼───┼──╳──────┤ x ├──┼───┼─»\n",
      "╰───╯  │   │         ╰───╯  │   │         ╰───╯  │   │         ╰───╯  │   │ »\n",
      "───────┼───┼────────────────┼───┼────────────────┼───┼────────────────┼───┼─»\n",
      "     ╭─┴─╮ │              ╭─┴─╮ │              ╭─┴─╮ │              ╭─┴─╮ │ »\n",
      "─────┤ x ├─╳──────────────┤ x ├─╳──────────────┤ x ├─╳──────────────┤ x ├─╳─»\n",
      "     ╰───╯                ╰───╯                ╰───╯                ╰───╯   »\n",
      "\n",
      "################################################################################\n",
      "\n",
      "                           ╭─────────────╮╭────────────╮╭───╮\n",
      "───────────────────────────┤ r1(-0.7854) ├┤ r1(-1.571) ├┤ h ├\n",
      "        ╭────────────╮╭───╮╰──────┬──────╯╰─────┬──────╯╰───╯\n",
      "────────┤ r1(-1.571) ├┤ h ├───────┼─────────────●────────────\n",
      "   ╭───╮╰─────┬──────╯╰───╯       │                          \n",
      "─●─┤ h ├──────●───────────────────●──────────────────────────\n",
      " │ ╰───╯                                                     \n",
      "─╳───────────────────────────────────────────────────────────\n",
      " │                                                           \n",
      "─┼───────────────────────────────────────────────────────────\n",
      " │                                                           \n",
      "─╳───────────────────────────────────────────────────────────\n",
      "                                                             \n",
      "─────────────────────────────────────────────────────────────\n",
      "                                                             \n",
      "─────────────────────────────────────────────────────────────\n",
      "                                                             \n",
      "\n",
      "Measurement results for a=5 and N=21 with 3 qubits in the control register \n",
      "{ 000:2843 010:913 111:1850 001:1935 011:1830 100:2846 101:1861 110:922 }\n",
      "\n"
     ]
    }
   ],
   "source": [
    "control_register_size = 3\n",
    "work_register_size = 5\n",
    "values_for_a = [4, 5]\n",
    "idx = 1  # change to 1 to select 5\n",
    "N = 21\n",
    "shots = 15000\n",
    "\n",
    "print(\n",
    "    cudaq.draw(phase_kernel, control_register_size, work_register_size,\n",
    "               values_for_a[idx], N))\n",
    "\n",
    "results = cudaq.sample(phase_kernel,\n",
    "                       control_register_size,\n",
    "                       work_register_size,\n",
    "                       values_for_a[idx],\n",
    "                       N,\n",
    "                       shots_count=shots)\n",
    "print(\n",
    "    \"Measurement results for a={} and N={} with {} qubits in the control register \"\n",
    "    .format(values_for_a[idx], N, control_register_size))\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4718d04f",
   "metadata": {},
   "source": [
    "### Determining the order from the measurement results of the phase kernel\n",
    "\n",
    "The measurement results from sampling the `phase_kernel` can be used to estimate the phase of the operator $U|x\\rangle = ax\\mod N$. We will only be interested in the most probable non-zero outcomes of the sampling.  Therefore, we'll create a function `top_results` to extract those values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "9e96846f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def top_results(sample_results, zeros, threshold):\n",
    "    \"\"\"Function to output the non-zero results whose counts are above the given threshold\n",
    "    Returns\n",
    "    -------\n",
    "        dict[str, int]: keys are bit-strings and values are the respective counts  \n",
    "    \"\"\"\n",
    "    results_dictionary = {k: v for k, v in sample_results.items()}\n",
    "    if zeros in results_dictionary.keys():\n",
    "        results_dictionary.pop(zeros)\n",
    "    sorted_results = {\n",
    "        k: v for k, v in sorted(\n",
    "            results_dictionary.items(), key=lambda item: item[1], reverse=True)\n",
    "    }\n",
    "    top_key = next(iter(sorted_results))\n",
    "    max_value = sorted_results[top_key]\n",
    "    top_results_dictionary = {top_key: max_value}\n",
    "\n",
    "    for key in sorted_results:\n",
    "        if results_dictionary[key] > min(threshold, max_value):\n",
    "            top_results_dictionary[key] = results_dictionary[key]\n",
    "    return top_results_dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b55ab5a",
   "metadata": {},
   "source": [
    "Let's see how the `top_results` extracts the top results from the measurement results for a=4 and N=21 with 3 qubits in the control register that we computed above.  One of these top result is likely to be an estimate for the phase of $U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "id": "cb717231",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'100': 2846,\n",
       " '001': 1935,\n",
       " '101': 1861,\n",
       " '111': 1850,\n",
       " '011': 1830,\n",
       " '110': 922,\n",
       " '010': 913}"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example\n",
    "top_results(results, '000', 750)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cda8137a",
   "metadata": {},
   "source": [
    "The algorithm for translating the phase estimate of the operator $U|x\\rangle = ax\\mod N$ into the order of $a \\mod N$ involves continued fractions. The function `find_order_quantum` below carries out this computation.  To learn more about how the phase of $U$ relates to the period of $a \\mod N$, check out these three lectures by Scott Aaronson: lecture [19](https://www.scottaaronson.com/qclec/19.pdf), [20](https://www.scottaaronson.com/qclec/20.pdf), [21](https://www.scottaaronson.com/qclec/21.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "id": "a562fde9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_order_from_phase(phase, phase_nbits, a, N):\n",
    "    \"\"\"Uses continued fractions to find the order of a mod N  \n",
    "    Parameters\n",
    "    ----------\n",
    "    phase: int\n",
    "        Integer result from the phase estimate of U|x> = ax mod N\n",
    "    phase_nbits: int\n",
    "        Number of qubits used to estimate the phase\n",
    "    a: int\n",
    "        For this demonstration a is either 4 or 5\n",
    "    N: int\n",
    "        For this demonstration N = 21\n",
    "    Returns\n",
    "    -------\n",
    "    int: period of a mod N if found, otherwise returns None\n",
    "    \"\"\"\n",
    "\n",
    "    assert phase_nbits > 0\n",
    "    assert a > 0\n",
    "    assert N > 0\n",
    "\n",
    "    eigenphase = float(phase) / 2**(phase_nbits)\n",
    "\n",
    "    f = fractions.Fraction.from_float(eigenphase).limit_denominator(N)\n",
    "\n",
    "    if f.numerator == 1:\n",
    "        return None\n",
    "    eigenphase = float(f.numerator / f.denominator)\n",
    "    print('eigenphase is ', eigenphase)\n",
    "    coefficients_continued_fraction = list(\n",
    "        contfrac.continued_fraction(eigenphase))\n",
    "\n",
    "    convergents_continued_fraction = list(contfrac.convergents(eigenphase))\n",
    "    print('convergent sequence of fractions for this eigenphase is',\n",
    "          convergents_continued_fraction)\n",
    "    for r in convergents_continued_fraction:\n",
    "        print(\n",
    "            'using the denominators of the fractions in the convergent sequence, testing order =',\n",
    "            r[1])\n",
    "        if a**r[1] % N == 1:\n",
    "            print('Found order:', r[1])\n",
    "            return (r[1])\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47511fbe",
   "metadata": {},
   "source": [
    "We are now ready to combine all the elements above into a function to find the order of $a$ $\\mod N$ and test it in the factoring algoithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "id": "7ad0abcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_order_quantum(a, N):\n",
    "    \"\"\"The quantum algorithm to find the order of a mod N, when x = 4 or x =5 and N = 21\n",
    "    Parameters\n",
    "    ----------\n",
    "    a: int\n",
    "        For this demonstration a will be either 4 or 5\n",
    "    N: int\n",
    "        For this demonstration N will be 21\n",
    "    Returns\n",
    "    r: int the period if it is found, or None if no period is found\n",
    "    -------\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    if (a == 4 and N == 21) or (a == 5 and N == 21):\n",
    "        shots = 15000\n",
    "        if a == 4 and N == 21:\n",
    "            control_register_size = 3\n",
    "            work_register_size = 2\n",
    "        if a == 5 and N == 21:\n",
    "            control_register_size = 5\n",
    "            work_register_size = 5\n",
    "\n",
    "        #cudaq.set_random_seed(123)\n",
    "        results = cudaq.sample(phase_kernel,\n",
    "                               control_register_size,\n",
    "                               work_register_size,\n",
    "                               a,\n",
    "                               N,\n",
    "                               shots_count=shots)\n",
    "        print(\"Measurement results:\")\n",
    "        print(results)\n",
    "\n",
    "        # We will want to ignore the all zero result\n",
    "        zero_result = ''.join(\n",
    "            [str(elem) for elem in [0] * control_register_size])\n",
    "        # We'll only consider the top results from the sampling\n",
    "        threshold = shots * (.1)\n",
    "        most_probable_bitpatterns = top_results(results, zero_result, threshold)\n",
    "\n",
    "        for key in most_probable_bitpatterns:\n",
    "            # Convert the key bit string into an integer\n",
    "            # This integer divided by 8 is an estimate for the phase\n",
    "            reverse_result = key[::-1]\n",
    "            phase = int(reverse_result, 2)\n",
    "\n",
    "            print(\"Trying nonzero bitpattern from the phase estimation:\", key,\n",
    "                  \"=\", phase)\n",
    "            r = get_order_from_phase(phase, control_register_size, a, N)\n",
    "            if r == None:\n",
    "                print('No period found.')\n",
    "\n",
    "                continue\n",
    "\n",
    "            return r\n",
    "            break\n",
    "    else:\n",
    "        print(\n",
    "            \"A different quantum kernel is required for this choice of a and N.\"\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "id": "79927898",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trying a = 5\n",
      "Measurement results:\n",
      "{ 11010:452 11001:92 11101:47 00000:2485 01001:96 00100:131 01011:1747 00101:1654 10101:1736 00110:475 11111:23 01000:50 01101:38 01110:35 01010:420 11110:22 01111:26 00011:50 01100:120 10011:35 00010:25 00001:22 11011:1710 00111:79 11100:113 10001:14 10010:29 10100:129 10110:467 10000:2536 10111:81 11000:61 }\n",
      "\n",
      "Trying nonzero bitpattern from the phase estimation: 10000 = 1\n",
      "No period found.\n",
      "Trying nonzero bitpattern from the phase estimation: 01011 = 26\n",
      "eigenphase is  0.8125\n",
      "convergent sequence of fractions for this eigenphase is [(0, 1), (1, 1), (4, 5), (13, 16)]\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 5\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 16\n",
      "No period found.\n",
      "Trying nonzero bitpattern from the phase estimation: 10101 = 21\n",
      "eigenphase is  0.65\n",
      "convergent sequence of fractions for this eigenphase is [(0, 1), (1, 1), (1, 2), (2, 3), (13, 20)]\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 2\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 3\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 20\n",
      "No period found.\n",
      "Trying nonzero bitpattern from the phase estimation: 11011 = 27\n",
      "eigenphase is  0.8421052631578947\n",
      "convergent sequence of fractions for this eigenphase is [(0, 1), (1, 1), (5, 6), (16, 19)]\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 1\n",
      "using the denominators of the fractions in the convergent sequence, testing order = 6\n",
      "Found order: 6\n",
      "The order of a =5 is 6\n",
      "Found factors of N = 21: 7 and 3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(7, 3)"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_integer = 21\n",
    "initial_value_to_start = 5  # Try replacing 5 with 4\n",
    "quantum = True\n",
    "shors_algorithm(my_integer, initial_value_to_start, quantum)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c243a5ee",
   "metadata": {},
   "source": [
    "### Postscript\n",
    "Recent [work of Oded Regev](https://arxiv.org/abs/2308.06572) improves Shor's Algorithm by reducing the number of gates needed.  You can read more about it [here](https://www.quantamagazine.org/thirty-years-later-a-speed-boost-for-quantum-factoring-20231017/)."
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
