# How to factor 2048 bit RSA integers with less than a million noisy qubits Source: https://arxiv.org/html/2505.15917v1 Kiwi News: https://news.kiwistand.com/stories/How-to-factor-2048-bit-RSA-integers-with-less-than-a-million-noisy-qubits?index=0x68335606e179406b6b777570cd21ad2e9ba664495e7190b49cdbc013f5a293629ff01b99 Submitted by misha on 2025-05-25 UTC, 1 upvote, 1 comment ## Summary Craig Gidney estimates that a 2048 bit RSA integer could be factored in less than a week by a quantum computer with less than a million noisy qubits, down from the roughly 20 million qubits his 2019 co-authored estimate required. He attributes the qubit reduction mainly to approximate residue arithmetic from Chevignard, Fouque and Schrottenloher (2024), to yoked surface codes, and to magic state cultivation. The longer runtime comes mainly from more Toffoli gates and fewer magic state factories, though the paper reports a Toffoli count over 100 times lower than Chevignard, Fouque and Schrottenloher's 2024 method. The estimates assume a square grid of qubits, a 0.1% gate error rate, and a 1 microsecond surface code cycle time. ## Article Craig Gidney August 24, 2026 Abstract Planning the transition to quantum-safe cryptosystems requires understanding the cost of quantum attacks on vulnerable cryptosystems. In Gidney+Ekerå 2019, I co-published an estimate stating that 2048 bit RSA integers could be factored in eight hours by a quantum computer with 20 million noisy qubits. In this paper, I substantially reduce the number of qubits required. I estimate that a 2048 bit RSA integer could be factored in less than a week by a quantum computer with less than a million noisy qubits. I make the same assumptions as in 2019: a square grid of qubits with nearest neighbor connections, a uniform gate error rate of 0.1%0.1\%, a surface code cycle time of 1 microsecond, and a control system reaction time of 1010 microseconds. The qubit count reduction comes mainly from using approximate residue arithmetic (Chevignard+Fouque+Schrottenloher 2024), from storing idle logical qubits with yoked surface codes (Gidney+Newman+Brooks+Jones 2023), and from allocating less space to magic state distillation by using magic state cultivation (Gidney+Shutty+Jones 2024). The longer runtime is mainly due to performing more Toffoli gates and using fewer magic state factories compared to Gidney+Ekerå 2019. That said, I reduce the Toffoli count by over 100x compared to Chevignard+Fouque+Schrottenloher 2024. Data availability: Code and assets created for this paper are available on Zenodo [Gid25]. 1 Introduction 2 Methods 2.1 Approximate Residue Arithmetic 2.2 Approximate Period Finding 2.3 Ekerå-Håstad Period Finding 2.4 Arithmetic Optimizations 3 Results 3.1 Logical Costs 3.2 Physical Costs 4 Conclusion 5 Acknowledgments References A Detailed Mock-ups A.1 Reference Python Implementation of Algorithm A.2 Addition Operation A.3 Phaseup Operation A.4 Lookup Operation A.5 Frequency Basis Measurement Figure 1: Historical estimates, with comparable physical assumptions, of the physical qubit cost of factoring 2048 bit RSA integers. Includes overheads from fault tolerance, routing, and distillation. Results are from [Jon+12, Fow+12, OC17, GM19, GE21]. Results such as [Van+10] and [LN22] aren’t included because they target substantially different assumptions or cost models. Figure 2: Pareto frontiers achieved by this paper for the Toffoli and logical qubit cost of factoring nn bit RSA integers, for various values of nn. This paper uses notably fewer logical qubits than [GE21] (points labeled “GE21”) and notably fewer Toffolis than [CFS24] (points labeled “CFS24”). 1 Introduction In 1994, Peter Shor published a paper showing quantum computers could efficiently factor integers [Sho94], meaning the RSA cryptosystem [RSA78] wasn’t secure against quantum computers. Understanding the cost of quantum factoring is important for planning and coordinating the transition away from RSA, and other cryptosystems vulnerable to quantum computers. Correspondingly, since Shor’s paper, substantial effort has gone into understanding the cost of quantum factoring [Kni95, Bec+96, VBE96, Zal98, CW, Bea03, Zal06, Whi+09, Van+10, Fow+12, Jon+12, PG14, Eke16, HRS16, Gid17, OC17, GM19, Eke19, Eke21, GE21, LN22, MS22, Reg24, CFS24] (and more). A key metric is the number of qubits used by the algorithm, since this bounds the required size of quantum computer. Historically, there was no known way to factor an nn bit number using fewer than 1.5​n1.5n logical qubits [Zal06, Bea03, HRS16, Gid17]. Anecdotally, it was widely assumed that nn was the minimum possible because doing arithmetic modulo an nn bit number “required” an nn qubit register. May and Schlieper had shown in 2019 that in principle only a single output qubit was needed [MS22], but there was no known way to prepare a relevant output that didn’t involve intermediate values spanning nn qubits. In 2024, Chevignard and Fouque and Schrottenloher (CFS) solved this problem [CFS24]. They found a way to compute approximate modular exponentiations using only small intermediate values, destroying the anecdotal nn qubit arithmetic bottleneck. However, their method is incompatible with “qubit recycling” [PP00, ME99], reviving an old bottleneck on the number of input qubits. Shor’s original algorithm used 2​n2n input qubits [Sho94], but Ekerå and Håstad proved in 2017 that (0.5+ϵ)​n(0.5+\epsilon)n input qubits was sufficient for RSA integers [EH17]. So, by combining May et al’s result with Ekerå et al’s result, the CFS algorithm can factor nn bit RSA integers using only (0.5+ϵ)​n(0.5+\epsilon)n logical qubits [CFS24]. A notable downside of [CFS24] was its gate count. In [GE21], it was estimated that a 2048 bit RSA integer could be factored using 3 billion Toffoli gates and a bit more than 3​n3n logical qubits. Whereas [CFS24] uses 2 trillion Toffoli gates and a bit more than 0.5​n0.5n logical qubits. So [CFS24] is paying roughly 1000x more Toffolis for a 6x reduction in space. This is a strikingly inefficient spacetime trade-off. However, as I’ll show in this paper, the trade-off can be made orders of magnitude more forgiving. In addition to optimizing the Toffoli count of the algorithm, I’ll provide a physical cost estimate showing it should be possible to factor 2048 bit RSA integers in less than a week using less than a million physical qubits (under the assumptions mentioned in the abstract). The paper is structured as follows. In Section 2, I describe a streamlined version of the CFS algorithm. In Section 3, I estimate its cost. I first estimate Toffoli counts and logical qubit counts, and then convert these into physical qubit cost estimates accounting for the overhead of fault tolerance. Finally, in Section 4, I summarize my results. The paper also includes Appendix A, which shows more detailed mock-ups of the algorithm and its physical implementation. 2 Methods Table 1: Various bits of mostly-Python-inspired notation used in this paper. In this section, I present a variation of the CFS algorithm. The underlying ideas are the same, but many of the details are different. For example, I use fewer intermediate values and I extract the most significant bits of the result rather than the least significant bits. Beware that, as in [GE21], for simplicity, I will describe everything in terms of period finding against f⁡(e)=gemodNf(e)=g^{e}\bmod N (as in Shor’s original algorithm) despite actually intending to use Ekerå-Håstad-style period finding [EH17]. Ultimately everything decomposes into a series of quantum controlled multiplications, which fundamentally is what is actually being optimized, so what I describe will trivially translate. 2.1 Approximate Residue Arithmetic Consider an integer LL equal to the product of many ℓ\ell-bit primes from a set PP: L=∏p∈PpL=\prod_{p\in P}p (1) ∀p∈P:(len ​p=ℓ)∧IsPrime(p)\forall p\in P:(\text{len }p=\ell)\land\text{IsPrime(p)} (2) A simple way to do arithmetic on a value modulo LL is to store the value as a 2s complement integer in a len ​L\text{len }L bit register, and when operating on this register add or subtract multiples of LL as appropriate to keep it in the range [0,L)[0,L). This requires Ω⁡(log⁡L)\Omega(\log L) bits of storage. Residue arithmetic performs addition and multiplication modulo LL by separately performing the operations modulo each of the primes from PP. The final result can then be recovered using the Chinese remainder theorem. By the prime number theorem, PP can be chosen so that the bit length of primes in PP is ℓ=Θ⁡(log⁡log⁡L)\ell=\Theta(\log\log L). Therefore, using residue arithmetic introduces the possibility of arithmetic being exponentially more space efficient: |P|=Θ⁡(log⁡Llog⁡log⁡L)|P|=\Theta\left(\frac{\log L}{\log\log L}\right) (3) ℓ=Θ⁡(log⁡log⁡L)\ell=\Theta(\log\log L) (4) In the context of Shor’s algorithm, the key operation we want to compute is a modular exponentiation V=gemodNV=g^{e}\bmod N. Here gg is a randomly chosen classical constant in [2,N)[2,N), ee is a uniformly superposed m=O⁡(log⁡N)m=O(\log N) qubit integer, and NN is the number to factor. The computation of VV can be decomposed into a series of multiplications controlled by the qubits of ee: V=gemodN=(∏k=0m−1Mkek)modNV=g^{e}\bmod N=\left(\prod_{k=0}^{m-1}M_{k}^{e_{k}}\right)\bmod N (5) m=O⁡(log⁡N)m=O(\log N) (6) In the above equation, eke_{k} is the qubit at little-endian offset kk within the register ee and MkM_{k} is the classically precomputed constant Mk=g(2k)modNM_{k}=g^{(2^{k})}\bmod N (7) Because the factors of NN aren’t known, and because in practice those factors would be large, it’s not viable to perform residue arithmetic modulo NN. However, performing the arithmetic modulo any other modulus L≠NL\neq N creates an issue. If the accumulating product exceeds LL, then the register will wrap around. This shifts its value by a multiple of LL and, since LmodN≠0L\bmod N\neq 0, this offsets the tracked value relative to the true result mod NN. The following multiplications would then amplify this error out of control, ruining the computation. I fix this in a simple way, the same way Chevignard et al fix it, by picking LL to be larger than the largest possible product. This guarantees the wraparound issue never occurs: With this promise about the size of LL, we can rewrite the computation of VV to use residue arithmetic. Instead of computing the product directly, we can compute it using a dot product between each residue rjr_{j} (the product modulo a prime from LL) and its contribution factor uju_{j} (the smallest positive integer satisfying ujmodpi=δi,ju_{j}\bmod p_{i}=\delta_{i,j}): rj=(∏k=0m−1Mkek)modpjr_{j}=\left(\prod_{k=0}^{m-1}M_{k}^{e_{k}}\right)\bmod p_{j} (9) uj=(L/pj)⋅MultiplicativeInversepj​(L/pj)u_{j}=(L/p_{j})\cdot\text{MultiplicativeInverse}_{p_{j}}(L/p_{j}) (10) V\displaystyle V =(∏k=0m−1Mkek)modLmodN\displaystyle=\left(\prod_{k=0}^{m-1}M_{k}^{e_{k}}\right)\bmod L\bmod N (11) =(∑j=0|P|−1rj​uj)modLmodN\displaystyle=\left(\sum_{j=0}^{|P|-1}r_{j}u_{j}\right)\bmod L\bmod N To simplify later steps, I decompose rjr_{j} into its ℓ\ell bits rj,kr_{j,k}, and write VV in those terms: rj=∑k=0ℓ−1rj,k≪kr_{j}=\sum_{k=0}^{\ell-1}r_{j,k}\ll k (12) V=(∑j=0|P|−1∑k=0ℓ−1rj,k[uj≪k])modLmodNV=\left(\sum_{j=0}^{|P|-1}\sum_{k=0}^{\ell-1}r_{j,k}\Big[u_{j}\ll k\Big]\right)\bmod L\bmod N (13) The computation of VV is now a sum modulo LL then modulo NN, instead of a product modulo NN. Crucially, this fixes the issue where errors early in the process would be amplified by later operations. As a result, we’re in a position to start using approximations. Our goal now is to use approximations to extract the most significant bits of VV, without computing the entirety of VV. To quantify the error in these approximations, we’ll track its “modular deviation” ΔN​(a−b)\Delta_{N}(a-b). The modular deviation is the minimum number of increments or decrements needed to turn aa into bb modulo NN, divided by NN: ΔN​(a−b)=min⁡((a−b)modN,(b−a)modN)N\Delta_{N}(a-b)=\frac{\min((a-b)\bmod N,(b-a)\bmod N)}{N} (14) To approximate additions modulo NN, we’ll truncate all values down to ff bits. When asked to add an offset ss into an accumulator, we’ll instead add (smodN)≫t(s\bmod N)\gg t where t=len ​N−ft=\text{len }N-f. The accumulator will be truncated to ff bits, and will be operated modulo N≫tN\gg t instead of modulo NN. At any time, the approximate result can be extracted by left shifting the accumulator’s value by tt. Truncating introduces two sources of error. First, during additions, carries that would have propagated out of the removed part of the register into the kept part of the register will no longer happen. Second, the accumulating additions wrap around slightly too quickly due to (N≫t)≪t(N\gg t)\ll t being less than NN. For an individual truncated addition, both of these errors offset the approximate result by an amount O⁡(2t)O(2^{t}), which is a modular deviation of O⁡(2−f)O(2^{-f}). So the modular deviation is exponentially small in the number of kept bits, and will accumulate linearly with the number of operations, meaning a series of AA truncated additions has a modular deviation of at most O⁡(A⋅2−f)O(A\cdot 2^{-f}): t=(len ​N)−ft=(\text{len }N)-f (15) S\displaystyle S =(∑k=0A−1sk)modN\displaystyle=\left(\sum_{k=0}^{A-1}s_{k}\right)\bmod N (16) S~\displaystyle\widetilde{S} =(∑k=0A−1(skmodN)≫t)mod(N≫t)\displaystyle=\left(\sum_{k=0}^{A-1}(s_{k}\bmod N)\gg t\right)\bmod(N\gg t) ΔN​(S−(S~≪t))\displaystyle\Delta_{N}(S-(\tilde{S}\ll t)) ≤O⁡(A⋅2−f)\displaystyle\leq O(A\cdot 2^{-f}) This kind of approximated accumulation isn’t quite directly applicable to the current expression for VV. The issue is that VV isn’t being accumulated modulo NN, it’s being accumulated modulo LL and only at the end is the modulo NN performed. If we were to accumulate only modulo NN, instead of modulo LL then modulo NN, then each time the accumulator would have wrapped modulo LL we’d miss an offset of LmodNL\bmod N. However, LL is the product of the primes in our residue system, and we get to choose these primes, so we can pull a trick. We can optimize LL such that LmodNL\bmod N has negligible modular deviation. Numerically, it seems to be the case that picking random sets of small primes results in values of LmodNL\bmod N uniformly distributed over the range [1,N)[1,N). In cases I’ve tested, I’m consistently able to find an LL with deviation below 2−f2^{-f} with high probability by randomly sampling O⁡(2f)O(2^{f}) sets of small primes. I conjecture this is true in general (see Assumption 1). I’ll show later in Equation 23 that the required value of ff grows logarithmically with the problem size len ​N\text{len }N. Assuming the conjecture and the promise that ff will grow like O⁡(log⁡log⁡N)O(\log\log N), an LL with sufficiently small deviation can be found in polynomial time: ΔN​(L)<2−f(by brute force random search for a satisfying P)\Delta_{N}(L)<2^{-f}\;\;\;\;\text{(by brute force random search for a satisfying $P$)} (17) Assumption 1 (random search for small modular deviations is efficient): A set PP of small primes satisfying (∏P)≥Nm(\prod P)\geq N^{m} and ΔN​(∏P)<2−f\Delta_{N}(\prod P)<2^{-f} can be found in O⁡(2f⋅poly​(m⋅len ​N))O(2^{f}\cdot\text{poly}(m\cdot\text{len }N)) expected time by randomly varying the primes included in PP. As an example, here’s a set of 25000 primes I found in ten seconds using a 128 core machine. Each prime is 22 bits long, and the set achieves a modular deviation below 2−322^{-32} versus the RSA2048 challenge number [Wik24]: P1\displaystyle P_{1} =PrimesBetween​(3814620,222)\displaystyle=\text{PrimesBetween}(3814620,2^{22}) (18) P2\displaystyle P_{2} ={2097769,3484783,3814501,3814543,3814561,3814583,3814609}\displaystyle=\{2097769,3484783,3814501,3814543,3814561,3814583,3814609\} L\displaystyle L =∏(P1∪P2)\displaystyle=\prod(P_{1}\cup P_{2}) N\displaystyle N =RSA2048\displaystyle=\text{RSA}_{2048} ΔN​(L)\displaystyle\Delta_{N}(L) <2−32\displaystyle<2^{-32} Given the promise that the modular deviation of the wraparound error is at most 2−f2^{-f}, we can generalize Equation 16 from approximating arithmetic modulo NN to approximating arithmetic modulo LL then modulo NN: S\displaystyle S =(∑k=0A−1sk)modLmodN\displaystyle=\left(\sum_{k=0}^{A-1}s_{k}\right)\bmod L\bmod N (19) S~\displaystyle\widetilde{S} =(∑k=0A−1(skmodLmodN)≫t)mod(N≫t)\displaystyle=\left(\sum_{k=0}^{A-1}(s_{k}\bmod L\bmod N)\gg t\right)\bmod(N\gg t) ΔN​(S−(S~≪t))\displaystyle\Delta_{N}(S-(\tilde{S}\ll t)) ≤O⁡(A⋅2−f+A⋅ΔN​(L))\displaystyle\leq O(A\cdot 2^{-f}+A\cdot\Delta_{N}(L)) ≤O⁡(A⋅2−f)\displaystyle\leq O(A\cdot 2^{-f}) We can now apply Equation 19 to Equation 13, producing an expression for V~≈V≫t\widetilde{V}\approx V\gg t: V~\displaystyle\widetilde{V} =(∑j=0|P|−1∑k=0ℓ−1rj,k​([(uj≪k)modLmodN]≫t))mod(N≫t)\displaystyle=\left(\sum_{j=0}^{|P|-1}\sum_{k=0}^{\ell-1}r_{j,k}\left(\Big[(u_{j}\ll k)\bmod L\bmod N\Big]\gg t\right)\right)\bmod(N\gg t) (20) =(∑j=0|P|−1∑k=0ℓ−1rj,k​Cj,k)mod(N≫t)\displaystyle=\left(\sum_{j=0}^{|P|-1}\sum_{k=0}^{\ell-1}r_{j,k}C_{j,k}\right)\bmod(N\gg t) So, by using truncated residue arithmetic, the modular exponentiation has been approximated into an ff-bit dot product between the bits of the residues and a table of classically precomputable constants Cj,kC_{j,k}: Cj,k=([(uj≪k)modLmodN]≫t)mod(N≫t)C_{j,k}=\left(\Big[(u_{j}\ll k)\bmod L\bmod N\Big]\gg t\right)\bmod(N\gg t) (21) The computation of V~\widetilde{V} incurs O⁡(2−f)O(2^{-f}) modular deviation per truncated addition, and there are |P|⋅ℓ|P|\cdot\ell additions, so the total modular deviation of V~\widetilde{V} is at most ΔN​(V−(V~≪t))≤O⁡(|P|⋅ℓ⋅2−f)\Delta_{N}(V-(\widetilde{V}\ll t))\leq O(|P|\cdot\ell\cdot 2^{-f}) (22) It will be clear from the next subsection that, in order for period finding to work, it’s sufficient to achieve a constant total modular deviation (e.g. a modular deviation of 10%). Solving for ff in |P|⋅ℓ⋅2−f=O⁡(1)|P|\cdot\ell\cdot 2^{-f}=O(1) after expanding terms using Equation 3, Equation 4, Equation 6, and Equation 8 gives a bound on the size of ff: f=O⁡(log⁡log⁡N)f=O(\log\log N) (23) See Figure 3 for working Python code that computes V~≈gemodN\widetilde{V}\approx g^{e}\bmod N using this method. Figure 3: Example Python code that approximates V~≈gemodN\widetilde{V}\approx g^{e}\bmod N given a choice of PP and ff. Registers that would store quantum values during the quantum factoring are prefixed with “Q_”. All would-be-quantum registers are of size O⁡(log⁡log⁡N)O(\log\log N), except the input register Q_e. For simplicity, this code omits optimizations (like windowing) and crucial quantum details (like uncomputation). See Appendix A.1 for more detailed code. import math def approximate_modular_exponentiation( g: int, Q_e: int, N: int, P: list[int], f: int, ) -> int: """Computes␣an␣approximate␣modular␣exponentiation. ␣␣␣␣Args: ␣␣␣␣␣␣␣␣g:␣The␣base␣of␣the␣exponentiation. ␣␣␣␣␣␣␣␣Q_e:␣The␣exponent. ␣␣␣␣␣␣␣␣N:␣The␣modulus. ␣␣␣␣␣␣␣␣P:␣Small␣primes␣for␣the␣residue␣arithmetic. ␣␣␣␣␣␣␣␣f:␣Kept␣bits␣during␣approximate␣accumulation. ␣␣␣␣Returns: ␣␣␣␣␣␣␣␣An␣approximation␣of␣pow(g,␣Q_e,␣N). ␣␣␣␣␣␣␣␣The␣modular␣deviation␣of␣the␣approximation␣is ␣␣␣␣␣␣␣␣at␣most␣3␣*␣Q_e.bit_length()␣/␣2**f␣. ␣␣␣␣""" L = math.prod(P) ell = max(p.bit_length() for p in P) m = Q_e.bit_length() t = N.bit_length() - f assert L == math.lcm(*P) and L >= N ** m and (L % N) < (N >> f) Q_total = 0 for p in P: Q_residue = 1 for k in range(Q_e.bit_length()): precomputed = pow(g, 1 << k, N) % p if Q_e & (1 << k): Q_residue *= precomputed Q_residue %= p u = (L // p) * pow(L // p, -1, p) for k in range(ell): precomputed = (((u << k) % L % N) >> t) % (N >> t) if Q_residue & (1 << k): Q_total += precomputed Q_total %= N >> t return Q_total << t def test_approximate_modular_exponentiation(): import random import sympy N = int( "15226050279225333605356183781326374297180681149613" "80688657908494580122963258952897654000350692006139" ) g = random.randrange(2, N) Q_e = random.randrange(2**100) P = [ *sympy.primerange(239382, 2**18), 131101, 131111, 131113, 131129, 131143, 131149, 131947, 182341, 239333, 239347, ] f = 24 result = approximate_modular_exponentiation(g, Q_e, N, P, f) error = result - pow(g, Q_e, N) deviation = min(error % N, -error % N) / N ell = max(p.bit_length() for p in P) assert deviation <= 3 * len(P) * ell / 2**f 2.2 Approximate Period Finding Shor’s algorithm, as usually described, performs period finding against a periodic function f⁡(x)f(x). It relies crucially on the fact that f⁡(x)f(x) is exactly periodic, with f⁡(x+P)=f⁡(x)f(x+P)=f(x) for all values of xx. In this paper I’m computing an approximation f~≈f\widetilde{f}\approx f, so I can only guarantee approximate periodicity: ∀x,y∈ℤ:ΔN​(f~​(x+y​P)−f~​(x))≤ϵ.\forall x,y\in\mathbb{Z}:\Delta_{N}\left(\widetilde{f}(x+yP)-\widetilde{f}(x)\right)\leq\epsilon. (24) Exact period finding will fail if used on a function that only guarantees approximate periodicity, because the error in the approximation can make the function aperiodic. To fix this, we must avoid measuring the error. For example, instead of measuring the entire output of f~​(x)\widetilde{f}(x) you could measure a truncated output f~(x)//(100Nϵ)\widetilde{f}(x)//(100N\epsilon). This could work, but causes two problems. The first problem is that truncating the output means more inputs will be consistent with it. This will cause the system to collapse to a large superposition of periodic signals, instead of a simple periodic signal. I’ll address this by analyzing the behavior of period finding against superpositions of periodic signals, as was done in [MS22]. The second problem with truncating the output is that the parts of the output that aren’t measured would need to be uncomputed, by running the modular exponentiation process backwards. The modular exponentiation is already the most expensive part of the algorithm, so this would double the execution costs. I’ll address this by using superposition masking. Superposition masking uses sacrificial superpositions to control the information that can enter into a register or be revealed by measurements [Zal06, Gid19, JG20]. In the present context, we want to avoid learning the low bits of f~​(x)\widetilde{f}(x). We can achieve this by measuring f~​(x)+s\widetilde{f}(x)+s instead of f~​(x)\widetilde{f}(x), where ss is a uniform superposition over a contiguous range [0,⌈S​N⌉)[0,\lceil SN\rceil). I’ll address the value of SS later; for now note that larger values of SS will suppress error from the approximation but reduce the amount of information remaining for period finding. Because f~​(x)\widetilde{f}(x) happens to be computed with a series of additions into an output register, it’s possible to compute and measure f~​(x)+s\widetilde{f}(x)+s without needing to then uncompute f~​(x)\widetilde{f}(x) or ss. This is done by initializing the output register to ss (instead of 0), then adding the approximation f~​(x)\widetilde{f}(x) into the output register, and then measuring the entire output register. To understand this process in detail, let’s analyze the states that occur as the algorithm progresses. At the beginning of the algorithm, the output register is storing a uniform superposition (the mask). There is also the input register: an mm qubit register storing its own uniform superposition. Together these two registers form the initial state |ψ0⟩|\psi_{0}\rangle: |ψ0⟩\displaystyle|\psi_{0}\rangle =|0:2m:1⟩⊗|0:⌈SN⌉:1⟩\displaystyle=|0:2^{m}:1\rangle\otimes|0:\lceil SN\rceil:1\rangle (25) =12m​⌈S​N⌉​∑e=02m−1∑s=0⌈S​N⌉−1|e⟩⊗|s⟩\displaystyle=\frac{1}{\sqrt{2^{m}\lceil SN\rceil}}\sum_{e=0}^{2^{m}-1}\sum_{s=0}^{\lceil SN\rceil-1}|e\rangle\otimes|s\rangle The algorithm now adds the approximation f~​(e)\widetilde{f}(e) into the output register, modulo NN, forming the actual pre-measurement state |ψ1~⟩|\widetilde{\psi_{1}}\rangle: |ψ1~⟩\displaystyle|\widetilde{\psi_{1}}\rangle =1⌈S​N⌉​2m​∑e=02m−1∑s=0⌈S​N⌉−1|e⟩⊗|(s+f~​(e))modN⟩\displaystyle=\frac{1}{\sqrt{\lceil SN\rceil 2^{m}}}\sum_{e=0}^{2^{m}-1}\sum_{s=0}^{\lceil SN\rceil-1}|e\rangle\otimes|(s+\widetilde{f}(e))\bmod N\rangle (26) If we had used ff instead of f~\widetilde{f}, we would have produced the ideal pre-measurement state |ψ1⟩|\psi_{1}\rangle: |ψ1⟩\displaystyle|\psi_{1}\rangle =1⌈S​N⌉​2m​∑e=02m−1∑s=0⌈S​N⌉−1|e⟩⊗|(s+f⁡(e))modN⟩\displaystyle=\frac{1}{\sqrt{\lceil SN\rceil 2^{m}}}\sum_{e=0}^{2^{m}-1}\sum_{s=0}^{\lceil SN\rceil-1}|e\rangle\otimes|(s+f(e))\bmod N\rangle (27) Consider the infidelity 1−|⟨ψ1|ψ1~⟩|21-|\langle\psi_{1}|\widetilde{\psi_{1}}\rangle|^{2}. For each possible computational basis value ee of the input register, conditioning |ψ1~⟩|\widetilde{\psi_{1}}\rangle and |ψ1⟩|\psi_{1}\rangle on the input register storing ee produces a uniform superposition in the output register covering ⌈S​N⌉\lceil SN\rceil contiguous values. If the maximum modular deviation of f~\widetilde{f} is ϵ\epsilon, then the conditioned ranges in |ψ1~⟩|\widetilde{\psi_{1}}\rangle are offset from the conditioned ranges in |ψ1⟩|\psi_{1}\rangle by at most N​ϵN\epsilon. The width of the ranges is at least S/ϵS/\epsilon times wider than the offset from the deviation. Therefore the infidelity between the two conditioned states is a most ϵ/S\epsilon/S. But this is true for every condition, and so also bounds the total infidelity between the two states: 1−|⟨ψ1|ψ1~⟩|2≤ϵ/S1-|\langle\psi_{1}|\widetilde{\psi_{1}}\rangle|^{2}\leq\epsilon/S (28) We can treat this infidelity as contributing a chance of the algorithm failing, and thereby continue the analysis assuming we are working with the ideal state |ψ1⟩|\psi_{1}\rangle instead of the actual state |ψ1~⟩|\widetilde{\psi_{1}}\rangle. For example, if we set S=1000​ϵS=1000\epsilon then |ψ1~⟩|\widetilde{\psi_{1}}\rangle would have 99.9% overlap with |ψ1⟩|\psi_{1}\rangle and so any success rate analysis we did with |ψ1⟩|\psi_{1}\rangle would produce a result within 0.1% of the true success rate. The next step in the algorithm, which we’ll analyze with |ψ1⟩|\psi_{1}\rangle instead of |ψ1~⟩|\widetilde{\psi_{1}}\rangle, is to measure the output register. This produces some measurement result VV, and collapses the system into the post-measurement state |ψ2⟩|\psi_{2}\rangle. This state is a superposition of all the exponents whose masked output ranges overlapped with VV: |ψ2⟩\displaystyle|\psi_{2}\rangle ∝∑e=02m−1|e⟩⋅int​(0≤(V−f⁡(e))modN<⌈N​S⌉)\displaystyle\propto\sum_{e=0}^{2^{m}-1}|e\rangle\cdot\text{int}\left(0\leq(V-f(e))\bmod N<\lceil NS\rceil\right) (29) Let PP be the period of ff, and let RR be the set of exponents less than PP that are consistent with the measured value VV. We can rewrite |ψ2⟩|\psi_{2}\rangle in these terms: R={e∈ℕ|e0k>0. The sum over values of RR can be pushed rightward, so that it just counts how many values of RR are compatible with a choice of (r1,r2)(r_{1},r_{2}). When r1=r2r_{1}=r_{2}, there are (P−1w−1)\binom{P-1}{w-1} values of RR that satisfy r1∈R∧r2∈Rr_{1}\in R\land r_{2}\in R. When r1≠r2r_{1}\neq r_{2}, there are (P−2w−2)\binom{P-2}{w-2} satisfying values. This replaces the sum over values of RR with a choice of multiplier determined by r1=?r2r_{1}\stackrel{{\scriptstyle?}}{{=}}r_{2}: E⁡(|Bk>0|2)\displaystyle E(|B_{k>0}|^{2}) =1(Pw)​w​P∑r1=0P−1∑r2=0P−1ei​k​r1​2​π/Pe−ikr22π/P(∑R∈𝒫⁡([0,P))|R|=wint(r1∈R∧r2∈R))\displaystyle=\frac{1}{\binom{P}{w}wP}\sum_{r_{1}=0}^{P-1}\sum_{r_{2}=0}^{P-1}e^{ikr_{1}2\pi/P}e^{-ikr_{2}2\pi/P}\left(\sum_{\begin{subarray}{c}R\in\mathcal{P}([0,P))\\ |R|=w\end{subarray}}\text{int}(r_{1}\in R\land r_{2}\in R)\right) (39) =1(Pw)​w​P∑r1=0P−1∑r2=0P−1ei​r1​k​2​π/Pe−ir2k2π/P{r1=r2→(P−1w−1)r1≠r2→(P−2w−2)\displaystyle=\frac{1}{\binom{P}{w}wP}\sum_{r_{1}=0}^{P-1}\sum_{r_{2}=0}^{P-1}e^{ir_{1}k2\pi/P}e^{-ir_{2}k2\pi/P}\begin{cases}r_{1}=r_{2}&\rightarrow\binom{P-1}{w-1}\\ r_{1}\neq r_{2}&\rightarrow\binom{P-2}{w-2}\end{cases} Note that, for non-zero values of kk, the nested sum in the above equation would evaluate to 0 if the right side multiplier was replaced by a constant: 00|2)\displaystyle E(|B_{k>0}|^{2}) =(Pw)−11w​P∑r=0P−1ei​r​k​2​π/Pe−irk2π/P((P−1w−1)−(P−2w−2))\displaystyle=\binom{P}{w}^{-1}\frac{1}{wP}\sum_{r=0}^{P-1}e^{irk2\pi/P}e^{-irk2\pi/P}\left(\binom{P-1}{w-1}-\binom{P-2}{w-2}\right) (41) =(Pw)−1​1w​P​P​((P−1w−1)−(P−2w−2))\displaystyle=\binom{P}{w}^{-1}\frac{1}{wP}P\left(\binom{P-1}{w-1}-\binom{P-2}{w-2}\right) =P−wP⁡(P−1)\displaystyle=\frac{P-w}{P(P-1)} In other words, when ww randomly chosen simple periodic signals are superposed, the system behaves as if β0\beta_{0} takes a w/Pw/P cut and then β1,…,βP−1\beta_{1},\dots,\beta_{P-1} equally split whats left: E⁡(|βk|2)={k=0→wP0