{
"cells": [
{
"cell_type": "code",
"execution_count": 2,
"id": "8b60e095",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:19:30.859229Z",
"start_time": "2023-02-08T10:19:24.424011Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy.random import randint\n",
"from fractions import Fraction\n",
"from qiskit import assemble, Aer, transpile, QuantumCircuit, ClassicalRegister\n",
"from diskit import *"
]
},
{
"cell_type": "markdown",
"id": "928ff400",
"metadata": {},
"source": [
"# Distribute Shor's Algorithm"
]
},
{
"cell_type": "markdown",
"id": "9b5a290f",
"metadata": {},
"source": [
"The code for building the monolithic circuit for factoring 15 is taken from Qiskit Library. <br>\n",
"https://qiskit.org/textbook/ch-algorithms/shor.html\n",
"<br><br>\n",
"\n",
"In this example we will solve the period finding problem for $a=7$ and $N=15$. We provide the circuits for $U$ where:\n",
"\n",
"$$U|y\\rangle = |ay\\bmod 15\\rangle $$\n",
"\n",
"without explanation. To create $U^x$, we will simply repeat the circuit $x$ times. In the next section we will discuss a general method for creating these circuits efficiently. The function `c_amod15` returns the controlled-U gate for `a`, repeated `power` times."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b9c067ef",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:40:44.628790Z",
"start_time": "2023-02-08T10:40:44.608776Z"
},
"code_folding": []
},
"outputs": [],
"source": [
"def c_amod15(a, power):\n",
" \"\"\"Controlled multiplication by a mod 15\"\"\"\n",
" if a not in [2,4,7,8,11,13]:\n",
" raise ValueError(\"'a' must be 2,4,7,8,11 or 13\")\n",
" U = QuantumCircuit(4) \n",
" for iteration in range(power):\n",
" if a in [2,13]:\n",
" U.swap(2,3)\n",
" U.swap(1,2)\n",
" U.swap(0,1)\n",
" if a in [7,8]:\n",
" U.swap(0,1)\n",
" U.swap(1,2)\n",
" U.swap(2,3)\n",
" if a in [4, 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",
" c_U = U.control()\n",
" return c_U"
]
},
{
"cell_type": "markdown",
"id": "e8eccb80",
"metadata": {},
"source": [
"We also import the circuit for the QFT (you can read more about the QFT in the [quantum Fourier transform chapter](./quantum-fourier-transform.html#generalqft)):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c7701102",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:19:30.989315Z",
"start_time": "2023-02-08T10:19:30.959278Z"
}
},
"outputs": [],
"source": [
"def qft_dagger(n):\n",
" \"\"\"n-qubit QFTdagger the first n qubits in circ\"\"\"\n",
" qc = QuantumCircuit(n)\n",
" # Don't forget the Swaps!\n",
" for qubit in range(n//2):\n",
" qc.swap(qubit, n-qubit-1)\n",
" for j in range(n):\n",
" for m in range(j):\n",
" qc.cp(-np.pi/float(2**(j-m)), m, j)\n",
" qc.h(j)\n",
" qc.name = \"QFT†\"\n",
" return qc"
]
},
{
"cell_type": "markdown",
"id": "6a899281",
"metadata": {},
"source": [
"Factoring from Period Finding\n",
"\n",
"Not all factoring problems are difficult; we can spot an even number instantly and know that one of its factors is 2. In fact, there are [specific criteria](https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf#%5B%7B%22num%22%3A127%2C%22gen%22%3A0%7D%2C%7B%22name%22%3A%22XYZ%22%7D%2C70%2C223%2C0%5D) for choosing numbers that are difficult to factor, but the basic idea is to choose the product of two large prime numbers.\n",
"\n",
"A general factoring algorithm will first check to see if there is a shortcut to factoring the integer (is the number even? Is the number of the form $N = a^b$?), before using Shor’s period finding for the worst-case scenario. Since we aim to focus on the quantum part of the algorithm, we will jump straight to the case in which N is the product of two primes.\n",
"\n",
"To see an example of factoring on a small number of qubits, we will factor 15, which we all know is the product of the not-so-large prime numbers 3 and 5."
]
},
{
"cell_type": "markdown",
"id": "077b96db",
"metadata": {},
"source": [
"We will use 8 counting qubits, so in the topology the number of qubits must be $\\geq$ 8 + 4 "
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c671feed",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:40:38.470084Z",
"start_time": "2023-02-08T10:40:38.460067Z"
}
},
"outputs": [],
"source": [
"N = 15"
]
},
{
"cell_type": "markdown",
"id": "58b265ec",
"metadata": {},
"source": [
"The first step is to choose a random number, $a$, between $1$ and $N-1$:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "334404b0",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:40:39.771286Z",
"start_time": "2023-02-08T10:40:39.751277Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"7\n"
]
}
],
"source": [
"np.random.seed(1) # This is to make sure we get reproduceable results\n",
"a = randint(2, 15)\n",
"print(a)"
]
},
{
"cell_type": "markdown",
"id": "1db84ab7",
"metadata": {},
"source": [
"Next we quickly check it isn't already a non-trivial factor of $N$:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "b47f5cdc",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:40:43.598013Z",
"start_time": "2023-02-08T10:40:43.578002Z"
}
},
"outputs": [
{
"data": {
"text/plain": "1"
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from math import gcd # greatest common divisor\n",
"gcd(a, N)"
]
},
{
"cell_type": "markdown",
"id": "6027124c",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:42:42.422103Z",
"start_time": "2023-02-08T10:42:42.402095Z"
}
},
"source": [
"At First we see the monolithic implementation."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "b4833cd2",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T10:50:04.952182Z",
"start_time": "2023-02-08T10:50:04.932162Z"
}
},
"outputs": [],
"source": [
"def qpe_amod15(a):\n",
" n_count = 8\n",
" qc = QuantumCircuit(4+n_count, n_count)\n",
" for q in range(n_count):\n",
" qc.h(q) # Initialize counting qubits in state |+>\n",
" qc.x(3+n_count) # And auxiliary register in state |1>\n",
" for q in range(n_count): # Do controlled-U operations\n",
" qc.append(c_amod15(a, 2**q), \n",
" [q] + [i+n_count for i in range(4)])\n",
" qc.append(qft_dagger(n_count), range(n_count)) # Do inverse-QFT\n",
" qc.measure(range(n_count), range(n_count))\n",
" # Simulate Results\n",
" aer_sim = Aer.get_backend('aer_simulator')\n",
" # Setting memory=True below allows us to see a list of each sequential reading\n",
" t_qc = transpile(qc, aer_sim)\n",
" qobj = assemble(t_qc, shots=1)\n",
" result = aer_sim.run(qobj, memory=True).result()\n",
" readings = result.get_memory()\n",
" print(\"Register Reading: \" + readings[0])\n",
" phase = int(readings[0],2)/(2**n_count)\n",
" print(\"Corresponding Phase: %f\" % phase)\n",
" return phase, qc"
]
},
{
"cell_type": "markdown",
"id": "0eee5e8d",
"metadata": {},
"source": [
"We can modify the monolithic version to a distributed version of the function conveniently with minimal changes in structure and data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "d8a00508",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:01:06.316801Z",
"start_time": "2023-02-08T11:01:06.286797Z"
}
},
"outputs": [],
"source": [
"def qpe_amod15_dist(a, topology):\n",
" qregs = topology.get_regs()\n",
" n_count = 8\n",
" qc = QuantumCircuit(*qregs, ClassicalRegister(n_count))\n",
" for q in range(n_count):\n",
" qc.h(q) # Initialize counting qubits in state |+>\n",
" qc.x(3+n_count) # And auxiliary register in state |1>\n",
" for q in range(n_count): # Do controlled-U operations\n",
" qc.append(c_amod15(a, 2**q), \n",
" [q] + [i+n_count for i in range(4)])\n",
" qc.append(qft_dagger(n_count), range(n_count)) # Do inverse-QFT\n",
" qc.measure(range(n_count), range(n_count))\n",
" remapper = CircuitRemapper(circuit_topo)\n",
" dist_circ = remapper.remap_circuit(qc, decompose=True)\n",
"\n",
" # Simulate Results\n",
" aer_sim = Aer.get_backend('aer_simulator')\n",
" \n",
" # Setting memory=True below allows us to see a list of each sequential reading\n",
" t_qc = transpile(dist_circ, aer_sim)\n",
" qobj = assemble(t_qc, shots=1)\n",
" result = aer_sim.run(qobj, memory=True).result()\n",
" readings = result.get_memory() \n",
" \n",
" ## Need to manipulate readouts to exclude cat_measurement registers\n",
" final_readings = readings[0].split()\n",
" for reading in final_readings:\n",
" if len(reading) == n_count:\n",
" final_reading = reading\n",
" \n",
" print(\"Register Reading: \" + final_reading)\n",
" phase = int(final_reading,2)/(2**n_count)\n",
" print(\"Corresponding Phase: %f\" % phase)\n",
" return phase, dist_circ"
]
},
{
"cell_type": "markdown",
"id": "ea319dd4",
"metadata": {},
"source": [
"**Again, Topology qubits must be $\\geq$ 8 + 4**, we have taken 13 qubits with distribution over 2 qpus"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "4e49ab68",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:01:07.968911Z",
"start_time": "2023-02-08T11:01:07.948908Z"
}
},
"outputs": [
{
"data": {
"text/plain": "({'sys_shor0': [Qubit(QuantumRegister(6, 'sys_shor0'), 0),\n Qubit(QuantumRegister(6, 'sys_shor0'), 1),\n Qubit(QuantumRegister(6, 'sys_shor0'), 2),\n Qubit(QuantumRegister(6, 'sys_shor0'), 3),\n Qubit(QuantumRegister(6, 'sys_shor0'), 4),\n Qubit(QuantumRegister(6, 'sys_shor0'), 5)],\n 'sys_shor1': [Qubit(QuantumRegister(7, 'sys_shor1'), 0),\n Qubit(QuantumRegister(7, 'sys_shor1'), 1),\n Qubit(QuantumRegister(7, 'sys_shor1'), 2),\n Qubit(QuantumRegister(7, 'sys_shor1'), 3),\n Qubit(QuantumRegister(7, 'sys_shor1'), 4),\n Qubit(QuantumRegister(7, 'sys_shor1'), 5),\n Qubit(QuantumRegister(7, 'sys_shor1'), 6)]},\n {'sys_shor0': Qubit(QuantumRegister(1, 'com_sys_shor0'), 0),\n 'sys_shor1': Qubit(QuantumRegister(1, 'com_sys_shor1'), 0)})"
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"circuit_topo = Topology()\n",
"circuit_topo.create_qmap(2, [6, 7],\"sys_shor\")\n",
"circuit_topo.qmap, circuit_topo.emap"
]
},
{
"cell_type": "markdown",
"id": "6b70eec6",
"metadata": {},
"source": [
"From this phase, we can easily find a guess for $r$:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "1eb37bc1",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:08:47.046439Z",
"start_time": "2023-02-08T11:01:23.534506Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Single level decomposition cutoff of 1 minute reached. Performing transpilation with basis gates: cx, u1, u2, u3, id\n",
"Register Reading: 11000000\n",
"Corresponding Phase: 0.750000\n"
]
},
{
"data": {
"text/plain": "Fraction(3, 4)"
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phase, circ = qpe_amod15_dist(a, circuit_topo)\n",
"Fraction(phase).limit_denominator(15) # Denominator should (hopefully!) tell us r"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "6a9066aa",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:09:26.056248Z",
"start_time": "2023-02-08T11:09:26.036247Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4\n"
]
}
],
"source": [
"frac = Fraction(phase).limit_denominator(15)\n",
"s, r = frac.numerator, frac.denominator\n",
"print(r)"
]
},
{
"cell_type": "markdown",
"id": "ed53e24e",
"metadata": {},
"source": [
"Now we have $r$, we might be able to use this to find a factor of $N$. Since:\n",
"\n",
"$$a^r \\bmod N = 1 $$\n",
"\n",
"then:\n",
"\n",
"$$(a^r - 1) \\bmod N = 0 $$\n",
"\n",
"which means $N$ must divide $a^r-1$. And if $r$ is also even, then we can write:\n",
"\n",
"$$a^r -1 = (a^{r/2}-1)(a^{r/2}+1)$$\n",
"\n",
"(if $r$ is not even, we cannot go further and must try again with a different value for $a$). There is then a high probability that the greatest common divisor of $N$ and either $a^{r/2}-1$, or $a^{r/2}+1$ is a proper factor of $N$ [2]:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "76b9c8ca",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:09:30.032264Z",
"start_time": "2023-02-08T11:09:30.012277Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[3, 5]\n"
]
}
],
"source": [
"guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]\n",
"print(guesses)"
]
},
{
"cell_type": "markdown",
"id": "3797f81f",
"metadata": {},
"source": [
"The cell below repeats the algorithm until at least one factor of 15 is found. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "dc4b7fd3",
"metadata": {
"ExecuteTime": {
"end_time": "2023-02-08T11:17:59.247984Z",
"start_time": "2023-02-08T11:10:07.037441Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Attempt 1:\n",
"Single level decomposition cutoff of 1 minute reached. Performing transpilation with basis gates: cx, u1, u2, u3, id\n",
"Register Reading: 01000000\n",
"Corresponding Phase: 0.250000\n",
"Result: r = 4\n",
"Guessed Factors: 3 and 5\n",
"*** Non-trivial factor found: 3 ***\n",
"*** Non-trivial factor found: 5 ***\n"
]
}
],
"source": [
"a = 7\n",
"factor_found = False\n",
"attempt = 0\n",
"while not factor_found:\n",
" attempt += 1\n",
" print(\"\\nAttempt %i:\" % attempt)\n",
" phase,_ = qpe_amod15_dist(a, circuit_topo) # Phase = s/r\n",
" frac = Fraction(phase).limit_denominator(N) # Denominator should (hopefully!) tell us r\n",
" r = frac.denominator\n",
" print(\"Result: r = %i\" % r)\n",
" if phase != 0:\n",
" # Guesses for factors are gcd(x^{r/2} ±1 , 15)\n",
" guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]\n",
" print(\"Guessed Factors: %i and %i\" % (guesses[0], guesses[1]))\n",
" for guess in guesses:\n",
" if guess not in [1,N] and (N % guess) == 0: # Check to see if guess is a factor\n",
" print(\"*** Non-trivial factor found: %i ***\" % guess)\n",
" factor_found = True"
]
},
{
"cell_type": "markdown",
"id": "58d72c18",
"metadata": {},
"source": [
"References\n",
"\n",
"1. Stephane Beauregard, _Circuit for Shor's algorithm using 2n+3 qubits,_ [arXiv:quant-ph/0205095](https://arxiv.org/abs/quant-ph/0205095)\n",
"\n",
"2. M. Nielsen and I. Chuang, _Quantum Computation and Quantum Information,_ Cambridge Series on Information and the Natural Sciences (Cambridge University Press, Cambridge, 2000). (Page 633)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.7.10"
},
"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
}
},
"nbformat": 4,
"nbformat_minor": 5
}