Large-Scale Simulation of Shor’s Quantum Factoring Algorithm

Abstract

Shor’s factoring algorithm is one of the most anticipated applications of quantum computing. However, the limited capabilities of today’s quantum computers only permit a study of Shor’s algorithm for very small numbers. Here, we show how large GPU-based supercomputers can be used to assess the performance of Shor’s algorithm for numbers that are out of reach for current and near-term quantum hardware. First, we study Shor’s original factoring algorithm. While theoretical bounds suggest success probabilities of only 3–4%, we find average success probabilities above 50%, due to a high frequency of “lucky” cases, defined as successful factorizations despite unmet sufficient conditions. Second, we investigate a powerful post-processing procedure, by which the success probability can be brought arbitrarily close to one, with only a single run of Shor’s quantum algorithm. Finally, we study the effectiveness of this post-processing procedure in the presence of typical errors in quantum processing hardware. We find that the quantum factoring algorithm exhibits a particular form of universality and resilience against the different types of errors. The largest semiprime that we have factored by executing Shor’s algorithm on a GPU-based supercomputer, without exploiting prior knowledge of the solution, is 549,755,813,701 = 712,321 × 771,781. We put forward the challenge of factoring, without oversimplification, a non-trivial semiprime larger than this number on any quantum computing device.

1. Introduction

The challenge of factoring integers is one of the oldest problems in mathematics [1,2]. Famous mathematicians such as Fermat, Euler, and Gauss have made substantial contributions to the problem, and even algorithms discovered by the ancient Greeks—the Euclidean algorithm and the sieve of Eratosthenes—are still in use today. The state-of-the-art algorithms are based on the general number field sieve [3] and have recently achieved the factorization of RSA-250 from the famous RSA factoring challenge [4]. Still, all known algorithms exhibit at best subexponential time and space complexity [4,5]. The difficulty of solving this type of problem using classical computers is an integral aspect of modern data and communication security [6,7].

In 1994, Peter Shor proposed an algorithm to factor integers on quantum computers with an exponential speedup [8,9,10] over the best known classical algorithms. Factoring an L-bit integer N with the conventional Shor algorithm [10] requires at least $3L$ qubits: $L=\lfloor {log}_{2}N\rfloor +1$ qubits to represent N, and $t=\lceil 2{log}_{2}N\rceil \approx 2L$ qubits for the Quantum Fourier Transform (QFT), plus $O\left(L\right)$ qubits for the modular exponentiation [6,11]. Kitaev, Griffiths, and Niu realized that by replacing the QFT with a semiclassical Fourier transform, only a single qubit can be reused t times to obtain the same result [12,13] (also known as qubit recycling [14,15] or dynamic quantum computing [16]). It is thus possible to run Shor’s algorithm with only $L+1$ qubits to factor L-bit integers (which is less than required by the best adiabatic algorithm [17,18]). We refer to this variant as the iterative Shor algorithm.

The iterative Shor algorithm has been executed on real quantum computing devices to factor 15, 21, and 35 [15,19,20], without relying on oversimplification [21]. Implementing the algorithm for integers beyond 35 continues to pose substantial experimental challenges [6,22].

To study the performance of Shor’s algorithm for much larger integers than those feasible for testing on real quantum devices, we have developed a massively parallel simulator called shorgpu [23], specifically designed to execute the iterative Shor algorithm on multiple GPUs. Using shorgpu, we have examined over 60,000 factoring scenarios for integers up to $N_{\max }$ = 549,755,813,701, significantly surpassing previous achievements using statevector simulators [24,25,26], matrix product states [27,28] (in [28], the authors simulated 60 qubits to factor N = 961,307), and tensor networks [29,30]. Note that $N_{\max }$ is still “small” for cryptographic purposes. In order to handle integers of the size of $N_{\max }$, shorgpu uses a new technique (see Section 2) to perform the distributed memory communications.

To factor $N_{\max }$, the conventional Shor algorithm would need 117 qubits. The iterative Shor algorithm, however, needs only 40 qubits. It is important to note that the resulting quantum algorithm is still an honest implementation of Shor’s algorithm: it produces the same results, does not require exponentially large classical resources (given a large enough quantum computer) and, most importantly, does not exploit a priori knowledge of the factors [21].

We emphasize that for all our simulations, we do not require the solution of the factoring problem to be known. If one presumes knowledge of the solution, and one is not interested in simulating the effect of quantum errors, it is possible to study even larger, cryptographically relevant cases using Qunundrum [31].

The procedure used to factor integers is shown in Figure 1: first, a factoring problem is selected, consisting of a semiprime $N=p\times q$ to factor and a random integer $1<a<N$ comprime to N (i.e., $\gcd (a,N)=1$). Using this as input, shorgpu executes the iterative Shor algorithm with $n=L+1$ qubits to produce several bitstrings. Each bitstring j is processed using either Shor’s [8,9,10,32] or Ekerå’s [33,34] post-processing method, which may or may not produce a factor of N (see the yellow section in Figure 1). An important step on the way is to extract a candidate r for the order $\widehat{r}={\mathrm{ord}}_N\left(a\right)$. Here, ${\mathrm{ord}}_N\left(a\right)$ denotes the order of a modulo N, defined as the smallest exponent $\widehat{r}>0$ such that $a^{\widehat{r}}\mathrm{mod}N=1$.

Note that “success” is not guaranteed by Shor’s algorithm; in particular, the sampled bitstring might produce an $r\ne \widehat{r}$ that is not the order, or r might be odd, in which case Shor’s post-processing method is not guaranteed to work. However, if the blind application of Shor’s factoring procedure still yields at least one factor, we count this execution as a “lucky” case. As shown below, a “lucky” factor is found much more often than expected.

In principle, the green section in Figure 1 representing shorgpu can be completely replaced by a sufficiently large, error-corrected quantum computing device. With this in mind, we put forward the challenge of indirect quantum supremacy [35] (a.k.a. limited quantum speedup [36]) for a future quantum computer. Here, “indirect” means that the simulator (running on a conventional computer) is required to simulate an ideal quantum computational model that executes the same quantum algorithm as the quantum computer, without using any prior knowledge of the solution. More specifically, the challenge for a gate-based quantum computer would be to factor, using Shor’s algorithm without oversimplification [21], an “interesting” semiprime—where “interesting” means that the two distinct prime factors shall have the same number of digits—that is larger than the largest semiprime that can be factored by the quantum computer simulator.

1.1. Related Work

There is a large body of literature on Shor’s quantum factoring algorithm; the related work can be roughly classified into five main categories. In this section, we give a survey of their main goals and discuss several individual results.

1

Theory: the first class of articles focuses on theoretical perspectives such as algorithmic modifications and improved lower bounds on the success probability [11,14,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60], many of which consider the case that some parameters of Shor’s algorithm are modified and certain trade-offs are made. This class contains work that estimates the number of resources required when using different levels of quantum computer technology [6,22]. This line of work culminates in Ekerå’s post-processing algorithms [34], by which the success probability for a single run of the quantum part can be brought arbitrarily close to one (see below).

2

Simulation: second, Shor’s algorithm has been studied by using simulators running on conventional computers. Some use universal quantum computer simulators [24,25,26,61], sometimes also called Schrödinger simulators, since they propagate the full quantum statevector. Another approach is to use so-called Feynman simulators, which can only access certain amplitudes from the full statevector, but may require fewer computational resources; they are often based on tensor networks or matrix product states [27,28,29,30]. Finally, there is software designed to directly sample from the probability distributions generated by Shor’s algorithm (cf. Equation (26) below) and various extensions thereof. This class contains the suite of programs called Qunundrum [31,62], which can simulate distributions for large, cryptographically relevant cases. Note, however, that the solution to the factoring problem (i.e., the order or the discrete logarithm) must be known in advance, and the effect of errors in the quantum part cannot be simulated.

3

Alternative Algorithms: a third line of work studies alternative ways to use gate-based quantum computers to solve the factoring problem. Some of them use Shor’s discrete logarithm quantum algorithm [62,63,64], which is also an instance of the hidden subgroup problem [65]. In the Ekerå-Håstad scheme [63], the idea to factor a semiprime $N=pq$ is to pick a random $g\in {\mathbb{Z}}_N^{\ast }$, compute $y=g^{N+1}\mathrm{mod}N$ with unknown order r, and then obtain $d\equiv {log}_{g}y\equiv pq+1\equiv pq+1-\varphi \left(N\right)\equiv p+q\left(\mathrm{mod}r\right)$ (using that $r\mid \varphi \left(N\right)=(p-1)(q-1)$ [66]). If $r>p+q$ (which is the case for many g), we have $d=p+q$, and additionally knowing $N=pq$ allows one to compute p and q. Another alternative way to solve the factoring problem is given in [67] and is based on the classical number field sieve [3]. In particular, Bernstein et al. propose to use Grover’s quantum search algorithm [68] (and/or Shor’s algorithm for a much smaller subproblem) to accelerate the number field sieve. This proposal is asymptotically worse in time complexity than using Shor’s algorithm directly, but it requires fewer qubits and is therefore possibly easier to realize in near-term physical devices. Finally, Li et. al. [69] presented an algorithm with an exponential speedup (beyond the framework of the hidden subgroup problem [65]) that solves the square-free decomposition problem—a problem related to factoring in which the task is to find, for any integer $N>0$, the unique integers $N_r$ and $N_s^2$ of the square-free decomposition $N=N_r{N}_s^2$.

4

Gate-Based Experiments: fourth, there have been several experimental efforts to implement Shor’s factoring algorithm on existing gate-based quantum computer devices [15,19,20,70,71,72,73,74,75,76]. However, many of these have made use of prior knowledge about the factors to simplify the experimental setup [21]. In the extreme case (namely when a base $a\in {\mathbb{Z}}_N^{\ast }$ with order ${\mathrm{ord}}_N\left(a\right)=2$ is used), this reduces the computational problem to the equivalent of flipping coins. Experiments that have not used such an oversimplified method can be found in [15,19,20].

