Shor's algorithm

Usage estimate: Three seconds on an Eagle r3 processor (NOTE: This is an estimate only. Your runtime might vary.)

Shor's algorithm, developed by Peter Shor in 1994, is a groundbreaking quantum algorithm for factoring integers in polynomial time. Its significance lies in its ability to factor large integers exponentially faster than any known classical algorithm, threatening the security of widely used cryptographic systems like RSA, which rely on the difficulty of factoring large numbers. By efficiently solving this problem on a sufficiently powerful quantum computer, Shor's algorithm could revolutionize fields such as cryptography, cybersecurity, and computational mathematics, underscoring the transformative power of quantum computation.

This tutorial focuses on demonstrating Shor's algorithm by factoring 15 on a quantum computer.

First, we define the order finding problem and construct corresponding circuits from the quantum phase estimation protocol. Next, we run the order finding circuits on real hardware using shortest-depth circuits we can transpile. The last section completes Shor's algorithm by connecting the order finding problem to integer factorization.

We end the tutorial with a discussion on other demonstrations of Shor's algorithm on real hardware, focusing on both generic implementations and those tailored to factoring specific integers such as 15 and 21.

Note: This tutorial focuses more on the implementation and demonstration of the circuits concerning Shor's algorithm. For an in-depth educational resource on the material, please refer to the Fundamentals of quantum algorithms course by Dr. John Watrous, and papers in the References section.

Requirements

Before starting this tutorial, ensure that you have the following installed:

• Qiskit SDK v2.0 or later with visualization support (pip install 'qiskit[visualization]')
• Qiskit Runtime v0.40 or later (pip install qiskit-ibm-runtime)

Setup

import numpy as np
import pandas as pd
from fractions import Fraction
from math import floor, gcd, log
 
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFT, UnitaryGate
from qiskit.transpiler import CouplingMap, generate_preset_pass_manager
from qiskit.visualization import plot_histogram
 
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Step 1: Map classical inputs to a quantum problem

Background

Shor's algorithm for integer factorization utilizes an intermediary problem known as the order finding problem. In this section, we demonstrate how to solve the order finding problem using quantum phase estimation.

Phase estimation problem

In the phase estimation problem, we're given a quantum state $\ket{\psi}$ of $n$ qubits, along with a unitary quantum circuit that acts on $n$ qubits. We're promised that $\ket{\psi}$ is an eigenvector of the unitary matrix $U$ that describes the action of the circuit, and our goal is to compute or approximate the eigenvalue $\lambda = e^{2 \pi i \theta}$ to which $\ket{\psi}$ corresponds. In other words, the circuit should output an approximation to the number $\theta \in [0, 1)$ satisfying $U \ket{\psi}= e^{2 \pi i \theta} \ket{\psi}.$ The goal of the phase estimation circuit is to approximate $\theta$ in $m$ bits. Mathematically speaking, we would like to find $y$ such that $\theta \approx y / 2^m$, where $y \in {0, 1, 2, \dots, 2^{m-1}}$. The following image shows the quantum circuit that estimates $y$ in $m$ bits by making a measurement on $m$ qubits.

In the above circuit, top $m$ qubits are initiated in the $\ket{0^m}$ state, and the bottom $n$ qubits are initiated in $\ket{\psi}$, which is promised to be an eigenvector of $U$. The first ingredient in the phase estimation circuit are the controlled-unitary operations that are responsible for performing a phase kickback to their corresponding control qubit. These controlled unitaries are exponentiated in accordance to the position of the control qubit, ranging from the least significant bit to the most significant bit. Since $\ket{\psi}$ is an eigenvector of $U$, the state of the bottom $n$ qubits is not affected by this operation, but the phase information of the eigenvalue propagates to the top $m$ qubits.

It turns out that after the phase kickback operation via controlled-unitaries, all possible states of the top $m$ qubits are orthonormal to each other for each eigenvector $\ket{\psi}$ of the unitary $U$. Therefore, these states are perfectly distinguishable, and we can rotate the basis they form back to the computational basis to make a measurement. A mathematical analysis shows that this rotation matrix corresponds to the inverse quantum Fourier transform (QFT) in $2^m$-dimensional Hilbert space. The intuition behind this is that the periodic structure of the modular exponentiation operators is encoded in the quantum state, and the QFT converts this periodicity into measurable peaks in the frequency domain.

