shor_standard_qft.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This code has been adapted and modified from IBM Qiskit 2021 and also from https://github.com/ttlion/ShorAlgQiskit.\n",
    "# It uses the implementation as contained in the work of Stephane Beauregard (https://arxiv.org/abs/quant-ph/0205095)\n",
    "# Many thanks to IBM Qiskit team, Tiago Miguel (ttlion), Qubit by Qubit, Peter Shor and Stephane Beauregard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-18T07:04:42.950393Z",
     "start_time": "2021-05-18T07:04:42.931457Z"
    }
   },
   "outputs": [],
   "source": [
    "from typing import Optional, Union, Tuple, List\n",
    "import math\n",
    "import array\n",
    "import fractions\n",
    "import logging\n",
    "import numpy as np\n",
    "from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister, execute, IBMQ, transpile,BasicAer, Aer, assemble\n",
    "from qiskit.circuit import Gate, Instruction, ParameterVector\n",
    "from qiskit.circuit.library import QFT\n",
    "from qiskit.providers import BaseBackend, Backend\n",
    "from qiskit.quantum_info import partial_trace\n",
    "from qiskit.utils import summarize_circuits\n",
    "from qiskit.utils.arithmetic import is_power\n",
    "from qiskit.utils.validation import validate_min\n",
    "from qiskit.utils.quantum_instance import QuantumInstance\n",
    "import qiskit.visualization\n",
    "\n",
    "from qiskit.providers.aer import QasmSimulator\n",
    "from datetime import datetime\n",
    "import csv\n",
    "# provider = IBMQ.enable_account(\"PUT TOKEN HERE\")\n",
    "backend = QasmSimulator()\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\"  #\"last_expr\" or \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-18T07:53:55.033519Z",
     "start_time": "2021-05-18T07:53:54.942084Z"
    },
    "code_folding": [
     0,
     14,
     41,
     56,
     100,
     153,
     168,
     240,
     256
    ]
   },
   "outputs": [],
   "source": [
    "def get_angles(a: int, n) -> np.ndarray:\n",
    "        # \"\"\"Calculates the array of angles to be used in the addition in Fourier Space.\"\"\"\n",
    "        s = bin(int(a))[2:].zfill(n + 1)\n",
    "        angles = np.zeros([n + 1])\n",
    "        for i in range(0, n + 1):\n",
    "            for j in range(i, n + 1):\n",
    "                if s[j] == '1':\n",
    "                    angles[n - i] += math.pow(2, -(j - i))\n",
    "            angles[n - i] *= np.pi\n",
    "        return angles[::-1]\n",
    "        # This returns the angles in the opposite order\n",
    "\n",
    "\n",
    "\n",
    "def my_create_QFT(qft_num_qubits,approximation_degree: int = 0,do_swaps: bool = False,insert_barriers: bool = True, name: str = 'qft'):\n",
    "#     up_reg = QuantumRegister(size = qft_num_qubits, name=\"aux\")\n",
    "    circuit_qft = QuantumCircuit(qft_num_qubits)\n",
    "    i=qft_num_qubits-1\n",
    "    while i>=0:\n",
    "#         circuit_qft.h(up_reg[i])\n",
    "        circuit_qft.h(i)\n",
    "        j=i-1  \n",
    "        while j>=0:\n",
    "            if (np.pi)/(pow(2,(i-j))) > approximation_degree:\n",
    "#                 circuit_qft.cu1( (np.pi)/(pow(2,(i-j))) , up_reg[i] , up_reg[j] )\n",
    "                circuit_qft.cu1( (np.pi)/(pow(2,(i-j))) , i , j )\n",
    "                j=j-1   \n",
    "        if insert_barriers:\n",
    "            circuit_qft.barrier()\n",
    "        i=i-1  \n",
    "\n",
    "    \"\"\" If specified, apply the Swaps at the end \"\"\"\n",
    "    if do_swaps:\n",
    "        i=0\n",
    "        while i < ((qft_num_qubits-1)/2):\n",
    "#             circuit_qft.swap(up_reg[i], up_reg[qft_num_qubits-1-i])\n",
    "            circuit_qft.swap(i, qft_num_qubits-1-i)\n",
    "            i=i+1\n",
    "    circuit_qft.name = \"QFT\"\n",
    "    return circuit_qft\n",
    "\n",
    "def my_create_inverse_QFT(qft_num_qubits,approximation_degree: int = 0,do_swaps: bool = False,insert_barriers: bool = True, name: str = 'iqft'):\n",
    "    my_create_QFT_circuit = my_create_QFT(qft_num_qubits,approximation_degree,do_swaps,insert_barriers, name)\n",
    "    my_create_inverse_QFT_circuit = my_create_QFT_circuit.inverse()\n",
    "    my_create_inverse_QFT_circuit.name = \"QFT†\"\n",
    "    return my_create_inverse_QFT_circuit\n",
    "\n",
    "\n",
    "def phi_add_gate(size: int, angles: Union[np.ndarray, ParameterVector]) -> Gate:\n",
    "#     \"\"\"Gate that performs addition by a in Fourier Space.\"\"\"\n",
    "    circuit = QuantumCircuit(size, name=\"phi_add\")\n",
    "    for i, angle in enumerate(angles):\n",
    "        circuit.p(angle, i)\n",
    "    return circuit.to_gate()\n",
    "    \n",
    "\n",
    "def double_controlled_phi_add_mod_N(num_qubits: int, angles: Union[np.ndarray, ParameterVector],reg_size, a, N, n) -> QuantumCircuit:\n",
    "    # \"\"\"Creates a circuit which implements double-controlled modular addition by a.\"\"\"\n",
    "    circuit = QuantumCircuit(num_qubits, name=\"ccphi_add_mod_N\")\n",
    "\n",
    "    ctl_up = 0\n",
    "    ctl_down = 1\n",
    "    ctl_aux = 2\n",
    "\n",
    "    # get qubits from aux register, omitting the control qubit\n",
    "    qubits = range(3, num_qubits)\n",
    "\n",
    "    # store the gates representing addition/subtraction by a in Fourier Space\n",
    "    phi_add_a = phi_add_gate(len(qubits), angles)\n",
    "    iphi_add_a = phi_add_a.inverse()\n",
    "    \n",
    "    phi_add_N = phi_add_gate(reg_size - 1, get_angles(N, n))\n",
    "    iphi_add_N = phi_add_N.inverse()\n",
    "\n",
    "    circuit.append(phi_add_a.control(2), [ctl_up, ctl_down, *qubits])\n",
    "    circuit.append(iphi_add_N, qubits)\n",
    "    \n",
    "    qft = QFT(n + 1).to_instruction()\n",
    "#     qft = my_create_QFT(n + 1).to_instruction()\n",
    "    \n",
    "    iqft = QFT(n + 1).inverse().to_instruction()\n",
    "#     iqft = my_create_inverse_QFT(n + 1).to_instruction()\n",
    "    circuit.append(iqft, qubits)\n",
    "\n",
    "    circuit.cx(qubits[0], ctl_aux)\n",
    "    \n",
    "    circuit.append(qft, qubits)\n",
    "    circuit.append(phi_add_N, qubits)\n",
    "    circuit.append(iphi_add_a.control(2), [ctl_up, ctl_down, *qubits])\n",
    "    circuit.append(iqft, qubits)\n",
    "\n",
    "    circuit.x(qubits[0])\n",
    "    circuit.cx(qubits[0], ctl_aux)\n",
    "    circuit.x(qubits[0])\n",
    "\n",
    "    circuit.append(qft, qubits)\n",
    "    circuit.append(phi_add_a.control(2), [ctl_up, ctl_down, *qubits])\n",
    "    return circuit\n",
    "\n",
    "# \"\"\"Circuit that implements single controlled modular multiplication by a\"\"\"\n",
    "def controlled_multiple_mod_N(num_qubits: int, N: int, a: int, n, aux_reg_size):\n",
    "#     \"\"\"Implements modular multiplication by a as an instruction.\"\"\"\n",
    "    circuit = QuantumCircuit(\n",
    "        num_qubits, \n",
    "#         name=\"multiply_by_{}_mod_{}\".format(a % N, N),\n",
    "        name=r\"${0}^{{{1}^{{{2}}}}} mod{3}$\".format(2,2,int(math.log(math.log(a,2),2)), N)\n",
    "    )\n",
    "#     label = r\"${0}^{{{1}^{{{2}}}}} mod{3}$\".format(\"†\",\"y\")\n",
    "    down = circuit.qubits[1: n + 1]\n",
    "    aux = circuit.qubits[n + 1:]\n",
    "    qubits = [aux[i] for i in reversed(range(n + 1))]\n",
    "    ctl_up = 0\n",
    "    ctl_aux = aux[-1]\n",
    "    \n",
    "    angle_params = ParameterVector(\"angles\", length=len(aux) - 1)\n",
    "    double_controlled_phi_add = double_controlled_phi_add_mod_N(\n",
    "        len(aux) + 2, angle_params, aux_reg_size, a, N, n\n",
    "    )\n",
    "    \n",
    "    idouble_controlled_phi_add = double_controlled_phi_add.inverse()\n",
    "\n",
    "    qft_circuit = QFT(n + 1).to_instruction()\n",
    "#     qft_circuit = my_create_QFT(n + 1).to_instruction()\n",
    "\n",
    "    iqft_circuit = QFT(n + 1).inverse().to_instruction()\n",
    "#     iqft_circuit = my_create_inverse_QFT(n + 1).to_instruction()\n",
    "    \n",
    "    circuit.append(qft_circuit, qubits)\n",
    "    \n",
    "    # perform controlled addition by a on the aux register in Fourier space\n",
    "    for i, ctl_down in enumerate(down):\n",
    "        a_exp = (2 ** i) * a % N\n",
    "        angles = get_angles(a_exp, n)\n",
    "        bound = double_controlled_phi_add.assign_parameters({angle_params: angles})\n",
    "        circuit.append(bound, [ctl_up, ctl_down, ctl_aux, *qubits])\n",
    "\n",
    "    circuit.append(iqft_circuit, qubits)\n",
    "    \n",
    "    # perform controlled subtraction by a in Fourier space on both the aux and down register\n",
    "    for j in range(n):\n",
    "        circuit.cswap(ctl_up, down[j], aux[j])\n",
    "    circuit.append(qft_circuit, qubits)\n",
    "    \n",
    "    a_inv = modinv(a, N)\n",
    "    for i in reversed(range(len(down))):\n",
    "        a_exp = (2 ** i) * a_inv % N\n",
    "        angles = get_angles(a_exp, n)\n",
    "        bound = idouble_controlled_phi_add.assign_parameters({angle_params: angles})\n",
    "        circuit.append(bound, [ctl_up, down[i], ctl_aux, *qubits])\n",
    "\n",
    "    circuit.append(iqft_circuit, qubits)\n",
    "    \n",
    "    return circuit\n",
    "\n",
    "def modinv(a: int, m: int) -> int:\n",
    "    # \"\"\"Returns the modular multiplicative inverse of a with respect to the modulus m.\"\"\"\n",
    "    def egcd(a: int, b: int) -> Tuple[int, int, int]:\n",
    "        if a == 0:\n",
    "            return b, 0, 1\n",
    "        else:\n",
    "            g, y, x = egcd(b % a, a)\n",
    "            return g, x - (b // a) * y, y\n",
    "\n",
    "    g, x, _ = egcd(a, m)\n",
    "    if g != 1:\n",
    "        raise ValueError(\"The greatest common divisor of {} and {} is {}, so the \"\n",
    "                         \"modular inverse does not exist.\".format(a, m, g))\n",
    "    return x % m\n",
    "\n",
    "def get_factors(N: int, a: int, measurement: str) -> Optional[List[int]]:\n",
    "    # \"\"\"Apply the continued fractions to find r and the gcd to find the desired factors.\"\"\"\n",
    "    x_final = int(measurement, 2)\n",
    "    #print('In decimal, x_final value for this result is: {}.'.format(x_final))\n",
    "\n",
    "    if x_final <= 0:\n",
    "        fail_reason = 'x_final value is <= 0, there are no continued fractions.'\n",
    "    else:\n",
    "        fail_reason = None\n",
    "        #print('Running continued fractions for this case.')\n",
    "\n",
    "    # Calculate T and x/T\n",
    "    T_upper = len(measurement)\n",
    "    T = pow(2, T_upper)\n",
    "    x_over_T = x_final / T  ## this is our theta\n",
    "\n",
    "    # Cycle in which each iteration corresponds to putting one more term in the\n",
    "    # calculation of the Continued Fraction (CF) of x/T\n",
    "\n",
    "    # Initialize the first values according to CF rule\n",
    "    i = 0\n",
    "    b = array.array('i')\n",
    "    t = array.array('f')\n",
    "\n",
    "    b.append(math.floor(x_over_T))\n",
    "    t.append(x_over_T - b[i])\n",
    "\n",
    "    exponential = 0.0\n",
    "    while i < N and fail_reason is None:\n",
    "        # From the 2nd iteration onwards, calculate the new terms of the CF based on the previous terms as the rule suggests\n",
    "        if i > 0:\n",
    "            try:\n",
    "                b_temp = math.floor(1 / t[i - 1])\n",
    "            except ZeroDivisionError as err:\n",
    "                b_temp = 0\n",
    "            b.append(b_temp)\n",
    "            try:\n",
    "                t_temp = (1 / t[i - 1]) - b[i]\n",
    "            except ZeroDivisionError as err:\n",
    "                t_temp = 0\n",
    "            t.append(t_temp)  # type: ignore\n",
    "\n",
    "        # Calculate the denominator of the CF using the known terms\n",
    "        denominator = calculate_continued_fraction(b)\n",
    "\n",
    "        # Increment i for next iteration\n",
    "        i += 1\n",
    "\n",
    "        if denominator % 2 == 1:\n",
    "            #print('Odd denominator, will try next iteration of continued fractions.')\n",
    "            continue\n",
    "\n",
    "        # Denominator is even, try to get factors of N. Get the exponential a^(r/2)\n",
    "\n",
    "        if denominator < 1000:\n",
    "            try:\n",
    "                exponential = pow(a, denominator / 2)\n",
    "            except OverflowError as err:\n",
    "                exponential = 999999999\n",
    "\n",
    "        # Check if the value is too big or not\n",
    "        if exponential > 1000000:\n",
    "            if exponential == 999999999:\n",
    "                fail_reason = 'OverflowError'\n",
    "            else:\n",
    "                fail_reason = 'denominator of continued fraction is too big (> 10^9).'\n",
    "        else:\n",
    "            # The value is not too big, get the right values and do the proper gcd()\n",
    "            putting_plus = int(exponential + 1)\n",
    "            putting_minus = int(exponential - 1)\n",
    "            one_factor = math.gcd(putting_plus, N)\n",
    "            other_factor = math.gcd(putting_minus, N)\n",
    "\n",
    "            # Check if the factors found are trivial factors or are the desired factors\n",
    "            if any(factor in {1, N} for factor in (one_factor, other_factor)):\n",
    "                #print('Found just trivial factors, not good enough.')\n",
    "                # Check if the number has already been found, (use i - 1 because i was already incremented)\n",
    "                if t[i - 1] == 0:\n",
    "                    fail_reason = 'the continued fractions found exactly x_final/(2^(2n)).'\n",
    "            else:\n",
    "                return sorted((one_factor, other_factor))\n",
    "\n",
    "    # Search for factors failed, write the reason for failure to the debug logs\n",
    "    #print('Cannot find factors from measurement {0} because {1}'.format(measurement, fail_reason or 'it took too many attempts.'))\n",
    "    return None\n",
    "\n",
    "def calculate_continued_fraction(b: array.array) -> int:\n",
    "    # \"\"\"Calculate the continued fraction of x/T from the current terms of expansion b.\"\"\"\n",
    "\n",
    "    x_over_T = 0\n",
    "\n",
    "    for i in reversed(range(len(b) - 1)):\n",
    "        x_over_T = 1 / (b[i + 1] + x_over_T)\n",
    "\n",
    "    x_over_T += b[0]\n",
    "\n",
    "    frac = fractions.Fraction(x_over_T).limit_denominator()\n",
    "\n",
    "    #print('Approximation number %s of continued fractions:'.format(len(b)))\n",
    "    #print(\"Numerator:{0} \\t\\t Denominator: {1}.\".format(frac.numerator, frac.denominator))\n",
    "    return frac.denominator\n",
    "\n",
    "def process_results(sim_result, circuit, shots, N, a, n):\n",
    "    counts_result = sim_result.get_counts(circuit)\n",
    "    total_counts = len(counts_result)\n",
    "    counts_result_sorted = sorted(counts_result.items(), key=lambda x: x[1], reverse=True)\n",
    "\n",
    "\n",
    "    # \"\"\" Print info to user from the simulation results \"\"\"\n",
    "#     print('Printing the various results followed by how many times they happened (out of the {} cases):\\n'.format(shots))\n",
    "    \n",
    "    counts_result_keys = list(counts_result.keys())\n",
    "    counts_result_values = list(counts_result.values())\n",
    "    \n",
    "    #i=0\n",
    "    #while i < len(counts_result):\n",
    "        #print('Result \\\"{0}\\\" happened {1} times out of {2}\\n'.format(list(sim_result.get_counts().keys())[i],list(sim_result.get_counts().values())[i],shots))\n",
    "        #print('Result \\\"{0}\\\" happened {1} times out of {2}\\n'.format(counts_result_keys[i],counts_result_values[i],shots))\n",
    "        #i=i+1\n",
    "\n",
    "    prob_success=0\n",
    "    prob_failure=0\n",
    "    result_successful_counts = 0\n",
    "    result_failure_counts = 0\n",
    "#     len(counts_result_sorted)\n",
    "    \n",
    "    # For each simulation result, print proper info to user and try to calculate the factors of N\n",
    "    #for measurement in counts_result_keys:\n",
    "    for measurement, frequency in counts_result_sorted:\n",
    "        # Get the x_final value from the final state qubits\n",
    "        x_value = int(measurement, 2)\n",
    "        #prob_this_result = 100 * ( int(counts_result[measurement] ) ) / (shots)\n",
    "        prob_this_result = 100 * frequency/shots\n",
    "#         print(\"------> Analyzing result {0}. This result happened in {1:.4f} % of all cases\\nIn decimal, x_final value for this result is: {2}\".format(measurement,prob_this_result,x_value))\n",
    "        factors = get_factors(N, a, measurement) \n",
    "\n",
    "        if factors:\n",
    "            prob_success = prob_success + prob_this_result\n",
    "#             print('Found factors {0} from measurement {1} which is {2} in decimal.\\n'.format(factors, measurement, x_value))\n",
    "            result_successful_counts = result_successful_counts + 1\n",
    "\n",
    "            if factors not in result_factors:\n",
    "                result_factors.append(factors)\n",
    "        elif not factors:\n",
    "            prob_failure = prob_failure + prob_this_result\n",
    "            result_failure_counts = result_failure_counts + 1\n",
    "            \n",
    "        \n",
    "    return [result_factors, prob_success, prob_failure, total_counts, result_successful_counts,result_failure_counts]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-18T07:54:09.656997Z",
     "start_time": "2021-05-18T07:54:09.635367Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_shor(a,N,shots):\n",
    "    start_time_number = datetime.now()\n",
    "    start_time = start_time_number.strftime(\"%H:%M:%S\")\n",
    "    summary_result = dict()\n",
    "    validate_min('N', N, 3)\n",
    "    validate_min('a', a, 2)\n",
    "\n",
    "\n",
    "    if N < 1 or N % 2 == 0:\n",
    "        raise ValueError('The input needs to be an odd integer greater than 1.')\n",
    "\n",
    "    if a >= N or math.gcd(a, N) != 1:\n",
    "        raise ValueError('The integer a needs to satisfy a < N and gcd(a, N) = 1.')\n",
    "\n",
    "    n = math.ceil(math.log(N,2))\n",
    "    global result_factors\n",
    "    result_factors = []\n",
    "    tf, b, p = is_power(N, return_decomposition=True)\n",
    "    if tf:\n",
    "        print('The input integer is a power: {0}={1}^{2}.'.format(N, b, p))\n",
    "        result_factors.append(b)\n",
    "\n",
    "    # \"\"\"auxilliary quantum register used in addition and multiplication\"\"\"\n",
    "    aux_reg = QuantumRegister(size = n+2, name=\"aux_reg\")\n",
    "\n",
    "    up_reg = QuantumRegister(2*n, name = \"up_reg\")\n",
    "\n",
    "    # \"\"\"quantum register where the multiplications are made\"\"\"\n",
    "    down_reg = QuantumRegister(n, name = \"down_reg\")\n",
    "\n",
    "    # \"\"\"classical register where the measured values of the QFT are stored\"\"\"\n",
    "    up_classic = ClassicalRegister(2*n, name=\"up_classic\")\n",
    "\n",
    "    # \"\"\" Create Quantum Circuit \"\"\"\n",
    "    circuit = QuantumCircuit(up_reg ,down_reg ,aux_reg, up_classic, name=\"Shor circuit(N={}, a={})\".format(N, a))\n",
    "\n",
    "    # phi_add_N_gate = phiADD(circuit,q,a,N,inv)\n",
    "    phi_add_N_gate = phi_add_gate(aux_reg.size - 1, get_angles(N,n))\n",
    "    iphi_add_N_gate = phi_add_N_gate.inverse()\n",
    "\n",
    "    # \"\"\" Initialize down register to 1 and create maximal superposition in top register \"\"\"\n",
    "    circuit.h(up_reg)\n",
    "    circuit.x(down_reg[0])\n",
    "    # circuit.draw(filename = \"shor_standard_QFT\")\n",
    "\n",
    "    # \"\"\" Apply the multiplication gates as showed in the report in order to create the exponentiation \"\"\"\n",
    "    for i, ctl_up in enumerate(up_reg):  # type: ignore\n",
    "        a_aux = int(pow(a, pow(2, i)))\n",
    "        controlled_multiple_mod_N_circuit = controlled_multiple_mod_N(\n",
    "            len(down_reg) + len(aux_reg) + 1, N, a_aux,n,aux_reg.size\n",
    "        )\n",
    "        controlled_multiple_mod_N_result = controlled_multiple_mod_N_circuit.to_instruction()\n",
    "        circuit.append(\n",
    "            controlled_multiple_mod_N_result, [ctl_up, *down_reg, *aux_reg]\n",
    "        )\n",
    "#         circuit.draw()\n",
    "\n",
    "    iqft = QFT(len(up_reg)).inverse().to_instruction()\n",
    "    # iqft = my_create_inverse_QFT(len(up_reg)).to_instruction()\n",
    "    # iqft = my_create_inverse_QFT(len(up_reg), insert_barriers = False).to_gate(label = r\"$QFT^{{{0}}}$\".format(\"†\"))\n",
    "\n",
    "    circuit.append(iqft, up_reg)\n",
    "    circuit.measure(up_reg, up_classic)\n",
    "\n",
    "#     circuit.draw(filename = \"shor_standard_QFT_final_circuit\",fold = -1 )\n",
    "\n",
    "    #     print(summarize_circuits(circuit))\n",
    "    # circuit.draw()\n",
    "\n",
    "    print('Running with N={0} and a={1} with number of qubits n={2}'.format(N, a, 4*n + 2))\n",
    "    \n",
    "    qc_compiled = transpile(circuit, backend, optimization_level = 3)\n",
    "    job_sim_1 = backend.run(qc_compiled, shots=shots)\n",
    "    sim_result=job_sim_1.result()\n",
    "\n",
    "#     counts_result = sim_result.get_counts(circuit)\n",
    "#     len(counts_result)\n",
    "#     measurement_plot = qiskit.visualization.plot_histogram(counts_result,figsize=(20, 12) ,number_to_keep = 30,bar_labels=True, title = \"Measurement results from shor_standard_QFT circuit variant\" )\n",
    "\n",
    "    # measurement_plot.savefig(\"shor_standard_QFT_measurement\")\n",
    "#     measurement_plot\n",
    "\n",
    "    processed_result = process_results(sim_result, circuit, shots, N, a, n)\n",
    "    end_time_number = datetime.now()\n",
    "    end_time = end_time_number.strftime(\"%H:%M:%S\")\n",
    "    duration = end_time_number - start_time_number\n",
    "\n",
    "\n",
    "    print(\"Current Start Time =\", start_time)\n",
    "    print(processed_result)\n",
    "    print(\"Current End Time =\", end_time)\n",
    "\n",
    "    circuit_count_ops = circuit.count_ops()\n",
    "    circuit_decomposed = circuit.decompose()\n",
    "    circuit_decomposed_count_ops = circuit_decomposed.count_ops()\n",
    "    qc_compiled_count_ops = qc_compiled.count_ops()\n",
    "\n",
    "    summary_result[\"num_qubits\"] = n\n",
    "    summary_result[\"Number(N)\"] = N\n",
    "    summary_result[\"a\"] = a\n",
    "    summary_result[\"start_time\"] = start_time\n",
    "    summary_result[\"end_time\"] = end_time\n",
    "    summary_result[\"duration\"] = duration\n",
    "    summary_result[\"result_factors\"] = processed_result[0]\n",
    "    summary_result[\"prob_success\"] = processed_result[1]\n",
    "    summary_result[\"prob_failure\"] = processed_result[2]\n",
    "    summary_result[\"total_counts\"] = processed_result[3]\n",
    "    summary_result[\"result_successful_counts\"] = processed_result[4]\n",
    "    summary_result[\"result_failure_counts\"] = processed_result[5]\n",
    "\n",
    "\n",
    "    summary_result[\"circuit_width\"] = circuit.width()\n",
    "    summary_result[\"circuit_depth\"] = circuit.depth()\n",
    "    summary_result[\"circuit_size\"] = circuit.size() \n",
    "    summary_result[\"circuit_num_nonlocal_gates\"] = circuit.num_nonlocal_gates()\n",
    "    summary_result[\"circuit_num_ancillas\"] = circuit.num_ancillas\n",
    "    summary_result[\"circuit_num_clbits\"] = circuit.num_clbits\n",
    "    summary_result[\"circuit_num_qubits\"] = circuit.num_qubits\n",
    "    summary_result[\"circuit_num_ancillas\"] = circuit.num_ancillas\n",
    "    summary_result[\"circuit_num_of_count_ops\"] = len(circuit_count_ops)\n",
    "    summary_result[\"circuit_num_of_x\"] = circuit_count_ops.get('x')\n",
    "    summary_result[\"circuit_num_of_measure\"] = circuit_count_ops.get('measure')\n",
    "    summary_result[\"circuit_num_of_h\"] = circuit_count_ops.get('h')\n",
    "    summary_result[\"circuit_num_of_cswap\"] = circuit_count_ops.get('cswap')\n",
    "    summary_result[\"circuit_num_of_swap\"] = circuit_count_ops.get('swap')\n",
    "    summary_result[\"circuit_num_of_cx\"] = circuit_count_ops.get('cx')\n",
    "    summary_result[\"circuit_num_of_toffoli\"] = circuit_count_ops.get('toffoli')\n",
    "    summary_result[\"circuit_num_of_p\"] = circuit_count_ops.get('p')\n",
    "    summary_result[\"circuit_num_of_t\"] = circuit_count_ops.get('t')\n",
    "\n",
    "    summary_result[\"circuit_decomposed_width\"] = circuit_decomposed.width()\n",
    "    summary_result[\"circuit_decomposed_depth\"] = circuit_decomposed.depth()\n",
    "    summary_result[\"circuit_decomposed_size\"] = circuit_decomposed.size() \n",
    "    summary_result[\"circuit_decomposed_num_nonlocal_gates\"] = circuit_decomposed.num_nonlocal_gates()\n",
    "    summary_result[\"circuit_decomposed_num_ancillas\"] = circuit_decomposed.num_ancillas\n",
    "    summary_result[\"circuit_decomposed_num_clbits\"] = circuit_decomposed.num_clbits\n",
    "    summary_result[\"circuit_decomposed_num_qubits\"] = circuit_decomposed.num_qubits\n",
    "    summary_result[\"circuit_decomposed_num_ancillas\"] = circuit_decomposed.num_ancillas\n",
    "    summary_result[\"circuit_decomposed_num_of_count_ops\"] = len(circuit_decomposed_count_ops)\n",
    "    summary_result[\"circuit_decomposed_num_of_x\"] = circuit_decomposed_count_ops.get('x')\n",
    "    summary_result[\"circuit_decomposed_num_of_measure\"] = circuit_decomposed_count_ops.get('measure')\n",
    "    summary_result[\"circuit_decomposed_num_of_h\"] = circuit_decomposed_count_ops.get('h')\n",
    "    summary_result[\"circuit_decomposed_num_of_cswap\"] = circuit_decomposed_count_ops.get('cswap')\n",
    "    summary_result[\"circuit_decomposed_num_of_swap\"] = circuit_decomposed_count_ops.get('swap')\n",
    "    summary_result[\"circuit_decomposed_num_of_cx\"] = circuit_decomposed_count_ops.get('cx')\n",
    "    summary_result[\"circuit_decomposed_num_of_toffoli\"] = circuit_decomposed_count_ops.get('toffoli')\n",
    "    summary_result[\"circuit_decomposed_num_of_p\"] = circuit_decomposed_count_ops.get('p')\n",
    "    summary_result[\"circuit_decomposed_num_of_t\"] = circuit_decomposed_count_ops.get('t')\n",
    "\n",
    "    summary_result[\"qc_compiled_width\"] = qc_compiled.width()\n",
    "    summary_result[\"qc_compiled_depth\"] = qc_compiled.depth()\n",
    "    summary_result[\"qc_compiled_size\"] = qc_compiled.size() \n",
    "    summary_result[\"qc_compiled_num_nonlocal_gates\"] = qc_compiled.num_nonlocal_gates()\n",
    "    summary_result[\"qc_compiled_num_ancillas\"] = qc_compiled.num_ancillas\n",
    "    summary_result[\"qc_compiled_num_clbits\"] = qc_compiled.num_clbits\n",
    "    summary_result[\"qc_compiled_num_qubits\"] = qc_compiled.num_qubits\n",
    "    summary_result[\"qc_compiled_num_ancillas\"] = qc_compiled.num_ancillas\n",
    "    summary_result[\"qc_compiled_num_of_count_ops\"] = len(qc_compiled_count_ops)\n",
    "    summary_result[\"qc_compiled_num_of_x\"] = qc_compiled_count_ops.get('x')\n",
    "    summary_result[\"qc_compiled_num_of_measure\"] = qc_compiled_count_ops.get('measure')\n",
    "    summary_result[\"qc_compiled_num_of_h\"] = qc_compiled_count_ops.get('h')\n",
    "    summary_result[\"qc_compiled_num_of_cswap\"] = qc_compiled_count_ops.get('cswap')\n",
    "    summary_result[\"qc_compiled_num_of_swap\"] = qc_compiled_count_ops.get('swap')\n",
    "    summary_result[\"qc_compiled_num_of_cx\"] = qc_compiled_count_ops.get('cx')\n",
    "    summary_result[\"qc_compiled_num_of_toffoli\"] = qc_compiled_count_ops.get('toffoli')\n",
    "    summary_result[\"qc_compiled_num_of_p\"] = qc_compiled_count_ops.get('p')\n",
    "    summary_result[\"qc_compiled_num_of_t\"] = qc_compiled_count_ops.get('t')\n",
    "\n",
    "    return summary_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Process the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-18T07:11:58.167920Z",
     "start_time": "2021-05-18T07:11:57.932902Z"
    }
   },
   "outputs": [],
   "source": [
    "# Run for just a single number N\n",
    "%%time\n",
    "N = 33\n",
    "shots = 1024\n",
    "global result_factors\n",
    "\n",
    "all_summary_result_temp = []\n",
    "for random_a in range(16, 17):\n",
    "    if math.gcd(random_a,N) > 1:\n",
    "        continue\n",
    "    a = random_a\n",
    "    summary_result = my_shor(a,N,shots)\n",
    "    print(\"Finished running for a = {} and N = {}\\n\".format(a, N))\n",
    "    all_summary_result_temp.append(summary_result)\n",
    "    summary_result_list = []\n",
    "    for key, value in summary_result.items():\n",
    "        summary_result_list.append([key,value])\n",
    "    summary_result_list\n",
    "    \n",
    "    with open(\"a({0})_N({1})_standard.csv\".format(a, N), 'a') as myfile:\n",
    "        write = csv.writer(myfile)\n",
    "        #write.writerow(fields)\n",
    "        write.writerows(summary_result_list)\n",
    "\n",
    "all_summary_result_temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run for many numbers N.\n",
    "%%time\n",
    "shots = 1024\n",
    "global result_factors\n",
    "all_summary_result = []\n",
    "for N in [15, 21, 33, 35, 39, 51, 55, 57]:\n",
    "    for a in range(2, N):\n",
    "        if math.gcd(a,N) > 1:\n",
    "            continue\n",
    "        print(\"Beginning running for a = {} and N = {}\".format(a, N))\n",
    "        summary_result = my_shor(a,N,shots)\n",
    "        print(\"Finished running for a = {} and N = {}\\n\\n\".format(a, N))\n",
    "        all_summary_result.append(summary_result)\n",
    "\n",
    "all_summary_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%qiskit_copyright"
   ]
  }
 ],
 "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.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}