{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b01470a",
   "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": "code",
   "execution_count": null,
   "id": "76cd6637",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-19T10:16:06.644552Z",
     "start_time": "2021-05-19T10:15:51.905687Z"
    }
   },
   "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",
    "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",
   "id": "02bd623e",
   "metadata": {},
   "source": [
    "### Define functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db9328a2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-19T10:16:06.753809Z",
     "start_time": "2021-05-19T10:16:06.692206Z"
    },
    "code_folding": [
     1,
     30,
     82,
     193,
     199
    ]
   },
   "outputs": [],
   "source": [
    "# \"\"\" Function to check if N is of type q^p\"\"\"\n",
    "def check_if_power(N):\n",
    "#     \"\"\" Check if N is a perfect power in O(n^3) time, n=ceil(logN) \"\"\"\n",
    "    b=2\n",
    "    while (2**b) <= N:\n",
    "        a = 1\n",
    "        c = N\n",
    "        while (c-a) >= 2:\n",
    "            m = int( (a+c)/2 )\n",
    "\n",
    "            if (m**b) < (N+1):\n",
    "                p = int( (m**b) )\n",
    "            else:\n",
    "                p = int(N+1)\n",
    "\n",
    "            if int(p) == int(N):\n",
    "                print('N is {0}^{1}'.format(int(m),int(b)) )\n",
    "                return True\n",
    "\n",
    "            if p<N:\n",
    "                a = int(m)\n",
    "            else:\n",
    "                c = int(m)\n",
    "        b=b+1\n",
    "\n",
    "    return False\n",
    "\n",
    "def egcd(a, b):\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",
    "def modinv(a, m):\n",
    "    g, x, y = egcd(a, m)\n",
    "    if g != 1:\n",
    "        raise Exception('modular inverse does not exist')\n",
    "    else:\n",
    "        return x % m\n",
    "\n",
    "def create_QFT(circuit,up_reg,n,with_swaps):\n",
    "    i=n-1\n",
    "    while i>=0:\n",
    "        circuit.h(up_reg[i])        \n",
    "        j=i-1  \n",
    "        while j>=0:\n",
    "            if (np.pi)/(pow(2,(i-j))) > 0:\n",
    "                circuit.cu1( (np.pi)/(pow(2,(i-j))) , up_reg[i] , up_reg[j] )\n",
    "                j=j-1   \n",
    "        i=i-1  \n",
    "\n",
    "    if with_swaps==1:\n",
    "        i=0\n",
    "        while i < ((n-1)/2):\n",
    "            circuit.swap(up_reg[i], up_reg[n-1-i])\n",
    "            i=i+1\n",
    "\n",
    "def create_inverse_QFT(circuit,up_reg,n,with_swaps):\n",
    "    if with_swaps==1:\n",
    "        i=0\n",
    "        while i < ((n-1)/2):\n",
    "            circuit.swap(up_reg[i], up_reg[n-1-i])\n",
    "            i=i+1\n",
    "    i=0\n",
    "    while i<n:\n",
    "        circuit.h(up_reg[i])\n",
    "        if i != n-1:\n",
    "            j=i+1\n",
    "            y=i\n",
    "            while y>=0:\n",
    "                 if (np.pi)/(pow(2,(j-y))) > 0:\n",
    "                    circuit.cu1( - (np.pi)/(pow(2,(j-y))) , up_reg[j] , up_reg[y] )\n",
    "                    y=y-1   \n",
    "        i=i+1\n",
    "\n",
    "def getAngle(a, N):\n",
    "    s=bin(int(a))[2:].zfill(N) \n",
    "    angle = 0\n",
    "    for i in range(0, N):\n",
    "        if s[N-1-i] == '1': \n",
    "            angle += math.pow(2, -(N-i))\n",
    "    angle *= np.pi\n",
    "    return angle\n",
    "\n",
    "def getAngles(a,N):\n",
    "    s=bin(int(a))[2:].zfill(N) \n",
    "    angles=np.zeros([N])\n",
    "    for i in range(0, N):\n",
    "        for j in range(i,N):\n",
    "            if s[j]=='1':\n",
    "                angles[N-i-1]+=math.pow(2, -(j-i))\n",
    "        angles[N-i-1]*=np.pi\n",
    "    return angles\n",
    "\n",
    "def ccphase(circuit, angle, ctl1, ctl2, tgt):\n",
    "    circuit.cu1(angle/2,ctl1,tgt)\n",
    "    circuit.cx(ctl2,ctl1)\n",
    "    circuit.cu1(-angle/2,ctl1,tgt)\n",
    "    circuit.cx(ctl2,ctl1)\n",
    "    circuit.cu1(angle/2,ctl2,tgt)\n",
    "\n",
    "def phiADD(circuit, q, a, N, inv):\n",
    "    angle=getAngles(a,N)\n",
    "    for i in range(0,N):\n",
    "        if inv==0:\n",
    "            circuit.u1(angle[i],q[i])\n",
    "        else:\n",
    "            circuit.u1(-angle[i],q[i])\n",
    "\n",
    "def cphiADD(circuit, q, ctl, a, n, inv):\n",
    "    angle=getAngles(a,n)\n",
    "    for i in range(0,n):\n",
    "        if inv==0:\n",
    "            circuit.cu1(angle[i],ctl,q[i])\n",
    "        else:\n",
    "            circuit.cu1(-angle[i],ctl,q[i])\n",
    "        \n",
    "def ccphiADD(circuit,q,ctl1,ctl2,a,n,inv):\n",
    "    angle=getAngles(a,n)\n",
    "    for i in range(0,n):\n",
    "        if inv==0:\n",
    "            ccphase(circuit,angle[i],ctl1,ctl2,q[i])\n",
    "        else:\n",
    "            ccphase(circuit,-angle[i],ctl1,ctl2,q[i])\n",
    "        \n",
    "def ccphiADDmodN(circuit, q, ctl1, ctl2, aux, a, N, n):\n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 0)\n",
    "    phiADD(circuit, q, N, n, 1)\n",
    "#     phiADD(circuit, q, a,N, 1)\n",
    "    create_inverse_QFT(circuit, q, n, 0)\n",
    "    circuit.cx(q[n-1],aux)\n",
    "    create_QFT(circuit,q,n,0)\n",
    "    cphiADD(circuit, q, aux, N, n, 0)\n",
    "#     cphiADD(circuit, q, aux, a, n, 0)\n",
    "    \n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 1)\n",
    "    create_inverse_QFT(circuit, q, n, 0)\n",
    "    circuit.x(q[n-1])\n",
    "    circuit.cx(q[n-1], aux)\n",
    "    circuit.x(q[n-1])\n",
    "    create_QFT(circuit,q,n,0)\n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 0)\n",
    "\n",
    "def ccphiADDmodN_inv(circuit, q, ctl1, ctl2, aux, a, N, n):\n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 1)\n",
    "    create_inverse_QFT(circuit, q, n, 0)\n",
    "    circuit.x(q[n-1])\n",
    "    circuit.cx(q[n-1],aux)\n",
    "    circuit.x(q[n-1])\n",
    "    create_QFT(circuit, q, n, 0)\n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 0)\n",
    "    cphiADD(circuit, q, aux, N, n, 1)\n",
    "#     cphiADD(circuit, q, aux, a, n, 1)\n",
    "    create_inverse_QFT(circuit, q, n, 0)\n",
    "    circuit.cx(q[n-1], aux)\n",
    "    create_QFT(circuit, q, n, 0)\n",
    "    phiADD(circuit, q, N, n, 0)\n",
    "#     phiADD(circuit, q, a, N, 0)\n",
    "    ccphiADD(circuit, q, ctl1, ctl2, a, n, 1)\n",
    "\n",
    "def cMULTmodN(circuit, ctl, q, aux, a, N, n):\n",
    "    \n",
    "\n",
    "#     up_reg = QuantumRegister(1, name = \"up_reg\")\n",
    "#     down_reg = QuantumRegister(n, name = \"down_reg\")\n",
    "#     up_classic = ClassicalRegister(2*n, name=\"up_classic\")\n",
    "#     c_aux = ClassicalRegister(1, name = \"aux_classic\")\n",
    "#     cMULTmodN_circuit = QuantumCircuit(\n",
    "#         up_reg ,down_reg , aux,up_classic, c_aux, \n",
    "#         name=r\"${0}^{{{1}^{{{2}}}}} mod{3}$\".format(2,2,int(math.log(math.log(a,2),2)), N)\n",
    "#     )\n",
    "# create_QFT(cMULTmodN_circuit,aux,n+1,0)\n",
    "#     for i in range(0, n):\n",
    "#         ccphiADDmodN(cMULTmodN_circuit, aux, q[i], ctl, aux[n+1], (2**i)*a % N, N, n+1)\n",
    "# create_inverse_QFT(cMULTmodN_circuit, aux, n+1, 0)\n",
    "#     for i in range(0, n):\n",
    "#         circuit.cswap(ctl,q[i],aux[i])\n",
    "#         cMULTmodN_circuit.cswap(ctl,q[i],aux[i])\n",
    "#     create_QFT(cMULTmodN_circuit, aux, n+1, 0)\n",
    "#         ccphiADDmodN_inv(cMULTmodN_circuit, aux, q[i], ctl, aux[n+1], math.pow(2,i)*a_inv % N, N, n+1)\n",
    "#     create_inverse_QFT(cMULTmodN_circuit, aux, n+1, 0)\n",
    "#     cMULTmodN_circuit_instruction = cMULTmodN_circuit.to_instruction()\n",
    "#     circuit.append(cMULTmodN_circuit_instruction, [ctl, *down_reg, *aux])\n",
    "        \n",
    "    create_QFT(circuit,aux,n+1,0)\n",
    "    for i in range(0, n):\n",
    "        ccphiADDmodN(circuit, aux, q[i], ctl, aux[n+1], (2**i)*a % N, N, n+1)\n",
    "        \n",
    "    create_inverse_QFT(circuit, aux, n+1, 0)\n",
    "    \n",
    "    for i in range(0, n):\n",
    "        circuit.cswap(ctl,q[i],aux[i])\n",
    "\n",
    "    a_inv = modinv(a, N)\n",
    "    create_QFT(circuit, aux, n+1, 0)\n",
    "\n",
    "    \n",
    "    i = n-1\n",
    "    while i >= 0:\n",
    "        ccphiADDmodN_inv(circuit, aux, q[i], ctl, aux[n+1], math.pow(2,i)*a_inv % N, N, n+1)\n",
    "        i -= 1\n",
    "    create_inverse_QFT(circuit, aux, n+1, 0)\n",
    "    \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",
    "    return frac.denominator\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^3).'\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",
    "    return None\n",
    "\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",
    "    counts_result_keys = list(counts_result.keys())\n",
    "    counts_result_values = list(counts_result.values())\n",
    "\n",
    "    prob_success=0\n",
    "    prob_failure=0\n",
    "    result_successful_counts = 0\n",
    "    result_failure_counts = 0\n",
    "    \n",
    "\n",
    "    for initial_undesired_measurement, frequency in counts_result_sorted:\n",
    "        measurement = initial_undesired_measurement.split(\" \")[1]\n",
    "        x_value = int(measurement, 2)\n",
    "        prob_this_result = 100 * frequency/shots\n",
    "        factors = get_factors(N, a, measurement) \n",
    "\n",
    "        if factors:\n",
    "            prob_success = prob_success + prob_this_result\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",
    "    return [result_factors, prob_success, prob_failure, total_counts, result_successful_counts,result_failure_counts]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc3a086e",
   "metadata": {},
   "source": [
    "### Initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8668d24",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-19T10:16:16.739385Z",
     "start_time": "2021-05-19T10:16:16.733334Z"
    },
    "code_folding": [
     0
    ]
   },
   "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",
    "    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",
    "\n",
    "\n",
    "    # \"\"\"auxilliary quantum register used in addition and multiplication\"\"\"\n",
    "    aux = QuantumRegister(size = n+2, name=\"aux_reg\")\n",
    "\n",
    "    # \"\"\"single qubit where the sequential QFT is performed\"\"\"\n",
    "    up_reg = QuantumRegister(1, name = \"up_reg\")\n",
    "\n",
    "    down_reg = QuantumRegister(n, name = \"down_reg\")\n",
    "\n",
    "    # \"\"\"classical register where the measured values of the sequential QFT are stored\"\"\"\n",
    "    up_classic = ClassicalRegister(2*n, name=\"up_classic\")\n",
    "\n",
    "    # \"\"\"classical bit used to reset the state of the top qubit to 0 if the previous measurement was 1\"\"\"\n",
    "    c_aux = ClassicalRegister(1, name = \"aux_classic\")\n",
    "\n",
    "    # \"\"\" Create Quantum Circuit \"\"\"\n",
    "    circuit = QuantumCircuit(up_reg ,down_reg , aux,up_classic, c_aux)\n",
    "\n",
    "    circuit.x(down_reg[0])\n",
    "#     circuit.draw(filename = \"shor_semiclassical_QFT_initialization\")\n",
    "\n",
    "\n",
    "    for i in range(0, 2*n):\n",
    "        circuit.x(up_reg).c_if(c_aux, 1)\n",
    "        circuit.h(up_reg)\n",
    "        cMULTmodN(circuit, up_reg[0], down_reg, aux, a**(2**(2*n-1-i)), N, n)\n",
    "    #     later confirm if this should be up_reg[i] instead of up_reg[0]\n",
    "\n",
    "        for j in range(0, 2**i):\n",
    "            circuit.u1(getAngle(j, i), up_reg[0]).c_if(up_classic, j)\n",
    "            circuit.h(up_reg)\n",
    "            circuit.measure(up_reg[0], up_classic[i])\n",
    "            circuit.measure(up_reg[0], c_aux[0])\n",
    "\n",
    "    # circuit.draw(filename = \"shor_semiclassical_QFT_final_circuit\")\n",
    "    circuit.draw()\n",
    "\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",
    "\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_semiclassical_QFT_measurement_result\")\n",
    "# measurement_plot\n",
    "\n",
    "    processed_result = process_results(sim_result, circuit, shots, N, a, n)\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": "code",
   "execution_count": null,
   "id": "3bbfd6d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run for just a single number N\n",
    "%%time\n",
    "N = 21\n",
    "shots = 1024\n",
    "global result_factors\n",
    "\n",
    "all_summary_result_temp = []\n",
    "for random_a in range(2, N):\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})_semiclassical.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,
   "id": "b1256231",
   "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\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff52fdc",
   "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": 5
}
