Fourier Methods

Kickback from controlled powers will leave a register of dimension \(N\) in a state like

\[\frac{1}{\sqrt N}\sum_{k=0}^{N-1}e^{2\pi ijk/N}\ket k.\]

Each basis state carries a phase that advances by the same step \(2\pi j/N\), like the hand of a clock. Suppose \(j\) is an unknown integer between \(0\) and \(N-1\). Can a measurement tell us \(j\)? In the computational basis every outcome has probability \(1/N\), whatever \(j\) is. The information sits entirely in the phases.

From the readout of Deutsch’s circuit we know what would make \(j\) readable: the states for different \(j\) must be orthogonal, and then a change of basis maps each one to \(\ket j\). Are they orthogonal? Write \(\omega=e^{2\pi i/N}\). The inner product of the states for \(j\) and \(j'\) is \(\frac1N\sum_k\omega^{(j'-j)k}\), so we need to know when such sums vanish.

The powers \(1,\omega,\omega^2,\ldots,\omega^{N-1}\) are \(N\) equally spaced points on the unit circle, the \(N\)th roots of unity, and \(\omega^N=1\). For an integer \(k\) not divisible by \(N\), they cancel:

\[\sum_{x=0}^{N-1}\omega^{kx}=0.\]

Here is a proof that needs no geometry. Call the sum \(S\). Multiplying by \(\omega^k\) shifts every term to the next one, and the last term \(\omega^{k(N-1)}\) becomes \(\omega^{kN}=1\), the first. So \(\omega^kS=S\). Since \(\omega^k\neq1\), \(S=0\). Geometrically, \(N\) arrows equally spaced around a circle (or wound around it several times) sum to nothing. When \(N\) divides \(k\), every term is \(1\) and the sum is \(N\).

For \(j\neq j'\) the difference \(j'-j\) lies strictly between \(-N\) and \(N\) and is not zero, so \(N\) does not divide it, and the inner product vanishes. Each state has length one. So the \(N\) clock states are orthonormal, and they form a basis. The unitary that maps each basis state to its clock state is the quantum Fourier transform:

\[F_N\ket x=\frac{1}{\sqrt N}\sum_{y=0}^{N-1}\omega^{xy}\ket y.\]

It sends an orthonormal basis to an orthonormal basis, so it is unitary. Its inverse \(F_N^\dagger\ket y=\frac{1}{\sqrt N}\sum_x\omega^{-xy}\ket x\) maps the clock state for \(j\) back to \(\ket j\), which is the readout we wanted.

For \(N=2\), \(\omega=-1\) and \(F_2=H\). For \(N=2^n\), \(F_N\) is not \(H^{\otimes n}\): the Hadamard transform uses signs \((-1)^{x\cdot z}\) with the bitwise product, while \(F_N\) uses \(\omega^{xy}\) with ordinary multiplication of integers. The two are Fourier transforms over the two groups we have met: bit strings under XOR, and integers under addition modulo \(N\).

What does \(F_N\) do to a state whose amplitudes repeat? Take \(N=8\) and the state supported on every second position,

\[\ket{\psi_0}=\tfrac12(\ket0+\ket2+\ket4+\ket6).\]

What to observe

Two outcomes, \(0\) and \(4\), each with probability \(\frac12\). A pattern of period 2 in a register of size 8 produces peaks at multiples of \(8/2=4\).

The amplitude of \(\ket y\) is \(\frac{1}{2\sqrt8}\sum_{k=0}^{3}\omega^{2ky}\). With \(\omega^2=e^{2\pi i/4}\), this is a sum of fourth roots of unity raised to the power \(y\): it is \(4\) when \(4\mid y\) and \(0\) otherwise. So

\[F_8\ket{\psi_0}=\tfrac{1}{\sqrt2}(\ket0+\ket4).\]

Now shift the pattern by one: \(\ket{\psi_1}=\tfrac12(\ket1+\ket3+\ket5+\ket7)\). Every term gains an extra factor \(\omega^{y}\), so \(F_8\ket{\psi_1}=\tfrac{1}{\sqrt2}(\ket0-\ket4)\). The probabilities are unchanged. The Fourier transform converts a shift into a phase, exactly as the Hadamard transform did in Simon’s algorithm, and a phase is invisible to measurement.

In general, if \(r\) divides \(N\) and a state is uniform on \(x_0,x_0+r,x_0+2r,\ldots\), its Fourier transform is uniform over the \(r\) multiples of \(N/r\), the frequencies \(0,N/r,2N/r,\ldots\), with phases depending on \(x_0\). For \(N=8\) and \(r=2\) these are \(0\) and \(4\), as we found. Measuring gives a random frequency \(kN/r\), which carries information about \(r\). When \(r\) does not divide \(N\), the peaks are no longer exact: they spread over nearby integers, and recovering \(r\) then needs more work.

A unitary on \(n=\log_2N\) qubits is a \(2^n\times2^n\) matrix; a general one needs exponentially many gates. The Fourier transform is special. Write \(x=\sum_jx_j2^j\) and \(y=\sum_ky_k2^k\). The phase \(\omega^{xy}=e^{2\pi ixy/2^n}\) depends on \(xy/2^n\) only modulo one, and

\[\frac{xy}{2^n}=\sum_k y_k\,\frac{x\,2^k}{2^n}\equiv\sum_ky_k\,(0.x_{n-k-1}\cdots x_1x_0)\pmod 1,\]

because multiplying \(x\) by \(2^k\) shifts its binary point, and the integer part contributes nothing. Each output bit \(y_k\) receives a phase determined by a binary fraction of the input bits. The sum over \(y\) therefore factors into a product over output bits:

\[F_N\ket x=\bigotimes_{k=n-1}^{0}\frac{\ket0+e^{2\pi i\,(0.x_{n-k-1}\cdots x_0)}\ket1}{\sqrt2},\]

with the highest output bit written first. This product representation says the output is unentangled for every basis input, and each factor can be built separately. The highest output bit gets \(e^{2\pi i\,0.x_0}=(-1)^{x_0}\), a Hadamard on \(x_0\). The lowest gets the full fraction \(0.x_{n-1}\cdots x_0\).

The circuit builds the factors in place. Apply \(H\) to the top input bit \(x_{n-1}\), giving phase \(e^{2\pi i\,0.x_{n-1}}\); then controlled phase gates \(P(\pi/2)\) from \(x_{n-2}\), \(P(\pi/4)\) from \(x_{n-3}\), and so on, add \(0.0x_{n-2}\), \(0.00x_{n-3}\), …, completing the factor for \(y_0\). Repeat on the next wire with one fewer rotation. The factors end in reversed order, so swaps restore the bit order.

The circuit uses \(n\) Hadamards, \(n(n-1)/2\) controlled phase gates, and \(\lfloor n/2\rfloor\) swaps: \(O(n^2)\) gates for a transform on a space of dimension \(2^n\). The inverse \(F_N^\dagger\) is the same circuit run backwards with negated angles. The construction, and its approximate version below, are due to Coppersmith (1994).

The smallest rotations in the QFT have angles \(2\pi/2^n\). For large \(n\) they are physically tiny and, in fault-tolerant gate sets, expensive to synthesize. Do we need them?

A controlled phase \(\mathrm{C}P(\theta)\) differs from the identity only on \(\ket{11}\), by \(e^{i\theta}-1\), whose magnitude is at most \(|\theta|\). So omitting it changes the operation by at most \(|\theta|\) in operator norm. By the error-composition bound for products of unitaries, omitting several gates changes the whole circuit by at most the sum of their angles. If we drop every rotation with angle below \(2\pi/2^b\), each wire loses rotations \(2\pi/2^{b+1},2\pi/2^{b+2},\ldots\), totalling less than \(2\pi/2^b\). Over \(n\) wires, the approximate QFT has error less than \(2\pi n/2^b\).

To reach error \(\epsilon\), choose \(b\approx\log_2(2\pi n/\epsilon)\). Each wire then keeps at most \(b\) rotations, so the gate count falls from \(O(n^2)\) to \(O(n\log(n/\epsilon))\). In phase estimation a modest constant error in the QFT only slightly lowers the success probability, so \(b\) of order \(\log n\) suffices. Each remaining rotation must still be synthesized, and its synthesis error is added to the total in the same way. When a resource estimate quotes the QFT as “\(O(n^2)\) gates”, ask whether it is exact or approximate, and in which gate set.

How do we write a clock pattern whose step is an unknown eigenphase? Let \(U\) be a unitary with eigenstate \(\ket u\) and eigenvalue \(e^{2\pi i\phi}\), where \(0\le\phi<1\). Phase estimation is the task of finding \(\phi\). Assume we can prepare \(\ket u\) and apply controlled-\(U^{2^j}\) for \(j=0,\ldots,m-1\).

Take \(m\) control qubits, each in \(\ket+\). Control qubit \(j\) (counting the lowest as \(j=0\)) controls \(U^{2^j}\), which kicks back \(e^{2\pi i\,2^j\phi}\) onto its \(\ket1\) branch, because \(U^k\ket u=e^{2\pi ik\phi}\ket u\). The target stays \(\ket u\). The control register becomes

\[\bigotimes_{j=m-1}^{0}\frac{\ket0+e^{2\pi i\,2^j\phi}\ket1}{\sqrt2}=\frac{1}{\sqrt{2^m}}\sum_{k=0}^{2^m-1}e^{2\pi ik\phi}\ket k,\]

where the equality expands the product: a basis state \(\ket k\) with bits \(k_j\) collects the phase \(\prod_je^{2\pi i\,k_j2^j\phi}=e^{2\pi ik\phi}\). This state is the steady phase ramp the Fourier transform was built to read.

The controlled powers are the expensive part. Controlled-\(U^{2^j}\) built by repeating controlled-\(U\) costs \(2^j\) uses, so all powers together cost \(1+2+\cdots+2^{m-1}=2^m-1\) uses of controlled-\(U\). Each extra bit of precision doubles the cost. For some unitaries, \(U^{2^j}\) has a shortcut. Modular multiplication by \(a\), raised to the power \(2^j\), is multiplication by \(a^{2^j}\bmod L\), which costs the same as multiplication by \(a\) once \(a^{2^j}\) is precomputed classically. That shortcut is what makes Shor’s algorithm efficient. For many other unitaries, including those that simulate physical systems, there is generally no such shortcut, and the doubling is real.

When does the inverse transform read the phase perfectly? Compare the control register with the Fourier transform of a basis state. With \(M=2^m\),

\[F_M\ket j=\frac{1}{\sqrt M}\sum_ke^{2\pi ijk/M}\ket k.\]

If \(\phi=j/M\) exactly, for an integer \(j\), the control register is \(F_M\ket j\). Applying the inverse quantum Fourier transform \(F_M^\dagger\) returns \(\ket j\), and measuring gives \(j\) with certainty. In binary, \(\phi=0.j_{m-1}\cdots j_0\): phase estimation reads the first \(m\) binary digits of the phase.

Here is an \(m=3\) example. Let \(\phi=3/8=0.011_2\). The controls receive phases \(e^{2\pi i\cdot3/8}\) (lowest, from \(U\)), \(e^{2\pi i\cdot6/8}\) (from \(U^2\)), and \(e^{2\pi i\cdot12/8}=e^{2\pi i\cdot4/8}=-1\) (highest, from \(U^4\)). The register is \(F_8\ket3\), and \(F_8^\dagger\) maps it to \(\ket{011}\).

The product form explains why the digits come out separately. Control qubit \(j\) carries phase \(e^{2\pi i\,2^j\phi}\), and multiplying \(\phi\) by \(2^j\) shifts its binary point, so modulo one this phase depends only on the digits of \(\phi\) from position \(j+1\) onward. The highest control, with \(U^{2^{m-1}}\), carries \(e^{2\pi i\,0.j_0}=\pm1\): it holds the last digit as a sign, readable by a single Hadamard. The others need the controlled phase corrections inside \(F_M^\dagger\), which subtract the contributions of digits already known. If we measure the highest control first, those corrections can be applied classically, one digit at a time, so a single control qubit suffices; we build that version once measurement becomes a technique of its own.

Most phases are not exact multiples of \(1/M\). Then the inverse transform cannot concentrate all amplitude on one outcome. How much lands near the right answer?

