factoring_integers_with_Shor's_Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factoring Integers With Shor’s Algorithm\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/bin/bash: /home/ryukijano/miniconda3/envs/cuda-quantum/lib/libtinfo.so.6: no version information available (required by /bin/bash)\n"
     ]
    }
   ],
   "source": [
    "!pip install contfrac==1.0.0 -q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shor’s Algorithm\n",
    "\n",
    "The factoring problem for an integer \\( N \\) aims to find integers \\( p \\) and \\( q \\) which are factors of \\( N \\). This problem is particularly interesting when \\( N \\) is not prime and the solutions \\( p \\) and \\( q \\) are non-trivial. In other words, we attempt to find integers \\( p \\) and \\( q \\) satisfying \\( N = p \\cdot q \\) with \\( 1 < p < N \\) and \\( 1 < q < N \\).\n",
    "\n",
    "The algorithm consists of two main ideas:\n",
    "\n",
    "1. Reduce the problem of factoring the integer \\( N \\) to the order-finding problem.\n",
    "2. Solve the order-finding problem: Given integers \\( x \\) and \\( N \\) such that \\( x \\) and \\( N \\) share no common factors (i.e., the greatest common divisor of \\( x \\) and \\( N \\) is 1), find the smallest positive integer \\( r \\) which satisfies \\( x^r \\equiv 1 \\ (\\text{mod} \\ N) \\). This value \\( r \\) is referred to as the order of \\( x \\).\n",
    "\n",
    "These two ideas are combined in the following steps in Shor’s algorithm:\n",
    "\n",
    "1. Rule out the easy case that \\( N \\) is even.\n",
    "2. Select a random integer \\( x \\) between 1 and \\( N-1 \\).\n",
    "3. Check to see if \\( x \\) is a factor of \\( N \\) (if so, we’re done!).\n",
    "4. Find the order of \\( x \\) (i.e., find \\( r \\) such that \\( x^r \\equiv 1 \\ (\\text{mod} \\ N) \\)).\n",
    "5. Check to see if \\( \\gcd(x^{r/2} - 1, N) \\) or \\( \\gcd(x^{r/2} + 1, N) \\) are integers and share common, non-trivial divisors with \\( N \\) (if so, we’re done!).\n",
    "6. If no factor is found, repeat steps 1 through 5.\n",
    "\n",
    "The function `shors_algorithm` below outlines these steps for both the classical and the quantum implementation of Shor’s algorithm. For demonstration purposes, we will also control the initial random integer selected in step 2 so that we can investigate cases in which the selected integer produces common divisors of \\( N \\) in step 5 and others in which step 5 produces no common factors.\n",
    "\n",
    "**Note about terminology:** Some literature refers to the “period-finding problem” in Shor’s algorithm. The order-finding problem, as described above, can be recast as finding the period of the function \\( f(x) = x^a \\ (\\text{mod} \\ N) \\), by noticing that the period of \\( f \\) is one more than the order of \\( x \\). 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
   "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\\). With this information we know that \\(a^r \\equiv 1 \\ (\\text{mod} \\ N)\\). Rewritten in another way:\n",
    "\n",
    "\\[a^r - 1 \\equiv 0 \\ (\\text{mod} \\ N)\\]\n",
    "\n",
    "If \\(r\\) is even, we can rewrite \\(a^r - 1\\) as \\((a^{r/2})^2 - 1\\) where \\(a^{r/2} \\not\\equiv -1 \\ (\\text{mod} \\ N)\\). Next, we can factor equation \\(a^r - 1\\) using the identity \\(x^2 - y^2 = (x - y)(x + y)\\):\n",
    "\n",
    "\\[(a^{r/2} - 1)(a^{r/2} + 1) \\equiv 0 \\ (\\text{mod} \\ N)\\]\n",
    "\n",
    "If, in addition, the equation\n",
    "\n",
    "\\[a^{r/2} \\not\\equiv -1 \\ (\\text{mod} \\ N)\\]\n",
    "\n",
    "is satisfied, then at least one of the terms \\(a^{r/2} - 1\\) or \\(a^{r/2} + 1\\) must share a common factor with \\(N\\). Peter Shor 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 \\(a^{r/2} \\not\\equiv -1 \\ (\\text{mod} \\ N)\\) is satisfied and searches for a non-trivial factor of \\(N\\)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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",
   "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 of \\(a\\). While a straightforward classical algorithm can solve this problem, it is notably inefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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",
   "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": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trying a = 100\n",
      "The order of a = 100 is 315\n",
      "No non-trivial factor found\n",
      "retrying...\n",
      "Trying a = 1257\n",
      "The order of a = 1257 is 630\n"
     ]
    },
    {
     "ename": "TypeError",
     "evalue": "cannot unpack non-iterable NoneType object",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 7\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;66;03m# Edit the value in the line below to try out different initial guesses for a.\u001b[39;00m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;66;03m# What happens when you choose a = 42 for the integer 123?\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# What happens when you choose a = 100 for the integer 123?\u001b[39;00m\n\u001b[1;32m      5\u001b[0m initial_value_to_start \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m100\u001b[39m  \u001b[38;5;66;03m# edit this value; it should be less than my_integer\u001b[39;00m\n\u001b[0;32m----> 7\u001b[0m \u001b[43mshors_algorithm\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmy_integer\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minitial_value_to_start\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[2], line 56\u001b[0m, in \u001b[0;36mshors_algorithm\u001b[0;34m(N, initial, quantum)\u001b[0m\n\u001b[1;32m     49\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe order of a = \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m is \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(a, r))\n\u001b[1;32m     51\u001b[0m \u001b[38;5;66;03m# 4. If the order of a is found and it is\u001b[39;00m\n\u001b[1;32m     52\u001b[0m \u001b[38;5;66;03m# * even and\u001b[39;00m\n\u001b[1;32m     53\u001b[0m \u001b[38;5;66;03m# * not a^(r/2) = -1 (mod N),\u001b[39;00m\n\u001b[1;32m     54\u001b[0m \u001b[38;5;66;03m# then test a^(r/2)-1 and a^(r/2)+1 to see if they share factors with N.\u001b[39;00m\n\u001b[1;32m     55\u001b[0m \u001b[38;5;66;03m# We also want to rule out the case of finding the trivial factors: 1 and N.\u001b[39;00m\n\u001b[0;32m---> 56\u001b[0m divisor1, divisor2 \u001b[38;5;241m=\u001b[39m test_order(a, r, N)\n\u001b[1;32m     57\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m (divisor1 \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m0\u001b[39m):  \u001b[38;5;66;03m# test_order will return a 0 if no factor is found\u001b[39;00m\n\u001b[1;32m     58\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFound factors of N = \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m and \u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\n\u001b[1;32m     59\u001b[0m         N, divisor1, divisor2))\n",
      "\u001b[0;31mTypeError\u001b[0m: cannot unpack non-iterable NoneType object"
     ]
    }
   ],
   "source": [
    "my_integer = 1349  #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 = 100  # edit this value; it should be less than my_integer\n",
    "\n",
    "shors_algorithm(my_integer, initial_value_to_start, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the Order-Finding Problem with a Quantum Algorithm\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 \\(a\\). 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\\) 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.\n",
    "\n",
    "\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 \\(a\\). 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`.*\n",
    "\n",
    "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: \\(U|x\\rangle = |a^x \\mod N\\rangle\\), which can be repeatedly applied \\(r\\)-times to \\(a\\) to carry out modular exponentiation \\(a^r \\mod N\\).\n",
    "\n",
    "### Inverse Quantum Fourier Transform\n",
    "\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.\n",
    "\n",
    "```python\n",
    "def qft(qubits):\n",
    "    \"\"\"Quantum Fourier Transform on the given qubits.\"\"\"\n",
    "    n = len(qubits)\n",
    "    for i in range(n):\n",
    "        for j in range(i):\n",
    "            cudaq.cp(qubits[j], qubits[i], np.pi / 2**(i-j))\n",
    "        cudaq.h(qubits[i])\n",
    "\n",
    "# Create the inverse Quantum Fourier Transform kernel\n",
    "iqft = cudaq.adjoint(qft)\n",
    "```"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAADXCAIAAADTBBLDAAAgAElEQVR4Ae2deVQUV/bHmRkzE0cnJ3pOcqJZTNQIatyio8GggHaLBlQERUcIiAKSIOpBkUUJGGXTIASEuCASV1YVkQioICoqIKJRFkEBw9LgQhQa2br7/fxZM0VTDU11dVV1dfd9f+h7VW+593tffc571NI6CBIoAAqAAhqqgI6G+gVugQKgACiAAHAwCUABbVegs7MzJCQkiMW0Y8eO9vZ2FnQHwLEgMgwBCnBagdTU1PXr1yeymHg83qFDh1gQBQDHgsgwBCjAaQVSU1P9/PzYNDEqKgoAx6bgMBYooL0KAOC0N/bgOSig8QoA4DQ+xOAgKKC9CgDgtDf24DkooPEKAOA4FOLU1NStW7f6KJL27t0rEok45AOYAgpwSQEuA+7Fixdr165dRymtXbtW/e6iGhkZpaSk/KZIGj16dHV1NZdmFNgCCnBIAS4DLiMjY/369TcopejoaPUDnImJybNnzxSaHZaWlgA4hRSDylqlAMcB5+3tTS0cCQkJADhq0kErUEBzFNAkwF18k7DYqBhwTU1NFOYIrOAoiAZNQAE5CmgM4CZMmKCjoxMZGYk5q0rAPXz4UEdH58mTJ3J07/WULOBu3bplYGBga2tbWlqalZU1YcIEV1fXuro6vDlsUXEpIAMKyCqgGYDr7Ow0NTUdMGCAQCDAfFQl4DZu3BgfHy+rdb9HZAGHENLV1X3w4AFCqKura8iQIUKhULofrgHu8uXL//jHP2bMmFFcXCxtJ+RBAZUooHaAa2pqSktLCwoK2rx584YNG3bu3HnixAknJ6cdO3b8/vvvuIaqBNzSpUsfP36Mm1JRUVFOLhkZGRFuMggEguHDh2Nd5eXl6evr491imcWLF2dnZ5PrnvFaV65ceeedd3T+l27fvs34kNwbICcnh3tGsWTRzZs379+/z9Jg5IbZv3//Dz/8QLhqGC1GRkYGBASQsS4mJsbT0xM3pra2dufOnXPmzHF3d09OTs7Kyrpx48apU6cWLVr04YcfOjs737t3D6+sSsCZm5vX1NTgptjY2PyHXBo/fjwBcHFxcdOnT494k8zMzLy8vPBusQyPx1u0aBG57hmvZWJi8vbbb2N8+9e//rVw4ULGh+TYACtWrBg/fjzHjGLPnBkzZpiZmbE3HomRZs+eTflOJeFaI1ncu3fv9OnTSZj2HyMjo02bNmHd3r59m8/nh4WFdXR0SA+Ulpa2ePHiP//8Mzk5ef78+XFxcdhZDgFO2lz5edktqrOz84EDB7BWCxYsyMjIIPTAtS3q2LFjBw0aNGDAAFNTU4KpUAQF2FdAXbaoPB4P+0uUtER3797l8XhVVVXYwZcvX1pYWDQ0NCCENARwenp65eXlCCGRSDR06FDCH+AQQlwDXFdXl62t7caNG6XjBHlQQFUKqAvgZP9mXV9fb2JicvPmTWnpmpqaysrKNAFwQqFw+/btAwcOjI+PFwgEbm5uw4YNS09Pl/aWg4BDCO3Zs+f48eMEO6EICqhEAXUBnKw4L1++lPNlYA1Zwcm6TTjCtRUcAI4QICiqVgH1BZx83QBw8vVh8Cys4BgUF7pWUAEAnIKCkahOuItKosV/q8jeZOi3Lazg+pUIKmizAgA4+qNPGXA8Ho/P55srkj799FPCkyX0+6Ngj7CCU1AwqM6gAlwGXGFh4fDhw+0ppdGjR6vsZXvKgBMKhTUKpufPnzM4Oyh1DYCjJBs0YkQBLgMOIfTgwYPfqSb1AxwjEWa9UwAc65LDgH0qwHHA9Wk3iRMAOBIiMVAFAMeAqNAlRQUAcBSFk9OM8hZVTp9qdAoAp0bB0nhTAXD0hxgABw/60j+roEdKCgDgKMkmtxEADgAnd4LASfYUAMDRrzUADgBH/6yCHikpAICjJJvcRgA4AJzcCQIn2VMAAEe/1gA4ABz9swp6pKQAAI6SbHIbAeAAcHInCJxkTwEAHP1aA+AAcPTPKuiRkgIAOEqyyW0EgAPAyZ0gcJI9BQBw9GsNgAPA0T+roEdKCgDgKMkmtxEADgAnd4LASfYUAMDRrzUADgBH/6yCHikpAICjJJvcRgA4AJzcCQIn2VMAAEe/1gA4ABz9swp6pKQAAI6SbP9rJJFI/pft/h8AB4Drng2QU6kCALje5RcIBL2fkDq6atUqV1dXqQP/zQLgAHCyswKOqEQBAuDi4+OnTZtmb29fWlrKkD1RUVGHDh1iqHPpbql/8PL1zxH6+vpK99Vr/tatWwA4WWXge3CymsARVSlAABxCyMzMLDU1lTl7VAa4e/fuhYaGnjp1SiKRdHV17d+/Pzw8vKysrL6+Pigo6OrVq/v3709NTRUKhTweb968eQcPHhSLxSEhIWlpaeHh4TU1NdKtXitFO+CSk5M9PDy8FEmhoaFdXV3MRYtCzwA4CqJBE4YU6AtwhKu+tbX14MGDBw4c6OzsjIuLO3nypEgkOnfuXHR0tFAozM3NDQgIuHbt2rFjx9ra2lJSUvbt25eSkoIQysvL++WXX6qrq3H7yQPuxYsX9vb2zpTSmjVreqzgWlpa5s6d29bW5uTkdPv2bQ8Pj4yMjJaWFjMzs0ePHm3ZssXKyqqzs9PAwKC5ufn8+fNbtmwRiUQIoYSEhIkTJ8bFxd25c0e6VVVVFe2AMzIySktLy1Qk6erqSouLq6zCDABOheLD0AQF+gKcWCwmXPWnTp2ysrJCCHl5ed29e9fV1bWwsPDp06fLli1rb2+fNGlSbGxscHDw4cOHU1NTOzo6AgMDk5OTQ0JCJBIJn8/HcIEQIg+4jIyMDRs2FFBKhw4d6gG4I0eOeHp6IoQ6OzslEsknn3yC3R/w8fGJiIgIDg6OjY1FCC1atKi2tjYzM9PDwwNTKjs7e926dQgh2Va0Aw5+F5UwO6EICiipQF+AQwgRrnqJRDJx4sSXL1+GhoZKJJL3338/6U3CLn9DQ0Ps9zmrqqp0dXXnzp2bk5Njbm7u7++flJS0YcOGhoYGzFSFAOft7U3NQeIv26ekpDg4OCCERCJRWVnZZ5999uLFC4SQi4vLlStXgoOD4+LiMMDV1NRcvHjR3d09Jyenvr4+OzsbN4LQCgDXa2xgBderLHBQJQrIB5z0VY8hz9bW9saNG6+3buPGjWtubkYI3b9/HyFkaGjY2dmJbelaW1tLSkrmzZvn6Oh4/vx5hFBdXV1TUxPmIEOAk0gk6W8SNgoRcCKRyMnJKSIiws/P79GjR+np6atWrTp58uTWrVvr6urmzZtnZWWVk5Ojq6sbGBj47NkzQ0PDH3/8USQSOTg4zJo1q6CgACEk3Qoh5OrqOmPGjPLyckLkKN9FhRUcQUkoggJKKkAAXHx8/JgxY6ytrXNycghXPULoyZMn06ZNw0a8ePHihg0bjh8/HhYWdvXq1REjRvj7+yOEAgMD3d3dk5OTIyIiamtr7e3t4+Pj3d3dxWIx1pAhwOnp6eno6ERFRWGjEAGHHW1ubsbtkEgkra2t2HHZf9va2jo6OmSPy2+F1be0tKypqZFt2+8RWcDdunXLwMDA1ta2tLQ0KytrwoQJrq6udXV1eFeWlpbwNzhcDciAAgQFCIAjnJUtSl/1EokEW8RJV+t8k1paWvCD2F4QLzIBuI6ODjMzs7feequxsREbqHfA4UYwmpk6dSrBZ5LDyQIOIaSrq/vgwQOEUFdX15AhQ4RCoXRvADhpNVSbf30n7sMPP/zoo4/c3NxUawmMjiugKODwhpQzSgLu6dOnKSkpO3bs2Lhxo4uLi5+fX2xs7Jo1awICArDNMmaYygBnZ2c3fvx4aXVKSkqKySUjIyPsD5l4c4FAMHz4cKyYl5enr6+Pn8IyixcvvnDhArnuWarl6uq6fft2lgbj0jAWFhZ/+9vfdHR0/vrXv3p6enLJNJZsuXz5cmFhIUuDkRsmMjLyhx9+IFw1jBYjIyN37NhBxroDBw54vrnzidlTU1Pj5+c3Z84cLy+vlJSUK1eu5Ofnnz17dvHixR9//LGjo+OdO3dwy1UDOIlEoqOjc/LkSdwOhJC9vb0tufTFF18QABcXFzd9+vSIN8nMzMzLy0u6Z4QQn8+3sLAg1z1LtT7//PNJkyaxNBiXhnnnnXd03qS//OUvEydO5JJpLNny9ddfL1myhKXByA1jbGyM3yQkXDsMFffu3auvr0/GOh6Pt2nTJsyMgoICPp+/d+9ewmOtZ8+eNTc3b25uTklJWbBgwbFjx7D6qgEcQqihocHIyOj1uxAU5JPdojo7Ox84cADrasGCBRkZGYRuYYtKEESFxczMzAEDBrz99ts6OjpVVVUqtASGxhXoa4va1NQk/ec2vL7yGWpbVD6fX1FRQRi9qKiIx+M9fvwYOy4UCpcuXYq9RUoRcO7u7kVFRYRhGhsb8VsThFO9FhctWoQ/FNNrhb4OygJOT08Pu0srEomGDh1K+AMcQggA15eYKjleUFAwatSoP//8UyWjw6CyCsgCrqKiYtWqVYmJiYcPH3Zzc8NWTMHBwR9//LG/v//PP//83ZuUlZW1ZMmSkJCQWbNm2dvbBwQEzJw5U7Z/2SPUAFdWVkboqra2dt68efn5+dLHX7x4gdWkCLjr16/Lct3Pz6+trU16GPl5Wh4TEQqF27dvHzhwYHx8vEAgcHNzGzZsWHp6OmFoABxBENUWxWLx9OnTVWsDjC6tAAFw7e3t48ePx9dKMTExTk5OWP1Bgwbh90bd3NxiY2OxW5bLly/HHpfbt2/fy5cvpTvvNU8NcLJdCYVCwnZVuk4PwEm/SkZ4fayxsTE6OvrEiRPY66jBwcElJSXSr5slJSUNHjw4PDwc+wsf/urZ48ePAwMDc3Jy9uzZQ7CDFsBJOyMnD4CTIw77pwBw7Gsuf0QC4E6fPi29EHv16tWgQYOwNc3gwYNbWlpu3759//79wsJC/C/6OOAqKyvJPB1BF+Dk+9UDcNKvkhFeHzM2Nm5ubl65cmVWVtbrbbmXl1dycrL062ZisVhfX18oFEokEum2nZ2d69atW716dUREBOFhGQCc1n4uCQAn/7Jk/ywBcLt27VqxYoW0Ge+99x7GssGDBwcFBY0bNy4vL0+6Ag446YNy8ioAHPamBXaDkvD6GI/HEwgEDg4OmJM//fRTcnKy9OtmCKGZM2diW1RCW19f36SkJFlXAXAAONlZAUdUogABcOnp6Xw+H7ekq6tr6NCh0iu469evFxYW4hUQQmoDOOxVMsLrYwEBAVFRUbm5uZhLGOASEhLw180QQgYGBq2trfHx8YS2vr6+mZmZ0lpgecqA4/P5xsbGZoqkkSNHEp4skbWH5SPa/C4qrOBYnmz9DkcAnFAo1NXVbWhoOHXq1KVLl/bt22djY4N1gm1RZTtkDnBFRUUffPDBt5TSyJEju78mIv0qGeH1MUdHx/nz55uYmOzevbuhoWH+/PnW1tYBAQH462YIIT8/v61bt545c0a6bWVlpbGxsZ2dney3fykD7tWrVwIFE5k/CsjGjNEjADhG5YXOFVKAADiE0IMHD5YvX469J7BixYqWlpbm5uaQkJCBAwdu27atsrIS718oFIaGho4bN87CwiIxMRE/Lj9DfouKEHr06FEJ1dQNOFmDMC6Ul5e7u7sjhMRi8dq1a+/evYvVlH3dTPruSb9MoQw4WTvV8QgATh2jpqk2ywIO87S6ujonJ8fCwiIgIAD7fAhdCigEOGUGlQc4rN/29nYXF5fIyMiYmBgfH59ef0GGggUAOPgbHIVpA02YUKAvwGFjrV692sjISPbJMGUs4RDgMDdEIpFCz/H26zwADgDX7ySBCuwoIB9wTNjAOcDR7iQADgBH+6SCDqkpAICjppu8VgA4AJy8+QHnWFQAAEe/2AA4ABz9swp6pKQAAI6SbHIbAeAAcHInCJxkTwEAHP1aA+AAcPTPKuiRkgIAOEqyyW0EgAPAyZ0gcJI9BQBw9GsNgAPA0T+roEdKCgDgKMkmtxEADgAnd4LASfYUAMDRrzUADgBH/6yCHikpAICjJJvcRgA4AJzcCQIn2VMAAEe/1gA4ABz9swp6pKQAAI6SbHIbAeAAcHInCJxkTwEAHP1aA+AAcPTPKuiRkgIAOEqyyW0EgAPAyZ0gcJI9BQBw9GsNgAPA0T+roEdKCgDgKMkmtxEADgAnd4LASfYUAMDRrzVlwMXHx2/atGmLIik4OJjwq6z0+6Ngj/DJcgUFg+oMKsBlwDU1NVlbWztQSnZ2dv1/spwhXSkDzsjIKD09/ZIiaezYsdXV1Qw5Qq1bABw13aAVEwpwGXAZGRkbN24sopRiYmLUD3AmJiaK/gYg/LI9E1cF5T7hZwMpS8dQQ44Dztvbm5rjPX7ZnloXlFtRXsEB4ChrzpGGADiOBAI3Q2MAJxaLz5w5c+7cOcw1ABweYrYzsEVlW3EYr28FNAZwo0eP1tHROXDgAOarhgDu1q1bBgYGtra2paWlWVlZEyZMcHV1raurwwMKW1RcCi5kYAXHhShI26AZgGtvbzc1Nf373//e2NiIeachgEMI6erqPnjwACHU1dU1ZMgQoVAoHT8AnLQaKs9rM+DCw8Pfeuutr776qra2VuWBwA1QO8A1NDQkJyf7+fm5uro6Oztv27YtOjrazs5u165dJSUluF8cAtzdu3fvkEuGhoaEmwwCgWD48OGYV3l5efr6+riHWGbRokW//fYbue5ZqvXdd9/5+PiwNBjHhrl9+/bnn3/OMaPYMOfXX38dNGiQzv8SG0OSGyM8PPyHH34gXDWMFiMjI/38/MhYFxUV5enpiRtTVVW1bdu2uXPn+vj4pKWl5ebm3r59+/XVbW5uPmLEiNWrVxcWFuKVVQm4JUuWSO8inZyc1pBLEyZMIAAuLi5u+vTpEW+SmZmZl5cX7iGW4fP5VlZW5LpnqZaent6XX37J0mAcG8be3v6f//wnx4xiwxwejzdw4ECMb++88w4bQ5Ibg8/nU75TSbjWSBb37t1rYGBAxjoTE5NNmzZh3ebl5fH5/H379olEIumBTp8+bWlp2dLSkpaW9s0338TGxmJnVQm4ZcuWFRcXS1tJMi97F9XZ2Rn/s+KCBQsyMjIIXcEWlSCIaotau0UViUTvvvsutojbtWuXaqMgPbq6bFF5PF5VVZW05QihW7du8Xg8fMvf2tq6fPny+vp6hJAqAffzzz8PGTKEYCuZoizg9PT0ysvLEUIikWjo0KGEP8AhhABwZIRlrY7WAg4h9OrVKxMTk8jISNbUJjOQugCuoqKC4E5NTQ2fz79165b08ZcvX5aVlakYcAihqKgoabNI5qUBJxQKt2/fPnDgwPj4eIFA4ObmNmzYsPT0dEJXADiCIKotajPgEEIuLi43b95UbQgIo6sL4AhmI4Ta29sJ21XpOqpcwUnboVBeGnAkGwLgSArFTjUAHAAuKirq0KFDZOZbRkYG5b8PAuDIKMxIHXjQlxFZ1aFTWMFhuzcAXO+zFVZwveuiPkdhBQcrOFjB9Xm9mpiYzJo1a74iaeTIkc+fP++zR1WcgBWcKlTnxJiwglNoBXfnzp1hw4b9h1IaMWKEyr4mQnmutbe3P1UwtbS0UB6OoYYAOIaE5X63ADiFAIcQqq6uLqea1A9w3J/BZCwEwJFRSSPrAOAUBZwy0wAAp4x61NsC4Khrp+YtAXAAODWfwiTMB8CREEkzqwDgAHCaObOlvQLASauhVXkAHABO8yc8AE7zY9yHhwA4AFwfU0ODDgPgNCiYirnCQcCVlJQMGTJkFovp7bffzs/PV0w4SrXhJgMl2ZRuBIBTWkJ17YCDgFNXKUnYDYAjIRIDVQBwDIiqHl0C4NiMEwCOTbW7xwLAdWuhZTkAHJsBB8CxqXb3WAC4bi20LAeAYzPgADg21e4eCwDXrYWW5QBwbAYcAMem2t1jAeC6tdCyHACOzYAD4NhUu3ssAFy3FlqWA8CxGXAAHJtqd48FgOvWQstyADg2Aw6AY1Pt7rEAcN1aaFkOAMdmwAFwbKrdPRYArlsLLcsB4NgMOACOTbW7xwLAdWuhZTkAHJsBB8CxqXb3WAC4bi20LAeAYzPgADg21e4eCwDXrYWW5QBwbAYcANdD7ezs7EWLFq1iPhkbG5uYmDA/zqqVK1feu3evh5N9FFpaWuzs7Fgwyc7ObuzYsSwMtGrVKg8Pjz7cJR4+cuSIlZUVC1bNmDHD3NychYFMTU3r6+uJfmpZGQDXI+CbN28+e/bsPQ1KsbGx/v7+PZzso1BUVGRjY6NBrv+/K5999lkf7hIPr1ixIi8vT5Pcd3JyOnr0KNFPLSurH+D8/PxMTU0tFEnLly9va2sjE9ktW7ZcunSJTE11qZOYmBgUFETG2qKiIgcHBzI11ajOjBkzSFprbW1dVlZGsrJaVAsLCwPAqR/gFixYUFFR8ViR9PXXX5NcqwPg1OLSJW8kAI68VhpZU/0AZ2ZmJhAIFArGwoULAXD9KgYrOFjB9TtJ1K4CAK5HyGAF10MO9S/ACk79Y6iUB6oE3O3bt1+8eKGo+bCCU0gx+BscSbngb3AkhVKvaioD3KVLl3R0dNrb2xXVSxZwBQUF8+bNMzAwqKurKyoqGjlypIeHh0QiwXuGLSouhZwMbFFhiypneqjpKZUBbu3atVevXqWgmizgEEIXL17U09NDCO3Zs6e4uJjQLQCOIEivRQAcAK7XiaHWB1UGOAsLi9raWly748ePHyWXZs6c2etNhrFjx65Zs6awsBDvE8/weLyIiAgy3VtYWGjeYyIrV64k4/vOnTu//fZbXDTNyEycOJGM70ePHuXxeBoGuD179jg7O5N0n3K1kpISLk8VlQHO3Ny8pqYGlyY8PDyMXJoxY0avgFu3bp2pqSneoXRm7ty527dvJ9O9qamp5gFuyZIlZHzfvHmzjY2NtG4akP/iiy/I+B4WFmZoaKh5gLOxsSHpPuVqvS4puDNzuAI48or0ukU9c+ZMXFzcu++++8cff8h2BVtUWU1kj8AWVcMABw/6IoQ0AXCZmZmJiYkIIRcXFy8vL9lLFwAnq4nsEQAcAE52Vqj7EfUGnEAgWL9+/ahRo7DHTezs7N57773jx48TogKAIwjSaxEAB4DrdWKo9UH1BhxJ6QFwZIQCwAHgyMwT9aoDgOsRL3iToYcc6l+ANxnUP4ZKeaB+gDM1Na2pqXmhSCL/YSwAnFKziXuNAXDciwmrFqkf4NavX6+vrz9XkWRoaNja2kpGVxYAFxERQcYSuupw6lWtxsbG+Ph4ulwj0w+nAMdy6OEuqlreRSUzrSnXYQFwLi4u27dvP3ToUGlpKWU7yTfkFOAaGhq+/vrr8PDw8+fPt7S0kPeCck1OAY7l0APgAHDEC4cFwF2/fh0h1NrampmZGRERsXnz5n379v35559EU2gqcwpwLS0tv//+O0KooqLi8OHDu3fv9vDwuHbtGk2+9tINpwDHcugBcAA44iXBAuCwIcVicUpKCo/HGzFiREBAwMuXL4mm0FTmFOBwn2pqary9vYcNG2Zubn7lyhX8OO0ZTgEO84610APgAHDEC4oFwD1//jw4OPjTTz+dM2fOqVOnRCIR0Qhay1wD3OXLly0tLT/44ANPT8/Hjx/T6msvnXEKcCyHHgAHgCNeEiwAzsHBYcqUKUVFRcSxmSlzCnC1tbUDBgzw9/cn+RMZykvCKcCxHHoAHACOeAWxALh9+/bl5+dHRkaGhYXt3bv3+vXr9+7do/blKKL1vZU5BbjGxsaDBw8mJCRgr3afPHmyrq4uOjq6N8PpOcYpwLEcegAcAI54FbEAOOkhRSJRfn6+jY3Nb7/9Jn2cxjynAEfwq7Gx8fjx49988w3hOI1FTgFO2i8WQg+AA8BJT7n/z7MMOGz4xsZG7N4i0Ro6ylwGHObfxYsX6XC09z44CzjMXEZDD4ADwBGvCpUAjmgErWXuA45Wd4mdcRxwRHNpLQPgAHDECQWAIyqi5mUAnJoHUFnz1e9VLWU9ltseACdXHvU7CYBTv5jRajEAroecKgFcc3Mzc0+EcX+Lev/+/R4xoLXAccAxGnrYosIWlXgxsQC4yMjI73omMzOz2NhYoik0lTkFuOrq6p6uf+fs7Kyvr0+Tr710wynAsRx6ABwAjnhJsAA4Dw+PgoKC/Px8Ly+vmzdv5uXlff/99/n5+URTaCpzCnAVFRW+vr537twJDg7+7bff8vLyLl26ZG1tTZOvvXTDKcCxHHoAHACOeEmwADhsyOLi4ry8PCxfU1Nz4sQJoik0lTkFONynmJgYPL9r1y48T3uGU4DDvGMt9AA4DgFuxowZU5VLurq6H3/8sXJ9TNXV1WXnZwPFYrGVldXRo0djY2MtLCwYfdl+1KhRZGQZO3aslZUV7YjptcMTJ054enqeOnXK398/ODi41zq0HBw9ejQZ36dOnaqnp8fOJ8tZC31ISMinn35K0n3K1Y4ePUpLpBjqBG4y9BCWtRUcQqi+vn7nzp1eXl4PHz7sYQStBW6u4BBCZ8+edXZ2TkhI6OzspNXjHp1xcAXHWuhhBcehFVyPWam6ApuAKy4u/u6778zMzFJTU5nzmJuA6+zsjI6O5vP5bm5u2C+iMaQANwHHTugBcAA44mXFGuDEYvG6desuX75848aNTZs2FRcXE02hqcxNwB07diw2NraoqOj06dOenp40+dpLNxwEHGuhB8DRDziJRNLLLOvtkLm5eU1NTW9nVHmMNcDl5+dfvnwZc1UgEMj+litdKnATcLt378YdDAkJwfO0ZzgIONZCD4CjE3BVVVVTp07NzMwkOUe1HHBCoXDZsmVJSUnJyckrVqxg9CZDUFAQmaCw+buooaGhu3fvPn/+fFRUlL+/PxnzqNXhIOBYCz0Ajk7AvZ5/vr6+ADjy12FFRcWWLVtcXFwYfZSfmys4iURy5MiRlStXxsTEMPrxSw4CDvtJChZCD4CTBzihUHjwTUpPT4+Ojq6qqoqNjU1OTohOvKoAAA2mSURBVO7q6tq/f394eHhZWZlYLA4JCUlLSwsPD3+938QAV1BQcPDgwZKSEvmXupav4KTF6ejoYO4P7dwEnLT7jY2N0kV689wEHO4jo6EHwMkDHEIoOjra1dX12bNnY8aMEYvFbm5uxcXFHh4eGRkZLS0tZmZmVVVVCQkJEydOjIuLu3PnDga4gICAo0eP9nvvHwB38eLF8PDw0NDQzZs3p6Wl4ZOe3gw3Affs2bP9+/eHhYWFhoba29vT67J0b9wEHDuhB8D1A7hXr15NnDjx3r17kydPLiwsDA8Pl0gkn3zyCXYnwcfHJyIiIjs7e926ddiU8vX1Xbhw4bfffis9w/rKazngnj59am9v/+uvvx45cuTHH388d+5cX0IpeZybgAsICIiKijryJjk6Oirpo5zmHAQca6EHwPUDOISQra2tjY3NjRs3li5deubMGYTQZ599hu2nXFxcrly5kp2d7e3tjc0wX1/fCxcu2NjYkLlctRxwAoEAfz3r2bNn9+7dk3OVKnOKm4ALDw/v6urC/MrOzlbGQfltOQg41kIPgOsfcJmZmR4eHgihf//73x0dHQih9PT0VatWnTx5cuvWrWKx2MHBYdasWQUFBdXV1cbGxk5OTqmpqe+9997Zs2flzzwtBxxCiM/nGxgYLF261NDQUNt+k+Hq1avvv/++hYXF0qVLTU1N5U8VZc5yEHCshR4A1z/gEEIY17B/sakmkUhaW1uVmXYIoYULFzL612Vq5rH2HFxxcfFPP/308OHDR48epaenX7hwgZrB/bbi5grOx8fn7t27j96kgICAfr2gXIGDgGMt9AA4UoCjPLfkNCwvLx88eLA0NOVUZvMUa4ATCoXBwcGvXr3C3kwsKChgyE1uAi4mJga/z56SksKQ7wghDgKOtdAD4FQDuObmZh0dHewvevjMNjc3X6hcmj179pQpU5TrY+GUKVPY+ZpIaWnpBx98MGrUqC+//PKjjz5idIv6+h4RGVlmz569fPlyPCKMZlxcXIYMGTJ58uQpU6bMnDmTubHGjBlDxveFCxdOmjSJna+JsBb6kJAQ5a+IftU7ffo0c+FTvmfVfE1k48aNdnZ2YrEYd0AgENQrl6qrq8vLy5Xro37dunXsAK6xsTErKwtzv66ujtEV3NatW8nIkpmZuXr1ajwijGbi4uLw/hldwU2bNo2M7/X19UuXLmUHcKyFPjQ0NCIigqT7lKsJhUI8lBzMqAZw2N9Zm5qauKYIa1tU1hzn5haVNfc5uEVlzXfYolLcorq7uxcVFRHi1NjYKL0iI5yVLVpaWmrzy/aygjB0BABHUlhra2t2VnAk7VG+GgCOIuCuX78ue3/Az89PoZcK4TER5WcwmR4AcGRUev1+KACOpFDqVa17i5qbmxsQEHDt2rVjx44hhPLy8n755Zfq6mqEUGNjY3R09IkTJ4KCgq5evRocHFxSUtLW1paSkrJv376UlJSkpKTBgweHh4e/fmFLuu3jx48DAwNzcnL27NmDP9iJCQSAY2eiAOBI6gyAIymUelXrBlx7e/ukSZNiY2ODg4OTk5NDQkIkEgmfzxeJRMbGxs3NzStXrszKympqavLy8kpOTj58+HBqampHR0dgYKBYLNbX1xcKhRKJRLptZ2fnunXrVq9eHRER0dzcLC0NAE5aDebyADiS2gLgSAqlXtW6AYcQMjQ0fPbsGULI3Nzc398/KSlpw4YNDQ0NPB5PIBA4ODhgC7SffvopOTm5qqpKV1d37ty5OTk5CKGZM2diW1RCW19f36SkJFlRAHCymjBxBABHUlUAHEmh1KsaEXDYV0AcHR3Pnz+PEKqrq2tqasJejc7NzcV8wwCXkJDQ2tpaUlIyb948hJCBgUFra2t8fDyhbV8fiQPAsTNRAHAkdQbAkRRKvap1A+7q1asjRozAPq9aW1trb28fHx/v7u4uFosdHR3nz59vYmKye/fuhoaG+fPnW1tbBwQEuLu7JycnR0REvL5b4efnt3Xr1jNnzki3raysNDY2trOzEwgEBF20DXCtra0///yzn5+f7A1ogjL0FjkFuK6urufPn9ProPzeuPCYyB9//HH48OHKyspr167Jt5bes3AXtZ+7qNhXQ8rLy93d3RFCYrF47dq1d+/excLQ+Sa1tLTgUZH+7na/X3DUNsAlJCQ8efKkra0tISEBV4yFDBcAJxQKvby8tm7dumHDhps3b7LgNT4EFwC3f//+3Nzc2NhY7DrCbWM6A4DrB3BYANrb211cXCIjI2NiYnx8fMj/rIz8+Gkb4G7evLljxw5/f//IyEj5ytB7lguAQwhhXGtubmbu98N61Y0LgDt58uT169cfPXq0ZcuWXo1k6CAAjhTgMPVFIpFCz/H2GzNtAxxC6NWrVy9evMBuxSj0zGC/YsqpwBHAYRZ2dnZiE6m9vV2OzTSe4gLgsN0PQgj7Bo/yX+IhqQ8ATgHAkdSUfDUtBBxCKDIyks/nOzo6lpWVffPNN+TlolyTO4CrrKycP3/+nDlzCgsLHR0db9y4Qdkp8g05AjiE0IYNG3g8XlBQ0JUrVxwcHMi7QLkmAI5DgPPx8fFWLjk7O69YsUK5PrznzJnD6Mv2Bw4ceP2d5MrKytDQUEtLS8pzl3zDxMTEefPmkZFlzZo1K1euJN+zojWPHj3a1tYmEomCgoK8vLyuXLmiaA8U6o8bN46M797e3jNnzmTuVa2Ojo7Y2FiE0OXLl48dO2ZtbU3BF0Wb7NmzZ+HChSTdp1yN5TsniorQfRdV0ZZK1ies4C5cuJChXEpMTIyJiVGujwxra2tGAXfs2DHsh9yrq6vNzMyU1JBM88TEREdHRzKyREVF2dnZkemTWp2SkhInJ6eqqiqJRPL999+zA7jJkyeT8T0jI2PBggXMAQ4hFBQUhD1vkJub6+TkRE1DhVrt2bNny5YtJN2nXA172Ukhw9iszBXAsemznLFY+JoI/soa9jfNrq4uuu7b9OoXd7aoCCHcd/zPUrIvNffqBeWD3NmiIoREIhHmCDuhhy0qh7aolGcwvQ1ZABxu8K5duy5cuFBZWbl7925t+9GZa9euHTx48OHDh+fOncPWNbgs9GY4BTjcNXZCD4ADwOFT7r8Z1gDX0tJia2u7adOmr776CvttRqIpNJU5tYLDfQoLC9uxY4elpeWECRNcXFzw47RnOAg41kIPgAPAES8o1gDX0NCAP/Gblpb2+++/E02hqcxNwB08eBB7UKa1tVXbVnCshR4AB4AjUoQ1wL3+i8y2bdtsbW2PHz9+8uRJRn+TISgoiOhnb+WioiJ2Hl/A7iPPmTMnMDDw9SvPvr6+vZlDzzEOruBYCz0ADgBHvIrYBNzrT7BkZmZu2bIlJiaGufsM3FzBIYSePHkSFhbm7e1dWlpKDAN9ZW4Cjp3QA+AAcMQriWXAEYdnoMxZwDHgay9dchZwvdhK9yEAHACOOKcAcERF1LwMgFPzACprPjwH10NBAFwPOdS/AIBT/xgq5QEArod8rAFO+odQy8vLmftKGge3qC9fvpR+Z4DRDyhxEHCshR62qLBF7UG31x/8YA1wOTk52MdzSkpK1qxZw9wD/RwEXFdX19q1a7Fnmz09PS9evEgMA31lDgKOtdAD4ABwxCuJNcBh7107OTkxSjeEEAcBhxDq7Oxcu3atk5PThQsXiDGgtcxBwLEWegAcAI54MbEJuNLSUgsLC6a/8spNwGGLZQsLC6a/f8lNwLETegAcAE5lgCstLXV0dOzo6MjOzsa+L0I0haYyNwHn7e198eJFbB3HKOM4CDjWQg+AA8ARKcLaCu7cuXP4391yc3OfPHlCNIWmMgcB9/TpU/xbSZ2dnWlpaTT52ks3HAQca6EHwAHgiJcEa4AjDsxYmYOAY8zXXjrmIOB6sZKZQwA4ABxxZgHgiIqoeRkAp+YBVNZ8eA6uh4IAuB5yqH8BAKf+MVTKAwBcD/kAcD3kUP8CAE79Y6iUBwC4HvIB4HrIof4FAJz6x1ApDwBwPeQDwPWQQ/0LADj1j6FSHnAFcO3t7W3KpdbW1paWFuX6aHNzc2P0V7WUihWlxomJiTt27CAjy82bN9esWUNpEO42mj59Ohnf29raVqxYIf2GLHddIm1ZaGjooUOHSLpPuRr+Szqk7WK1osoAZ2FhUVdXh/tqZGRkoFyaNGnS6NGjlevDYPz48ZoHOD09PTKyTJ482crKCmlW+vzzz8n4bmBgMG7cOA0DXEhIyJgxY0i6T7laXFwcl6eMygC3bNmyBw8ecE0a2KJyLSJK2gNbVCUFVPfmKgOcu7v7tGnTuCYfAI5rEVHSHgCckgKqe3OVAU4ikXh6ejL3WwTUAgOAo6YbZ1sB4DgbGnYMUxng2HFP0VEAcIoqxvH6ADiOB4hp8wBwPRQGwPWQQ/0LADj1j6FSHgDgesgHgOshh/oXAHDqH0OlPADA9ZAPANdDDvUvAODUP4ZKeQCA6yEfAK6HHOpfAMCpfwyV8gAA10O+wMDA8ePHv36oXWPSiBEjjhw50sPJPgpVVVUffvihxjiOOTJy5Mg+3CUe/vbbb2fPnq1J7r///vs5OTlEP7WsDIDrEfDOzs4KjUtisbiHk30X6urqNMx78r/HKBQKNcz3R48e9R1qbTkDgNOWSIOfoIAWKgCA08Kgg8uggLYo8H8mADclStf0rgAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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",
   "metadata": {},
   "source": [
    "### Quantum Kernels for Modular Exponentiation\n",
    "\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 \\( n \\)-bit number could require a quantum kernel with \\( 2n + 3 \\) qubits to achieve accuracy for the continued fractions part of the algorithm and as many as \\( 72n^3 \\) quantum gates for the modular exponentiation (Beckman, Chari, Devabhaktuni, & Preskill, 1996). 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 \\( a \\) 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^x \\mod N \\).\n",
    "\n",
    "For the purposes of this demonstration, we will consider only the order-finding problem for the values of \\( N \\) to be either \\( 15 \\) or \\( 21 \\) with \\( a = 4 \\). We’ll be using the quantum circuits depicted in this [paper](https://arxiv.org/abs/quant-ph/9602016) and this [report](https://arxiv.org/abs/quant-ph/9503016), respectively.\n",
    "\n",
    "#### The Case \\( N = 21 \\) and \\( a = 5 \\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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": "code",
   "execution_count": 8,
   "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",
   "metadata": {},
   "source": [
    "#### The Case \\( N = 21 \\) and \\( a = 4 \\)\n",
    "\n",
    "The example below demonstrates an optimized work register that uses fewer gates and qubits compared to a straightforward implementation of modular multiplication, as seen in the previous case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": "code",
   "execution_count": 10,
   "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": 12,
   "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:2989 111:2036 010:1011 100:3082 001:1938 101:1962 011:1981 110:1001 }\n",
      "\n"
     ]
    }
   ],
   "source": [
    "control_register_size = 3\n",
    "work_register_size = 5\n",
    "values_for_a = [4, 5]\n",
    "idx = 1  # changing to 1  selects 5\n",
    "N = 21\n",
    "shots = 16000\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",
   "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. 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": 13,
   "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": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'100': 3082,\n",
       " '111': 2036,\n",
       " '011': 1981,\n",
       " '101': 1962,\n",
       " '001': 1938,\n",
       " '010': 1011,\n",
       " '110': 1001}"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example\n",
    "top_results(results, '000', 750)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": "code",
   "execution_count": 16,
   "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trying a = 5\n",
      "Measurement results:\n",
      "{ 10111:89 10000:2438 11010:441 10100:118 10110:415 00010:30 01101:40 01010:431 01110:15 10101:1685 11111:17 00110:460 10010:25 01000:53 11000:47 00100:122 00111:98 11100:127 11011:1780 11110:38 01111:24 00001:19 10001:13 10011:58 01100:108 00000:2510 01001:97 00101:1760 01011:1754 11101:51 11001:88 00011:49 }\n",
      "\n",
      "Trying nonzero bitpattern from the phase estimation: 10000 = 1\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": 17,
     "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cuda-quantum",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}