For a more in-depth understanding of why the QFT circuit is employed in Shor's algorithm, we refer the reader to the Fundamentals of quantum algorithms course.

We are now ready to use the phase estimation circuit for order finding.

Order finding problem

To define the order finding problem, we begin with some number theory concepts. First, for any given positive integer $N$, define the set $\mathbb{Z}_N$ as $\mathbb{Z}_N = \{0, 1, 2, \dots, N-1\}.$ All arithmetic operations in $\mathbb{Z}_N$ are performed modulo $N$. In particular, all elements $a \in \mathbb{Z}_n$ that are coprime with $N$ are special and constitute $\mathbb{Z}^*_N$ as $\mathbb{Z}^*_N = \{ a \in \mathbb{Z}_N : \mathrm{gcd}(a, N)=1 \}.$ For an element $a \in \mathbb{Z}^*_N$, the smallest positive integer $r$ such that $a^r \equiv 1 \; (\mathrm{mod} \; N)$ is defined as the order of $a$ modulo $N$. As we will see later, finding the order of an $a \in \mathbb{Z}^*_N$ will allow us factor $N$.

To construct the order finding circuit from the phase estimation circuit, we need two considerations. First, we need to define the unitary $U$ that will allow us to find the order $r$, and second, we need to define an eigenvector $\ket{\psi}$ of $U$ to prepare the initial state of the phase estimation circuit.

To connect the order finding problem to phase estimation, we consider the operation defined on a system whose classical states correspond to $\mathbb{Z}_N$, where we multiply by a fixed element $a \in \mathbb{Z}^*_N$. In particular, we define this multiplication operator $M_a$ such that $M_a \ket{x} = \ket{ax \; (\mathrm{mod} \; N)}$ for each $x \in \mathbb{Z}_N$. Note that it's implicit that we're taking the product modulo $N$ inside of the ket on the right-hand side of the equation. A mathematical analysis shows that $M_a$ is an unitary operator. Furthermore, it turns out that $M_a$ has eigenvector and eigenvalue pairs that allow us to connect the order $r$ of $a$ to the phase estimation problem. Specifically, for any choice of $j \in \{0, \dots, r-1\}$, we have that $\ket{\psi_j} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \omega^{-jk}_{r} \ket{a^k}$ is an eigenvector of $M_a$ whose corresponding eigenvalue is $\omega^{j}_{r}$, where $\omega^{j}_{r} = e^{2 \pi i \frac{j}{r}}.$

By observation, we see that a convenient eigenvector/eigenvalue pair is the state $\ket{\psi_1}$ with $\omega^{1}_{r} = e^{2 \pi i \frac{1}{r}}$. Therefore, if we could find the eigenvector $\ket{\psi_1}$, we could estimate the phase $\theta=1/r$ with our quantum circuit and therefore get an estimate of the order $r$. However, it's not easy to do so, and we need to consider an alternative.

Let's consider what the circuit would result in if we prepare the computational state $\ket{1}$ as the initial state. This is not an eigenstate of $M_a$, but it is the uniform superposition of the eigenstates we just described above. In other words, the following relation holds. $\ket{1} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} \ket{\psi_k}$ The implication of the above equation is that if we set the initial state to $\ket{1}$, we will obtain precisely the same measurement result as if we had chosen $k \in \{ 0, \dots, r-1\}$ uniformly at random and used $\ket{\psi_k}$ as an eigenvector in the phase estimation circuit. In other words, a measurement of the top $m$ qubits yields an approximation $y / 2^m$ to the value $k / r$ where $k \in \{ 0, \dots, r-1\}$ is chosen uniformly at random. This allows us to learn $r$ with a high degree of confidence after several independent runs, which was our goal.

Modular exponentiation operators

So far, we linked the phase estimation problem to the order finding problem by defining $U = M_a$ and $\ket{\psi} = \ket{1}$ in our quantum circuit. Therefore, the last remaining ingredient is to find an efficient way to define modular exponentials of $M_a$ as $M_a^k$ for $k = 1, 2, 4, \dots, 2^{m-1}$. To perform this computation, we find that for any power $k$ we choose, we can create a circuit for $M_a^k$ not by iterating $k$ times the circuit for $M_a$, but instead by computing $b = a^k \; \mathrm{mod} \; N$ and then using the circuit for $M_b$. Since we only need the powers that are powers of 2 themselves, we can do this classically efficiently by using iterative squaring.

Step 2: Optimize problem for quantum hardware execution

Specific example with $N = 15$ and $a=2$