5

Other Experiments: finally, quantum annealers and adiabatic quantum computers have been used to study alternative factoring algorithms [17,77,78,79,80,81,82]. The quantum annealing approach requires at most $O\left(L^2\right)$ qubits to factor an L-bit number. Quantum annealing and adiabatic quantum computation are technologically significantly ahead of gate-based quantum computing, in that larger quantum processing units with more than 5000 qubits exist and that they can solve much larger problems [83,84]. In particular, numbers up to and above 200,000 have been factored on the D-Wave 2000Q [78,79] and 1,005,973 has been factored using D-Wave hybrid [80]. Although significantly larger than the numbers factored on gate-based quantum computers (without oversimplification), these numbers are still much smaller than $N_{\max }$ = 549,755,813,701 factored in this work using shorgpu.

1.2. Outline

This paper is structured as follows. In Section 2, we describe the algorithmic details of shorgpu. In particular, we explain how to implement the modular multiplication as a systematic communication scheme between the compute nodes. In Section 3, we present our results from over 60,000 quantum computer simulations using up to 2048 GPUs. Section 4 contains our conclusions.

2. Materials and Methods

For almost all results reported in this work, we use shorgpu to simulate the iterative Shor algorithm (the source code is available online [23]). It propagates an $n=L+1$-qubit statevector $|\psi \rangle$ consisting of $2^{L+1}$ complex numbers through the quantum circuit for factoring an L-bit number shown in Figure 2. Each step in the quantum circuit corresponds to an operation on $|\psi \rangle$. In this section, we describe each of these operations in a linear algebraic context. The probability distribution generated by Shor’s algorithm is derived and visualized in Appendix A.

As the total memory is the bottleneck of such statevector simulations, the $2^{L+1}$ complex numbers $|\psi \rangle$ are distributed over the memory of up to 2048 GPUs (cf. Figure 1). Communication between the GPUs is managed using the Message Passing Interface (MPI) [85].

We use the Jülich Universal Quantum Computer Simulator (JUQCS) [24,25,86] for verification. JUQCS was previously used to simulate the conventional Shor algorithm for $N\le$ 65,531 with up to $n=48$ qubits on the Sunway TaihuLight and the K computer [25]. A few features had to be added to JUQCS to be able to also simulate the iterative Shor algorithm. The latter made it possible to simulate one bitstring for N = 4,194,293 with $n=23$ qubits in about 720 s (using four A100 GPUs). However, our JUQCS implementation of the oracle which performs the modular exponentiation becomes highly inefficient as the number n of qubits increases because it does not distribute well over many cores or GPUs. In contrast, the new, dedicated algorithm described below generates a bitstring in about 0.4 s for the same problem and the same number of GPUs. For the largest problem simulated (N = 549,755,813,701 with $n=40$ qubits), shorgpu generates a bitstring in about 200 s using 2048 GPUs. We verified that the iterative Shor algorithm simulated with shorgpu produces the same results as JUQCS for problems of the size that can be simulated with JUQCS.

2.1. Initialization

To simulate the iterative Shor algorithm for factoring an L-bit semiprime N, shorgpu simulates the full quantum circuit with $n=L+1$ qubits shown in Figure 2. This is conducted by computing all complex coefficients of the statevector

$$\begin{array}{c}\hfill |\psi \rangle =\sum _{k_L\cdots k_0=0,1}{\psi }_{k_L\cdots k_0}|k_L\cdots k_0\rangle =\left(\begin{array}{c}{\psi }_{0\cdots 00}\\ {\psi }_{0\cdots 01}\\ ⋮\\ {\psi }_{1\cdots 11}\end{array}\right).\end{array}$$

(1)

These $2^{L+1}$ complex double-precision numbers ${\psi }_{k_L\cdots k_0}$ are distributed over $N_{\mathrm{GPU}}\in \{2,4,8,\dots ,2048\}$ GPUs. The distributed memory communication between the GPUs uses CUDA-aware MPI. In our approach, each GPU is identified by its MPI rank, i.e., an $n_{\mathrm{global}}$-bit integer called $\mathtt{mpi}\_\mathtt{rank}=0,\dots ,N_{\mathrm{GPU}}-1$. Here, $n_{\mathrm{global}}$ denotes the number of so-called global qubits (see [24,25,86]). We have $N_{\mathrm{GPU}}=2^{n_{\mathrm{global}}}$ GPUs. The other $n_{\mathrm{local}}=n-n_{\mathrm{global}}$ qubits are called local qubits, since each GPU holds in its local memory all $2^{n_{\mathrm{local}}}$ complex coefficients

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\psi }_{\mathrm{bin}(\mathtt{mpi}\_\mathtt{rank})0\cdots 00}\\ {\psi }_{\mathrm{bin}(\mathtt{mpi}\_\mathtt{rank})0\cdots 01}\\ ⋮\\ {\psi }_{\mathrm{bin}(\mathtt{mpi}\_\mathtt{rank})1\cdots 11}\end{array}\right).\end{array}$$

(2)

The GPUs (i.e., the MPI processes) are further divided into two separate groups, identified by the most significant bit of the MPI rank,

$$\begin{array}{c}\hfill \mathrm{bin}(\mathtt{mpi}\_\mathtt{rank})=\mathtt{mpi}\_\mathtt{x}\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank}).\end{array}$$

(3)

Here, $\mathtt{mpi}\_\mathtt{x}=0,1$ identifies the group and $0\le \mathtt{mpi}\_\mathtt{xrank}<N_{\mathrm{GPU}}/2$ identifies the GPU within each group. Thus, shorgpu requires at least two GPUs to work (unless a single GPU is used with 2 MPI processes in overscheduling mode). The reason for the separation into two groups is that the implementation of the controlled modular multiplication gate (see below) requires an all-to-all communication between all GPUs with $\mathtt{mpi}\_\mathtt{x}=1$.

At the start of the simulation, the statevector $|\psi \rangle$ is initialized in the state $|+\rangle |0\cdots 01\rangle$, where $|+\rangle =(|0\rangle +|1\rangle )/\sqrt{2}$. This means that we set

$$\begin{array}{cc}\hfill {\psi }_{00\cdots 01}& =\frac{1}{\sqrt{2}},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\psi }_{10\cdots 01}& =\frac{1}{\sqrt{2}},\hfill \end{array}$$

(5)

and all other coefficients to zero. This type of initialization is always used unless shorgpu is used to assess the effect of quantum initialization errors (for information on this mode, see Section 2.6 below).

2.2. Controlled Modular Multiplication Gate

The first gate in each of the t stages in Figure 2 is the controlled modular multiplication gate (also called oracle gate), controlled by the first qubit. Mathematically, its operation is defined by

$$\begin{array}{c}\hfill \mathrm{C}{U}_{\mathtt{a}}|x\rangle |y\rangle =\left\{\begin{array}{cc}|0\rangle |y\rangle \hfill & (x=0)\hfill \\ |1\rangle |\mathtt{a}y\mathrm{mod}N\rangle \hfill & (x=1\mathrm{and}0\le y<N)\hfill \\ |1\rangle |y\rangle \hfill & (x=1\mathrm{and}N\le y)\hfill \end{array}\right.,\end{array}$$

(6)

where x denotes the first qubit, $y=0,\dots ,N_{\mathrm{GPU}}/2$ denotes the other qubits, and $\mathtt{a}\in \{a^{2^{t-1}}\mathrm{mod}N,a^{2^{t-2}}\mathrm{mod}N,\dots ,a\}$ stands for one of the powers of a in Figure 2. Note that each individual modular exponentiation is always precomputed for any realization of this circuit (we use the shift-and-multiply algorithm). This is independent of whether the circuit is executed by a quantum computer simulator or a real quantum computer (see also [6,19,20]).

Looking at Equation (6), we see that the oracle gate performs a permutation of all complex coefficients among the GPUs in the $\mathtt{mpi}\_\mathtt{x}=1$ group. shorgpu implements this unitary operation by computing, on each GPU, all indices of the coefficients that are sent to other GPUs (stored in a GPU buffer oracle_idxsend) and those that are received from other GPUs (stored in a GPU buffer oracle_idxrecv) using the precomputed modular inverse $\mathtt{ainv}={\mathtt{a}}^{-1}\mathrm{mod}N$, which is efficiently computable using the extended Euclidean algorithm. The MPI communication scheme for an example with 16 GPUs is shown in Figure 3.

The complexity of the permutation depends on the value of a. For instance, in the special case that the order of a is a small power of 2, we have $\mathtt{a}=1$ in many of the early stages, so the oracle gate would not require MPI communication between different GPUs. In general the communication scheme can be very complicated. Figure 3 shows a typical instance where each GPU sends (red arrows) and receives (blue arrows) some coefficients from other GPUs.

To implement this communication scheme between the GPUs, shorgpu uses non-blocking point-to-point communication in a circular fashion using MPI_Isend and MPI_Irecv. Additionally, before the send operations, each GPU first arranges all coefficients that are sent to a particular other GPU in a contiguous block of memory, schematically denoted by ${\psi }^{\left(\mathrm{cont}\right)}$. This is imperative since for large N, this part of the simulation takes a significant fraction of the total run time. Alternative implementations using one-sided communication such as MPI_Put, collective communication using MPI_Alltoallv, or communication based on custom MPI data types (see [85] for more information) performed significantly worse in our experiments.

2.3. Rotation Gate

After the oracle gate, each stage (except the first stage) of the quantum circuit in Figure 2 contains a sequence of rotation gates defined by

$$\begin{array}{c}\hfill R_l=\left(\begin{array}{cc}1& 0\\ 0& e^{2\pi i/2^l}\end{array}\right).\end{array}$$

(7)

These are controlled by bits resulting from previous measurements. Specifically, at the stage $\mathtt{cbit}=0,\dots ,t-1$, in which the classical bit $j_{\mathtt{cbit}}$ is being measured, the sequence of these controlled rotation gates reads

$$\begin{array}{c}\hfill \prod _{l=2}^{1+\mathtt{cbit}}\mathrm{C}{R}_l=\prod _{l=2}^{1+\mathtt{cbit}}\left(\begin{array}{cc}1& 0\\ 0& e^{2\pi i{j}_{1+\mathtt{cbit}-l}/2^l}\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& e^{i{\phi }_{\mathtt{cbit}}}\end{array}\right),\end{array}$$

(8)

where the phase ${\phi }_{\mathtt{cbit}}$ at stage $\mathtt{cbit}$ amounts to

$$\begin{array}{c}\hfill {\phi }_{\mathtt{cbit}}=2\pi \sum _{l=2}^{1+\mathtt{cbit}}\frac{j_{1+\mathtt{cbit}-l}}{2^l}=\frac{\pi j^{\left(\mathtt{cbit}\right)}}{2^{\mathtt{cbit}}},\end{array}$$

(9)

and $j^{\left(\mathtt{cbit}\right)}=j_{\mathtt{cbit}-1}{j}_{\mathtt{cbit}-2}\cdots j_1{j}_0$ is the integer assembled from all classical bits measured up to this point.

As the phase gate given by Equation (8) only affects coefficients ${\psi }_{k_L\cdots k_0}$ where the first qubit index $k_L=1$, this operation only needs to be implemented by the GPUs in the $\mathtt{mpi}\_\mathtt{x}=1$ group. This is performed directly after the implementation of the oracle gate, when moving all coefficients out of the contiguous memory blocks ${\psi }^{\left(\mathrm{cont}\right)}$, according to

$$\begin{array}{c}\mathrm{Re}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }\right)\hfill \\ \leftarrow cos\left({\phi }_{\mathtt{cbit}}\right)\mathrm{Re}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }^{\left(\mathrm{cont}\right)}\right)\hfill \\ -sin\left({\phi }_{\mathtt{cbit}}\right)\mathrm{Im}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }^{\left(\mathrm{cont}\right)}\right),\hfill \end{array}$$