The amplitude of outcome \(j\) is the overlap of the register with \(F_M\ket j\):

\[a_j=\frac1M\sum_{k=0}^{M-1}e^{2\pi ik(\phi-j/M)}.\]

Write \(\Delta=\phi-j/M\). This is a geometric series with ratio \(e^{2\pi i\Delta}\), so

\[|a_j|^2=\frac{1}{M^2}\left|\frac{1-e^{2\pi iM\Delta}}{1-e^{2\pi i\Delta}}\right|^2=\frac{\sin^2(\pi M\Delta)}{M^2\sin^2(\pi\Delta)}.\]

For the nearest grid point, \(|\Delta|\le\frac{1}{2M}\). Using \(|\sin(\pi M\Delta)|\ge2M|\Delta|\) (since \(\sin u\ge2u/\pi\) on \([0,\pi/2]\)) and \(|\sin\pi\Delta|\le\pi|\Delta|\), we get \(|a_j|^2\ge4/\pi^2\approx0.405\). The distribution is peaked around \(M\phi\), and falls off roughly as \(1/(M\Delta)^2\).

For \(M=8\) and \(\phi=0.3\), so that \(M\phi=2.4\), the outcomes \(2\) and \(3\) have probabilities about \(0.58\) and \(0.26\); together \(0.84\). The phase resolution is \(1/M\): outcomes within one grid spacing are likely, and outcomes \(e\) spacings away have probability falling like \(1/e^2\). How do we buy more confidence? A quoted bound (Cleve, Ekert, Macchiavello, and Mosca, 1998; Nielsen and Chuang, Section 5.2.1) says that to obtain \(\phi\) to within \(2^{-t}\) with probability at least \(1-\delta\), it suffices to use \(m=t+\lceil\log_2(2+1/(2\delta))\rceil\) control qubits. Repetition also works. Suppose each run lands within the target interval with probability above one half, and we take the median of \(R\) runs. The median falls outside the interval only if at least half the runs do, and by the majority-vote bound that probability falls exponentially in \(R\).

In terms of cost: precision \(\epsilon\) needs \(M\approx1/\epsilon\), hence about \(1/\epsilon\) uses of controlled-\(U\). Compare the Hadamard test, whose \(1/\epsilon^2\) runs each used \(U\) once. This quadratic improvement, \(1/\epsilon\) instead of \(1/\epsilon^2\), comes from coherence across the controlled powers, and amplitude estimation will achieve the same improvement for probabilities.

What if the target is not an eigenstate? We would like to write it as a combination of eigenstates and let linearity handle each one. Two facts are needed. First, eigenstates of a unitary with different eigenvalues are orthogonal. Suppose \(U\ket u=\lambda\ket u\) and \(U\ket v=\mu\ket v\) with \(\lambda\neq\mu\). Then

\[\braket uv=\braket{Uu}{Uv}=\lambda^*\mu\braket uv.\]

Since \(|\lambda|=1\), \(\lambda^*=1/\lambda\), so \(\lambda^*\mu=\mu/\lambda\neq1\), and therefore \(\braket uv=0\). Second, the eigenstates of a unitary can be chosen to form an orthonormal basis. This is the spectral theorem of linear algebra, which we assume. Together they let us write any state as \(\ket\psi=\sum_sc_s\ket{u_s}\), with \(U\ket{u_s}=e^{2\pi i\phi_s}\ket{u_s}\) and \(\sum_s|c_s|^2=1\). This eigenbasis expansion is the key: phase estimation is linear, so it acts on each term separately:

\[\ket{0^m}\ket\psi\longmapsto\sum_sc_s\,\ket{\tilde\phi_s}\ket{u_s},\]

where \(\ket{\tilde\phi_s}\) is the control register that phase \(\phi_s\) would produce. Measuring the controls returns an estimate of one \(\phi_s\), chosen with probability \(|c_s|^2\), and leaves the target close to \(\ket{u_s}\) (exactly, if the phases are distinct grid points).

