The Order-Finding Subroutine

Fix $N$ and an integer $a$ with $\gcd(a,N)=1$. The multiplicative order of $a$ modulo $N$, written $\mathrm{ord}_N(a)$, is the smallest positive integer $r$ such that $a^r \equiv 1 \pmod N$. Such an $r$ always exists (by Euler's theorem, $r \le \varphi(N)$), and the sequence $a^0, a^1, a^2, \dots \pmod N$ is periodic with period exactly $r$. Order finding is precisely the task of computing this period given $a$ and $N$. The key structural fact is this: consider the unitary operator $U_a$ acting on the $n$-qubit register of integers mod $N$ (padded up to dimension $2^n \ge N$) defined by

and acting as the identity (or arbitrarily, as long as it stays unitary) on basis states $\vert x\rangle$ with $x \ge N$. Because $a$ is invertible mod $N$, this map is a permutation of $\{0,\dots,N-1\}$, hence $U_a$ is a genuine unitary (a permutation matrix) on that subspace. Since $a^r \equiv 1$, we have $U_a^r = I$ on that subspace, so every eigenvalue of $U_a$ restricted there is an $r$-th root of unity: $e^{2\pi i k/r}$ for $k = 0, 1, \dots, r-1$. The problem of finding $r$ becomes the problem of reading off one of these eigenphases and recognizing the denominator $r$ hiding inside it. Phase estimation is built exactly for this: given a unitary $U$ and an eigenstate $\vert\psi\rangle$ with $U\vert\psi\rangle = e^{2\pi i \varphi}\vert\psi\rangle$, it produces an estimate of $\varphi$. The catch is that we do not have direct access to a single eigenstate of $U_a$ — we can only prepare basis states, which are superpositions of eigenstates. It turns out this is not a problem: preparing $\vert 1\rangle$ (or any convenient computational basis state coprime-compatible with the cycle structure) as the phase-estimation input state is an equal superposition over the eigenstates of $U_a$ with a specific, uniform distribution of eigenvalues $e^{2\pi i k/r}$, $k=0,\dots,r-1$. Phase estimation on this input therefore collapses, upon measurement, to a random $k$ and reports (an estimate of) $k/r$. The remaining work is to recover $r$ from noisy knowledge of $k/r$, and that is exactly what the continued fraction algorithm does: given a good rational approximation to $k/r$, it reconstructs $k/r$ in lowest terms, from which $r$ is read off as the denominator (after ruling out common-factor degeneracies).

The full algorithm has a quantum stage (phase estimation on modular exponentiation) and a classical stage (continued fractions plus a repeat-and-verify loop).

Setup. Let $N$ be the modulus, $n = \lceil \log_2 N \rceil$ so an $n$-qubit register can hold any residue mod $N$. Choose the phase-estimation precision register to have $t = 2n + 1$ qubits (enough bits to make continued-fraction recovery unique with high probability; the "+1" is a small safety margin used in the standard analysis).

Step 1 — Prepare registers. Initialize a $t$-qubit control register in $\vert 0\rangle^{\otimes t}$ and an $n$-qubit target register in $\vert 1\rangle$ (a convenient basis state, since $1$ generates a cyclic subgroup interaction with every $a$).

Step 2 — Superposition on the control register. Apply Hadamard to every control qubit, producing

Step 3 — Controlled modular exponentiation. Apply the unitary

implemented as a cascade of controlled-$U_a^{2^i}$ operations (the standard phase-estimation trick), each conditioned on one bit of $j$: qubit $i$ of the control register controls an application of $U_a^{2^i}$, i.e. "multiply the target register by $a^{2^i} \bmod N$ if that control bit is 1." This is exactly repeated modular squaring — the same fast-exponentiation trick used classically — realized reversibly as a quantum circuit. The result is

Step 4 — Inverse QFT on the control register. This is the phase-estimation readout step: apply $\mathrm{QFT}^{-1}$ to the $t$ control qubits and measure them in the computational basis, yielding a $t$-bit integer $y$.

Step 5 — Classical continued-fraction recovery. Compute $\varphi = y / 2^t$ and expand it as a continued fraction; take successive convergents $p_i/q_i$ until you find one with $q_i < N$ that is a strong candidate for $s/r$ in lowest terms. Output the denominator $q_i$ as a candidate order $r_{\text{cand}}$.

Step 6 — Verify and repeat. Check $a^{r_{\text{cand}}} \bmod N \stackrel{?}{=} 1$. If not (or if $\gcd(s, r) > 1$ made $q_i$ a proper divisor of $r$ rather than $r$ itself), repeat the whole procedure a small number of times and take the least common multiple of the candidate denominators obtained, then verify that LCM against $a$. ``` function ORDER-FIND(a, N): n ← ceil(log2(N)) t ← 2n + 1 # precision register size repeat up to M = O(log log N) times:

# --- quantum phase estimation stage --- prepare control register C of t qubits in |0...0> prepare target register T of n qubits in |1> apply H to every qubit of C for i in 0 .. t-1: apply controlled-U_a^(2^i) with control = C[i], target = T

# U_a^(2^i) |y> = |a^(2^i) * y mod N>, realized via

# repeated modular squaring, computed reversibly apply inverse-QFT to C measure C, obtaining integer y in [0, 2^t - 1]

# --- classical post-processing stage --- phi ← y / 2^t (p, q) ← BEST-CONVERGENT(phi, max_denominator = N)

# continued-fraction expansion of phi,

# returns the convergent p_i/q_i with q_i < N

# closest to phi candidates.append(q) r ← LCM(candidates) if a^r mod N == 1 and r > 0: return r else: return FAIL # retry ORDER-FIND from scratch ``` ``` function BEST-CONVERGENT(phi, max_denominator):

# standard continued fraction expansion of phi a0 ← floor(phi) terms ← [a0] x ← phi - a0 p_prev, q_prev ← 1, 0 p_cur, q_cur ← a0, 1 while q_cur <= max_denominator and x != 0: x ← 1 / x ai ← floor(x) p_next ← ai p_cur + p_prev q_next ← ai q_cur + q_prev if q_next > max_denominator: break p_prev, q_prev ← p_cur, q_cur p_cur, q_cur ← p_next, q_next x ← x - ai return (p_cur, q_cur) ``` The controlled-$U_a^{2^i}$ gates are the only "arithmetic" part of the circuit: each is a reversible modular multiplication by a fixed classical constant $a^{2^i} \bmod N$ (itself precomputed classically, once, by repeated squaring), implemented with reversible adders and reversible modular reduction circuits — standard building blocks with known polynomial-size constructions.

Correctness — why the phases are exact multiples of $1/r$. Because $U_a$ restricted to the cyclic set $\{a^0, a^1, \dots, a^{r-1}\} \bmod N$ is a cyclic permutation of order $r$, its eigenvalues on the span of these basis states are exactly the $r$-th roots of unity $e^{2\pi i k / r}$ for $k = 0, \dots, r-1$, with eigenvectors

A direct check confirms $U_a \vert u_k\rangle = e^{2\pi i k/r}\vert u_k\rangle$. The initial state $\vert 1 \rangle$ decomposes as the equal-weight superposition $\vert 1\rangle = \frac{1}{\sqrt r}\sum_{k=0}^{r-1} \vert u_k\rangle$. Running phase estimation on this input is, by linearity, equivalent to running it on a uniformly random $\vert u_k\rangle$: the measured control register collapses to an estimate of $k/r$ for a uniformly random $k \in \{0, \dots, r-1\}$.

Why continued fractions recover $r$. Phase estimation with $t$ control qubits produces $y$ such that $y/2^t$ approximates $k/r$ to within $1/2^t$ with probability at least $\approx 4/\pi^2$ (the standard phase estimation error bound), and with the choice $t = 2n+1 \ge 2\log_2 N + 1$, the approximation error satisfies

because $r < N \le 2^n$ implies $2r^2 < 2 \cdot 2^{2n} \le 2^t$. A classical theorem on continued fractions guarantees: if $|\varphi - p/q| < 1/(2q^2)$ for coprime integers $p, q$, then $p/q$ is necessarily one of the convergents of the continued-fraction expansion of $\varphi$. Since $k/r$ (in lowest terms, dividing out $\gcd(k,r)$) satisfies this bound, the convergent-search step is guaranteed to encounter it. The algorithm scans convergents with denominator $< N$ and returns the appropriate one — this recovers $s/r_0$ where $r_0 = r/\gcd(k,r)$, a divisor of $r$, not necessarily $r$ itself.

Why repetition fixes the $\gcd(k,r) > 1$ issue. If $k$ happens to share a common factor with $r$, the recovered denominator $r_0$ is a proper divisor of $r$, not $r$. However, $k$ is uniform on $\{0,\dots,r-1\}$, and it is a standard fact that a uniformly random $k$ satisfies $\gcd(k, r) = 1$ with probability at least $\Omega(1/\log\log r)$ (via bounds on the density of integers with few or large prime factors relative to $\varphi(r)/r$ type arguments — concretely, the probability that $k/r$ reduces to something with denominator exactly $r$ is bounded below by a constant fraction over $O(\log\log r)$ independent trials in expectation). Running the quantum stage $O(\log\log N)$ times and taking the LCM of the resulting candidate denominators drives the failure probability down to an arbitrarily small constant, and a final classical check $a^{r} \bmod N = 1$ certifies success — false accepts are impossible since verification is exact.

Time complexity.

• Quantum circuit size: $t = O(n)$ control qubits, and $t$ controlled modular-exponentiation gates, each a reversible modular multiplication circuit of size $O(n^2)$ (schoolbook) or $O(n \log n \log\log n)$ (fast multiplication), giving total gate count $O(n^3)$ (or $O(n^2 \log n \log \log n)$ with fast arithmetic).
• Inverse QFT: $O(t^2) = O(n^2)$ gates (or $O(n \log n)$ with the approximate QFT).
• Classical continued-fraction expansion: $O(n)$ Euclidean-algorithm-style steps, each with $O(n^2)$-time big-integer arithmetic (or $O(n \log n)$ with fast algorithms), so $O(n^3)$ total — negligible compared to the quantum circuit.
• Repetitions: $O(\log \log N) = O(\log n)$, an essentially negligible multiplicative overhead.

Overall, order finding runs in $O(n^3)$ elementary gates ($n = \log_2 N$), i.e. polynomial in $\log N$ — exponentially faster than the best known classical algorithms, which require time super-polynomial in $\log N$ (subexponential, roughly $L$-notation complexity like general number-field-sieve-style bounds for related problems, though order finding specifically has no known polynomial classical algorithm and no proof of hardness either — it is an unproven, empirically robust separation).

Space complexity. $O(n)$ qubits: $t = O(n)$ for the control register and $n$ for the target register, plus $O(n)$ ancilla qubits for the reversible arithmetic circuits (garbage collection of ancillas is possible via standard reversible-computation techniques, keeping the total workspace linear in $n$).

Honest tradeoffs. The polynomial gate count hides a large constant: fault-tolerant modular exponentiation on cryptographically-sized $N$ (2048-bit RSA moduli) requires millions of logical qubit-operations and, with current error-correction overheads, potentially billions of physical gate operations. The algorithm is probabilistic — repetition and verification are mandatory, not optional — and its correctness relies on classical arithmetic (GCD, LCM, verification) being cheap and exact, which they are, but errors in the quantum part (decoherence, gate infidelity) are not corrected by the classical post-processing and must be handled by Quantum Error Correction beneath this algorithm.

Take $N = 21$, $a = 2$. First check $\gcd(2, 21) = 1$: yes. The true order is found by direct computation:

So $r = 6$: the true order of $2$ mod $21$ is $6$. Now trace the quantum algorithm conceptually. We need $n = \lceil \log_2 21 \rceil = 5$ qubits for the target register, and $t = 2n+1 = 11$ control qubits, giving $2^{11} = 2048$ possible measurement outcomes for $y$. The eigenphases of $U_2$ on the cycle $\{1,2,4,8,16,11\}$ are $e^{2\pi i k/6}$ for $k = 0,\dots,5$. Suppose the measurement collapses to $k = 5$ (a case with $\gcd(5,6)=1$, the "clean" case). Phase estimation ideally would give $y/2^{11} \approx 5/6 = 0.8\overline{3}$. Concretely $5/6 \cdot 2048 = 1706.67$, so a likely measured value is $y = 1707$, giving

Run the continued fraction expansion of $\varphi = 1707/2048$:

The first convergent with denominator $< 21$ that closely matches is $5/6$, obtained at step 2, exactly the true $k/r = 5/6$. The algorithm outputs candidate denominator $q = 6$.

Verification. Check $2^6 \bmod 21 = 64 \bmod 21 = 1$. ✓ Success — $r = 6$ is confirmed on the first attempt, no repetition needed. Now consider a "unlucky" branch: suppose instead the measurement had collapsed to $k = 2$, where $\gcd(2,6) = 2$. Then $k/r = 2/6 = 1/3$ in lowest terms, and the continued-fraction search returns $q = 3$, a proper divisor of the true order. A single run reports the wrong order ($2^3 \bmod 21 = 8 \ne 1$, so the verification step would actually catch this directly!). This is exactly why Step 6 (verify, and if needed repeat and take LCMs) exists: with two independent runs, one landing on $k=5$ (denominator 6) and another on $k=4$ (giving $2/3$ in lowest terms, denominator 3), $\mathrm{lcm}(6,3) = 6$ recovers the correct order even if a single run alone was ambiguous.

Order finding is not deployed in production systems today, because no quantum computer yet has enough error-corrected qubits to run it on cryptographically relevant moduli — but its threat model already shapes real infrastructure decisions. Standards bodies such as NIST have run multi-year competitions to select post-quantum cryptographic algorithms (lattice-based schemes like Kyber and Dilithium) precisely because Shor's algorithm, built on this subroutine, would break RSA and elliptic-curve cryptography outright given a sufficiently large fault-tolerant quantum computer. Financial institutions, governments, and browser vendors are actively migrating TLS and VPN protocols to hybrid classical/post-quantum schemes today, years before such a machine might exist, as a direct consequence of this algorithm's existence. On the research side, order finding (and its close cousin, the abelian Hidden Subgroup Problem) is the canonical example used to teach the power of the quantum Fourier transform, and every subsequent quantum algorithm with exponential speedups over classical methods for algebraic problems — discrete logarithm, Pell's equation, certain hidden shift problems — reuses the same phase-estimation-plus-continued-fractions template. Experimental groups building superconducting, trapped-ion, and photonic quantum processors use small-scale factoring demonstrations (factoring 15, 21, and similar small numbers) as benchmark tasks that exercise controlled modular arithmetic, phase estimation, and measurement — making order finding one of the most-implemented "hard" quantum subroutines in the experimental literature, even at toy scale.

Ask the tutor: How does the choice of $a$ affect the probability that order finding succeeds on a given attempt, and why does averaging over random $a$ matter for the factoring reduction?