(10)

$$\begin{array}{c}\mathrm{Im}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }\right)\hfill \\ \leftarrow cos\left({\phi }_{\mathtt{cbit}}\right)\mathrm{Im}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }^{\left(\mathrm{cont}\right)}\right)\hfill \\ +sin\left({\phi }_{\mathtt{cbit}}\right)\mathrm{Re}\left({\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }^{\left(\mathrm{cont}\right)}\right).\hfill \end{array}$$

(11)

2.4. Hadamard Gate

The implementation of the Hadamard gate on the first qubit transforms the statevector coefficients as

$$\begin{array}{cc}\hfill & {\psi }_{0\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }\hfill \\ \hfill & \leftarrow \frac{{\psi }_{0\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }+{\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }}{\sqrt{2}},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & {\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }\hfill \\ \hfill & \leftarrow \frac{{\psi }_{0\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }-{\psi }_{1\mathrm{bin}(\mathtt{mpi}\_\mathtt{xrank})\ast \cdots \ast }}{\sqrt{2}}.\hfill \end{array}$$

(13)

For every GPU in the $\mathtt{mpi}\_\mathtt{x}=0$ group, this requires two-sided MPI communication with exactly one GPU in the $\mathtt{mpi}\_\mathtt{x}=1$ group.

2.5. Measurement Operation

At the end of each stage in Figure 2, the classical bit $j_{\mathtt{cbit}}$ is measured, where $\mathtt{cbit}=0,\dots ,t-1$ enumerates the stage. This amounts to adding up the probabilities

$$\begin{array}{c}\hfill p_1=\sum _{k_{L-1}\cdots k_0=0,1}{\left|{\psi }_{1k_{L-1}\cdots k_0}\right|}^2,\end{array}$$

(14)

which is an MPI reduction over all GPUs belonging to the $\mathtt{mpi}\_\mathtt{x}=1$ group. The probability to measure 1 (0) is then given by $p_1$ ($p_0=1-p_1$). This probability can be used to sample $j_{\mathtt{cbit}}$, which is performed by drawing a uniform random number $R\in [0,1)$, and assigning $j_{\mathtt{cbit}}=1$ if this $R<p_1$ and $j_{\mathtt{cbit}}=0$ otherwise.

2.6. Reset Operation

The reset operation performs both the von-Neumann projection of the statevector to the result of the measurement and the reinitialization of the first qubit in $|+\rangle$ at the same time. If the result of the measurement is given by $j_{\mathtt{cbit}}=0,1$, this operation is performed by transforming all coefficients according to

$$\begin{array}{c}\hfill \left(\begin{array}{c}{\psi }_{00\cdots 0}\\ ⋮\\ {\psi }_{01\cdots 1}\\ {\psi }_{10\cdots 0}\\ ⋮\\ {\psi }_{11\cdots 1}\end{array}\right)\leftarrow \left(\begin{array}{c}{\psi }_{j_{\mathtt{cbit}}0\cdots 0}/\sqrt{2p_{j_{\mathtt{cbit}}}}\\ ⋮\\ {\psi }_{j_{\mathtt{cbit}}1\cdots 1}/\sqrt{2p_{j_{\mathtt{cbit}}}}\\ {\psi }_{j_{\mathtt{cbit}}0\cdots 0}/\sqrt{2p_{j_{\mathtt{cbit}}}}\\ ⋮\\ {\psi }_{j_{\mathtt{cbit}}1\cdots 1}/\sqrt{2p_{j_{\mathtt{cbit}}}}\end{array}\right).\end{array}$$

(15)

Of course, in the case of quantum errors, the coefficients have to be replaced accordingly (cf. Equations (25a) and (25b)).

This operation requires an MPI transfer of all coefficients from the GPUs in the group $\mathtt{mpi}\_\mathtt{x}=j_{\mathtt{cbit}}$ to the GPUs in the group $\mathtt{mpi}\_\mathtt{x}=1-j_{\mathtt{cbit}}$.

2.7. Initialization Errors

There are two different types of initialization errors that shorgpu can simulate, namely an amplitude initialization error and a phase initialization error. In both cases, a slightly different initial state $|{+}^{\prime }\rangle$ is used instead of $|+\rangle$ for the first qubit in all stages $\mathtt{cbit}=0,\dots ,t-1$ of the circuit in Figure 2. The slightly erroneous state $|{+}^{\prime }\rangle$ is parameterized in terms of an error parameter $\delta \in [0,1]$. Our motivation to prioritize the recycled qubit for a study of initialization errors instead of the other “internal” qubits is that this qubit is measured and reinitialized successively in every stage of the iterative Shor algorithm.

2.7.1. Amplitude Initialization Error

We define an amplitude initialization error as the case in which, at the beginning of each stage in Figure 2, the quantum state is not initialized in the equal superposition $|+\rangle$ but the slightly unequal superposition

$$\begin{array}{c}\hfill |{+}_{\mathrm{ampl}}^{\prime }\left(\delta \right)\rangle =\sqrt{\frac{1+\delta }{2}}|0\rangle +\sqrt{\frac{1-\delta }{2}}|1\rangle .\end{array}$$

(16)

This expression is motivated by the observation that quantum computer prototypes from the NISQ era sometimes tend to prefer $|0\rangle$ over $|1\rangle$ when brought to a uniform superposition by multiple quantum gates [87]. Furthermore, one of the most prominent decoherence and noise processes in qubit systems is a decay from $|1\rangle$ to $|0\rangle$, a so-called $T_1$ relaxation process [88,89,90].

2.7.2. Phase Initialization Error

As a second type of initialization error, we consider a phase initialization error defined as

$$\begin{array}{c}\hfill |{+}_{\mathrm{phase}}^{\prime }\left(\delta \right)\rangle =\frac{1}{\sqrt{2}}|0\rangle +\frac{e^{i\pi \delta }}{\sqrt{2}}|1\rangle .\end{array}$$

(17)

This expression is motivated by the fact that in addition to the $T_1$ relaxation process, dephasing processes are other prominent consequences of decoherence and noise in quantum systems [88,89,91].

2.7.3. Effective Single-Qubit Error Probability

For both initialization errors, the error parameter $\delta \in [0,1]$ can be related to an effective, single-qubit error probability, defined as the probability that the erroneous state would correctly be observed as a $|+\rangle$ state when measured along the x axis:

$$\begin{array}{cc}\hfill p_{\mathrm{ampl}}^{\mathrm{error}}\left(\delta \right)& =1-{|\langle +|{+}_{\mathrm{ampl}}^{\prime }\left(\delta \right)\rangle |}^2=\frac{1-\sqrt{1-{\delta }^2}}{2},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill p_{\mathrm{phase}}^{\mathrm{error}}\left(\delta \right)& =1-{|\langle +|{+}_{\mathrm{phase}}^{\prime }\left(\delta \right)\rangle |}^2=\frac{1-cos\left(\pi \delta \right)}{2}.\hfill \end{array}$$

(19)

Note, however, that this interpretation is not unique; depending on the particular realization of the quantum circuit, there may be more reasonable, alternative interpretations of $\delta$ in relation to an effective error probability.

2.8. Measurement Errors

For quantum processors, a measurement is often a slow and susceptible process by which destructive influences from the environment can enter the quantum system [92,93,94,95,96]. Moreover, it is particularly challenging to implement quantum non-demolition readout required for midcircuit measurements [16,97]. We distinguish between two different types of measurement errors, namely a classical error corresponding to a misclassification of the quantum measurement result, and a quantum error that may occur during or before each measurement.

