Insights and Limitations of Shor's Factorization Algorithm on IBM Real Quantum Computers

Insights and Limitations of Shor’s Factorization Algorithm on IBM Real Quantum Computers Authors: Juan Manuel Villarreal D’Angelo, Gerard De Manresa, and Adriano Camps Tutor: Joseph Schindler Quantum Engineering Post-graduate Course (2024-2025) - Final project Fundaci´o Polit`ecnica de Catalunya, Barcelona, Spain [email protected], [email protected], [email protected] [email protected] Abstract—Quantum Shor’s algorithm offers an exponential speedup over classical methods for integer factorization, a central problem to modern cryptography, and the implications of Shor’s algorithm in factoring large numbers. The algorithm leverages on quantum parallelism to compute the period of a modular exponentiation function, employing the Inverse Quantum Fourier Transform (IQFT) to extract periodicity, and classical post-processing to derive the factors. This final report of the Quantum Engineering post-graduate course of the Fundaci´ o Polit` enica de Catalunya, presents the main findings encountered in the implementation and benchmarking of Schor’s algorithm in real IBM quantum computers. To accomplish this, we created a command line tool made in Python, that allows the integration of the circuits implementing the Shor’s algorithm with the overall set of configuration parameters needed to perform the benchmark. In this report, we first provide an overview of the quantum factoring problem, followed by a detailed description of Shor’s algorithm, and its implementation using the Qiskit Python library. The circuit depth and noise effects are analyzed across different backend types, also different values for the optimization level and approximation degrees are configured when performing the circuit transpilation. As mentioned earlier different backend types were utilized for performance comparison, namely backend class with or without a noise model (AerSim), a fake backend class provided by the Qiskit library (fake provider), and the real quantum hardware (IBM Quantum Processing Units in JanuaryJune 2025).The study follows by analyzing the different error mitigation techniques: dynamic decoupling, Pauli twirling, the SABRE Optimization (Stochastic Algorithm for Quantum Boolean Optimization and Routing), and any combination of them. Once the optimum configuration is found, the maximum number that can be reasonably factorized with success is explored. Finally, the main conclusions, findings, and possible future research lines are presented, including the factorization of a N≈1000 if appropriate parameters of the algorithm are selected. Index Terms—Quantum Computing, Shor’s Algorithm, Factorization, Quantum Fourier Transform, Modular Exponentiation, Cryptography I. INTRODUCTION Integer factorization is a computationally hard problem underpinning the security of widely used cryptographic systems like RSA. Classical algorithms, such as the quadratic sieve or general number field sieve, factor large integers in sub-exponential time, making them inefficient for large numbers. Shor’s algorithm, introduced by Peter Shor in 1994 [1], leverages on quantum computing to achieve polynomial-time factorization, posing a potential threat to RSA-based encryption. The first experimental demonstration of Shor’s algorithm dates from 2001[2], when N= 15 was factored. In 2007 Lu et al.[3] factorized N= 15 using a photonic implementation. As of today, the largest numbers factorized by true Shor’s algorithm is N= 21 on actual quantum hardware (2012)[4], and in an IBM quantum processor in 2021 [5]. In [6] it is estimated that factoring a security-standard 2,048-bit number would require a quantum computer with 20 million qubits. This work explores Shor’s algorithm, detailing its components—modular exponentiation, Quantum Fourier Transform (QFT), and classical post-processing—and their roles in achieving quantum speedup. We also address practical challenges, including circuit depth and noise in quantum hardware, which limit current implementations. Additionally, a Python-based simulation using Qiskit[10] is described to illustrate the algorithm’s mechanics for small integers. The significance of circuit depth and noise is analyzed to highlight barriers to scaling Shor’s algorithm on real quantum devices, emphasizing the need for advancements in quantum error correction. II. SHOR’SALGORITHM Before entering in depth in the description of the algorithm itself, its implementation, and the results obtained, we will explain Shor’s algorithm in a simple way. Shor’s algorithm is a way to break down large numbers into their smaller building blocks (called prime factors), for example, breaking 15 into 3 × 5. This is important because most of our internet security relies on the fact that factoring huge numbers (with hundreds of digits) is nearly impossible for regular computers, so if you can factor these numbers quickly, you can break encryption. A. Shor’s Algorithm for Dummies Shor’s algorithm relies on finding a hidden pattern. Let’s try to illustrate it with an example, trying to factorize the number 15. •Initialization: pick a random number smaller than 15 that doesn’t share factors with 15, for example 7. •Finding the repetition pattern: we look at powers of 7, but we will only care about the remainder after dividing it by 15. For example, 70= 1 (the remainder of 1 when divided by 15 is 1), 71= 7 (remainder 7), 72= 49 (remainder 4), 73= 343 (remainder 13), 74= 2401 (remainder 1), so we are back to 1. It is clear now that the pattern repeats every 4 steps: 1 →7→4→13 →1→7... •Factorization: Here, we use the pattern to find the factors. Since the pattern length is 4, we calculate: 74/2= 7 = 49, then we compute 49 −1 = 48 and 49 + 1 = 50. •Finding the greatest common divisor: The final step is to find the greatest common divisor (gcd) of the previous two numbers: gcd(48,15) = 3, gcd(50,15) = 5,(1) Finally, the result is 15 = 3 × 5. The quantum part relies on superposition to find the pattern length much faster than a classical computer by testing many possibilities simultaneously, as opposed to test each power one by one, which is much slower. Real encryption uses numbers with hundreds of digits. Finding their repeating patterns would take classical computers forever, but quantum computers could potentially do it in reasonable time. However, today’s quantum computers can only handle simple examples like this one, and we would need thousands of error-corrected qubits to threaten real encryption that uses 300+ digit numbers. Figure 1illustrates the implementation of Shor’s factorization algorithm for N= 15 and a= 4. Here is a description of the circuit architecture: •Initialization of the qubits and application of Hadamard gates (left section): Control Register are 9 qubits labeled control0 to control8, that are used to create a superposition of states and store the input xfor the function f(x). Note that the number of control qubits should be larger than n1= 2⌈log215⌉+ 1 = 9. Target Register are 4 qubits labeled target0 to target3. We need enough qubits to represent numbers up to 15. With 4 qubits, we can represent values from 0 to 15 (since 24= 16), which is sufficient. These hold the output of the function f(x) = axmod 15. Finally, output (lower line) is not a qubit, but represents the final Fig. 1: Quantum circuit implementing Shor’s algorithm for N= 15 and a= 4. Fig. 2: Sample periodic function axmod N,a= 4, N= 15 whose period we try to find. measurement of the control register after the Inverse Quantum Fourier Transform (IQFT), which gives the period. •Modular exponentiation (middle section): computes f(x) = axmod 15 for each x(Figure 2), and the result is stored in the control qubits. •Inverse Quantum Fourier Transform - IQFT (right Section): The IQFT converts the superposition of xvalues into a superposition of frequencies. Peaks in the frequency domain correspond to the period r, as illustrated in figure 3obtained in a real QPU for different values of a. Clean peaks at 0 and 0.5, i.e. separation equal to 1/2, corresponds to r= 2, and peaks at 0 and 0.25, 0.5, 0.75, i.e. separation equal to 1/4, correspond to r= 4. •Measurement (rightmost section): After the IQFT, measuring the control qubits (gray gauges) gives a value related to 29/r, which can be used to determine rin a classical way. This circuit demonstrates the core quantum part of Shor’s algorithm: the use of superposition and the IQFT to find the period of a function efficiently. B. Shor’s Algorithm: Detailed Explanation In this Section we will explain in detail the different operations that take place in the computation of the Shor’s algorithm. 2 (a) (b) (c) Fig. 3: Estimated phases from the order finding circuit for different values of ain a real QPU (same parameters in all three cases, irrelevant at this stage to understand the operation of the algorithm): a) Peak separation is 1/4. r= 4, b-c) Peak separation is 1/2. r= 2 1) Modular Exponentiation: Modular exponentiation is the computational core of Shor’s algorithm, enabling the quantum speedup that makes integer factorization exponentially faster than classical methods1. Given an integer Nto be factored and a randomly chosen integer asuch that gcd(a, N) = 1, Shor’s algorithm aims to find the order r, the smallest positive integer satisfying ar≡1 mod N. The function f(x) = axmod N is periodic with order r, and computing this function quantumly exploits superposition to evaluate it for many values of xsimultaneously, a key feature distinguishing quantum from classical approaches. The quantum circuit for modular exponentiation operates on two registers. The first register, with n1= 2⌈log2N⌉+ 1 qubits2, encodes the input x. The second register, with n2=⌈log2N⌉qubits, stores the output f(x). The algorithm begins by initializing the first register in a uniform superposition using Hadamard gates applied to each qubit: |ψ0⟩=1 √2n1 2n1−1 X x=0 |x⟩|+⟩.(2) The modular exponentiation is implemented via a unitary operator Uf, defined as: Uf|x⟩|y⟩=|x⟩|y·axmod N⟩,(3) where yis initially |1⟩in the second register to handle multiplicative operations. Applying Ufto the initial state produces: |ψ1⟩=1 √2n1 2n1−1 X x=0 |x⟩|axmod N⟩.(4) 1For a comprehensive introduction to Shor’s algorithm and its quantum foundations, see [7]. 2⌈ ⌉ denotes the ”ceil” function: round to the immediately larger integer. This entangled state encodes the periodic function f(x), with the order rembedded in the first register’s values (control0 to control8), setting the stage for the subsequent period finding via the Quantum Fourier Transform. Classically, modular exponentiation is computed efficiently using the square-and-multiply algorithm, which reduces the number of multiplications to O(log2x). Quantum implementations of such arithmetic operations are discussed in [8]. For an exponent xrepresented in binary as x=Pn1−1 j=0 xj2j, the exponentiation is: axmod N= n1−1 Y j=0 (a2jmod N)xj.(5) In the quantum circuit, this is implemented as a sequence of n1controlled modular multiplications. For each qubit |xj⟩in the first register, a controlled operation multiplies the second register by a2jmod N, if xj= 1. The values a2jmod Ncan be precomputed classically to reduce quantum gate count, though on-the-fly computation is possible for generality. Each modular multiplication is a reversible quantum operation. To multiply the second register by b=a2j mod N, a quantum circuit performs: |y⟩→|y·bmod N⟩,(6) conditioned on |xj⟩=|1⟩. Appendix 1 presents the implementation of the the exponential multiplication, and discusses ways to optimize it. 2) Inverse Quantum Fourier Transform: The Inverse Quantum Fourier Transform (IQFT) is a critical component of Shor’s algorithm, enabling the extraction of the order rof the function f(x) = axmod N, where Nis the integer to be factored and ais a randomly chosen integer such that gcd(a, N) = 1. The direct QFT transforms a quantum state encoding the input 3 register into the frequency domain, amplifying states corresponding to multiples of the period’s reciprocal, thus facilitating period finding. This is a fundamental step to achieve the exponential speedup of Shor’s algorithm over classical factorization methods, as it efficiently identifies the periodicity that leads to N’s factors. The QFT is the quantum analog of the classical discrete Fourier transform, operating on an n-qubit register with basis states |0⟩,...,|2n−1⟩. For an input state |x⟩, where x=Pn−1 j=0 xj2jis the binary representation, the QFT is defined as: QFT|x⟩=1 √2n 2n−1 X k=0 e2πixk/2n|k⟩.(7) In Shor’s algorithm, the IQFT is applied to the first register (with n1=⌈2 log2N⌉+1 qubits) after modular exponentiation, which produces the state given in eqn. (4). Since f(x) = axmod Nis periodic with order r, the first register’s state can be approximated (ignoring phase offsets for simplicity) as a superposition of states |x⟩where f(x)repeats. Because of the periodicity of f(x), the application of the IQFT to |ψ1⟩makes the amplitudes to have peaks at values of k≈2n1j/r, where jis an integer. Measuring the first register yields k, which is used to estimate rvia classical post-processing (e.g., continued fraction algorithm). Appendix 2 presents the QFT implementation and some ways to optimize it. 3) Classical Post-Processing: Classical post-processing is the final step in Shor’s algorithm, transforming the quantum measurement output into the factors of the composite integer N. After modular exponentiation and the Inverse Quantum Fourier Transform (IQFT), the algorithm measures the first register to obtain a value k, which encodes information about the order rof the function f(x) = axmod N, where ais a randomly chosen integer satisfying gcd(a, N)=13. The classical post-processing step uses the continued fraction algorithm to extract r, then computes the greatest common divisor (gcd) to derive non-trivial factors of N. This step is critical, as it connects the quantum computation’s probabilistic output with deterministic factor extraction, completing the factorization process with high probability. Measuring the first register yields k, an integer between 0 and 2n1−1. The goal is to estimate rfrom the fraction k/2n1, which approximates j/r for some j 3The selection of aplays a critical role in the success of Schor’s algorithm, as ahas to be a coprime number of N. If if is not, the transpilation of the code gives an error, as the matrices Ubare not unitary. For small values of N, as the values that can be factored in today’s QPU’s, we can know a priori the list of coprime factors smaller than N, but for very large numbers, as the ones used in cryptography, the values of aare selected randomly. coprime to r. The continued fraction algorithm is used to find the convergent j/r such that:  k 2n1−j r≤1 2n1+1 .(8) This ensures ris accurately recovered with high probability, provided n1is sufficiently large. The continued fraction algorithm expansion of a real number α= k/2n1is computed as: α=a0+1 a1+1 a2+1 ... ,(9) where aiare integers obtained via: •Set α0=α. •For each i, compute ai=⌊αi⌋,αi+1 = 1/(αi−ai), until αiis integer or the expansion terminates. The convergents of the continued fraction are fractions pm/qm, computed recursively: •Initialize: p0=a0,q0= 1;p1=a0a1+1,q1=a1. •For m≥2:pm=ampm−1+pm−2,qm= amqm−1+qm−2. The convergent pm/qmwith qm≤Nand gcd(pm, qm)=1is a candidate for j/r. The denominator qmis tested as the period r. Once ris found, it is checked for suitability: •rmust be even (since odd periods do not yield factors in this context). •rvalidates the equation: ar≡1 mod N. •ar/2≡ 1 mod N, to avoid trivial results. If these conditions hold, the factors of Nare computed as: f1= gcd(ar/2−1, N), f2= gcd(ar/2+1, N).(10) The gcd is efficiently computed using Euclid’s algorithm, which has complexity O((log2N)2)4If f1or f2are non-trivial (i.e., 1< fi< N), they are factors of N. If r is invalid or yields trivial factors, the algorithm repeats with a new a. The success of classical post-processing depends on the IQFT outputting ksuch that k/2n1≈j/r with gcd(j, r) = 1. Shor’s analysis shows this occurs with probability at least ϕ(r)/r ≥1/(2 log2log2r), where 4The Euclidean Algorithm is a method to find the greatest common divisor (gcd) of two non-negative integers. It is based on the principle that the gdc of two numbers also divides their difference. The procedure is as follows: given two non-negative integers aand b, where a≥b, divide aby bto obtain a quotient qand remainder r, such that a= b·q+r, 0≤r < b. If r= 0, then bis the gcd, but if r= 0, then set a=b,b=r, and repeat the previous step, The gcd is the last non-zero remainder. 4 ϕis Euler’s Totient function 5For large r, this probability is sufficiently high, and multiple measurements increase the likelihood of success. The choice of n1= ⌈2 log2N⌉+ 1 ensures the continued fraction algorithm converges to the correct rwith high probability. Classical post-processing is computationally efficient, with the continued fraction algorithm and gcd computation both running in polynomial time (O((log2N)3)). Unlike quantum steps, this phase is robust, as it relies on well-established classical algorithms. However, the quality of the quantum output kis critical. Noise in the quantum circuit (e.g., gate errors of 10−3to 10−2in 2025 hardware) can distort k, leading to incorrect convergents. Fault-tolerant quantum computing is thus essential for large N. For small N(e.g., N= 15), simulations in Qiskit demonstrate successful period finding (e.g., r= 4 for a= 7), but scaling to cryptographic sizes (e.g., 2048bit N) requires error-corrected qubits. The algorithm 1outlines the Shor’s algorithm implementation, and algorithm 2the classical postprocessing. III. COMMAND-LINE TOOL CREATED TO BENCHMARK SHOR’S ALGORITHM QUANTUM CIRCUIT As we introduced in the abstract, in order to accomplish the main objective of performing the benchmark over the Shor’s algorithm quantum circuit, we created a command line tool made in Python[11]. Among the main features are: •Account management (IBM cloud authentication). •Parametrization of the quantum circuit instantiated to perform the benchmark, providing the possibility to test different circuits implementations of Shor’s algorithm. •Parametrization of the different backend types instantiated, namely, backend class with or without a noise model (AerSim), a fake backend class provided by the Qiskit library (fake provider), and the real quantum hardware (IBM Quantum Processing Units). •Parametrization of the different availables QPUs to utilize when sending the jobs to sample to the IBM cloud. In this project the UPC provided 5Euler’s Totient Function, denoted ϕ(n), counts the number of positive integers up to a positive integer Nthat are coprime to N. Two integers are coprime if their greatest common divisor (gcd) is 1. Thus, ϕ(N)is the number of integers kin the range 1≤k≤Nsuch that gcd(k, N)=1. For a positive integer N, the totient function is defined as ϕ(N) = |{k∈Z|1≤k≤N, gcd(k, N)=1}|. If N has the prime factorization N=pe1 1pe2 2···pek k, where piare distinct primes and eiare their exponents, then ϕ(N) = NQp|N1−1 p,, where the product is over all distinct primes pdividing N. Cases of interest: for a prime p,ϕ(p) = p−1, for a prime power pk, ϕ(pk) = pk−pk−1=pk−1(p−1), and for coprime integers m and n,ϕ(m·n) = ϕ(m)·ϕ(n). Algorithm 1 Shor’s Algorithm Implementation Input: Number N, arguments (args), a-coefficients Output: Factors of N Check if Nis even or a prime power; return factors if true Compute gcd(a, N); return factors if gcd = 1 Initialize Shor with backend, sampler, transpiler parameters for each ain a-coefficients do Create circuit: Create control (n1), target (n2=) registers Apply Hadamard gates to control register Apply controlled modular exponentiation for a2jmod N Apply inverse QFT to control register Measure control register to obtain k Transpile circuit (inlcuding SABRE optimization only if set) Run sampler with shots (dynamical decoupling and Pauli twirling are included on the sampling only if set) Process results: compute rvia continued fractions Compute factors: gcd(ar/2±1, N) if factors are non-trivial then Return: Factors end if end for Return: Failure (repeat with new a) free and payed access to the following IBM QPUs: ibm_aachen,ibm_strasbourg, ibm_brussels for payed access, and ibm_sherbrooke,ibm_brisbane, ibm_torino for free access, all of them under the FPC account on IBM cloud. It should be noted that before this access was provided we had access to ibm_kyiv with our personal accounts before his decommission on April 18th 2025, this is worth to mention because one of our best results in terms of circuit transpilation and results came up from this QPU. •Parametrization of the main transpilation features when using Qiskit transpiler. •Parametrization of the SABRE optimization feature when using Qiskit transpiler. •Parametrization of the main sampling features when using Qiskit sampler. •Transpile using a Fake provider (ibm kyiv in this case) and then setting the initial layout of the transpiled circuit into another ibm qpu transpilation process for circuit depth comparison. The key point is that, the QPU fake provider class and the target QPU need to belong to the same family of processor 5 Algorithm 2 Classical Post-Processing in Shor’s Algorithm Input: Measured value k, number of qubits n1, integer N, base a Output: Factors of Nor failure Compute α=k/2n1 Initialize continued fraction: α0=α,a0=⌊α0⌋, p0=a0,q0= 1 Initialize: p1=a0a1+ 1,q1=a1, where a1= ⌊1/(α0−a0)⌋ for each m= 2,3, . . . until qm> N do Compute αm= 1/(αm−1−am−1),am=⌊αm⌋ Update: pm=ampm−1+pm−2, qm=amqm−1 + qm−2 if qm≤Nand gcd(pm, qm) = 1 then Set r=qm if ris even and ar≡1 mod Nthen Compute f1= gcd(ar/2−1, N), f2= gcd(ar/2+ 1, N) if 1< f1< N or 1< f2< N then Return:f1, f2 end if end if end if end for Return: Failure (repeat with new a) type. •Save transpiled circuits to disk in order to avoid transpiling them again when sending them to the IBM QPU for sampling. •Obtain results from sampling jobs submitted previously. •Generate all the visual outputs that the current work is presenting: quantum circuit, quantum circuit physical layout, probability distribution for normalized order count, with the circuit statistics and noise percentage included. •Calculate a set of optimal values of abased on the number that we want to factor. •Full integration with Qiskit SDK, Qiskit runtime environment, Qiskit visualization tools. IV. DESCRIPTION OF THE IMPLEMENTED CODE The implementation of Shor’s algorithm provided leverages Qiskit, a Python-based quantum computing framework, to simulate or execute the algorithm for factoring small (TBD) integers on quantum hardware or simulators. The code orchestrates the algorithm’s key components: modular exponentiation, Inverse Quantum Fourier Transform (IQFT), and classical post-processing. It supports both classical pre-checks and quantum computation, with optimizations for real quantum hardware, including transpilation and error mitigation. This section details the implementation’s structure, functionality, and practical considerations, highlighting its alignment with Shor’s algorithm’s theoretical framework. The entry point is main.py, which initializes the program by loading settings and parsing commandline arguments using CommandLineParser and ConfigurationIni classes. It supports two modes: processing results from a previous job (via job_id) or running Shor’s algorithm to factor a number N(via Factorizer). The Factorizer class in arithmetics.py coordinates the factorization process, performing classical checks (e.g., if Nis even or a prime power) before invoking the quantum algorithm via the Shor class in phase_estimation.py. The implementation is modular, with separate modules for circuit construction (circuit.py), sampling (sampling.py), and transpilation (transpilation.py). In arithmetics.py, the Factorizer class first attempts classical factorization to avoid unnecessary quantum computation. For an input number N, it checks: •If Nis even, it returns factors 2and N/2. •If Nis a prime power (N=db), it identifies factors dand db−1. •For a chosen a, it computes gcd(a, N). If gcd(a, N)= 1, it yields the factors directly. If these fail, it selects the coefficients aeither randomly, user-specified, or as coprimes of N), using get_circuits_coefficients method, then passes them to the Shor class for quantum period finding. The circuit.py module defines a QCircuit abstract base class and the RegisterQC that inherits from the QCircuit class and implements the create_circuit method. The create_circuit method in RegisterQC builds the quantum circuit for a given Nand a: •Registers: Initializes a control register with n1= 2⌈log2N⌉+ 1 qubits, a target register with n2= ⌈log2N⌉qubits, and a classical register for measurements. •Initialization: Applies an X gate to the target register’s least significant qubit (setting it to |1⟩) and Hadamard gates to the control register to create a superposition (eqn.(2)). •Modular Exponentiation: For each control qubit j, computes b=a2jmod Nclassically and applies a controlled unitary gate Uf, where: Uf|x⟩|y⟩= |x⟩|y·axmod N⟩.(eqn. (3)). The unitary is constructed as a permutation matrix in b_mod_n, mapping |y⟩→|b·ymod N⟩, which produces |ψ1⟩. 6 •IQFT: Applies an inverse QFT to the control register using Qiskit’s QFT library, transforming the state to amplify amplitudes at k≈2n1j/r. •Measurement: Measures the control register, yielding k. The transpilation.py module’s Transpiler class adapts the circuit for specific quantum hardware. It uses Qiskit’s pass_manager method with parameters like optimization level, basis gates, and routing methods (e.g., SABRE). The SabreLayout class optimizes qubit mapping and routing, with options for multiple iterations and trials to minimize circuit depth and gate count. The get_isa_statistics method collects metrics (e.g., two-qubit gate count, circuit depth), critical for assessing hardware feasibility. For example, a circuit for N= 15 may require hundreds of gates, challenging for NISQ devices with error rates of 10−3to 10−2. The sampling.py module defines Sampler classes (AerSampler for simulators, IBMRuntimeSampler for IBM Quantum hardware). The Shor class in phase_estimation.py orchestrates circuit execution: •Creates circuits for multiple a-coefficients. •Transpiles circuits using Transpiler class. •Submits jobs to the sampler with configurable shots (e.g., 1024). •Optionally waits for results, processes measurement counts, and extracts candidate orders rvia get_candidate_rs method. The _set_sampler_params include error mitigation techniques like dynamical decoupling and Pauli twirling, reducing decoherence and gate errors on real hardware. Post-processing occurs in job_results.py via find_nontrivial_factors method, which: •Uses the measured kto compute α=k/2n1. •Applies the continued fraction algorithm to find r. •Computes factors as gcd(ar/2±1, N). The plot_results_distribution function visualizes measurement outcomes, aiding analysis of period candidates. The implementation is optimized for small (TBD) N (e.g., N= 15), where a= 7 yields r= 4, producing factors 3 and 5. For larger N, the circuit depth (thousands of gates for 2048-bit N) and noise (gate errors ∼10−3) make execution infeasible on current hardware, as it will be shown later. The code mitigates this with transpilation optimizations and error mitigation, but fault-tolerant quantum computing is required for cryptographic applications. This implementation effectively simulates Shor’s algorithm for small N, benefiting from Qiskit’s modularity and hardware-aware optimizations. Scaling to cryptographic sizes would require advances in quantum hardFig. 4: Periodicity of axmod N for N= 15,a= 2,3,4...14. Note that for a= 6,10 the plot is flat, and for a= 4,a= 11 and a= 14 there is an integer number of periods. ware and error correction techniques. In the last section of this report, the authors will explore the maximum achievable Nthat can be factorized in the available IBM quantum computers. In order to gain a better understanding on the performance of Shor’s algorithm, Figure 4represents the periodicity of the function axmod N for N= 15, and a= 2,3,4...14. Valid values of aare those that are coprimes6of N, namely 2, 4, 7, 8, 11, 13 and 14. Values of athat are not coprimes of Nwill lead to nonunitary matrices and transpilation errors, but in the figure we included them all, as in the original Shor’s algorithm the values of aare randomly chosen, and there are high chances that the selected ais not actually a coprime of N. Additionally, to find the period of the function axmod N,amust be large enough so that at least two periods of the function fit in [0, N −1], that is: aN−1 min ≥2N. For N≥4,amin = 2, so the condition is always satisfied. Table Ilists all the semiprime numbers smaller than 77, and some of their coprimes. On the other hand, an interesting case occurs when r= 2, then for a=amax =N−1, then the gcd(ar/2±1, N) = gcd(a±1, N) = gcd(N−1±1, N)is gcd(N−2, N) = 1, and gcd(N, N) = N, which are the trivial factors of N. V. MATERIALS AND METHODS In this Section we first analyze the Circuit Depth and the Total Number of Gates for the thre QPUs available in the free IBM plan, namely ibm_Kyiv,ibm_Brisbane, ibm_Sherbrooke,ibm_Torino vs. those in the pay per use QPUs, namely ibm_Aachen, ibm_Brussels,ibm_Strasbourg... Then, we revise the quantum hardware specifications, and we finalize with an analysis of the different error mitigation 6Coprime numbers are two or more integers that have no common positive divisor other than 1, i.e. their greatest common divisor (gcd) is 1. 7 TABLE I: Semiprime numbers smaller or equal than 145, prime factors, coprime numbers smaller than N. NPrime factors Coprimes < N 42×2 1,3 62×3 1,5 93×3 1,2,4,5,7,8 10 2×5 1,3,7,9 14 2×7 1,3,5,9,11,13 15 3×5 1,2,4,7,8,11,13,14 21 3×7 1,2,4,5,8,10,11,...,20 22 2×11 1,3,5,7,9,13,15,...,21 25 5×5 1,2,3,4,6,7,8,...,24 26 2×13 1,3,5,7,9,11,15,...,25 33 3×11 1,2,4,5,7,8,10,...,32 34 2×17 1,3,5,7,9,11,13,...,33 35 5×7 1,2,3,4,6,8,9,...,34 38 2×19 1,3,5,7,9,11,13,...,37 39 3×13 1,2,4,5,7,8,10,...,38 46 2×23 1,3,5,7,9,11,13,...,45 49 7×7 1,2,3,4,5,6,8,...,48 51 3×17 1,2,4,5,7,8,10,...,50 55 5×11 1,2,3,4,6,7,8,...,54 57 3×19 1,2,4,5,7,8,10,...,56 58 2×29 1,3,5,7,9,11,13,...,57 62 2×31 1,3,5,7,9,11,13,...,61 65 5×13 1,2,3,4,6,7,8,...,64 69 3×23 1,2,4,5,7,8,10,...,68 74 2×37 1,3,5,7,9,11,13,...,73 77 7×11 1,2,3,4,5,6,8,...,76 82 2×41 1, 3, 5, 7, 9, 11, 13, ... 81 85 5×17 1, 2, 3, 4, 6, 7, 8, ... 84 86 2×43 1, 3, 5, 7, 9, 11, 13, ... 85 87 3×29 1, 2, 4, 5, 7, 8, 10, ... 86 91 7×13 1, 2, 3, 4, 5, 6, 8, ... 90 93 3×31 1, 2, 4, 5, 7, 8, 10, ... 92 94 2×47 1, 3, 5, 7, 9, 11, 13, ... 93 95 5×19 1, 2, 3, 4, 6, 7, 8, ... 94 106 2×53 1, 3, 5, 7, 9, 11, 13, ... 105 111 3×37 1, 2, 4, 5, 7, 8, 10, ... 110 115 5×23 1, 2, 3, 4, 6, 7, 8, ... 114 118 2×59 1, 3, 5, 7, 9, 11, 13, ... 117 119 7×17 1, 2, 3, 4, 5, 6, 8, ... 118 121 11 ×11 1, 2, 3, 4, 5, 6, 7, ... 120 122 2×61 1, 3, 5, 7, 9, 11, 13, ... 121 123 3×41 1, 2, 4, 5, 7, 8, 10, ... 122 129 3×43 1, 2, 4, 5, 7, 8, 10, ... 128 133 7×19 1, 2, 3, 4, 5, 6, 8, ... 132 134 2×67 1, 3, 5, 7, 9, 11, 13, ... 133 141 3×47 1, 2, 4, 5, 7, 8, 10, ... 140 142 2×71 1, 3, 5, 7, 9, 11, 13, ... 141 143 11 ×13 1, 2, 3, 4, 5, 6, 7, ... 142 145 5×29 1, 2, 3, 4, 6, 7, 8, ... 144 and optimization techniques, namely the Dynamical Decoupling, the Pauli Twirling, and the SABRE Optimization, respectively. To do this, the tool presented in Section III was used. A. Quantum Hardware IBM quantum computing resources available for this work are categorized into free and paid tiers, each offering distinct access levels and capabilities for implementing algorithms like Shor’s. Free-tier access, typically provided through platforms like IBM Quantum Experience, allows researchers to experiment with quantum hardware and simulators at no cost, but with significant limitations. These include restricted access to a subset of quantum processors (often smaller systems with fewer qubits, e.g., 127-qubit systems like ibm_kyiv,ibm_brisbane and ibm_sherbrooke), limited monthly runtime (e.g., 10 minutes), and lower priority in job queues, leading to longer wait times. For Shor’s algorithm, which requires deep circuits with thousands of gates for factoring even small numbers (e.g., N= 15), free-tier access often results in incomplete executions due to time constraints and noise accumulation on NISQ devices. Paid-tier access, such as IBM Quantum’s Premium or Enterprise plans, unlocks advanced hardware like ibm_aachen,ibm_brussels, and ibm_strasbourg, with higher qubit counts (156 qubits) and improved error rates (e.g., 6.28 ×10−4to 3.84 ×10−3for two-qubit gates). These plans offer dedicated runtime, and higher queue priority. For Shor’s algorithm, paid hardware enables execution of larger circuits with better fidelity, which is crucial for period finding and factorization of numbers beyond toy examples. However, even paid-tier systems struggle with cryptographic-scale N(e.g., 2048-bit), requiring fault-tolerant quantum computing not yet available in 2025. The choice between free and paid hardware thus balances cost, accessibility, and the feasibility of executing Shor’s algorithm effectively. In this study we have used only the following IBM Quantum processors to implement Shor’s algorithm, specifically ibm_aachen,ibm_torino,ibm_brisbane, and at the beginning of the semester also ibm_Kyiv which exhibit outstanding performances. Their specifications, including processor type, number of qubits, two-qubit gate error rates (best, layered, and median), coherence times, CNOT-length operations per second and basis gates 7, are summarized in Table II. Additionally, highlights for ibm_kyiv are provided for context, though only used at the beginning of this study while it was active, but not with the final version of the code. The specifications in Table II reveal the capabilities and limitations of each system for implementing Shor’s algorithm. Heron processors (ibm_aachen, 72 qubit gates: ECR gate (Echoed Cross-Resonance Gate): a native entangling gate in IBM’s superconducting QPUs, that uses crossresonance interactions with an echo sequence to reduce noise, modulating the target qubit based on the control qubit’s state via microwave pulses. CZ (Controlled-Z Gate): applies a Z (phase flip) operation to the target qubit if the control qubit is in the 1state, introducing a π phase shift to the 11 state, RZZ (Rotation around ZZ Axis): applies a rotation by an angle θaround the ZZ axis of the two-qubit Pauli operator (Z⊗Z), introducing a phase to the 11 state, used in variational algorithms or error correction, reducing circuit depth. 1 qubit gates: ID (Identity Gate), RX (Rotation around X Axis of the Bloch sphere by an angle θ), RZ (Rotation around Z Axis) of the Bloch sphere by an angle θ, applying a phase shift. SX (Square-Root of X gate): applies aπ/2rotation around the X axis, equivalent to the square root of the X gate, and X (Pauli X gate) used for bit flips 8 ibm_torino) feature lower error rates (e.g., 0.628 × 10−3for ibm_aachen) and higher qubit counts (up to 156), supporting deeper circuits with up to 250K CLOPS (Circuit Layer Operations Per Second or number or single-qubit rotations and single set of random two-qubit gates applied across a subset of qubits per second). Eagle r3 systems (ibm_brussels,ibm_strasbourg, ibm_brisbane,ibm_sherbrooke) offer consistent performance with 127 qubits, but higher error rates (e.g., 3.84 ×10−3for ibm_strasbourg) and shorter coherence times (100–150 µs) limit their ability to execute large-scale Shor’s circuits without significant error mitigation. The ibm_kyiv system, with moderate error rates and 200K CLOPS, aligns with other Eagle r3 systems in performance. Detailed system specifications for these QPUs are available from IBM Quantum [12]. B. Error Mitigation and Optimization Techniques 1) Error Mitigation - Dynamical Decoupling: Dynamical decoupling (DD) is an error mitigation technique designed to extend qubit coherence times by suppressing decoherence and noise in quantum systems, critical for executing deep circuits like those in Shor’s algorithm[13] [16]. This technique is particularly relevant for NISQ devices, such as IBM Quantum’s Heron and Eagle processors, which exhibit coherence times of 100–300 µs and gate error rates of 10−3to 10−2. In Shor’s algorithm, where modular exponentiation and the Inverse Quantum Fourier Transform (IQFT) require thousands of gates, decoherence from environmental interactions (e.g., magnetic fluctuations) can corrupt the quantum state. DD applies a sequence of fast, periodic singlequbit gates (e.g., X,Y, or XX pulses) to average out noise, effectively refocusing the qubit state. The implementation in the provided code (phase_estimation.py) configures DD via the DynamicalDecouplingOptions class within _set_sampler_params. Parameters include enable=True,sequence_type (e.g., XX or XY 4), extra_slack_distribution, scheduling_method, and skip_reset_qubits, allowing customization for specific hardware. For example, the XY 4sequence applies X,Y,−X,−Y pulses, mitigating low-frequency noise with a cycle time shorter than the coherence time. The benefits are important: DD can extend coherence times by 20–50% on devices like ibm_aachen (Heron r2, 200–300 µs), enabling deeper circuits (e.g., 250K CLOPS) to execute before decoherence dominates [17] [16]. The effectiveness depends on precise timing and hardware calibration, and it is less effective against gatespecific errors, necessitating complementary techniques like Pauli Twirling. 2) Error Mitigation - Pauli Twirling: Pauli Twirling (PT) is an error mitigation strategy that randomizes gate errors to approximate them as stochastic Pauli channels, simplifying error correction and improving measurement fidelity[14]. This technique is vital for Shor’s algorithm, where the QFT and modular exponentiation circuits are very sensitive to coherent errors that accumulate over thousands of gates on NISQ devices. Current IBM Quantum hardware (e.g., ibm_brussels, Eagle r3) with error rates of 2.86 ×10−3benefits from PT to mitigate systematic errors. In the code (phase_estimation.py), PT is implemented via TwirlingOptions in _set_sampler_params, with parameters like enable_gates=True, enable_measure,num_randomizations, shots_per_randomization, and strategy. For each gate (e.g., CNOT, Rz), random Pauli operators (I,X,Y,Z) are inserted before and after, transforming coherent errors into depolarizing noise. This randomization averages out phase errors, aligning the error model with zero-noise extrapolation or mitigation protocols. The advantage is improved fidelity for expectation values, critical for extracting the period rfrom QFT measurements. On ibm_strasbourg (Eagle r3, 3.84 × 10−3error), twirling can reduce effective error rates by 10–30%, enhancing the success probability of Shor’s algorithm for N= 15. However, it increases runtime due to additional shots (e.g., 1024 per randomization) and may not fully mitigate high-frequency noise or crosstalk, limiting its impact on large-scale circuits. Combining DD and twirling, as supported by the code, offers a robust mitigation strategy for NISQ constraints. 3) Optimization Techniques (Depth Reduction): SABRE Transpilation: SABRE is a transpilation optimization technique that reduces circuit depth by optimizing qubit mapping and routing, essential for executing Shor’s algorithm on constrained quantum hardware[15]. Shor’s algorithm circuits, with O((log2N)3)gate complexity, face depth limitations on devices like ibm_torino (Heron r1, 210K CLOPS, 150–200 µs coherence), where gate errors (1.03 ×10−3) accumulate rapidly. The transpilation.py module implements SABRE via the SabreLayout pass within Qiskit’s pass_manager. Configured in _set_transpile_optimization_params (phase_estimation.py), SABRE uses parameters like sabre_optimization,max_iterations, layout_trials, and swap_trials. It dynamically maps logical qubits to physical qubits based on a coupling map (e.g., ibm_aachen’s 156-qubit topology), minimizing the number of SWAP gates required for two-qubit interactions (e.g., CNOTs in 9 N= 15 and a= 4 N= 15 and a= 14 (a) (b) (c) (d) (e) (f) (g) Fig. 6: Summary Table of Optimization and Error Mitigation Results (see text for details) 16 TABLE IX: Circuit depth, total gates, histogram noise, and period-finding outcomes for N= 15,a= 4,14 on ibm_aachen QPU, with error mitigation, approximation degree 0.7, and optimization level 3. Techniques include Dynamical Decoupling (DD), Pauli Twirling (PT), SABRE optimization (SO), and all their combinations. Noise indicates histogram noise level Table 6aOptimization Technique Circuit Depth Total Gates Histogram Noise Period Outcome a) 4 DD 157 1121 50.3% Period found b) 4 PT 157 1121 10.5% Period found c) 4 SO 201 1306 7.7% Period found d) 4 DD + PT 157 1121 99.5% Period found e) 4 DD + SO 201 1306 41.9% Period found f) 4 PT + SO 201 1306 29.0% Period found g) 4 DD + PT + SO 201 1306 98.6% Period found a) 14 DD 99 694 33.4% Period found b) 14 PT 99 694 11.6% Period found c) 14 SO 108 776 6.2% Period found d) 14 DD + PT 99 694 96.9% Period found (wrong) e) 14 DD + SO 108 776 46% Period found f) 14 PT + SO 88 686 15.3% Period found g) 14 DD + PT + SO 88 686 97.6% Period found (wrong) we successfully factored up to N= 26 transpiling and executing the circuits for all values of acoprime of N, and finding the correct period rfor at least one value of a11. The use of identical optimization parameters (optimization level 3, approximation factor 0.7, and Pauli Twirling) allows us to isolate the effects of problem size on quantum circuit complexity and performance. Using these parameters, we also successfully factored N= 51, N= 55,N= 57, and N= 65, by searching at least a value of aamong all the coprimes of N, for which the correct period ris found, as suggested in Shor’s algorithm. In the next sections a more efficient approach is presented that allows to factor even larger N’s. For N= 51 all values of acoprime of Nlead to feasible circuits. Results are shown in Table X. We believe this is due to the fact that all values of rare powers of 2. For larger values of Nthe fraction of the total number of circuits that can be transpiled rapidly decreases. Results for for N= 55,N= 57, and N= 65 are presented in Tables XI,XII, and XIII. Figures 7,8, and 9present the actual circuit layouts and their probability distributions, respectively. •Circuit Feasibility and Scaling: The transition 11Note: Shor’s algorithm doesn’t need to work for all values of ait only needs to work for a sufficient fraction of them to be practical. When factoring an integer N, a random value of acoprime to Nhas to be selected and we have to try to find the period of the function f(x) = axmod N. The algorithm succeeds in factoring Nwhen this period ris even and ar/21(modN), because then gcd(ar/21, N) gives a non-trivial factor. The crucial theoretical result is that for any odd composite Nthat is not a prime power, at least half of the values of acoprime to Nwill yield a period that leads to successful factoring. This means that there are at least a 50% chance of success on each attempt. In practice, if one choice of a fails, another random ais selected and the process is repeated. The expected number of attempts is at most 2, and with high probability it will succeed in a few tries. from N= 55 to N= 65 reveals significant scaling challenges in quantum factoring. For N= 55, 7 out of 40 coprime values (a= 12,21,23,32,34,43,54) produced feasible circuits that could be transpiled and executed successfully, requiring between 14,819 and 81,641 total gates and circuit depths between 2,708 and 14,708. For N= 57, only 3 out of 36 coprimes (a= 20,37,56) were feasible, with gate counts between 29,938 and 40,428 and depths between 5,723 and 7,494. For N= 65, only 2 out of 47 coprimes (a= 14,51) were feasible, with gate counts of 161,286 and 166,981, and circuit depths of 29,280 and 30,129, respectively. Non-feasible cases for N > 65 required over 2 million gates, a ∼13–15×increase over feasible cases, highlighting a bimodal distribution and extreme dependence on the choice of a. •Gate Distribution Patterns: For N= 55, RZ and CZ gates dominated complexity, with no RX or RZZ gates. For N= 57, SX gates were prominent alongside CZ gates. For N= 65, a stark shift occurred: no RZ or SX gates were used, with CZ gates (e.g., 62348 for a= 14, 64816 for a= 51) and RX gates (e.g., 29627 for a= 14, 30524 for a= 51) dominating, indicating a change in circuit synthesis strategy or basis gate decomposition on the ibm_aachen QPU. •Circuit Depth Analysis: Circuit depths for N= 55 ranged from 2,708 to 14,708 (average ∼11,600), while N= 57 feasible cases had lower depths (5,723–7,494). For N= 65, feasible cases had significantly higher depths (29,280 for a= 14, 30,129 for a= 51), reflecting increased circuit complexity despite fewer feasible cases, possibly 17 due to larger modular exponentiation requirements. •Noise Performance Comparison: Histogram noise levels for N= 55 varied widely (10.3%–55.2%), while N= 57 feasible cases showed better consistency (13.7%–19.3%). For N= 65, noise levels were higher (56.3% for a= 14, 55.3% for a= 51), suggesting that larger Nvalues increase susceptibility to decoherence and gate errors, even in feasible cases, likely due to deeper circuits. •Success Rate: Success rates decreased with increasing N: 17.5% (7/40) for N= 55, 8.3% (3/36) for N= 57, and only 4.3% (2/47) for N= 65. The consistent r= 2 period finding in successful N= 57 and N= 65 cases, compared to mixed r= 2 and r= 4 for N= 55, suggests that shorter periods produce more favorable quantum interference patterns, which will be used later to optimize the overall algorithm, although the reduced number of feasible avalues for N= 65 indicates tighter constraints. These results highlight the exponential scaling challenges in near-term quantum factoring. The drastic reduction in feasible cases from N= 55 to N= 57 and N= 65 confirms current hardware limitations, with N= 65 pushing the boundaries of what is achievable running Shor’s algorithm12. The sensitivity to acould be leveraged by developing classical pre-processing to identify tractable cases before transpilation, as explored in later sections. This fact and further tunning of the circuit optimizations will show how to beat this record much further. 3) Circuit Performance Dependence on the particular QPU: We started performing most of our studies with IBM Kyiv QPU, having achieved already the capability to factorize N= 51 (see Fig. 10). However, on March 18th, 2025, the machine was no longer available, and we had to use other machines. Unfortunately, despite having similar specifications and theoretical performance, it was not possible to reproduce the results. It was hypothesized that this was due to less noisy qubits, or mapping into the circuit layout, but the coherence times are very similar. To address this issue, ibm_Kyiv’s circuit layout was literally ”implanted” in four other computers (ibm_Brisbane,ibm_Sherbrooke, ibm_Brussels and ibm_Strasbourg), which means that the circuit topology was exactly the same. Despite this, results were dramatically different (Table XIV ). ibm_Kyiv’s depth of 5,309 is approximately five times smaller than the others (25,325–25,326), which are nearly identical. This suggests that ibm_Kyiv’s circuit 12Recall that the largest value of Nthat has been factorized in an IBM NISQ is N= 21 [5]. (a) N= 55, a= 12 (a) N= 55, a= 21 (a) N= 55, a= 23 (a) N= 55, a= 32 (a) N= 55, a= 34 (b) N= 55, a= 43 (c) N= 55, a= 54 Fig. 7: Summary Table of Probability Distributions for N= 55 and a) a= 12, b) a= 21, c) a= 23, d) a= 32, e) a= 34, f) a= 43, and g) a= 54, the only 7 a’s for which the total number of gates is 50,000 gates. r= 4 in all cases, except r= 2 for a= 21 and a= 54. For a= 54 the correct non-trivial factors were found 18 TABLE X: Quantum Circuit Statistics for N= 51 and all aCoefficients, transpiled and successfully executed in an ibm_aachen QPU with Optimization Degree = 3, Approximation Factor = 0.7, and Pauli Twirling. Number of RX, RZ, SX, X, CZ, and RZZ gates not shown for the sake of simplificy. a r Total Gates Circuit Depth Histogram Noise (%) 2 4 120545 21999 56.0% 4 4 7993 14377 54.1% 5 16 159739 28752 54.0% 7 16 167981 30384 10.4% 8 8 125261 22862 12.3% 10 16 165820 29740 50.8% 11 16 162097 29364 8.5% 13 4 83114 15222 55.7% 14 16 165516 30046 17.4% 16 2 40912 22464 15.6% 20 8 124846 22762 20.5% 22 16 168841 30405 11.0% 23 16 165424 30063 6.8% 25 8 121014 21991 53.1% 26 8 120188 21737 13.8% 28 16 168554 30523 12.3% 29 16 163706 29340 7.9% 31 16 167359 30516 15.3% 32 8 119874 21498 17.3% 35 2 36800 6641 19.8% 37 16 163522 29621 8.6% 38 4 78433 14041 13.8% 40 16 164321 29704 15.9% 41 16 159830 28623 13.7% 43 8 124186 22325 9.8% 44 16 169469 30758 16.0% 46 16 162817 29426 5.0% 47 4 81351 14835 16.0% 49 8 122173 21871 9.4% 50 2 22227 4001 22.4% was optimized through a far more efficient transpilation, reduced swap gate insertion, or a topology-aligned qubit mapping. The uniformity among the other QPUs indicates a more standardized compilation approach, with minor variations possibly due to QPU-specific connectivity. In terms of the total number of gates ibm_Kyiv’s 46,447 gates are significantly lower than ibm_Strasbourg’s 189,416 to ibm_Brisbane’s 209,487. The 10–15% variation among the other QPUs suggests differences in transpilation or routing overhead. RZ gates dominate (˜ 40% of the total), followed by SX ( 25–30%), ECR ( 12–15%), and X ( 5%) gates, reflecting the structure of Shor’s algorithm, which relies heavily on quantum Fourier transforms and modular arithmetic. All QPUs use the Eagle r3 architecture with the basis gate set {ecr, id, rz, sx, x}. Thus, differences arise from transpilation, and are not due to the set of native gates. In terms of the error rates ibm_Sherbrooke and ibm_Brisbane have the lowest error rates, while (a) N= 57, a= 20 (b) N= 57, a= 37 (c) N= 57, a= 56 Fig. 8: Circuit layout and Probability Distributions for N= 57 and a) a= 20, b) a= 37, and c) a= 56. (a) N= 65, a= 14 (b) N= 65, a= 51 Fig. 9: Circuit layout and Probability Distributions for N= 65 and a) a= 14, and b) a= 51. ibm_Strasbourg has the highest. ibm_Kyiv’s moderate errors do not explain either its compact circuit. Among the others, ibm_Sherbrooke’s lower layered error may reduce its gate count by minimizing swap gates. ibm_Strasbourg’s higher errors may increase routing overhead, but its has a lower gate count which suggests that transpilation dominates. Coherence times and number of qubits are pretty uniform 19 TABLE XI: Quantum Circuit Statistics for N= 55 and Different aCoefficients, transpiled and successfully executed in an ibm_aachen QPU with Optimization Degree = 3, Approximation Factor = 0.7, and Pauli Twirling ar RX Gates RZ Gates SX Gates X Gates CZ Gates RZZ Gates Total Gates Circuit Depth Histogram Noise 2 20 0 157208 207123 1520 99718 0 540867 97577 – 3 20 0 128448 168759 1241 82130 0 443273 80240 – 4 10 0 157776 207743 1470 99736 0 541719 97416 – 6 10 0 160202 210171 1365 100747 0 548277 98428 – 7 20 0 157059 205616 1432 98474 0 536858 96030 – 8 20 0 160987 211291 1505 101461 0 551608 99152 – 9 10 0 160953 209666 1505 100006 0 547208 97585 – 12 4 0 22424 29764 220 14317 0 77707 13983 55.2 % 13 20 0 160229 210390 1473 101030 0 549068 98727 – 14 10 0 158188 208378 1448 100295 0 543846 97913 – 16 5 0 158489 209443 1436 100980 0 546469 98659 – 17 20 0 157520 207896 1534 100343 0 543092 98110 – 18 20 0 161023 211001 1547 101239 0 550867 98830 – 19 10 0 161848 212553 1484 102212 0 554871 99942 – 21 2 0 11768 15836 95 7710 0 41350 7548 11.5 % 23 4 0 23708 31171 216 15029 0 81641 14708 10.3 % 24 10 0 159993 210273 1403 100882 0 548393 98516 – 26 5 0 161906 211698 1444 101325 0 552441 98839 – 27 20 0 157804 207533 1484 99893 0 542094 97536 – 28 20 0 163940 214697 1360 102870 0 559992 100490 – 29 10 0 161899 212874 1442 102197 0 555055 99831 – 31 5 0 161352 210506 1560 100813 0 549877 98352 – 32 4 0 23998 31624 225 15157 0 82603 14827 11.5 % 34 2 0 10919 14591 139 7086 0 38299 6943 17.4 % 36 5 0 160071 211031 1521 101569 0 550640 99207 – 37 20 0 158598 207544 1456 99212 0 541350 96694 – 38 20 0 159748 210176 1505 101262 0 548950 98929 – 39 10 0 158917 207965 1533 99688 0 543025 97373 – 41 10 0 158216 209086 1391 100759 0 545424 98472 – 42 20 0 158271 207320 1432 99373 0 541137 96949 – 43 4 0 23479 31035 213 14947 0 81098 14658 55.0 % 46 10 0 162943 214242 1395 102824 0 558556 100490 – 47 20 0 161462 212062 1455 101386 0 552435 98950 – 48 20 0 159982 210329 1403 100979 0 548605 98605 – 49 10 0 163255 213794 1479 102408 0 557557 99896 – 51 10 0 165191 215721 1427 103214 0 562688 100643 – 52 20 0 158870 209726 1443 100878 0 546876 98600 – 53 20 0 163249 213981 1459 102752 0 558234 100354 – 54 2 0 4170 5612 31 2793 0 14819 2708 17.3 % (100–150 µs) and number of qubit is the same (127), so they cannot contribute to observed differences either. CLOPS ranges from 150K to 220K, with Kyiv’s having a moderate CLOPS (200K). Therefore it does not explain either its efficiency. After all these analyses, we discard all possible causes except that ibm_Kyiv’s more compact circuit may be likely due to a very optimized transpilation. This highlights the importance of transpilation with NISQ and shows directions to enhance circuit efficiency. 4) Optimizing of Shor’s Algorithm for Real Quantum Hardware: Classical Preprocessing for Maximum Factorable N, the Optima aValues: In this section, our primary objective is to investigate how far we could extend the factorization capabilities of Shor’s algorithm beyond the brute force approach of trying avalues randomly. To do this we use the most powerful QPU available to us ibm_aachen’s quantum processor, which features 156 qubits and employs a Heron r2-type QPU. Ideally, when (a) N= 51,a= 4 Fig. 10: Circuit layout and Probability Distributions for N= 51 and a) a= 4 in Kyiv QPU. Same circuit layout was implanted in Brisbane, Sherbrooke, Brussels and Strasbourg QPUs, but it did not work. executed on quantum hardware, Shor’s algorithm would randomly test multiple values of auntil it identifies one that yields a useful period, ideally, the smallest possible period (r= 2). However, due to current hardware 20 TABLE XII: Quantum Circuit Statistics for N= 57 and Different aCoefficients, transpiled and successfully executed in an ibm_aachen QPU with Optimization Degree = 3, Approximation Factor = 0.7, and Pauli Twirling ar RX Gates RZ Gates SX Gates X Gates CZ Gates RZZ Gates Total Gates Circuit Depth Histogram Noise 2 18 0 150058 210092 1511 100605 0 537899 98288 – 4 9 0 155039 216480 1468 103706 0 554437 101234 – 5 18 0 154044 214057 1415 102035 0 547789 99487 – 7 3 0 125055 174430 1836 84400 0 450081 82001 – 8 6 0 147436 207160 1639 99530 0 529853 97120 – 10 18 0 145709 202668 1556 96995 0 520263 94657 – 11 6 0 145428 205146 1480 99044 0 525190 96693 – 13 18 0 153616 214271 1421 102386 0 548300 99927 – 14 18 0 154590 215664 1461 103191 0 552086 100785 – 16 9 0 155217 216668 1494 103521 0 554360 101081 – 17 18 0 151090 211148 1372 100937 0 540289 98420 – 20 2 0 9929 14399 96 7117 0 37146 6992 19.3% 22 18 0 153095 213698 1552 102085 0 546836 99608 – 23 18 0 153045 213893 1470 102406 0 547592 99954 – 25 9 0 151890 210513 1431 100419 0 539544 97885 – 26 6 0 146254 205214 1466 98969 0 525805 96445 – 28 9 0 150352 210133 1413 100494 0 537867 98127 – 29 18 0 151741 211981 1372 101196 0 542085 98696 – 31 6 0 148426 209586 1595 101553 0 536955 99212 – 32 18 0 152732 214500 1487 102834 0 548719 100417 – 34 18 0 152184 211835 1351 100990 0 542202 98520 – 35 18 0 135616 188128 1701 90164 0 484400 87693 – 37 2 0 11007 15716 112 7649 0 40428 7494 17.0% 40 18 0 150210 211526 1509 101733 0 541489 99411 – 41 18 0 153221 213952 1483 102328 0 547624 99900 – 43 9 0 150764 211621 1456 101718 0 542089 99494 – 44 18 0 149231 209170 1478 100117 0 535248 97709 – 46 6 0 144019 203294 1592 98177 0 520511 95809 – 47 18 0 153145 214086 1469 102691 0 548423 100290 – 49 3 0 149566 210474 1469 101475 0 538484 99012 – 50 6 0 149295 209577 1558 101115 0 536870 98706 – 52 18 0 151506 212377 1533 101885 0 543779 99671 – 53 18 0 152997 214160 1495 102465 0 547775 99938 – 55 9 0 149761 211135 1626 101400 0 540333 99115 – 56 2 0 7840 11695 61 5839 0 29938 5723 13.7% limitations, this approach is highly time-consuming and impractical. In the previous sections, we adopted a methodological approach: for each value of Nto be factored, we generated all coprime values of N, and transpiled each of them in the corresponding quantum circuit. Once we identified the athat produced the lowest circuit depth number, we executed that circuit on the real ibm_aachen quantum processor. However, as Nincreased, the process became extremely timeconsuming. For example, when factoring N= 57, we evaluated 35 coprime values of a, specifically (2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 20, 22, 23, 25, 26, 28, 29, 31, 32, 34, 35, 37, 40, 41, 43, 44, 46, 47, 49, 50, 52, 53, 55, 56). On average, each transpilation required approximately two hours, resulting in a total of nearly 70 hours to identify the optimal a. This exhaustive process becomes increasingly unsustainable for larger numbers, particularly given current hardware constraints. Despite working with a 156-qubit machine, the practical limitations of circuit depth and gate fidelity meant that we remained restricted to relatively small integers, far below those relevant to real-world cryptographic applications. To overcome this bottleneck, based on our previous experience of transpiling all the admissible values of a, we realized that when searching for the athat produces the lowest circuit depth a pattern appeared. As you can observe in tables XI,XII, and XIII, the bolded lines mark the transpiled circuits that possess the lowest circuit depth, and for those the corresponding rvalue is always 2. So our objective was come up with a way to know in advance those avalues that produce periods (r) equal to 2, and that as a consequence generate circuits with the lowest circuit depth. To achieve this, we developed a classical pre-processing algorithm designed to identify values of athat are more likely to yield the smallest possible period (r= 2). Specifically, the script computes the order r, the smallest positive integer such that ar≡1 mod N, for all values of a that are coprime to a given N. For example, when factoring N= 111, the algorithm identified three values of a (38, 73, and 110) that yielded the minimal period r= 2.These values were then prioritized for circuit transpilation and execution. The pseudocode 4 outlines the steps of the Order F inding Algorithm. 21 TABLE XIII: Quantum Circuit Statistics for N= 65 and Different aCoefficients, transpiled and successfully executed in an ibm_aachen QPU with Optimization Degree = 3, Approximation Factor = 0.7, and Pauli Twirling ar RX Gates RZ Gates SX Gates X Gates CZ Gates RZZ Gates Total Gates Circuit Depth Histogram Noise 2 12 501779 0 0 801352 1067092 3814 2745207 494956 – 3 12 493606 0 0 790288 1051119 3712 2702897 486471 – 4 6 500675 0 0 801224 1064951 3557 2740934 493856 – 6 12 494715 0 0 790745 1052929 3690 2707979 488012 – 7 12 388444 0 0 622867 821824 5455 2133977 380830 – 8 4 62192 0 0 98891 131754 448 339705 61398 – 9 6 494455 0 0 790234 1051727 3722 2706468 487558 – 11 12 493703 0 0 786674 1048440 3528 2698068 486934 – 12 4 60169 0 0 94070 126971 472 326722 59410 – 14 2 29627 0 0 46570 62348 255 161286 29280 56.3 % 16 3 499119 0 0 801622 1063634 3707 2737053 492187 – 17 12 490198 0 0 785547 1043214 3586 2684924 483500 – 18 4 59563 0 0 94725 126267 481 325607 58805 – 19 12 492384 0 0 788748 1047422 3739 2697183 485807 – 21 4 62024 0 0 98082 131503 502 338739 61248 – 22 12 487345 0 0 774785 1020404 4780 2648215 480529 – 23 12 483761 0 0 765402 1010165 4824 2623290 477198 – 24 12 493646 0 0 778364 1030721 4825 2673210 486868 – 27 4 59479 0 0 91739 122869 586 319407 58768 – 28 12 488873 0 0 771352 1021096 4715 2648383 482237 – 29 6 491283 0 0 779478 1028541 4639 2667700 484345 – 31 4 60512 0 0 93209 125114 528 324487 59770 – 32 12 408765 0 0 650431 853605 5862 2227898 401711 – 33 12 488508 0 0 772836 1020430 4684 2648578 481685 – 34 4 60920 0 0 96113 126928 596 329928 60047 – 36 6 498448 0 0 791834 1043665 5019 2707576 491257 – 37 12 502659 0 0 799816 1053434 4853 2732195 495348 – 38 4 59175 0 0 92340 122867 613 319439 58476 – 41 12 500031 0 0 794740 1048339 4698 2717560 492778 – 42 12 493457 0 0 784074 1033628 4836 2680857 486533 – 43 12 498278 0 0 794292 1059431 3758 2724257 491618 – 44 4 61510 0 0 97481 130095 463 335401 60707 – 46 12 488687 0 0 782823 1040982 3638 2678352 481904 – 47 4 58237 0 0 92234 123235 516 317657 57527 – 48 12 496332 0 0 794592 1057297 3639 2718189 489201 – 49 6 489532 0 0 783312 1042643 3556 2681457 482625 – 51 2 30524 0 0 48562 64816 211 166981 30129 55.3% 53 4 59830 0 0 94526 126313 549 326108 59108 – 54 12 490829 0 0 776372 1025957 4816 2661871 484122 – 56 6 496250 0 0 785113 1037350 4851 2691157 489460 – 57 4 60684 0 0 93657 125149 515 325409 59921 – 58 12 489736 0 0 768783 1020456 4908 2647082 483212 – 59 12 490126 0 0 777376 1026004 4789 2661668 483385 – 61 3 498438 0 0 795187 1045814 4761 2712417 491412 – 62 12 493043 0 0 781927 1031655 4775 2676686 486280 – 63 12 497617 0 0 791111 1042877 4781 2704992 490742 – 64 2 29012 0 0 45589 60240 262 157026 28587 – This optimization significantly reduced the number of quantum circuits that needed to be transpiled and executed. As a result, we were able to factor larger numbers, such as 111, 115, 123, directly on the hardware, demonstrating a more scalable and efficient approach to implement Shor’s algorithm within current technological limits. It is worth mentioning that for these values of N, SABRE optimization had to be activated in order to potentially reduce the circuit depth. Using the above technique, Fig. 11 shows the quantum circuit, circuit layout, and the probability distributions for the maximum value of Nwe have been able to factor N= 123 for a= 40 and a= 80 with and without Pauli Twirling and with and without SABRE optimization. An interesting fact that has been observed in all the optimum circuits analyzed is that all the exponentiation blocks except the first one are equal to [1modN](e.g. Fig. 11 (a) ) which is probably the reason of the significantly smaller number of gates, and significantly improved performance. VII. DISCUSSION The implementation and evaluation of Shor’s algorithm presented in this report illustrates its theoretical power and its current practical limitations. While Shor’s algorithm promises exponential speedup over classical algo22 (a) N= 123, a= 40 (b) N= 123,a= 83 (c) Quantum circuit N= 123,a= 40 (for a= 83 only 1st exponentiation block changes) (d) a= 40, Pauli Twirling, SABRE Optimization (e) a= 83, Pauli Twirling, SABRE Optimization (f) a= 40, No Pauli Twirling, No SABRE Optimization (g) a= 83, No Pauli Twirling, No SABRE Optimization (h) a= 40, Pauli Twirling, No SABRE Optimization (i) a= 83, Pauli Twirling, No SABRE Optimization (j) a= 40, No Pauli Twirling, SABRE Optimization (k) a= 83, No Pauli Twirling, SABRE Optimization Fig. 11: Quantum circuit (top), and probability distributions for N= 123 and a= 40 (left column) and a= 83 (right column): a-b) with PT and SO, c-d) without PT, without SO, e-f) with PT and without SO, and g-h) without PT and with SO. 23 TABLE XIV: Comparison of Circuit Parameters Across Different QPUs For N= 51, and a= 4 QPU Two-qubit ECR Gates Single-qubit RZ Gates Single-qubit SX Gates Single-qubit X Gates Total Gates Circuit Depth Kyiv 5502 21 554 14 155 621 46 447 5309 Brisbane 25 741 95 992 61 864 4710 209 487 25 326 Sherbrooke 25 743 91 728 59 566 5897 204 124 25 326 Brussels 25 740 96 312 59 777 6296 209 298 25 325 Strasbourg 25 743 81 894 49 707 10 896 189 416 25 326 Algorithm 4 Order Finding Algorithm Input: Integer N, base a Output: Order ror None for afrom 2 to Ndo if gcd(a, N) = 1 (check if a,Ncoprime) then v= 1 for rfrom 1 to Ndo v= (v×a) mod N(update vwith mod mult) if v= 1 (detect order) then Return:r else(no order yet) Return: None end if if r=None (process valid r)then if r < rmin (find smallest order) then rmin =r p= [(a, r)] else(handle other orders) r=rmin p.append((a, r)) end if end if end for end if end for Return: Failure (no order found) rithms for integer factorization —in principle a direct threat to RSA-based cryptographic systems— our study shows how far current quantum hardware remains from making this promise true. A. Practical Feasibility on Current Quantum Hardware Factoring ”small” numbers, such as N= 15,N= 51, and up to N= 123, serves as a practical benchmark for assessing quantum circuit performance on NISQ devices. Our experiments across IBM quantum processors (ibm_brisbane,ibm_torino,ibm_aachen, and ibm_kyiv) reveal significant variations in circuit performance due to differences in hardware architecture, qubit connectivity, and transpilation efficiency. Notably, the Heron-based ibm_aachen processor consistently outperformed Eagle-based processors (e.g. ibm_brisbane,ibm_torino) due to lower gate error rates (0.375% for CZ gates), improved coherence times (up to 450 µs with dynamical decoupling), and enhanced qubit connectivity. For instance, ibm_aachen achieved successful factorization of N= 123 with a= 40 and a= 83, as shown in Fig. 11, with circuit depths around 100,000 and total gate counts exceeding 500,000, even with optimizations like SABRE and error mitigation technique like Pauli Twirling. In terms of noise, it is worth noting the the best performance are achieved for a= 40 when neither Pauli Twirling, nor SABRE optimization are used (2.6%), and for a= 83 when only SABRE optimization is used (4.1%). The combined effect of both leads to high noise levels (22.2% and 45.6%). SABRE Optimization does not seem to show a consistent correlation with lower noise, i.e. the lowest noise (2.6%) occurs without SABRE (Fig. 11f), while higher noise levels (80.6%) also appear without it (Fig. 11i), suggesting SABRE’s impact on noise is context-dependent. With respect to the case when Pauli Twirling is used, the lowest noise -although not the lowest absolute valueis 13.5% (Fig. 11h), but also the highest one 80.6% (Fig. 11i), so its potential noise mitigation effects may be arguable in this case. This highlights that it is difficult to infer general rules for circuit optimization in terms of noise. Circuit depths remain a critical limitation, and even for small N, circuits often exceed 50,000 in depth and include tens of thousands of gates, pushing the coherence limits of current hardware (100–150 µs). Error mitigation techniques, such as Pauli Twirling and Dynamical Decoupling, improved success rates (e.g., 17.5% for N= 55, 8.3% for N= 57, 4.3% for N= 65), but were insufficient for fault-tolerant execution. The ibm_kyiv QPU, before its retirement on March 18, 2025, demonstrated superior performance for N= 51 with a circuit depth of 1,383 and 14,256 gates, likely due to highly optimized transpilation and topologyaligned qubit mapping. Replicating this circuit layout on other QPUs (ibm_brisbane,ibm_sherbrooke, ibm_strasbourg, and ibm_brussels) resulted in significantly higher depths (around 25,326) and gate counts (189,416–209,487), underscoring the critical role of transpilation in performance. The best configurations on ibm_aachen QPU (e.g., a= 14 for N= 15, although it leads to the trivial factorization) still yielded more than 1000 gates, indicating that fault-tolerant quan24 tum computing is essential to progress beyond small N. B. Optimization Trade-Offs and Algorithm Tuning Our transpilation benchmarking reveals that optimization level and approximation degree significantly impact circuit complexity. Studies on dynamical decoupling and pulse-level optimizations on IBM quantum computers, as explored in [17], provide further insights into these effects. Higher optimization levels (e.g., level 3) reduced circuit depth by up to 40%, while an approximation degree of 0.7 decreased gate counts by approximately 25%. However, these optimizations often compromise fidelity, as seen in cases where extreme approximations led to unrealistic circuit depths (e.g., circuit depth=0). The choice of the base aalso affects significantly the performance. For N= 123, values like a= 40 and a= 83 yielded circuits with minimal periods (r= 2), reducing modular exponentiation complexity and improving success rates. Our classical pre-processing algorithm, outlined in Algorithm 4, identifies such optimal avalues by computing the order rfor coprime a, significantly reducing the number of quantum circuits needed. This enabled successful factorization of larger numbers like N= 111,N= 115, and even N= 123 on 156-qubit hardware. SABRE optimization proved particularly effective for larger N, reducing histogram noise (e.g., 10.7% with SABRE vs. 69.0% without for N= 123,a= 40). Pauli Twirling was essential for increasing maximum circuit depth, though, for smaller circuits, SABRE Optimization introduced a minor overhead in terms of number of gates. These findings highlight the delicate balance between circuit depth, gate fidelity, and success probability, necessitating careful tuning for NISQ devices. C. Scalability Limitations and Requirements For cryptographically significant integers (e.g., 2048-bit RSA keys), Shor’s algorithm would require millions of gates, thousands of qubits, and long coherence times, all of which exceed current NISQ capabilities by far. Even with the best IBM QPU (ibm_aachen), we are orders of magnitude away from such requirements. The exponential models fitted to our transpilation benchmarking data predict that even a small increase in Nor less favorable choices of awill dramatically increase circuit requirements. This reinforces the need for significant improvements in fault-tolerant QC, quantum error correction (e.g., surface codes), hardware improvements in qubit connectivity and fidelity, etc. Our theoretical analysis in Appendix III estimates maximum circuit depths of 97–119 for IBM QPUs under optimistic error tolerances (50% over 100 layers), yet observed depths were much higher. This suggests undocumented error mitigation or higher tolerances in IBM systems. The modular exponentiation step remains the primary bottleneck, contributing the majority of gates and depth due to its O((log2N)3)complexity. Optimizations like precomputing a2jmod Nand approximate QFT reduce complexity, but cannot bridge the gap to cryptographic scales without fault-tolerant quantum computing. Surface codes, requiring thousands of physical qubits per logical qubit, are currently infeasible on 127–156-qubit systems. VIII. CONCLUSIONS AND FUTURE RESEARCH LINES This study has successfully implemented Shor’s algorithm on IBM quantum hardware for small integers up to N= 123, demonstrating its theoretical feasibility on NISQ devices while highlighting significant practical limitations. The Heron-based ibm_aachen QPU outperformed Eagle-based processors due to lower error rates, longer coherence times, and better connectivity, achieving factorizations like N= 123 with optimized a values (a= 40,83). However, the algorithm’s scalability is severely constrained by circuit depth, gate count, and noise, with modular exponentiation identified as the primary bottleneck due to its cubic complexity in log2N. Key findings include the critical role of transpilation and optimization techniques. SABRE optimization and Pauli Twirling significantly reduced noise and improved success rates, particularly for larger N, while classical pre-processing to select optimal avalues minimized circuit complexity, enabling factorization of N= 111, N= 115, and N= 123. Despite these advancements, current hardware remains orders of magnitude away from factoring cryptographically relevant numbers (e.g., 2048-bit RSA keys), which would require millions of gates and thousands of error-corrected qubits. The ibm_kyiv QPU’s exceptional performance for N= 51 suggests that highly optimized transpilation and topology alignment are crucial for maximizing NISQ performance. The implementation and benchmarking of Shor’s algorithm on current quantum hardware have revealed critical insights into its practical execution and limitations. The following recommendations and key observations are derived from our study to guide future implementations and optimizations of Shor’s algorithm on NISQ devices: •Circuit Depth and Gate Count Constraints: Empirical analysis indicates that quantum circuits for Shor’s algorithm perform reliably only when the circuit depth is ≤100,000 and the total gate count is ≤500,000. These thresholds reflect the coherence time and gate fidelity limitations of current NISQ hardware, necessitating careful circuit optimization to remain within these bounds. •Error Mitigation and Optimization Strategies: Pauli Twirling is essential to extend the maximum achievable circuit depth by mitigating coher25