Two consequences shape what follows. First, phase estimation samples eigenvalues; it does not list them. To learn the eigenvalue of one particular eigenstate \(\ket{u_0}\), we need an initial state with substantial overlap \(|c_0|^2\) with it, and the expected number of repetitions is about \(1/|c_0|^2\). Preparing such a state can be the hard part. Second, phase estimation followed by measurement prepares an eigenstate. That makes it a filter. Polynomial transformations, developed later, refine this filter so that it acts on all eigenvalues coherently without measuring.

Order finding, next, exploits the first consequence in an unusual way: there, the easy-to-prepare state \(\ket1\) is an equal superposition of eigenstates, and the random eigenphase it yields is exactly what we need.

Which unitary carries the question? Let \(L\) have \(n\) bits and \(\gcd(a,L)=1\). On an \(n\)-qubit register define

\[U_a\ket y=\ket{ay\bmod L}\quad(0\le y<L),\qquad U_a\ket y=\ket y\quad(L\le y<2^n).\]

Because \(a\) has an inverse modulo \(L\), multiplication by \(a\) permutes the residues \(0,\ldots,L-1\), so \(U_a\) is a permutation and hence unitary. It can be built from reversible modular addition, using \(O(n^2)\) gates with schoolbook arithmetic.

The return time is visible in the orbit of \(1\): \(1\to a\to a^2\to\cdots\to a^{r-1}\to1\). Take \(L=15\), \(a=7\):

\[1\to7\to49\bmod15=4\to28\bmod15=13\to91\bmod15=1,\]

The smallest \(r>0\) with \(a^r\equiv1\pmod L\) is the order of \(a\) modulo \(L\), so here \(r=4\). Classically, finding \(r\) by stepping along the orbit takes up to \(r\) steps, and \(r\) can be nearly as large as \(L\), exponential in \(n\).

For phase estimation we need controlled-\(U_a^{2^j}\). Since \(U_a^{2^j}=U_{a^{2^j}\bmod L}\), compute \(a^{2^j}\bmod L\) classically by repeated squaring, \(m\) squarings in all, and build a controlled multiplier by that constant. All the controlled powers together form modular exponentiation, \(\ket k\ket y\mapsto\ket k\ket{a^ky\bmod L}\), which costs \(m\) controlled multiplications, \(O(mn^2)\) gates. This avoids the exponential cost of building each power by repeating controlled-\(U_a\).

What are the eigenstates of \(U_a\)? Since \(U_a\) cycles \(1\to a\to\cdots\to a^{r-1}\to1\), it acts on the orbit like a shift on \(r\) positions. We saw with periodic states that the Fourier transform turns a shift into a phase, so Fourier-type states on the orbit should be eigenstates. For \(s=0,1,\ldots,r-1\), define

\[\ket{u_s}=\frac{1}{\sqrt r}\sum_{k=0}^{r-1}e^{-2\pi isk/r}\ket{a^k\bmod L}.\]

Applying \(U_a\) shifts \(k\) to \(k+1\), and the same substitution we used for the phase-gradient state gives \(U_a\ket{u_s}=e^{2\pi is/r}\ket{u_s}\). The eigenphases are \(s/r\): the order is in the denominator.

We cannot prepare \(\ket{u_s}\) without knowing \(r\). But we do not need to. Sum over \(s\):

\[\frac{1}{\sqrt r}\sum_{s=0}^{r-1}\ket{u_s}=\frac1r\sum_k\Big(\sum_se^{-2\pi isk/r}\Big)\ket{a^k}=\ket1,\]

because the inner sum is \(r\) for \(k=0\) and \(0\) otherwise. The state \(\ket1\), trivial to prepare, is an equal superposition of all \(r\) eigenstates. This is the heart of order finding: since phase estimation acts on a superposition of eigenstates term by term, phase estimation on \(\ket1\) returns an estimate of \(s/r\) for a uniformly random \(s\in\{0,\ldots,r-1\}\).

For \(L=15\), \(a=7\), \(r=4\), and \(m=3\) control qubits, the phases \(s/4\in\{0,\frac14,\frac12,\frac34\}\) are exact grid points \(j/8\) with \(j\in\{0,2,4,6\}\). Each outcome occurs with probability \(\frac14\). The outcome \(6\) gives \(6/8=3/4\), so \(r=4\). The outcome \(4\) gives \(1/2\), which suggests \(r=2\); checking \(7^2\bmod15=4\neq1\) rejects it, and another run is needed. The outcome \(0\) gives nothing.