2.8.1. Classical Measurement Error

We define a classical measurement error as a misclassification that occurs right after the quantum measurement process with a given, constant error probability $\delta$. It is defined by flipping only the resulting bit $j_{\mathtt{cbit}}$, while leaving the internal quantum state unchanged.

Simulating a classical measurement error requires a second sampling step, by drawing another uniform random number $R_2\in [0,1)$ and flipping the bit if $R_2<\delta$. In case of a misclassification error, we simply use $j_{\mathtt{cbit}}^{\left(\mathrm{observed}\right)}=1-j_{\mathtt{cbit}}$ for the classical bitstring in Figure 2. The quantum state, however, is left in its original state with the first qubit projected on $|j_{\mathtt{cbit}}\rangle$.

Note that even such a single misclassification error can have non-trivial consequences, since this error affects the angles of all subsequent rotation gates (see Figure 2). This has an influence on the measurements of the following bits $\mathtt{cbit}+1,\dots ,t-1$. Therefore, a single bit flip error can induce a change in more than one classical bit of the output bitstring j.

2.8.2. Quantum Measurement Error

Quantum errors are conventionally modeled as operations $\rho \mapsto \mathcal{E}\left(\rho \right)$ on the system’s density matrix $\rho =|\psi \rangle \langle \psi |$. If such an operation is a completely positive, trace-preserving map, it is called a quantum channel or error channel (see [32,98,99] for more information).

We model a quantum measurement error by applying a depolarizing error channel in every measurement process (which, on quantum computer hardware, is a time evolution that can take a significant amount of time [96]). The depolarizing error channel is defined by the quantum operation

$$\begin{array}{cc}{\mathcal{E}}_{\mathrm{dep}}\left(\tilde{\rho }\right)\hfill & =(1-p_x-p_y-p_z)\rho \hfill \\ \hfill & +p_x{\sigma }^x\tilde{\rho }{\sigma }^x+p_y{\sigma }^y\tilde{\rho }{\sigma }^y+p_z{\sigma }^z\tilde{\rho }{\sigma }^z,\hfill \end{array}$$

(20)

where $\tilde{\rho }$ is a single-qubit density matrix, $({\sigma }^x,{\sigma }^y,{\sigma }^z)$ are the Pauli matrices, and $(p_x,p_y,p_z)$ represent the error probabilities for the respective Pauli errors. After the application of ${\mathcal{E}}_{\mathrm{dep}}$ to the first qubit, the density matrix that describes the state of the full quantum computer reads

$$\begin{array}{cc}\rho \hfill & =({\mathcal{E}}_{\mathrm{dep}}\otimes I)\left(|\psi \rangle \langle \psi |\right)\hfill \\ \hfill & =\sum _{\overline{k}{\overline{k}}^{\prime }}{\mathcal{E}}_{\mathrm{dep}}({\psi }_{0\overline{k}}{\psi }_{0{\overline{k}}^{\prime }}^{\ast }|0\rangle \langle 0|+{\psi }_{0\overline{k}}{\psi }_{1{\overline{k}}^{\prime }}^{\ast }|0\rangle \langle 1|\hfill \\ \hfill & +{\psi }_{1\overline{k}}{\psi }_{0{\overline{k}}^{\prime }}^{\ast }|1\rangle \langle 0|+{\psi }_{1\overline{k}}{\psi }_{1{\overline{k}}^{\prime }}^{\ast }|1\rangle \langle 1|)\otimes |\overline{k}\rangle \langle {\overline{k}}^{\prime }|,\hfill \end{array}$$

(21)

where I is an identity operation on the remaining L qubits, and $\overline{k},{\overline{k}}^{\prime }=k_{L-1}\cdots k_0$ enumerate their $2^L$ different indices.

A measurement of the first qubit is quantum-mechanically described by the measurement operators ${\mathcal{M}}_0=|0\rangle \langle 0|\otimes I$ and ${\mathcal{M}}_1=|1\rangle \langle 1|\otimes I$. Using Equation (20), we find the probability to measure $j_{\mathtt{cbit}}=0,1$ as

$$\begin{array}{cc}{p}_{j_{\mathtt{cbit}}}^{\prime }\hfill & =\mathrm{Tr}{\mathcal{M}}_{j_{\mathtt{cbit}}}\rho {\mathcal{M}}_{j_{\mathtt{cbit}}}^{†}\hfill \\ \hfill & =(1-p_x-p_y)p_{j_{\mathtt{cbit}}}+(p_x+p_y)p_{1-j_{\mathtt{cbit}}},\hfill \end{array}$$

(22)

where $p_{j_{\mathtt{cbit}}}={\sum }_{\overline{k}}{|{\psi }_{j_{\mathtt{cbit}}\overline{k}}|}^2$ is computed according to Equation (14). A calculation of the post-measurement state ${\rho }_{j_{\mathtt{cbit}}}^{\prime }$ yields

$$\begin{array}{cc}{\rho }_{j_{\mathtt{cbit}}}^{\prime }\hfill & =\frac{M_{j_{\mathtt{cbit}}}\rho {\mathcal{M}}_{j_{\mathtt{cbit}}}^{†}}{\mathrm{Tr}{M}_{j_{\mathtt{cbit}}}\rho {\mathcal{M}}_{j_{\mathtt{cbit}}}^{†}}\hfill \\ \hfill & =p_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}\rangle \langle {\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}|\hfill \\ \hfill & +p_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}\rangle \langle {\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}|,\hfill \end{array}$$

(23)

where

$$\begin{array}{cc}\hfill p_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}& =\frac{(1-p_x-p_y)p_{j_{\mathtt{cbit}}}}{(1-p_x-p_y)p_{j_{\mathtt{cbit}}}+(p_x+p_y)p_{1-j_{\mathtt{cbit}}}},\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill p_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}& =\frac{(p_x+p_y)p_{1-j_{\mathtt{cbit}}}}{(1-p_x-p_y)p_{j_{\mathtt{cbit}}}+(p_x+p_y)p_{1-j_{\mathtt{cbit}}}},\hfill \end{array}$$

(24b)

and

$$\begin{array}{cc}\hfill |{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}\rangle & =\sum _{\overline{k}}\frac{{\psi }_{j_{\mathtt{cbit}}\overline{k}}}{\sqrt{p_{j_{\mathtt{cbit}}}}}|j_{\mathtt{cbit}}\overline{k}\rangle ,\hfill \end{array}$$

(25a)

$$\begin{array}{cc}\hfill |{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}\rangle & =\sum _{\overline{k}}\frac{{\psi }_{(1-j_{\mathtt{cbit}})\overline{k}}}{\sqrt{p_{1-j_{\mathtt{cbit}}}}}|j_{\mathtt{cbit}}\overline{k}\rangle .\hfill \end{array}$$

(25b)

Here, the superscript “correct” (“error”) refers to the probability and the state in the case that no error (an error) has occurred. Furthermore, the expressions show that both Pauli x and y errors only occur in combination, so we define the joint quantum error probability $\delta =p_x+p_y$, by analogy with the classical case.

As in the classical case, a simulation of the quantum error process requires two sampling operations: first, a random number $R\in [0,1)$ is sampled to assign the measurement result with probability $p_{j_{\mathtt{cbit}}}^{\prime }$ given by Equation (22), i.e., we assign $j_{\mathtt{cbit}}=1$ if this $R<p_1^{\prime }$ and $j_{\mathtt{cbit}}=0$ otherwise.

Second, a random number $R_2\in [0,1)$ is sampled to determine whether an error has happened or not. If $R_2<p_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}$, an error has happened while measuring $j_{\mathtt{cbit}}$, and the simulation continues with the state $|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}\rangle$ given by Equation (25b). Otherwise, the simulation continues with the state $|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}\rangle$.

Note that the quantum error has a more direct influence on the quantum state than the classical error, since the projection to either $|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{correct}\right)}\rangle$ or $|{\psi }_{j_{\mathtt{cbit}}}^{\left(\mathrm{error}\right)}\rangle$ directly affects the quantum state, not only implicitly through the angles of subsequent rotation gates.

2.9. Details on Memory and Computing Time

The largest part of the memory needed by shorgpu is taken by the coefficients of the statevector $|\psi \rangle$ (see Equation (1)). For a 40-qubit iterative Shor circuit, which can be used to factor 39-bit integers, the statevector needs $2^{40}$ complex double-precision floating point numbers, so $16\times 2^{40}\mathrm{B}=16\mathrm{TiB}$. For performance reasons, two statevector buffers are used in the implementation of the oracle gate and the following single-qubit gates. In addition to the two statevector buffers, shorgpu requires two 32-bit integer buffers for the implementation of the oracle gate, called oracle_idxrecv and oracle_idxsend (see above). Each of these takes another $4\times 2^{40}\mathrm{B}=4\mathrm{TiB}$. The total GPU memory required is thus slightly larger than $40\mathrm{TiB}$. When using $N_{\mathrm{GPU}}=2048$ GPUs, the required memory per GPU is slightly larger than $20\mathrm{GiB}$.

We performed all simulations on JUWELS Booster [100,101], a GPU cluster with 3744 NVIDIA A100 Tensor Core GPUs [102], each of which has $40\mathrm{GiB}$ of GPU memory. Note that the implementation of the algorithm requires the number of GPUs to be a power of two (cf. Section 2.1), so the maximum number of NVIDIA A100 GPUs that we can use on JUWELS Booster is 2048. The total computing time used to perform the simulations amounts to 594 core years (corresponding to 49.5 GPU years since each node contains four A100 GPUs and 48 physical CPU cores). We note that the total computing time is $22\%$ of the 2700 core years used for the recent factoring record of RSA-250—a number with 829 binary digits from the famous RSA factoring challenge [4].