We can pause here to discuss a specific example and construct the order finding circuit for $N=15$. Note that possible nontrivial $a \in \mathbb{Z}_N^*$ for $N=15$ are $a \in \{2, 4, 7, 8, 11, 13, 14 \}$. For this example, we choose $a=2$. We will construct the $M_2$ operator and the modular exponentiation operators $M_2^k$.

The action of $M_2$ on the computational basis states is as follows. $M_2 \ket{0} = \ket{0} \quad M_2 \ket{5} = \ket{10} \quad M_2 \ket{10} = \ket{5}$ $M_2 \ket{1} = \ket{2} \quad M_2 \ket{6} = \ket{12} \quad M_2 \ket{11} = \ket{7}$ $M_2 \ket{2} = \ket{4} \quad M_2 \ket{7} = \ket{14} \quad M_2 \ket{12} = \ket{9}$ $M_2 \ket{3} = \ket{6} \quad M_2 \ket{8} = \ket{1} \quad M_2 \ket{13} = \ket{11}$ $M_2 \ket{4} = \ket{8} \quad M_2 \ket{9} = \ket{3} \quad M_2 \ket{14} = \ket{13}$ By observation, we can see that the basis states are shuffled, so we have a permutation matrix. We can construct this operation on four qubits with swap gates. Below, we construct the $M_2$ and the controlled-$M_2$ operations.

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)
 
    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)
 
    U = U.to_gate()
    U.name = f"M_{b}"
 
    return U

# Get the M2 operator
M2 = M2mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)
 
    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)
 
    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()
 
    return c_U

# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Gates acting on more than two qubits will be further decomposed into two-qubit gates.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

Now we need to construct the modular exponentiation operators. To obtain enough precision in the phase estimation, we will use eight qubits for the estimation measurement. Therefore, we need to construct $M_b$ with $b = a^{2^k} \; (\mathrm{mod} \; N)$ for each $k = 0, 1, \dots, 7$.

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a

k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]
 
print(b_list)

Output:

[2, 4, 1, 1, 1, 1, 1, 1]

As we can see from the list of $b$ values, in addition to $M_2$ that we previously constructed, we also need to build $M_4$ and $M_1$. Note that $M_1$ acts trivially on the computational basis states, so it is simply the identity operator.

$M_4$ acts on the computational basis states as follows. $M_4 \ket{0} = \ket{0} \quad M_4 \ket{5} = \ket{5} \quad M_4 \ket{10} = \ket{10}$ $M_4 \ket{1} = \ket{4} \quad M_4 \ket{6} = \ket{9} \quad M_4 \ket{11} = \ket{14}$ $M_4 \ket{2} = \ket{8} \quad M_4 \ket{7} = \ket{13} \quad M_4 \ket{12} = \ket{3}$ $M_4 \ket{3} = \ket{12} \quad M_4 \ket{8} = \ket{2} \quad M_4 \ket{13} = \ket{7}$ $M_4 \ket{4} = \ket{1} \quad M_4 \ket{9} = \ket{6} \quad M_4 \ket{14} = \ket{11}$

Therefore, this permutation can be constructed with the following swap operation.

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)
 
    U.swap(1, 3)
    U.swap(0, 2)
 
    U = U.to_gate()
    U.name = f"M_{b}"
 
    return U

# Get the M4 operator
M4 = M4mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)
 
    U.swap(1, 3)
    U.swap(0, 2)
 
    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()
 
    return c_U

# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Gates acting on more than two qubits will be further decomposed into two-qubit gates.

circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

We saw that $M_b$ operators for a given $b \in \mathbb{Z}^*_N$ are permutation operations. Due to the relatively small size of the permutation problem that we have here, since $N=15$ requires only four qubits, we were able to synthesize these operations directly with SWAP gates by inspection. In general, this might not be a scalable approach. Instead, we might need to construct the permutation matrix explicitly, and use Qiskit's UnitaryGate class and transpilation methods to synthesize this permutation matrix. However, this can result in significantly deeper circuits. An example follows.

def mod_mult_gate(b, N):
    """
    Modular multiplication gate from permutation matrix.
    """
    if gcd(b, N) > 1:
        print(f"Error: gcd({b},{N}) > 1")
    else:
        n = floor(log(N - 1, 2)) + 1
        U = np.full((2**n, 2**n), 0)
        for x in range(N):
            U[b * x % N][x] = 1
        for x in range(N, 2**n):
            U[x][x] = 1
        G = UnitaryGate(U)
        G.name = f"M_{b}"
        return G

