26. Shor's Algorithm
26.1Integer factorization
Given an odd composite N, find a nontrivial factor (then recurse). This is the problem behind RSA: multiplying two large primes is instant, factoring the product is (as far as anyone knows) hard, and every internet key exchange of the last four decades has bet on that asymmetry. Shor's 1994 algorithm factors N in polynomial time — O(n³)-class gates for n-bit N on a quantum computer — and remains the single most consequential result in the field: it converted quantum computing from a curiosity into a cryptographic deadline. Note what Shor does not do: it does not search for factors; it finds periods, and factors fall out.
26.2Classical difficulty
The best classical algorithm, the general number field sieve, runs in exp(((64/9)^{1/3} + o(1))·(ln N)^{1/3}·(ln ln N)^{2/3}) — subexponential, superpolynomial. RSA-2048 sits beyond it by an enormous margin (records: RSA-250 in 2020 took ~2700 core-years). No polynomial classical factoring algorithm is known despite centuries of work, which is why the problem is a trustworthy hardness assumption — and why Shor's polynomial-time quantum algorithm is a genuine complexity-theoretic event, not an incremental speedup. The number-theoretic detour through periods is what makes it work: factoring and period finding are, it turns out, the same problem wearing different clothes.
26.3Modular arithmetic
Pick a random a with gcd(a, N) = 1 (if the gcd is bigger, you already found a factor — lucky draw). Work in the multiplicative group Z*_N of residues coprime to N, which has size φ(N) (Euler's totient). By Lagrange's theorem, every element has a finite order: the smallest r with a^r ≡ 1 (mod N), and r divides φ(N). The powers a^0, a^1, a^2, … therefore cycle with period r. For N = pq with distinct odd primes, a random a has even order with a^{r/2} ≢ −1 (mod N) with probability ≥ 1/2 — the condition that makes the final gcd step split N.
26.4Period finding
The reduction: if you can find the order r of a random a mod N, you can factor. When r is even and a^{r/2} ≢ −1 (mod N), write x = a^{r/2}; then x² ≡ 1 (mod N) while x ≢ ±1, so N divides x²−1 = (x−1)(x+1) without dividing either factor — and gcd(x−1, N) and gcd(x+1, N) are nontrivial factors. Classically, finding r means computing powers until one hits 1: O(r) ≤ O(N) multiplications, or O(√r) with baby-step giant-step. For RSA-scale N, r is astronomically large — period finding is the exponential wall, and it is exactly the wall the QFT dismantles.
26.5Quantum period finding
The quantum core. Choose Q = 2^t with N² ≤ Q < 2N² (t ≈ 2n+1). Prepare (1/√Q)·Σ_{x=0}^{Q−1}|x⟩|1⟩, and apply the modular-exponentiation unitary |x⟩|y⟩ ↦ |x⟩|y·a^x mod N⟩ — built from repeated squaring along the binary digits of x, the O(n³)-gate behemoth:
reg 1 (t ≈ 2n qubits): |0> ──H^⊗t────────[QFT over Q]──measure y
│
reg 2 (n qubits): |1> ──────U_pow── U_pow: |y> ↦ |y·a^x mod N>Now register 2 holds a^x mod N — periodic with period r — entangled with register 1. The entanglement is the point: register 1 is a coherent superposition of all x, tagged by the period structure of their images.
26.6QFT connection
Apply the QFT over Q to register 1. The periodic train in the entangled state diffracts: amplitude concentrates on integers y with |y/Q − k/r| ≤ 1/(2Q) for some k — multiples of Q/r, up to rounding. Measure: y/Q is a fraction within 1/(2Q) of k/r. Now the classical step — the continued fraction expansion of y/Q with denominators bounded by N recovers k/r whenever gcd(k, r) = 1, which happens for a constant fraction of shots (the peaks carry probability ≥ 4/π²·φ(r)/r in total). Verify a^r ≡ 1 (mod N); if verification fails (you caught r/gcd(k,r), or an odd r), resample — a constant expected number of rounds. This is Simon's skeleton with (Z, +) in place of (Z₂)ⁿ.
26.7The complete algorithm
Everything assembled — the quantum part is one boxed line:
import math, random
def shor(N): # classical driver
if N % 2 == 0:
return 2
while True:
a = random.randrange(2, N - 1)
g = math.gcd(a, N)
if g > 1:
return g # lucky draw: factor for free
r = quantum_order(a, N) # QFT period finding (26.5-26.6)
if r % 2 != 0:
continue # odd order: retry with new a
x = pow(a, r // 2, N)
if x == N - 1:
continue # a^(r/2) = -1: retry
f1, f2 = math.gcd(x - 1, N), math.gcd(x + 1, N)
if 1 < f1 < N:
return f1 # recurse on f1 and N // f1quantum_order(a, N) runs the 26.5–26.6 circuit a constant number of times. Total cost: O(n³) gates naively (O(n²·polylog n) with fast arithmetic), O(1) expected repetitions — each repetition factors N with probability ≥ 1/2.
26.8Toy implementation
N = 15, a = 7 — the standard end-to-end example. The powers cycle with period 4; the toy uses Q = 16 (smaller than the required Q ≥ N², but exact here because r = 4 divides 16):
x : 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
7^x mod 15 : 1 7 4 13 1 7 4 13 1 7 4 13 1 7 4 13
└── period r = 4 ──┘
QFT over Q = 16 → peaks at y ∈ {0, 4, 8, 12} (each with probability 1/4)
y/Q = 4/16, 8/16, 12/16 = 1/4, 1/2, 3/4 → continued fractions → r = 4
7^(r/2) = 7^2 = 49 ≡ 4 (mod 15), 4 ≢ −1 ≡ 14 → condition satisfied
gcd(4 − 1, 15) = 3 and gcd(4 + 1, 15) = 5 → 15 = 3 × 5And the full simulation — oracle, QFT, measurement, continued fractions, gcd:
import numpy as np
from fractions import Fraction
rng = np.random.default_rng(3)
N, a, Q = 15, 7, 16
f = np.array([pow(a, x, N) for x in range(Q)]) # 7^x mod 15
psi = np.zeros(Q * N, dtype=complex)
for x in range(Q):
psi[x * N + f[x]] = 1 # |x>|a^x mod N>
H = np.exp(2j * np.pi * np.outer(np.arange(Q), np.arange(Q)) / Q) / np.sqrt(Q)
psi = np.kron(H, np.eye(N)) @ psi # QFT over Q on register 1
p = (np.abs(psi.reshape(Q, N)) ** 2).sum(axis=1)
y = rng.choice(Q, p=p) # measure register 1
r = Fraction(y, Q).limit_denominator(N).denominator # continued fractions
print("y =", y, " r candidate =", r)
if r % 2 == 0 and pow(a, r, N) == 1:
print("factors:", np.gcd(pow(a, r // 2, N) - 1, N), np.gcd(pow(a, r // 2, N) + 1, N))
else:
print("unusable shot (y=0 or k shares a factor with r) — rerun")The else branch is not decoration: y = 0 (probability ¼) and y = 8 (which yields 1/2, i.e. r/gcd(k,r) = 2) both land there, exactly as 26.6 promised.
from fractions import Fraction
for y in (4, 8, 12): # measured peaks, Q = 16, true r = 4
fr = Fraction(y, 16).limit_denominator(15)
print(f"{y}/16 -> {fr} r-candidate = {fr.denominator}")26.9Resource requirements
Logical circuit: t ≈ 2n+1 counting qubits, n-qubit arithmetic registers, QFT O(n²) gates, and modular exponentiation dominating at O(n³) gates naively — O(n²·polylog n) with advanced multiplication circuits. End-to-end estimates for RSA-2048: Gidney–Ekerå (2019) — about 20 million noisy physical qubits running 8 hours; Gidney (2025) — under 1 million noisy qubits running under a week, thanks to better arithmetic and error-correction layouts. Both assume surface-code fault tolerance at physical error rates near 10⁻³. Against today's machines — hundreds to a few thousand physical qubits — the gap is roughly three orders of magnitude in qubit count. That is the honest planning horizon: reachable, but not imminent, engineering.
26.10Implications for cryptography
Shor breaks RSA, finite-field Diffie–Hellman, and elliptic-curve cryptography once fault-tolerant machines at the 26.9 scale exist. Symmetric primitives survive: Grover costs only a factor of 2 in effective key length, so AES-256 remains sound. The operative threat is harvest now, decrypt later — adversaries recording encrypted traffic today to decrypt after a CRQC exists — which makes migration urgent for long-lived secrets. NIST standardized the replacements in 2024: ML-KEM (FIPS 203, lattices), ML-DSA (FIPS 204), SLH-DSA (FIPS 205, hash-based). The migration is a decade-scale program touching every protocol you have ever shipped — and, closer to home, it is where cryptanalysis careers meet quantum computing (51.9). The algorithm is 30 years old; its consequences are still being built.