3. Results

In this section, we describe and interpret the results obtained from simulating Shor’s algorithm according to Figure 1. For our analysis, we generated 61,362 factoring problems $(N,a)$, 52,077 of which were chosen to have uniformly distributed prime factors to ensure unbiased results, and the rest comprise individual factoring problems for large semiprimes (see Appendix B). We consider Shor’s original post-processing in Section 3.1 and Ekerå’s post-processing in Section 3.2.

3.1. Using Shor’s Post-Processing

A given factoring problem for Shor’s algorithm consists of a semiprime N and a random integer $1<a<N$ coprime to N. For each such factoring problem $(N,a)$, Shor’s algorithm produces a sample of M bitstrings (we typically consider $M=1024$ samples). Each bitstring j is analyzed using the so-called standard procedure (see Appendix C). If all checks on j from the standard procedure pass, the algorithm was successful and we count the bitstring j as “success”. However, if certain checks on j fail, we still evaluate j and test if it yields a factor of N. If it does, we count this factoring attempt as “lucky”. Figure 4a shows a scatter plot of all “success” and “success+lucky” probabilities for all uniformly distributed problems ($n<30$ qubits) and the individual large problems ($30\le n\le 40$ qubits).

Surprisingly, “lucky” occurs much more often than expected. In Figure 4b, we see that on average only $25\%$ of all bitstrings yield “success”. Including the “lucky” cases, however, a factor of the semiprime N can be extracted from over $50\%$ of all bitstrings on average. Additionally, all average success probabilities are significantly larger than the theoretical bound of 3–4% (see Appendix A.2). We hypothesize that asymptotically, the average success probability for “success + lucky” approaches $50\%$ from above (further evidence is given in Section 3.1.1 below, where we give a classification of the different “lucky” scenarios and show that the main contribution saturates around $25\%$). This observation is remarkable, as it shows that factoring a semiprime with Shor’s algorithm is often successful, even though the order-finding procedure actually fails.

Simulating Shor’s algorithm for semiprimes N between 536,870,861 and 549,755,813,701 requires substantial computational resources. Therefore, only individual cases are shown in Figure 4a. These cases correspond to the largest “interesting” semiprimes for a given number of qubits $n=30,\dots ,40$. A noteworthy case is the factoring problem for $(N,a)$ = (8,589,933,181, 3,974,323,683) ($n=34$ qubits). Here, the “lucky” cases raise the success probability from $56.25\%$ to $100\%$ (yellow square) among all M bitstrings. Furthermore, the factoring problem for $(N,a)$ = (274,877,906,893, 226,009,433,972) ($n=39$ qubits, second from the right) has a success probability of $0\%$ when the sufficient conditions for Shor’s algorithm are presupposed (green circle). However, when ignoring the violations of these conditions, we find that Shor’s algorithm can indeed factor N with a “lucky” success probability of $12.5\%$ (yellow squares).

The unexpectedly large success probabilities when the lucky cases are included prompt the question “how many bitstrings do we need to sample until a factor is found?” This is a relevant question, since for large problems, computing time on both classical and quantum computers is an essential resource. Figure 5a demonstrates that, for more than half of all factoring problems examined, the first sampled bitstring already yields a factor of N. Furthermore, in only $7.7\%$ of all factoring problems $(N,a)$, none of the 1024 bitstrings produced a factor. In this case, the reason is usually that the choice of a was bad, which can be estimated to happen with probability 50% (see Proposition A3 in Appendix A.2). Clearly, the failure probability of $7.7\%$ is much smaller than the theoretical estimate would suggest.

Figure 5b further reveals that, even when the order-finding procedure in Shor’s algorithm fails, the first bitstring often still produces a factor. In $38\%$ of all cases, the first bitstring yields the order of a modulo N (leftmost bar). From Figure 4b, we know that on average, $75\%$ of all factoring problems can be solved after the order is known (blue triangles). Thus we expect approximately $38\%\times 75\%\approx 29\%$ of all factoring problems to be solved by the correct order after the first bitstring. However, in Figure 5a, we see that $56\%$ of all problems are solved by processing the first bitstring. This percentage obviously is much larger than $29\%$, implying that it is easier to find a factor with Shor’s algorithm than to solve the underlying order-finding problem.

Another interesting result is observed when reducing the number of bits t in the sampled bitstring below the recommended $\lceil 2{log}_{2}N\rceil$ (cf. Appendix C). This saves resources in both versions of Shor’s algorithm. For the conventional Shor algorithm, it reduces the required number of qubits and gates required for the QFT. For the iterative Shor algorithm, it linearly reduces the number of quantum gates and thus the execution time.

Surprisingly, in almost all cases, reducing the number of bits t still allows for a successful factorization. Three representative cases are shown in Figure 6. First, Figure 6a shows that reducing t may even increase the frequency of “lucky” factorizations to over $99\%$, as it does for $10\le t\le 13$ in this case. Second, in Figure 6b, we see that even though the success probabilities decrease with t, at half of the recommended number of classical bits, that is at $t=14$, there are still “lucky" cases, allowing for successful factorization. Finally, in the case shown in Figure 6c, the success and lucky probabilities are essentially constant for $4\le t\le 28$.

Although it is known that reducing t may still allow for non-zero success probabilities [39,44,46,61], the surprising robustness (or even increase) of the “lucky” success probabilities has not been appreciated. In conclusion, Shor’s algorithm can still be successful (sometimes even more successful) if much less classical bits t are sampled than the recommended $t=\lceil 2{log}_{2}N\rceil$.

3.1.1. Classification of the “Lucky” Scenarios

If a sampled bitstring j does not pass the standard tests required by Shor’s algorithm, but still produces a factor with the procedure shown in Figure 1, we call this a “lucky” case. As shown above, this happens much more often than expected. In this section, we explain and classify the different scenarios that can happen.

For a given factoring problem with an L-bit semiprime N and an integer a coprime to N, let $\widehat{r}$ be the multiplicative order of a modulo N. Furthermore, let r be the denominator and k be the numerator extracted from the convergent $k/r$ to $j/2^t$ using the continued fractions algorithm from the standard procedure (i.e., the largest $r<N$ such that $k/r$ is a convergent to $j/2^t$ with $|k/r-j/2^t|\le 1/2r^2$ [10]). We distinguish between three scenarios in which the standard checks for Shor’s algorithm fail:

• (n,e) $r\ne \widehat{r}$ is not the order but r is even,
• (n,o) $r\ne \widehat{r}$ is not the order and r is odd,
• (o,o) $r=\widehat{r}$ is the order but the order is odd.

In Figure 7, we show a breakdown of the average success probability for these scenarios. We see that the (n,e) scenario makes up a large fraction of the successful factorizations and its relevance grows for larger integers. In contrast, the (n,o) and (o,o) scenarios, where the extracted r is odd, only matter for smaller integers. We explain the reasons for this in the discussions of each individual scenario below.

Figure 8 shows the number of cases for each scenario among the 52,077 uniformly drawn factoring problems. We see that the cases where bitstrings yield a factor in the (n,e) scenario are responsible for a significant fraction of all successful factorizations. Indeed, as Figure 7 suggests, the relevance of this scenario also grows on average and tends to saturate above $25\%$ for larger L. We expect that this contribution persists for even larger semiprimes.

Contributions from the (n,o) and (o,o) scenarios seem to be responsible only for a small number of successful factorizations, and mostly only for small semiprimes up to $L=10$. It is remarkable, however, that for many factoring problems that can be factored in the (o,o) scenario, 50–100% of all bitstrings yield a factor (see Figure 8c).

To understand these effects, we discuss the different scenarios individually. The goal is to obtain an understanding for why the different scenarios occur. The number-theoretic ideas are similar to the algorithm used in [33], which can be traced back to Miller’s algorithm ([103], Lemma 5).

3.1.2. The (n,e) Scenario

From the quantum circuit of Shor’s algorithm, one can compute the probability distribution for the bitstrings j that are sampled at the measurement [104] (see Appendix A for the derivation),

$$\begin{array}{cc}\hfill p_{\widehat{r},t}\left(j\right)& =\frac{\widehat{r}}{2^{2t}}{\left(\frac{sin(\pi \widehat{r}j\lfloor \frac{2^t}{\widehat{r}}\rfloor /2^t)}{sin(\pi \widehat{r}j/2^t)}\right)}^2+\frac{2^t-\widehat{r}\lfloor \frac{2^t}{\widehat{r}}\rfloor }{2^{2t}}\frac{sin(\pi \widehat{r}j[2\lfloor \frac{2^t}{\widehat{r}}\rfloor +1]/2^t)}{sin(\pi \widehat{r}j/2^t)},\hfill \end{array}$$

(26)

where $\lfloor 2^t/\widehat{r}\rfloor$ denotes the integral number of times that $\widehat{r}$ fits into $2^t$. Note that for a given factoring problem $(N,a)$ with t classical bits per bitstring, this distribution only depends on the order $\widehat{r}$ of a modulo N. This is a consequence of the fact that the QFT in Shor’s algorithm is used to determine the period of the function $f\left(k\right)=a^k\mathrm{mod}N$, which is exactly $\widehat{r}$.

The distribution in Equation (26) is shown for a few representative cases in Figure A1 in Appendix A. It is strongly peaked at $\widehat{r}$ bitstrings (see also [51])

$$\begin{array}{c}\hfill j\in \{\mathrm{round}(\widehat{k}\times 2^t/\widehat{r}):\widehat{k}=0,\dots ,\widehat{r}-1\},\end{array}$$

(27)