When \(r\) does not divide \(2^m\), phase estimation returns an integer \(j\) with \(j/2^m\approx s/r\). How do we recover \(r\) from an approximation? Two facts help. The fraction \(s/r\) has a small denominator, \(r<L\), and the approximation can be made very good. We need a way to find the fractions with small denominators that lie close to a given number.

Euclid’s algorithm supplies one. Take a worked case. Let \(L=21\) and \(a=2\). The powers of \(2\) modulo \(21\) are \(2,4,8,16,11,1\), so \(r=6\). Take \(m=9\), so \(2^m=512\). Suppose \(s=1\); then \(512/6\approx85.3\), and the most likely outcome is \(j=85\). Run Euclid on \(85/512\):

\[512=6\cdot85+2,\qquad85=42\cdot2+1,\qquad2=2\cdot1.\]

The quotients say that

\[\frac{85}{512}=\cfrac{1}{6+\cfrac{1}{42+\cfrac12}},\]

written \([0;6,42,2]\). An expression \(a_0+1/(a_1+1/(a_2+\cdots))\) of this kind is a continued fraction. Stopping it early gives the convergents, here \(0\), \(\tfrac16\), \(\tfrac{42}{253}\), \(\tfrac{85}{512}\). Each convergent is a fraction with a small denominator that approximates the number well. The last convergent with denominator below \(L=21\) is \(1/6\). Candidate \(r=6\); check \(2^6=64\equiv1\pmod{21}\). Confirmed.

When is the fraction we want guaranteed to be a convergent? Number theory supplies the answer, which we quote (Legendre’s theorem; see Hardy and Wright, An Introduction to the Theory of Numbers, Theorem 184): if

\[\left|\frac{j}{2^m}-\frac sr\right|<\frac{1}{2r^2},\]

then \(s/r\), in lowest terms, is one of the convergents of \(j/2^m\). How many control qubits make this hold? When phase estimation returns the nearest grid point, \(|j/2^m-s/r|\le1/2^{m+1}\). Since \(r<L\), it suffices that \(1/2^{m+1}\le1/(2L^2)\), that is, \(2^m\ge L^2\), or about \(m=2n\). In the example, \(512\ge441\), and indeed \(|85/512-1/6|=1/1536<1/72\).

The overall order-finding procedure is then: run phase estimation on \(\ket1\) with \(m\approx2n+1\); compute the continued-fraction convergents of \(j/2^m\); for each denominator \(q<L\), test whether \(a^q\equiv1\); combine denominators from a few runs by least common multiples if needed. Why might a run fail? The nearest grid point appears with probability at least \(4/\pi^2\). And if \(s\) and \(r\) share a factor, the convergent is \(s/r\) in lowest terms, whose denominator is only a divisor of \(r\). The chance that a random \(s\) shares no factor with \(r\) is at least of order \(1/\log\log r\), a quoted number-theory estimate (Shor, 1997), and in practice usually much larger; a constant number of repetitions, plus classical checks, gives high confidence. The quantum part does one thing: it produces a sample of \(s/r\) to high precision. Everything else is classical number theory.

How can an order reveal a factor? Suppose \(L\) is odd, not a prime power, and we want a nontrivial factor. Pick \(a\) at random in \(\{2,\ldots,L-1\}\). If \(\gcd(a,L)>1\), we have a factor already. Otherwise find the order \(r\) of \(a\). If \(r\) is even, then

\[a^r-1=(a^{r/2}-1)(a^{r/2}+1)\equiv0\pmod L.\]

So \(L\) divides the product. It does not divide \(a^{r/2}-1\), since \(r\) is the smallest exponent with \(a^r\equiv1\). If also \(a^{r/2}\not\equiv-1\pmod L\), then \(L\) does not divide \(a^{r/2}+1\) either, and the prime factors of \(L\) must be split between the two factors. Then \(\gcd(a^{r/2}-1,L)\) is a nontrivial factor, computed by Euclid’s algorithm.