# Let's build M2 using the permutation matrix definition
M2_other = mod_mult_gate(2, 15)
 
# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2_other, inplace=True)
circ = circ.decompose()
 
# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)
 
print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.decompose().draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

qubits: 4
2q-depth: 94
2q-size: 96
Operator counts: OrderedDict({'cx': 45, 'swap': 32, 'u': 24, 'u1': 7, 'u3': 4, 'unitary': 3, 'circuit-335': 1, 'circuit-338': 1, 'circuit-341': 1, 'circuit-344': 1, 'circuit-347': 1, 'circuit-350': 1, 'circuit-353': 1, 'circuit-356': 1, 'circuit-359': 1, 'circuit-362': 1, 'circuit-365': 1, 'circuit-368': 1, 'circuit-371': 1, 'circuit-374': 1, 'circuit-377': 1, 'circuit-380': 1})

Let's compare these counts with the compiled circuit depth of our manual implementation of the $M_2$ gate.

# Get the M2 operator from our manual construction
M2 = M2mod15()
 
# Add it to a circuit
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ = circ.decompose(reps=3)
 
# Transpile the circuit and get the depth
coupling_map = CouplingMap.from_line(4)
pm = generate_preset_pass_manager(coupling_map=coupling_map)
transpiled_circ = pm.run(circ)
 
print(f"qubits: {circ.num_qubits}")
print(
    f"2q-depth: {transpiled_circ.depth(lambda x: x.operation.num_qubits==2)}"
)
print(f"2q-size: {transpiled_circ.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circ.count_ops()}")
transpiled_circ.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

qubits: 4
2q-depth: 9
2q-size: 9
Operator counts: OrderedDict({'cx': 9})

As we can see, the permutation matrix approach resulted in a significantly deep circuit even for a single $M_2$ gate compared to our manual implementation of it. Therefore, we will continue with our previous implementation of the $M_b$ operations.

Now, we are ready to construct the full order finding circuit using our previously defined controlled modular exponentiation operators. In the following code, we also import the QFT circuit from the Qiskit Circuit library, which uses Hadamard gates on each qubit, a series of controlled-U1 (or Z, depending on the phase) gates, and a layer of swap gates.

# Order finding problem for N = 15 with a = 2
N = 15
a = 2
 
# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation
 
# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]
 
# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)
 
# Initialize the target register to the state |1>
circuit.x(num_control)
 
# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator
 
# Apply the inverse QFT to the control register
circuit.compose(QFT(num_control, inverse=True), qubits=control, inplace=True)
 
# Measure the control register
circuit.measure(control, output)
 
circuit.draw("mpl", fold=-1)

Output:

Note that we omitted the controlled modular exponentiation operations from the remaining control qubits because $M_1$ is the identity operator.

Note that later in this tutorial, we will run this circuit on the ibm_marrakesh backend. To do this, we transpile the circuit according to this specific backend and report the circuit depth and gate counts.

service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")
pm = generate_preset_pass_manager(optimization_level=2, backend=backend)
 
transpiled_circuit = pm.run(circuit)
 
print(
    f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}"
)
print(
    f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}"
)
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(
    output="mpl", fold=-1, style="clifford", idle_wires=False
)

Output:

2q-depth: 187
2q-size: 260
Operator counts: OrderedDict({'sx': 521, 'rz': 354, 'cz': 260, 'measure': 8, 'x': 4})

Step 3: Execute using Qiskit primitives

First, we discuss what we would theoretically obtain if we ran this circuit on an ideal simulator. Below, we have a set of simulation results of the above circuit using 1024 shots. As we can see, we get an approximately uniform distribution over four bitstrings over the control qubits.

# Obtained from the simulator
counts = {"00000000": 264, "01000000": 268, "10000000": 249, "11000000": 243}

plot_histogram(counts)

Output:

By measuring the control qubits, we obtain an eight-bit phase estimation of the $M_a$ operator. We can convert this binary representation to decimal to find the measured phase. As we can see from the above histogram, four different bitstrings were measured, and each of them corresponds to a phase value as follows.

# Rows to be displayed in table
rows = []
# Corresponding phase of each bitstring
measured_phases = []
 