where $\widehat{k}$ enumerates the $\widehat{r}$ peaks. Given j, the continued fractions algorithm yields a convergent $k/r=\widehat{k}/\widehat{r}$ to $j/2^t$. However, if $\widehat{k}$ and $\widehat{r}$ have a common factor, the denominator r from the extracted convergent will not be equal to the order $\widehat{r}$. For instance, this is the reason that the “success” cases in Figure 4a are typically below $50\%$, since every second peak corresponds to an even $\widehat{k}$, and an even order $\widehat{r}$ is a sufficient condition for success; hence, at least a factor of two is lost in $\widehat{k}/\widehat{r}$. We note that with a very small probability, this procedure may also yield an $r>\widehat{r}$ if the sampled bitstring j is not at one of the $\widehat{r}$ peaks of $p_{\widehat{r},t}\left(j\right)$.

To understand why $r\ne \widehat{r}$ may still yield a factor of N, we consider the case that the order $\widehat{r}$ yields a factor of N (as Figure 4b shows, this case occurs with approximately $75\%$ frequency). In this case, we have

$$\begin{array}{cc}\hfill \gcd (a^{\widehat{r}/2}-1,N)& =p,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \gcd (a^{\widehat{r}/2}+1,N)& =q,\hfill \end{array}$$

(29)

where p and q are the two prime factors of N. Let $2^d$ be the largest power of 2 in $\widehat{r}$ such that $2\nmid \widehat{r}/2^d$ (meaning 2 does no divide $\widehat{r}/2^d$). Note that often, $2^d\ge 4$ since the multiplicative order of the whole group ${\mathbb{Z}}_N^{\ast }$ is $\varphi \left(N\right)=(p-1)(q-1)$, which is at least divisible by 4 (here, $\varphi \left(N\right)=|Z_N^{\ast }|$ is Euler’s totient function). In this case, Equation (28) can be written as

$$\begin{array}{cc}\hfill \gcd ((a^{\widehat{r}/4}-1)(a^{\widehat{r}/4}+1),N)& =p,\hfill \end{array}$$

(30)

so either $(a^{\widehat{r}/4}-1)$ or $(a^{\widehat{r}/4}+1)$ contain the prime factor p. Since every second $\widehat{k}$ in Equation (27) is even, it is likely that $r\in \{\widehat{r}/2,\widehat{r}/4,\dots ,\widehat{r}/2^{d-1}\}$ (each case decreasing in likelihood). If $r=\widehat{r}/2$, Equation (30) shows that by testing both $\gcd (a^{r/2}\pm 1,N)$, a factor will be found. If $r=\widehat{r}/4$, knowing that r is even, we can further write Equation (30) as

$$\begin{array}{cc}\hfill \gcd ((a^{\widehat{r}/8}-1)(a^{\widehat{r}/8}+1)(a^{\widehat{r}/4}+1),N)& =p,\hfill \end{array}$$

(31)

so if p does not happen to be in $(a^{\widehat{r}/4}+1)$, also a factor will be found. This reasoning can be iterated up to the unlikely case that $r=\widehat{r}/2^{d-1}$, where Equation (28) becomes

$$\begin{array}{cc}\gcd (\hfill & (a^{\widehat{r}/2^d}-1)(a^{\widehat{r}/2^d}+1)\hfill \\ \hfill & \times (a^{\widehat{r}/2^{d-1}}+1)\cdots (a^{\widehat{r}/4}+1),N)=p.\hfill \end{array}$$

(32)

Similarly, if $r=\widehat{r}/3$, we can write Equation (28) as

$$\begin{array}{cc}\hfill \gcd ((a^{r/2}-1)(a^r+a^{r/2}+1),N)& =p,\hfill \end{array}$$

(33)

in which case we also find a factor if p happens to be in the first part of the product. Finally, if $r=\widehat{r}/l$ with $l>4$, we can write Equation (28) as

$$\begin{array}{cc}\gcd (\hfill & (a^{r/2}-1)\hfill \\ \hfill & \times (a^{r/2\times (l-1)}+a^{r/2\times (l-2)}+\cdots +a^{r/2}+1),N)=p.\hfill \end{array}$$

(34)

Note that as l grows, this case becomes increasingly unlikely since l would need to be a factor of $\widehat{k}$ already. However, also in this case, there is a small chance that when evaluating $\gcd (a^{r/2}-1,N)$, a factor can be found. We remark that when $r/2=\widehat{r}/2l$ is prime, the decomposition in Equation (34) is irreducible [66], such that no further polynomial in $a^{r/2}$ including p can be factored out.

The distribution of the fractions $r/\widehat{r}$ (extracted from the data generated from the uniformly distributed factoring problems) is shown in Figure 9a. Indeed, we see that very often a small multiple of r is equal to the order $\widehat{r}$. In particular, the significance of fractions up to $r/\widehat{r}=1/11$ seems to increase with L, i.e., with increasingly large semiprimes N. This observation agrees well with the argument given in [33].

Interesting examples for the lucky (n,e) scenario are the individual problems for $n=34$ and $n=39$ discussed in Section 3.1. In particular, the $n=39$ case with N = 274,877,906,893 = 364,303 × 754,531, a = 226,009,433,972 and order $\widehat{r}$ = 45,812,798,010 violates the condition $a^{\widehat{r}/2}\not\equiv -1\left(\mathrm{mod}N\right)$. As this is one of the sufficient conditions for Shor’s algorithm to guarantee successful factorization, the corresponding “success” probability is zero (green circle). However, $12.5\%$ of the sampled bitstrings yield even integers $r\in \{\widehat{r}/3,\widehat{r}/5,\widehat{r}/111\}$, which still allow for a successful “lucky” factorization of N (the corresponding quadratic residues $a^r$ do not have a trivial square root, i.e., $a^{r/2}\not\equiv -1\left(\mathrm{mod}N\right)$).

For large N, the (n,e) scenario makes up the majority of all “lucky” factorizations (see Figure 7). We hypothesize that on average, the probability of “success+lucky” factorizations asymptotically approaches $50\%$ due to the (n,e) scenario.

3.1.3. The (n,o) Scenario

In the (n,o) scenario, the bitstring j yields an integer r that is neither the order $\widehat{r}$ nor even. Using the reasoning from the (n,e) scenario, this happens only if $r\mid \widehat{r}/2^d$, so all d powers of 2 must have been in $\widehat{k}$ from Equation (27) already (or $d=0$, in which case $\widehat{r}$ is odd). This is an unlikely scenario given that all $\widehat{k}\in \{0,\dots ,\widehat{r}-1\}$ occur with roughly the same probability, see Figure A1 in Appendix A. Moreover, as Figure 8 shows, the frequency of this scenario tends to zero for larger semiprimes N. However, it does occur for smaller semiprimes with up to $25\%$ frequency on average (see the red area in Figure 7), so it is instructive to understand how a factor can be found in this case. In what follows, we exclude the irrelevant case $r=1$ as it will never yield a factor.

Let $r=\widehat{r}/l$ where $2^d\mid l$ and $2\nmid r$ (note that $d=0$ is possible if the order $\widehat{r}$ is odd, but we always have $l\ge 2$ since $r\ne \widehat{r}$). In this case, using the procedure depicted in Figure 1, we test

$$\begin{array}{c}\hfill \gcd (a^{\lfloor r/2\rfloor }\pm 1,N)=\gcd (a^{(\widehat{r}/l-1)/2}\pm 1,N)\end{array}$$

(35)

to find a factor of N.

Since Equation (35) seems somewhat arbitrary from the perspective of the original theory behind Shor’s algorithm, one might think that a factor may only be found by coincidence. For instance, when $N=p\times q$ has a very small prime factor (say 3, 5, or 7), then whatever number is computed by Equation (35) might have a chance of including the small prime factor.

That this reasoning does not always hold is shown in Figure 10a, where we list the percentage of bitstrings that yield a factor in the (n,o) scenario as a function of the smallest prime factor of N. Indeed, we see that also larger prime factors can be found in certain cases. The most important of these are the cases in which a has a small order with respect to either p or q (purple squares, black crosses, and red pluses in Figure 10a). Indeed, one can prove

Theorem 1.

If

$$\begin{array}{c}\hfill \min \{{\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right)\}\in \{1,2\},\end{array}$$

(36)

then $\gcd (a^{\lfloor r/2\rfloor }\pm 1,N)$ with $r\ne \widehat{r}$ odd yields a factor of N.

Proof. 

Without loss of generality, we assume ${\mathrm{ord}}_p\left(a\right)\in \{1,2\}$. This implies that $a\equiv \pm 1\left(\mathrm{mod}p\right)$ (because if ${\mathrm{ord}}_p\left(a\right)=1$, we have $a\equiv 1\left(\mathrm{mod}p\right)$, and if ${\mathrm{ord}}_p\left(a\right)=2$, we have $a\equiv -1\left(\mathrm{mod}p\right)$). Hence,

$$\begin{array}{c}\hfill x_p:=a^{\lfloor r/2\rfloor }\mathrm{mod}p={(\pm 1)}^{\lfloor r/2\rfloor }\mathrm{mod}p=\pm 1.\end{array}$$

(37)

Moreover, since $\widehat{r}=\mathrm{lcm}({\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right))$, where $\mathrm{lcm}$ denotes the least common multiple, ${\mathrm{ord}}_p\left(a\right)\in \{1,2\}$ implies that ${\mathrm{ord}}_q\left(a\right)\in \{\widehat{r},\widehat{r}/2\}$.

Thus we have ${\mathrm{ord}}_q\left(a\right)\ge \widehat{r}/2\ge \widehat{r}/l>\widehat{r}/l-1$ (where $l\ge 2$ is defined above Equation (35)), and therefore

$$\begin{array}{c}\hfill x_q:=a^{\lfloor r/2\rfloor }\mathrm{mod}q=a^{(\widehat{r}/l-1)/2}\mathrm{mod}q\ne \pm 1\end{array}$$