For \(L=15\), \(a=7\), \(r=4\): \(7^2=49\equiv4\), and \(\gcd(3,15)=3\), \(\gcd(5,15)=5\). For \(L=21\), \(a=2\), \(r=6\): \(2^3=8\), \(\gcd(7,21)=7\), \(\gcd(9,21)=3\).

This factoring reduction can fail: \(r\) may be odd, or \(a^{r/2}\equiv-1\). For \(L=15\), \(a=14\) has \(r=2\) and \(14\equiv-1\), so the gcds are trivial. How often does it fail? For \(L=15\) the coprime choices are \(a=2,4,7,8,11,13,14\), and only \(a=14\) fails. A quoted number-theory result (Nielsen and Chuang, Appendix 4) says that for odd \(L\) with at least two distinct prime factors, a random coprime \(a\) succeeds with probability at least \(1/2\). So a few random choices suffice. This is Shor’s algorithm (1994/1997). Prime powers \(L=p^k\) are detected and factored classically.

Where does the cost of Shor’s algorithm go? Take an \(n\)-bit \(L\) and follow the circuit from left to right.

Phase estimation needs about \(2n\) control qubits, and each controls one modular multiplication by a precomputed constant. With schoolbook arithmetic a multiplication costs \(O(n^2)\) gates, so modular exponentiation costs \(O(n^3)\) gates and dominates. Faster multiplication algorithms lower the exponent asymptotically. The arithmetic acts on an \(n\)-qubit target with \(O(n)\) workspace.

The inverse QFT on \(2n\) qubits then costs \(O(n^2)\) gates exactly, or \(O(n\log n)\) approximately. If its corrections are applied after measurement by classical feedback, one recycled control qubit replaces all \(2n\). What remains is classical: continued fractions, gcds and modular checks, all polynomial in \(n\), repeated over a few runs.

Asymptotically this is polynomial, \(O(n^3)\), against the best known classical algorithm, the number field sieve, which runs in time roughly \(\exp(c\,n^{1/3}(\log n)^{2/3})\). The speedup is superpolynomial, but it is a comparison with the best known classical method: no proof rules out a fast classical factoring algorithm.

What does this mean for a real key? A concrete resource estimate must fix the error-correction scheme and hardware assumptions. Factoring a 2048-bit RSA modulus requires billions of logical operations, each protected by error correction. A careful 2019 estimate found about 20 million noisy physical qubits running for about eight hours (Gidney and Ekerå); a 2025 estimate by Gidney reduced this to under a million noisy qubits running for under a week, using improved arithmetic and error-correction techniques (Gidney 2025). Such estimates depend on assumed physical error rates and cycle times, and they change as methods improve.

What do Simon’s problem and order finding have in common? In each, a function \(f\) on a group is constant on the cosets of a hidden subgroup and distinct across cosets. In Simon’s problem the group is \(\{0,1\}^n\) under XOR and the subgroup is \(\{0,s\}\). In order finding the group is the integers under addition and the subgroup is the multiples of \(r\), since \(a^x\) is constant on each coset \(x_0+r\mathbb Z\) and differs between cosets. The task in the hidden subgroup problem is to find the subgroup.

The algorithm is the same in both: evaluate \(f\) on a uniform superposition, measure the output to collapse onto a random coset, apply the Fourier transform of the group, and measure. The Fourier transform turns the random coset shift into a phase and leaves a sample that constrains the subgroup. Classical processing of a few samples recovers it.

Does it work for every group? For every finite abelian group it works efficiently, as a quoted result. Discrete logarithms, the basis of elliptic-curve cryptography, are an abelian hidden subgroup problem over \(\mathbb Z_N\times\mathbb Z_N\), solved in Shor’s same paper. For some nonabelian groups, where the order of combination matters, the same Fourier sampling is known to give too little information per sample. A polynomial number of queries still suffices in principle, but no efficient way to process them is known. The table lists two such cases, the dihedral and symmetric groups; for the dihedral group the best known algorithm (Kuperberg) is subexponential, not polynomial.

Fourier methods resolve periodicity, but they return useful outcomes only with some probability, and estimating a quantity to precision \(\epsilon\) still costs on the order of \(1/\epsilon\) controlled operations. When a useful outcome has small probability, can we make it large without knowing the outcome in advance?