{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:59:57.837259Z",
     "start_time": "2021-05-20T16:59:57.815321Z"
    }
   },
   "outputs": [],
   "source": [
    "# This code has been adapted and modified from IBM Qiskit 2021 and \n",
    "# uses the constant optimized modular exponentiation circuit for mod 15 as contained\n",
    "# in https://arxiv.org/abs/1202.6614.\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "from decimal import *\n",
    "import matplotlib.pyplot as plt\n",
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qiskit import Aer, transpile, assemble\n",
    "import pandas as pd\n",
    "from fractions import Fraction\n",
    "# \n",
    "# import math\n",
    "# from math import gcd\n",
    "# from numpy.random import randint\n",
    "# \n",
    "# \n",
    "# from decimal import *\n",
    "print(\"Imports Successful\")\n",
    "\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:21:39.075443Z",
     "start_time": "2021-05-20T16:21:39.056495Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_mod(a,n):\n",
    "    getcontext().prec = 27\n",
    "    return round((Decimal(a)/Decimal(n) - Decimal(a)//Decimal(n) ) * n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:51:57.551427Z",
     "start_time": "2021-05-20T16:51:57.527722Z"
    }
   },
   "outputs": [],
   "source": [
    "def constant_optimized_modular_exponentation_modulus15(a, power):\n",
    "    if a not in [2,7,8,11,13]:\n",
    "        raise ValueError(\"'a' must be 2,7,8,11 or 13\")\n",
    "    U = QuantumCircuit(4)        \n",
    "    for iteration in range(power):\n",
    "        if a in [2,13]:\n",
    "            U.swap(0,1)\n",
    "            U.swap(1,2)\n",
    "            U.swap(2,3)\n",
    "        if a in [7,8]:\n",
    "            U.swap(2,3)\n",
    "            U.swap(1,2)\n",
    "            U.swap(0,1)\n",
    "        if a == 11:\n",
    "            U.swap(1,3)\n",
    "            U.swap(0,2)\n",
    "        if a in [7,11,13]:\n",
    "            for q in range(4):\n",
    "                U.x(q)\n",
    "    U = U.to_gate()\n",
    "    U.name = \"%i^%i mod 15\" % (a, power)\n",
    "    control_U = U.control()\n",
    "    return control_U\n",
    "\n",
    "def inverse_qft(n):\n",
    "    circuit = QuantumCircuit(n)\n",
    "    \n",
    "    for i in range(n//2):\n",
    "        circuit.swap(i, n-1-i)\n",
    "        \n",
    "    for j in range(n):\n",
    "        for m in range(j):\n",
    "            circuit.cp(-np.pi/float(2**(j-m)), m, j)\n",
    "        circuit.h(j)\n",
    "    circuit.name = \"QFT†\"\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:49:58.246094Z",
     "start_time": "2021-05-20T16:49:54.646660Z"
    }
   },
   "outputs": [],
   "source": [
    "N = 15\n",
    "a = 7\n",
    "n_count = 8\n",
    "counting_register = QuantumRegister(size = n_count, name = \"counting_register\")\n",
    "acting_register = QuantumRegister(size = 4, name=\"acting_register\")\n",
    "classic_register = ClassicalRegister(size = n_count, name=\"classic_register\")\n",
    "qc = QuantumCircuit(counting_register, acting_register ,classic_register)\n",
    "initial_state = [1,0]\n",
    "\n",
    "for q in range(8):\n",
    "    qc.initialize(initial_state, q)\n",
    "qc.draw(output = 'mpl', filename = \"Step0\")\n",
    "\n",
    "for q in range(n_count):\n",
    "    qc.h(q)\n",
    "\n",
    "qc.draw(output = 'mpl', filename = \"Step1\")\n",
    "\n",
    "qc.x(3+n_count)\n",
    "\n",
    "qc.draw(output = 'mpl', filename = \"Step1b\")\n",
    "\n",
    "for q in range(n_count):\n",
    "    qc.append(constant_optimized_modular_exponentation_modulus15(a, 2**q), \n",
    "             [q] + [i+n_count for i in range(4)])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:47:24.112263Z",
     "start_time": "2021-05-20T16:47:22.598362Z"
    }
   },
   "outputs": [],
   "source": [
    "qc.measure(range(n_count,n_count + 4), range(4))\n",
    "qc.barrier()\n",
    "qc.draw(output = 'mpl', filename = \"Step2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:52:12.264039Z",
     "start_time": "2021-05-20T16:52:03.364918Z"
    }
   },
   "outputs": [],
   "source": [
    "qc.append(inverse_qft(n_count), range(n_count))\n",
    "qc.draw(output = 'mpl', filename = \"Step3\")\n",
    "\n",
    "# Measure circuit\n",
    "qc.measure(range(n_count), range(n_count))\n",
    "qc.draw(output = 'mpl', filename = \"Step4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:53:15.945286Z",
     "start_time": "2021-05-20T16:52:57.138633Z"
    }
   },
   "outputs": [],
   "source": [
    "qasm_sim = Aer.get_backend('qasm_simulator')\n",
    "t_qc = transpile(qc, qasm_sim)\n",
    "qobj = assemble(t_qc)\n",
    "results = qasm_sim.run(qobj).result()\n",
    "counts = results.get_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:54:08.214833Z",
     "start_time": "2021-05-20T16:54:07.941591Z"
    }
   },
   "outputs": [],
   "source": [
    "plot_histogram(counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T16:57:32.296109Z",
     "start_time": "2021-05-20T16:57:32.228154Z"
    }
   },
   "outputs": [],
   "source": [
    "rows, measured_phases = [], []\n",
    "for output in counts:\n",
    "    decimal = int(output, 2)  \n",
    "    phase = decimal/(2**n_count)  \n",
    "    measured_phases.append(phase)\n",
    "    \n",
    "    rows.append([f\"{output}(bin) = {decimal:>3}(dec)\", \n",
    "                 f\"{decimal}/{2**n_count} = {phase:.2f}\"])\n",
    "headers=[\"Register Output\", \"Phase\"]\n",
    "df = pd.DataFrame(rows, columns=headers)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T17:00:09.389291Z",
     "start_time": "2021-05-20T17:00:09.370463Z"
    }
   },
   "outputs": [],
   "source": [
    "rows = []\n",
    "for phase in measured_phases:\n",
    "    frac = Fraction(phase).limit_denominator(15)\n",
    "    rows.append([phase, f\"{frac.numerator}/{frac.denominator}\", frac.denominator])\n",
    "\n",
    "headers=[\"Phase\", \"Fraction\", \"Guess for r\"]\n",
    "df = pd.DataFrame(rows, columns=headers)\n",
    "my_period_r = max(df[\"Guess for r\"])\n",
    "print(\"My period (r) is %i\" % my_period_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T17:00:21.088572Z",
     "start_time": "2021-05-20T17:00:20.749883Z"
    }
   },
   "outputs": [],
   "source": [
    "# Confirm that the period is 4\n",
    "xvals = np.arange(N)\n",
    "xvals = [x.item() for x in xvals]\n",
    "yvals = [my_mod(a**x, N) for x in xvals]\n",
    "\n",
    "fig, ax = plt.subplots();\n",
    "ax.plot(xvals, yvals, linewidth=1, linestyle='dotted', marker='x');\n",
    "ax.set(xlabel='$x$', ylabel='$%i^x$ mod $%i$' % (a, N),\n",
    "       title=\"Example of Periodic Function in Shor's Algorithm\");\n",
    "try: \n",
    "    r = yvals[1:].index(1) +1 \n",
    "    plt.annotate(s = '', xy=(0,1), xytext=(r,1), arrowprops=dict(arrowstyle='<->'));\n",
    "    plt.annotate(s = '$r=%i$' % r, xy=(r/3,1.5));\n",
    "except ValueError:\n",
    "    print('Could not find a period')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-20T17:08:12.197591Z",
     "start_time": "2021-05-20T17:08:12.185056Z"
    }
   },
   "outputs": [],
   "source": [
    "first_shared_factor = math.gcd((7**(int(my_period_r/2)) + 1), 15)\n",
    "first_shared_factor\n",
    "second_shared_factor = math.gcd((7**(int(my_period_r/2)) - 1), 15)\n",
    "second_shared_factor"
   ]
  },
  {
   "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.3"
  },
  "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
}