(38)

(since otherwise $\widehat{r}/l-1$ would be a multiple of the order of a modulo q). Applying the Chinese remainder theorem [66] to Equations (37) and (38), we obtain

$$\begin{array}{c}\hfill a^{\lfloor r/2\rfloor }\equiv x_{q}p{p}_q+x_{p}q{q}_p\left(\mathrm{mod}N\right),\end{array}$$

(39)

where $p_q$ is the inverse of $p\left(\mathrm{mod}q\right)$ and $q_p$ is the inverse of $q\left(\mathrm{mod}p\right)$. Using that $p{p}_q+q{q}_p\equiv 1\left(\mathrm{mod}N\right)$, we finally have

$$\begin{array}{cc}\hfill \gcd (a^{\lfloor r/2\rfloor }\pm 1,N)& =\gcd ((x_q\pm 1)p{p}_q+(x_p\pm 1)q{q}_p,N).\hfill \end{array}$$

(40)

Thus, since $x_q\ne \pm 1$, the case where $x_p\pm 1=0$ will yield the factor p.    □

One can estimate how often the situation given by Equation (36) happens. Choosing a uniformly at random from $\{2,\dots ,N-1\}$ with $\gcd (a,N)=1$ is equivalent to choosing $a\mathrm{mod}p$ ($a\mathrm{mod}q$) uniformly at random from $\{1,\dots ,p-1\}$ ($\{1,\dots ,q-1\}$) with the exception of $a\mathrm{mod}p=a\mathrm{mod}q=1$. Thus there are $(p-2)(q-2)-1\sim pq$ choices for a, $2p+2q-13\sim p+q$ of these satisfying either $a\mathrm{mod}p=\pm 1$ or $a\mathrm{mod}q=\pm 1$. Therefore, the contribution of cases with $\min \{{\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right)\}=1,2$ becomes negligible for large p and q. We remark that individual cases with $\min \{{\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right)\}=3,4,\dots ,$ may further contribute to lucky (n,o) factorizations, as Figure 10a shows.

3.1.4. The (o,o) Scenario

In the (o,o) scenario, the bitstring j yields the order $r=\widehat{r}$ and the order is odd (see also [15,54,55,56]). Interestingly, as Figure 10b suggests, this case can be fully classified:

Theorem 2.

$$\begin{array}{c}\hfill \gcd (a^{\lfloor r/2\rfloor }-1,N)\mathit{yields}a\mathit{factor}\iff 1\in \{a\mathrm{mod}p,a\mathrm{mod}q\}.\end{array}$$

(41)

Note that $a>1$ by construction, so one of $\{a\mathrm{mod}p,a\mathrm{mod}q\}$ is larger than 1, and in particular $r>1$.

Proof. 

The “⇐” case is a special case of Theorem 1 from the (n,o) scenario. If $a\mathrm{mod}p=1$, we have $x_p=1$ in Equation (37) and $x_q\ne \pm 1$ in Equation (38), so $\gcd (a^{\lfloor r/2\rfloor }-1,N)=p$ by Equation (40).

We therefore only need to show the “⇒” case. Without loss of generality, let $\gcd (a^{\lfloor r/2\rfloor }-1,N)=p$, so $p\mid (a^{\lfloor r/2\rfloor }-1)$. From this follows that $a^{\lfloor r/2\rfloor }\equiv 1\left(\mathrm{mod}p\right)$, so ${\mathrm{ord}}_p\left(a\right)\mid \lfloor r/2\rfloor =(r-1)/2$. However, we also have ${\mathrm{ord}}_p\left(a\right)\mid r$ because $r=\mathrm{lcm}({\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right))$. Since $\gcd \left(\right(r-1)/2,r)=1$ (using the Euclidean algorithm), this is only possible if ${\mathrm{ord}}_p\left(a\right)=1$, which means $a\mathrm{mod}p=1$.    □

Theorem 3.

The “+” case, i.e., $\gcd (a^{\lfloor r/2\rfloor }+1,N)$, never gives a factor in the case $r={\mathrm{ord}}_N\left(a\right)$ odd.

Proof. 

Assume $\gcd (a^{(r-1)/2}+1,N)=p$. This means that $a^{(r-1)/2}\equiv -1\left(\mathrm{mod}p\right)$, and thus $a^{r-1}\equiv 1\left(\mathrm{mod}p\right)$. The former implies that ${\mathrm{ord}}_p\left(a\right)\nmid (r-1)/2$ and the latter implies that ${\mathrm{ord}}_p\left(a\right)\mid r-1$. Therefore, $2\mid {\mathrm{ord}}_p\left(a\right)$. From $r=\mathrm{lcm}({\mathrm{ord}}_p\left(a\right),{\mathrm{ord}}_q\left(a\right))$ then follows that $2\mid r$, but this is a contradiction because r was assumed to be odd.    □

Finally, we can show that it does not matter whether we round $r/2$ up or down, i.e., whether we take $\lfloor r/2\rfloor =(r-1)/2$ or $\lceil r/2\rceil =(r+1)/2$ in Figure 1, since one of them yields a factor whenever the other one also yields a factor:

$$\begin{array}{cc}\gcd (a^{\lfloor r/2\rfloor }-1,N)\hfill & =\gcd (a^{\lfloor r/2\rfloor }(a^{\lceil r/2\rceil }-1),N)\hfill \\ \hfill & =\gcd (a^{\lceil r/2\rceil }-1,N),\hfill \end{array}$$

(42)

where we used that $a^{\lfloor r/2\rfloor }{a}^{\lceil r/2\rceil }\equiv a^{(r-1)/2}{a}^{(r+1)/2}\equiv a^r\equiv 1\left(\mathrm{mod}N\right)$ and furthermore that $a^{\lfloor r/2\rfloor }$ cannot have a common factor with N (since a was chosen coprime to N),

A special, additional condition that yields a success probability of almost $100\%$ with the (o,o) scenario is when $\widehat{r}$ is prime, as indicated in Figure 10b. In this case, almost all bitstrings j are sampled at the peaks of Shor’s bitstring distribution $p_{\widehat{r},t}\left(j\right)$ given by Equation (27) and directly yield the order $r=\widehat{r}$, because $\widehat{k}$ and $\widehat{r}$ are always coprime. The only exceptions are either when j belongs to the first peak corresponding to $\widehat{k}=0$, or when j lies in the neighborhood of one of the peaks of Equation (26) where the probability is small.

3.2. Using Ekerå’s Post-Processing

In 2022, Ekerå has proven a lower bound for the success probability that takes into account additional, efficient classical post-processing procedures [34] (his implementation of the procedures can be found in [105]). While with Shor’s post-processing, a factor is found with more than $50\%$ probability on average after a single run (see Figure 5a), with Ekerå’s post-processing, it is possible to increase this probability arbitrarily close to unity. The bound reads

$$\begin{array}{c}\hfill p(\mathrm{success}\mid N)\ge \underset{j\mathrm{sampled}\pm B\mathrm{bitstrings}\mathrm{around}\mathrm{peak}}{\underbrace{\left(1-\frac{1}{{\pi }^2}\left(\frac{2}{B}+\frac{1}{B^2}+\frac{1}{3B^2}\right)-\frac{{\pi }^2(2B+1)}{\sqrt{2^{m+\ell }}}\right)}}\underset{\mathrm{peak}\mathrm{yields}\mathrm{order}\widehat{r}}{\underbrace{\left(1-\frac{1}{clogcm}\right)}}\underset{\mathrm{order}\widehat{r}\mathrm{yields}\mathrm{factors}\mathrm{of}N}{\underbrace{\left(1-2^{-k}\left(\genfrac{}{}{0pt}{}{n_F}{2}\right)-\frac{1}{2{\varsigma }^2{log}^2\varsigma L}\right)}},\end{array}$$

(43)

where L is the bit length of N, $m+\ell =t$ is the number of classical bits obtained from Shor’s algorithm, $n_F$ is the number of distinct prime factors of N (i.e., $n_F=2$ for semiprimes), and $B,c,k,\varsigma \ge 1$ are constants of the post-processing algorithms that can be freely selected. We choose $m=L$ and $\ell =t-L$ so that the results are in line with the analysis presented above. Note that the only technical requirement is $2^m>\widehat{r}$ and $2^{m+\ell }>{\widehat{r}}^2$, so $m=L-1$ is possible [34]; this does not make a difference for the results that follow. We remark that the three factors in Equation (43) are directly related to the three Propositions A1–A3 discussed in Appendix A.2. We discuss each factor in turn.

The first factor in Equation (43) comes from the idea that whenever the bitstring j is not sampled at one of the $\widehat{r}$ peaks (see Figure A1 in Appendix A), it is often sampled very close to a peak. Thus, one can try out all bitstrings in the range $\{j-B,\dots ,j+B\}$ for some small B. The probability to find the peak among these bitstrings can be estimated from the distribution $p_{\widehat{r},t}\left(j\right)$ in Equation (26). Instead of $4/{\pi }^2$ (see Equation (A20); this would correspond to $B=0$), we then obtain a larger probability, given by the first factor. Here, we choose $B=L$, i.e., the number of bits in N.

The second factor in Equation (43) stems from the idea that, when the continued fraction method does not yield the order $\widehat{r}$, it will often yield a large divisor $r=\widehat{r}/D$ for some small D (see Figure 9). Starting from r, Ekerå gives several classical algorithms in [34] to efficiently recover the real order $\widehat{r}=r\times D$. The corresponding success probability is given by $1-1/clogcm$, where $c\ge 1$ is a parameter that is free to choose. Its derivation is based on the probability that D is $cm$-smooth, meaning that $D>0$ is not divisible by any prime power larger than $cm$. For our numerical work we choose $c=1$.