for output in counts:
    decimal = int(output, 2)  # Convert bitstring to decimal
    phase = decimal / (2**num_control)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append(
        [
            f"{output}(bin) = {decimal:>3}(dec)",
            f"{decimal}/{2 ** num_control} = {phase:.2f}",
        ]
    )
 
# Print the rows in a table
headers = ["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

            Register Output           Phase
0  00000000(bin) =   0(dec)    0/256 = 0.00
1  01000000(bin) =  64(dec)   64/256 = 0.25
2  10000000(bin) = 128(dec)  128/256 = 0.50
3  11000000(bin) = 192(dec)  192/256 = 0.75

Recall that the any measured phase corresponds to $\theta = k / r$ where $k$ is sampled uniformly at random from $\{0, 1, \dots, r-1 \}$. Therefore, we can use the continued fractions algorithm to attempt to find $k$ and the order $r$. Python has this functionality built in. We can use the fractions module to turn a float into a Fraction object, for example:

Fraction(0.666)

Output:

Fraction(5998794703657501, 9007199254740992)

Because this gives fractions that return the result exactly (in this case, 0.6660000...), this can give gnarly results like the one above. We can use the .limit_denominator() method to get the fraction that most closely resembles our float, with a denominator below a certain value:

# Get fraction that most closely resembles 0.666
# with denominator < 15
Fraction(0.666).limit_denominator(15)

Output:

Fraction(2, 3)

This is much nicer. The order (r) must be less than N, so we will set the maximum denominator to be 15:

# Rows to be displayed in a table
rows = []
 
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append(
        [phase, f"{frac.numerator}/{frac.denominator}", frac.denominator]
    )
 
# Print the rows in a table
headers = ["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Output:

   Phase Fraction  Guess for r
0   0.00      0/1            1
1   0.25      1/4            4
2   0.50      1/2            2
3   0.75      3/4            4

We can see that two of the measured eigenvalues provided us with the correct result: $r=4$, and we can see that Shor’s algorithm for order finding has a chance of failing. These bad results are because $k = 0$, or because $k$ and $r$ are not coprime - and instead of $r$, we are given a factor of $r$. The easiest solution to this is to simply repeat the experiment until we get a satisfying result for $r$.

So far, we implemented the order finding problem for $N=15$ with $a=2$ using the phase estimation circuit on a simulator. The last step of Shor's algorithm will be to relate the order finding problem to the integer factorization problem. This last part of the algorithm is purely classical and can be solved on a classical computer after the phase measurements have been obtained from a quantum computer. Therefore, we defer the last part of the algorithm until after we demonstrate how we can run the order finding circuit on real hardware.

Hardware runs

Now we can run the order finding circuit that we previously transpiled for ibm_marrakesh. Here we turn to dynamical decoupling (DD) for error suppression, and gate twirling for error mitigation purposes. DD involves applying sequences of precisely timed control pulses to a quantum device, effectively averaging out unwanted environmental interactions and decoherence. Gate twirling, on the other hand, randomizes specific quantum gates to transform coherent errors into Pauli errors, which accumulate linearly rather than quadratically. Both techniques are often combined to enhance the coherence and fidelity of quantum computations.

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)
 
# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True
 
pub = transpiled_circuit
job = sampler.run([pub], shots=1024)

result = job.result()[0]
counts = result.data["out"].get_counts()

plot_histogram(counts, figsize=(35, 5))

Output:

As we can see, we obtained the same bitstrings with highest counts. Since quantum hardware has noise, there is some leakage to other bitstrings, which we can filter out statistically.

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2
 
for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value
 
print(counts_keep)

Output:

{'00000000': 58, '01000000': 41, '11000000': 42, '10000000': 40}

Step 4: Post-process and return result in desired classical format

Integer Factorization

So far, we discussed how we can implement the order finding problem using a phase estimation circuit. Now, we connect the order finding problem to integer factorization, which completes Shor's algorithm. Note that this part of the algorithm is classical.

We now demonstrate this using our example of $N = 15$ and $a = 2$. Recall that the phase we measured is $k / r$, where $a^r \; (\textrm{mod} \; N) = 1$ and $k$ is a random integer between $0$ and $r - 1$. From this equation, we have $(a^r - 1) \; (\textrm{mod} \; N) = 0,$ which means $N$ must divide $a^r-1$. If $r$ is also even, then we can write $a^r -1 = (a^{r/2}-1)(a^{r/2}+1).$ If $r$ is not even, we cannot go further and must try again with a different value for $a$; otherwise, there is 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$.

Since some runs of the algorithm will statistically fail, we will repeat this algorithm until at least one factor of $N$ is found.

The cell below repeats the algorithm until at least one factor of $N=15$ is found. We will use the results of the hardware run above to guess the phase and the corresponding factor in each iteration.

a = 2
N = 15
 
FACTOR_FOUND = False
num_attempt = 0
 
while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")
 
    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")
 
    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

Output:

ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.25
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

Discussion

Related work

In this section, we discuss other milestone work that has demonstrated Shor's algorithm on real hardware.

The seminal work [3] from IBM® demonstrated Shor's algorithm for the first time, factoring the number 15 into its prime factors 3 and 5 using a seven-qubit nuclear magnetic resonance (NMR) quantum computer. Another experiment [4] factored 15 using photonic qubits. By employing a single qubit recycled multiple times and encoding the work register in higher-dimensional states, the researchers reduced the required number of qubits to one-third of that in the standard protocol, utilizing a two-photon compiled algorithm. A significant paper in the demonstration of Shor's algorithm is [5], which uses Kitaev's iterative phase estimation [8] technique to reduce the qubit requirement of the algorithm. Authors used seven control qubits and four cache qubits, together with the implementation of modular multipliers. This implementation, however, requires mid-circuit measurements with feed-forward operations and qubit recycling with reset operations. This demonstration was done on an ion-trap quantum computer.

More recent work [6] focused on factoring 15, 21, and 35 on IBM Quantum® hardware. Similar to previous work, researchers used a compiled version of the algorithm that employed a semi-classical quantum Fourier transform as proposed by Kitaev to minimize the number of physical qubits and gates. A most recent work [7] also performed a proof-of-concept demonstration for factoring the integer 21. This demonstration also involved the use of a compiled version of the quantum phase estimation routine, and built upon the previous demonstration by [4]. Authors went beyond this work by using a configuration of approximate Toffoli gates with residual phase shifts. The algorithm was implemented on IBM quantum processors using only five qubits, and the presence of entanglement between the control and register qubits was verified successfully.

Scaling of the algorithm

We note that RSA encryption typically involves key sizes on the order of 2048 to 4096 bits. Attempting to factor a 2048-bit number with Shor's algorithm will result in a quantum circuit with millions of qubits, including the error correction overhead and a circuit depth on the order of a billion, which is beyond the limits of current quantum hardware to execute. Therefore, Shor's algorithm will require either optimized circuit construction methods or robust quantum error correction to be practically viable for breaking modern cryptographic systems. We refer you to [9] for a more detailed discussion on resource estimation for Shor's algorithm.

Challenge

Congratulations for finishing the tutorial! Now is a great time to test your understanding. Could you try to construct the circuit for factoring 21? You can select an $a$ of your own choice. You will need to decide on the bit accuracy of the algorithm to choose the number of qubits, and you will need to design the modular exponentiation operators $M_a$. We encourage you to try this out yourself, and then read about the methodologies shown in Fig. 9 of [6] and Fig. 2 of [7].

def M_a_mod21():
    """
    M_a (mod 21)
    """
 
    # Your code here
    pass

References

1. Shor, Peter W. "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer." SIAM review 41.2 (1999): 303-332.
2. IBM Quantum Fundamentals of Quantum Algorithms course by Dr. John Watrous.
3. Vandersypen, Lieven MK, et al. "Experimental realization of Shor's quantum factoring algorithm using nuclear magnetic resonance." Nature 414.6866 (2001): 883-887.
4. Martin-Lopez, Enrique, et al. "Experimental realization of Shor's quantum factoring algorithm using qubit recycling." Nature photonics 6.11 (2012): 773-776.
5. Monz, Thomas, et al. "Realization of a scalable Shor algorithm." Science 351.6277 (2016): 1068-1070.
6. Amico, Mirko, Zain H. Saleem, and Muir Kumph. "Experimental study of Shor's factoring algorithm using the IBM Q Experience." Physical Review A 100.1 (2019): 012305.
7. Skosana, Unathi, and Mark Tame. "Demonstration of Shor’s factoring algorithm for N=21 on IBM quantum processors." Scientific reports 11.1 (2021): 16599.
8. Kitaev, A. Yu. "Quantum measurements and the Abelian stabilizer problem." arXiv preprint quant-ph/9511026 (1995).
9. Gidney, Craig, and Martin Ekerå. "How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits." Quantum 5 (2021): 433.