# The original code is sourced from:
# https://github.com/munich-quantum-toolkit/bench/blob/main/src/mqt/bench/benchmarks/shor.py
# https://quantum.cloud.ibm.com/docs/en/tutorials/shors-algorithm.
#
# This program is adapted for use as a benchmark in controlled software testing experiments.
# Modifications made to the original code include (for Apache License 2.0):
# - Modify the code in the `QuantumCircuit` class
# - Decouple the dependency:
# - Use the local QFT implementation instead of Qiskit's built-in QFT
# - Move the synthesis function from `qft_decompose_full.py` to here
#
# Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
# Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
# All rights reserved.
#
# SPDX-License-Identifier: MIT
#
# Licensed under the MIT License
"""Shor benchmark definition."""
# Code is based on Qiskit, therefore we have to add the following information:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2019, 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
from __future__ import annotations
from typing import TYPE_CHECKING
if TYPE_CHECKING: # pragma: no cover
from numpy.typing import NDArray
import math
import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit, QuantumRegister
# from qiskit.synthesis import synth_qft_full
# Local QFT
from ..quantum_fourier_transform.qft import QFT
_SIZE_TO_PARAMS = {
18: (15, 4), # "small"
42: (821, 4), # "medium"
58: (11777, 4), # "large"
74: (201209, 4), # "xlarge"
}
# ---------- benchmark registration ----------
from ..benchmark_registry import register_benchmark
from pathlib import Path
@register_benchmark(
Path(__file__).stem,
family=Path(__file__).resolve().parent.name,
description="Shor's algorithm for integer factorization",
class_name="Shor",
source={
"repo": "https://github.com/munich-quantum-toolkit/bench/tree/main",
"file": "src/mqt/bench/benchmarks/shor.py",
"sdk": "Qiskit",
"available_doc": False
},
testability_refactoring=[
"Structure reorganization",
"Dependency decoupling"
]
)
def create_shor(to_be_factored_number: int, a: int):
return Shor(to_be_factored_number=to_be_factored_number, a=a)
# def create_circuit(circuit_size: int) -> QuantumCircuit:
# """Construct Shor's circuit based on total qubit count.
# Arguments:
# circuit_size: one of 18, 42, 58, 74.
# Raises:
# ValueError: if the size is not available.
# Returns:
# QuantumCircuit implementing Shor's algorithm for the chosen size.
# """
# try:
# n, a = _SIZE_TO_PARAMS[circuit_size]
# except KeyError as exc:
# valid = ", ".join(map(str, sorted(_SIZE_TO_PARAMS)))
# msg = f"No Shor instance for circuit_size={circuit_size}. Available: {valid}."
# raise ValueError(msg) from exc
# return create_circuit_from_num_and_coprime(n, a)
# def create_circuit_from_num_and_coprime(num_to_be_factorized: int, a: int = 2) -> QuantumCircuit:
# """Returns a quantum circuit implementing the Shor's algorithm.
# Arguments:
# num_to_be_factorized: number which shall be factorized
# a: any integer that satisfies 1 < a < num_to_be_factorized and gcd(a, num_to_be_factorized) = 1
# """
# qc = Shor().construct_circuit(num_to_be_factorized, a)
# qc.measure_all()
# qc.name = "shor"
# return qc
[docs]
def create_instance(choice: str) -> list[int]:
"""Returns the number to be factorized and the integer a for the Shor's algorithm."""
instances = {
"xsmall": [9, 4], # 18 qubits
"small": [15, 4], # 18 qubits
"medium": [821, 4], # 42 qubits
"large": [11777, 4], # 58 qubits
"xlarge": [201209, 4], # 74 qubits
}
return instances[choice]
[docs]
class Shor(QuantumCircuit):
"""Shor's algorithm implementation."""
def __init__(self, to_be_factored_number: int, a: int):
r"""
Args:
to_be_factored_number: The integer to be factored, denoted as :math:`N`.
a: An integer coprime to :math:`N`, denoted as :math:`a`.
"""
self.validate_input(to_be_factored_number, a)
self._N = to_be_factored_number
self._a = a
num_bits_necessary = to_be_factored_number.bit_length()
up_qreg = QuantumRegister(2 * num_bits_necessary, name="up")
down_qreg = QuantumRegister(num_bits_necessary, name="down")
aux_qreg = QuantumRegister(num_bits_necessary + 2, name="aux")
super().__init__(
up_qreg,
down_qreg,
aux_qreg,
name=f"Shor(N={to_be_factored_number}, a={a})",
)
self._build()
[docs]
@staticmethod
def get_angles(a: int, n: int) -> NDArray[np.float64]:
"""Calculates the array of angles to be used in the addition in Fourier Space."""
bits_little_endian = (f"{a:b}".zfill(n))[::-1]
angles = np.zeros(n)
for i in range(n):
for j in range(i + 1):
k = i - j
if bits_little_endian[j] == "1":
angles[i] += pow(2, -k)
return angles * np.pi
[docs]
@staticmethod
def phi_add_gate(angles: NDArray[np.float64] | ParameterVector) -> QuantumCircuit:
"""Gate that performs addition by a in Fourier Space."""
circuit = QuantumCircuit(len(angles), name="phi_add_a")
for i, angle in enumerate(angles):
circuit.p(angle, i)
return circuit
[docs]
def synth_qft_full(
self,
num_qubits: int,
do_swaps: bool = True,
approximation_degree: int = 0,
insert_barriers: bool = False,
inverse: bool = False,
name: str | None = None,
) -> QuantumCircuit:
"""Construct a circuit for the Quantum Fourier Transform using all-to-all connectivity.
.. note::
With the default value of ``do_swaps = True``, this synthesis algorithm creates a
circuit that faithfully implements the QFT operation. This circuit contains a sequence
of swap gates at the end, corresponding to reversing the order of its output qubits.
In some applications this reversal permutation can be avoided. Setting ``do_swaps = False``
creates a circuit without this reversal permutation, at the expense that this circuit
implements the "QFT-with-reversal" instead of QFT. Alternatively, the
:class:`~.ElidePermutations` transpiler pass is able to remove these swap gates.
Args:
num_qubits: The number of qubits on which the Quantum Fourier Transform acts.
do_swaps: Whether to synthesize the "QFT" or the "QFT-with-reversal" operation.
approximation_degree: The degree of approximation (0 for no approximation).
It is possible to implement the QFT approximately by ignoring
controlled-phase rotations with the angle beneath a threshold. This is discussed
in more detail in https://arxiv.org/abs/quant-ph/9601018 or
https://arxiv.org/abs/quant-ph/0403071.
insert_barriers: If ``True``, barriers are inserted for improved visualization.
inverse: If ``True``, the inverse Quantum Fourier Transform is constructed.
name: The name of the circuit.
Returns:
A circuit implementing the QFT operation.
"""
# _warn_if_precision_loss(num_qubits - approximation_degree - 1)
circuit = QuantumCircuit(num_qubits)
for j in reversed(range(num_qubits)):
circuit.h(j)
num_entanglements = max(0, j - max(0, approximation_degree - (num_qubits - j - 1)))
for k in reversed(range(j - num_entanglements, j)):
# Use negative exponents so that the angle safely underflows to zero, rather than
# using a temporary variable that overflows to infinity in the worst case.
lam = np.pi * (2.0 ** (k - j))
circuit.cp(lam, j, k)
if insert_barriers:
circuit.barrier()
if do_swaps:
for i in range(num_qubits // 2):
circuit.swap(i, num_qubits - i - 1)
if inverse:
circuit = circuit.inverse()
# It is important to set the name afte the circuit's generic "inverse" is called,
# since that will add ``_dg`` to the name.
if name is not None:
circuit.name = name
return circuit
[docs]
def double_controlled_phi_add_mod_n(
self,
angles: NDArray[np.float64] | ParameterVector,
c_phi_add_n: QuantumCircuit,
iphi_add_n: QuantumCircuit,
qft: QuantumCircuit,
iqft: QuantumCircuit,
) -> QuantumCircuit:
"""
Creates a circuit which implements double-controlled modular addition by a.
"""
ctrl_qreg = QuantumRegister(2, "ctrl")
b_qreg = QuantumRegister(len(angles), "b")
flag_qreg = QuantumRegister(1, "flag")
circuit = QuantumCircuit(ctrl_qreg, b_qreg, flag_qreg, name="ccphi_add_a_mod_N")
cc_phi_add_a = self.phi_add_gate(angles).control(2)
cc_iphi_add_a = cc_phi_add_a.inverse()
circuit.compose(cc_phi_add_a, [*ctrl_qreg, *b_qreg], inplace=True)
circuit.compose(iphi_add_n, b_qreg, inplace=True)
circuit.compose(iqft, b_qreg, inplace=True)
circuit.cx(b_qreg[-1], flag_qreg[0])
circuit.compose(qft, b_qreg, inplace=True)
circuit.compose(c_phi_add_n, [*flag_qreg, *b_qreg], inplace=True)
circuit.compose(cc_iphi_add_a, [*ctrl_qreg, *b_qreg], inplace=True)
circuit.compose(iqft, b_qreg, inplace=True)
circuit.x(b_qreg[-1])
circuit.cx(b_qreg[-1], flag_qreg[0])
circuit.x(b_qreg[-1])
circuit.compose(qft, b_qreg, inplace=True)
circuit.compose(cc_phi_add_a, [*ctrl_qreg, *b_qreg], inplace=True)
return circuit
[docs]
def controlled_multiple_mod_n(
self,
num_bits_necessary: int,
to_be_factored_number: int,
a: int,
c_phi_add_n: QuantumCircuit,
iphi_add_n: QuantumCircuit,
qft: QuantumCircuit,
iqft: QuantumCircuit,
) -> QuantumCircuit:
"""Implements modular multiplication by a as an instruction."""
ctrl_qreg = QuantumRegister(1, "ctrl")
x_qreg = QuantumRegister(num_bits_necessary, "x")
b_qreg = QuantumRegister(num_bits_necessary + 1, "b")
flag_qreg = QuantumRegister(1, "flag")
circuit = QuantumCircuit(ctrl_qreg, x_qreg, b_qreg, flag_qreg, name="cmult_a_mod_N")
angle_params = ParameterVector("angles", length=num_bits_necessary + 1)
modulo_adder = self.double_controlled_phi_add_mod_n(angle_params, c_phi_add_n, iphi_add_n, qft, iqft)
def append_adder(adder: QuantumCircuit, constant: int, idx: int) -> None:
partial_constant = (pow(2, idx, to_be_factored_number) * constant) % to_be_factored_number
angles = self.get_angles(partial_constant, num_bits_necessary + 1)
bound = adder.assign_parameters({angle_params: angles})
circuit.append(bound, [*ctrl_qreg, x_qreg[idx], *b_qreg, *flag_qreg])
circuit.compose(qft, b_qreg, inplace=True)
# perform controlled addition by a on the aux register in Fourier space
for i in range(num_bits_necessary):
append_adder(modulo_adder, a, i)
circuit.compose(iqft, b_qreg, inplace=True)
# perform controlled subtraction by a in Fourier space on both the aux and down register
for i in range(num_bits_necessary):
circuit.cswap(ctrl_qreg, x_qreg[i], b_qreg[i])
circuit.compose(qft, b_qreg, inplace=True)
a_inv = pow(a, -1, mod=to_be_factored_number)
modulo_adder_inv = modulo_adder.inverse()
for i in reversed(range(num_bits_necessary)):
append_adder(modulo_adder_inv, a_inv, i)
circuit.compose(iqft, b_qreg, inplace=True)
return circuit
[docs]
def power_mod_n(self, num_bits_necessary: int, to_be_factored_number: int, a: int) -> QuantumCircuit:
"""Implements modular exponentiation a^x as a quantum circuit."""
up_qreg = QuantumRegister(2 * num_bits_necessary, name="up")
down_qreg = QuantumRegister(num_bits_necessary, name="down")
aux_qreg = QuantumRegister(num_bits_necessary + 2, name="aux")
circuit = QuantumCircuit(up_qreg, down_qreg, aux_qreg, name=f"{a}^x mod {to_be_factored_number}")
qft = self.synth_qft_full(num_bits_necessary + 1, do_swaps=False)
iqft = qft.inverse()
# Create gates to perform addition/subtraction by N in Fourier Space
phi_add_n = self.phi_add_gate(self.get_angles(to_be_factored_number, num_bits_necessary + 1))
iphi_add_n = phi_add_n.inverse()
c_phi_add_n = phi_add_n.control(1)
# Apply the multiplication gates as showed in
# the report in order to create the exponentiation
for i in range(2 * num_bits_necessary):
partial_a = pow(self._a, pow(2, i), self._N)
modulo_multiplier = self.controlled_multiple_mod_n(
num_bits_necessary, self._N, partial_a, c_phi_add_n, iphi_add_n, qft, iqft
)
circuit.compose(modulo_multiplier, [up_qreg[i], *down_qreg, *aux_qreg], inplace=True)
return circuit
def _build(self):
"""Construct quantum part of the algorithm.
Arguments:
to_be_factored_number: The odd integer to be factored, has a min. value of 3.
a: Any integer that satisfies 1 < a < N and gcd(a, N) = 1.
Returns:
Quantum circuit.
"""
to_be_factored_number = self._N
a = self._a
self.validate_input(to_be_factored_number, a)
# Get n value used in Shor's algorithm, to know how many qubits are used
num_bits_necessary = to_be_factored_number.bit_length()
# quantum register where the sequential QFT is performed
up_qreg = QuantumRegister(2 * num_bits_necessary, name="up")
# quantum register where the multiplications are made
down_qreg = QuantumRegister(num_bits_necessary, name="down")
# auxiliary quantum register used in addition and multiplication
aux_qreg = QuantumRegister(num_bits_necessary + 2, name="aux")
# Create Quantum Circuit
circuit = QuantumCircuit(up_qreg, down_qreg, aux_qreg, name=f"Shor(N={to_be_factored_number}, a={a})")
# Create maximal superposition in top register
circuit.h(up_qreg)
# Initialize down register to 1
circuit.x(down_qreg[0])
# Apply modulo exponentiation
modulo_power = self.power_mod_n(num_bits_necessary, to_be_factored_number, a)
circuit.compose(modulo_power.decompose(reps=4), circuit.qubits, inplace=True)
# Apply inverse QFT
iqft = QFT(len(up_qreg)).inverse()
circuit.compose(iqft.to_gate(), up_qreg, inplace=True)
# Final composition into self
self.compose(circuit, self.qubits, inplace=True)