Finally, the third factor in Equation (43) follows from the algorithm presented in [33] (see also [54]). This algorithm describes the factoring of an arbitrary composite integer N (with $n_F\ge 2$ distinct prime factors) given the order $\widehat{r}$ of a single element $a\in {\mathbb{Z}}_N^{\ast }$ selected uniformly at random. The corresponding success probability depends on two parameters $k,\varsigma \ge 1$ ($\varsigma$ is called c in [33]) that can be freely selected. For our numerical work we choose $k=100$ and $\varsigma =1$.

Probably the most important consequence of Ekerå’s result given by Equation (43) is that, as the size of the factoring problem becomes very large, i.e., $N,\widehat{r}\to \infty$, the success probability approaches one. This trend can already be seen in Figure 11a (gray line), which shows that the bound is increasing—even though it is already quite large for our modest choice of $(B,c,k,\varsigma )=(L,1,100,1)$. Furthermore, the gray diamonds in Figure 11a represent the actual success probabilities, obtained from applying Ekerå’s post-processing to the largest scenarios studied above. They are all larger than $93\%$ and thus even closer to unity than expected (this potential underestimation was noted in [34]).

3.3. Errors during the Execution of Shor’s Algorithm

With Ekerå’s post-processing [33,34], the expected success probability using only a single run of the quantum part of the algorithm can be brought arbitrarily close to $100\%$ by properly selecting the constants $(B,c,k,\varsigma )$ (cf. Equation (43)). However, these probabilistic estimates still require a successful execution of the quantum part of Shor’s algorithm. Since fully error-corrected, fault-tolerant quantum computers will probably not become available for several years to come [106,107,108], it is an interesting, relevant question to study how the performance of the post-processing algorithms is affected by errors during the execution of the quantum algorithm.

In this section, we consider five different models for errors arising in the quantum part of Shor’s algorithm. Each of these is shown in the inset of Figure 11b, which schematically marks the places in the iterative Shor algorithm (cf. Figure 2) at which the respective errors may occur.

1

Classical measurement errors (blue squares) are defined as misclassifications occurring directly after each quantum measurement process with a constant error probability $p_{\mathrm{error}}\left(\delta \right)=\delta$ (see Section 2.8.1).

2

Quantum measurement errors (yellow circles) are modeled as depolarizing quantum noise during the measurement process with effective error probability $p_{\mathrm{error}}\left(\delta \right)=\delta =p_x+p_y$ (see Section 2.8.2).

3

Amplitude initialization errors (green upward-pointing triangles) are modeled by initializing the recycled qubit not in the uniform superposition $|+\rangle =(|0\rangle +|1\rangle )/\sqrt{2}$, but by increasing the amplitude of $|0\rangle$ as a function of $\delta$. The effective error probability is $p_{\mathrm{error}}\left(\delta \right)=(1-\sqrt{1-{\delta }^2})/2$ (see Section 2.7.1).

4

Phase initialization errors (red down-pointing triangles) are defined by introducing a relative phase $e^{i\pi \delta }$ between the states $|0\rangle$ and $|1\rangle$ in the initialization. The effective error probability is $p_{\mathrm{error}}\left(\delta \right)=(1-cos\left(\pi \delta \right))/2$ (see Section 2.7.2).

5

Bit flip errors (purple stars) are defined by flipping each bit in the final bitstring j with probability $p_{\mathrm{error}}\left(\delta \right)=\delta$. This error model, in contrast to the others, does not affect the execution of the quantum part of the iterative Shor algorithm. While such an error (e.g., a fault in the classical computer memory) may be considered unlikely, it is still interesting to compare its consequences to the errors in the quantum part.

We consider the case that for each of these errors, Ekerå’s post-processing algorithm is applied to the resulting bitstrings, without the user being aware that one or more errors may have occurred.

Figure 11a shows the success probabilities for the different errors and problem sizes. We see that errors with $\delta =0.01$ (which corresponds to $1\%$ error probability for the blue squares, yellow circles, and purple stars) can decrease the success probability below the bound Equation (43) indicated by the solid gray line. Furthermore, the success probabilities show a decrease as a function of problem size that rivals the increasing success probability from the bound Equation (43). Nevertheless, Ekerå’s post-processing algorithm still produces correct factors even in the presence of errors.

Figure 11b shows the scaling of the success probability as a function of the effective single-qubit error probability for the 30-qubit case N = 536,870,903. For all errors, we see that the performance of the factoring algorithm including Ekerå’s post-processing scales similarly, despite the fundamental differences in the error models. This type of universal behavior is an interesting and unexpected observation.

It is also instructive to compare the simulation results to ${(1-p_{\mathrm{error}})}^t$ (black dash-dotted line in Figure 11b), which represents the probability to obtain a bitstring for which no error occurred (under the assumption that errors for individual bits are independent). Since we know from Figure 11a that the considered 30-qubit case is solved with $100\%$ success, ${(1-p_{\mathrm{error}})}^t$ represents the assumption that errors are destructive, i.e., an error in one of the bits prevents a successful factorization. Hence, the systematic gap between the dash-dotted line and the other dashed lines in Figure 11b shows that the quantum factoring problem can still be solved with Ekerå’s post-processing in the presence of errors. The fact that the success probability is systematically larger than the one for independent errors is encouraging, because an error at one stage in the iterative Shor algorithm affects the operation of all subsequent gates that depend on previous measurement results (see inset of Figure 11b). Such an error can thus propagate through the quantum algorithm and induce further, correlated errors. Our simulation results reveal a certain resilience when using Ekerå’s post-processing in combination with the iterative Shor algorithm for factoring integers.

3.4. Discussion of Limitations and Future Directions

Some of our design choices, made to achieve a large-scale simulation of Shor’s algorithm for many factoring scenarios, result in certain practical limitations to what we can simulate. In this section, we list these design choices, state the accompanying limitations, and discuss interesting future research directions that alternative choices could offer.

1

In a practical realization of Shor’s quantum circuit shown in Figure 2, most of the work is expected to be in the implementation of the exponentiation in terms of the controlled modular multiplications (see Section 2.2). Our choice to simulate the multiplications using direct permutations with no extra qubits, while allowing the large-scale MPI scheme sketched in Figure 3, prevents the direct simulation of quantum errors during the multiplications (which is why, in Figure 11, essentially only initialization errors before and measurement errors after the multiplications are shown). The alternative would be to implement a general multiplication circuit using standard quantum gates and additional workspace qubits. To pursue this research direction to allow the study of errors during the multiplications, an informative exposition to start from is the construction by Gidney and Ekerå [6], which combines and optimizes many techniques discovered over the past decades to implement the modular multiplications.

2

Although Shor’s order-finding algorithm is the most prominent quantum algorithm for factoring, a practical solution of the factoring problem on gate-based quantum computers might rather use the Ekerå-Håstad factoring scheme [63] based on the discrete logarithm quantum algorithm (see point 3 in Section 1.1). Instead of $t\approx 2L$ stages in the iterative quantum circuit (cf. Figure 2) using the semiclassical Fourier transform, this algorithm requires at most $t\approx 3L/2$ stages, with a systematic option to reduce t further at the cost of reducing the success probability below 99% [64]. In the context of quantum circuit simulation, the Ekerå-Håstad scheme would save valuable execution time (cf. Section 2.9), allowing for more statistics to be gathered for larger factoring scenarios.

3

The shorgpu implementation used for this work maintains two full statevector buffers psi and psibuf, which reduce the simulation time by enabling contiguous memory transfer through the MPI network (see [23]). However, the total amount of memory fixes the maximum number of qubits that can be simulated, which puts a limit on the size of simulatable factoring problems. An alternative choice would be to use only a single statevector buffer, thereby having to replace the contiguous memory transfer with interleaved communication and computation. This choice (potentially combined with reducing computing time by switching to the Ekerå-Håstad scheme, see previous item) would allow the simulation of yet another qubit, to push the boundary of simulatable factoring problems and the threshold of the proposed challenge one step further.

4. Conclusions

In this paper, we have introduced a method to simulate the iterative Shor algorithm on supercomputers with thousands of GPUs. The simulation software [23] allowed us to push the size of factoring problems far beyond what has been achieved previously. We have used the simulation software to perform an in-depth analysis of the iterative Shor algorithm.

Using Shor’s original post-processing [8,9,10,32], we have shown that a significant amount of “lucky” factorizations raises the expected success probability from 3–$4\%$ to above $50\%$. We have given number-theoretic arguments for the existence of the lucky cases, and we hypothesize that they continue to contribute with approximately $25\%$ beyond the size of integer factoring problems investigated in this paper.

Using Ekerå’s post-processing [33,34], the success probability for a factoring scenario can be brought close to unity using only a single bitstring obtained by executing the iterative Shor algorithm. However, Ekerå’s post-processing method assumes that the quantum part has been executed without errors, an assumption which is unlikely to hold for quantum processors in the near future. Therefore, we have studied how additional classical and quantum errors, as present in today’s quantum information processing hardware [106], influence the performance of the post-processing procedure. Remarkably, we find that Ekerå’s post-processing procedure exhibits a particular form of universality and resilience. Here, “universality” means that the decrease in success probability is roughly independent of the particular type of error and “resilience” means that the success probability is systematically larger than the success probability expected from independent bit flip errors.

Although these results might inspire confidence in the quantum factoring procedure, the first successful factorization of a cryptographically relevant number—say RSA-2048 from the famous RSA factoring challenge—is still out of reach [6,22]. Therefore, a more modest challenge towards true quantum supremacy might be to demonstrate that a real quantum computing device can factorize an interesting semiprime which is larger than $N_{\max }$ = 549,755,813,701. In fact, since gate-based quantum computers might already require full error correction for this purpose, it is conceivable that this challenge is first met by a quantum annealer [17,77,78,79,80,81,82].