Chapter 15: Shor's Algorithm — Factoring Large Numbers in Polynomial Time (and Why It Breaks RSA Encryption)

Learning Objectives

By the end of this chapter, you will be able to:

• Explain the reduction from integer factoring to period finding.
• Construct the modular exponentiation circuit for a given $N$ and $a$.
• Apply the quantum Fourier transform to extract the period $r$ from a quantum state.
• Perform the continued fraction algorithm to recover $r$ from a phase estimate.
• Implement Shor's algorithm in Qiskit to factor $N = 15$.
• Analyze the complexity of Shor's algorithm and its implications for RSA cryptography.
• Discuss the state of experimental demonstrations of Shor's algorithm.
• Derive the probability distribution of measurement outcomes in period finding.
• Understand why Shor's algorithm works from both linear algebra and number theory perspectives.
• Evaluate the resource requirements for factoring cryptographically relevant integers.

15.1 The Factoring Problem and RSA

15.1.1 Why Factoring Matters

The security of the RSA public-key cryptosystem — used to secure HTTPS connections, email, and digital signatures for billions of people — rests entirely on a single empirical fact: multiplying two large primes is easy; factoring their product is hard. Given two 1024-bit primes $p$ and $q$, a classical computer can compute $N = pq$ in microseconds. But given only $N$, the best known classical algorithm — the general number field sieve — runs in sub-exponential time:

$$T_{\text{GNFS}}(N) = \exp\left( \left( \sqrt[3]{\frac{64}{9}} + o(1) \right) (\ln N)^{1/3} (\ln \ln N)^{2/3} \right)$$

For a 2048-bit modulus, this translates to billions of core-years. The hardness of factoring is not proven — no one has shown that factoring lies outside $\mathbf{P}$ — but decades of effort by the world's best number theorists have failed to find a polynomial-time classical algorithm.

Historical Context. The RSA cryptosystem was publicly proposed in 1977 by Ron Rivest, Adi Shamir, and Leonard Adleman, building on earlier classified work by Clifford Cocks at GCHQ (1973). Its security assumption — that factoring is hard — had been studied since antiquity. Eratosthenes devised the sieve algorithm around 200 BCE. The modern era of factoring began with Morrison and Brillhart's continued fraction method (1975), which factored the 39-digit Fermat number $F_7$. Each subsequent algorithm — quadratic sieve, elliptic curve method, number field sieve — pushed the boundary further, but the sub-exponential scaling remained. Shor's 1994 result changed everything by showing that a quantum computer could factor in polynomial time, an exponential speedup over the best known classical algorithms.

Why should you care? Every time you visit an HTTPS website, your browser performs an RSA (or ECC) key exchange. If someone built a large-scale quantum computer tomorrow, they could decrypt the vast majority of encrypted internet traffic. This isn't a hypothetical concern — intelligence agencies are reportedly storing encrypted data today in anticipation of future quantum capabilities, a strategy called "harvest now, decrypt later."

Common Misconception: "Shor's algorithm breaks all cryptography."

Shor's algorithm breaks public-key cryptography based on factoring (RSA) or discrete logarithms (ECC, Diffie-Hellman). It does not break symmetric-key cryptography (AES, ChaCha20). Grover's algorithm provides only a quadratic speedup against symmetric ciphers, which can be neutralized by doubling the key size. The real threat is to public-key infrastructure — the system of certificate authorities, key exchange, and digital signatures that underpins internet security.

15.1.2 The RSA Trapdoor

RSA works as follows. Alice chooses two large primes $p, q$, computes $N = pq$, and selects a public exponent $e$ coprime to $\phi(N) = (p-1)(q-1)$. She publishes $(N, e)$. Bob encrypts a message $m$ as $c = m^e \bmod N$. Alice decrypts using her private key $d = e^{-1} \bmod \phi(N)$ by computing $c^d \bmod N = m$.

Why does decryption work? By Euler's theorem, $m^{\phi(N)} \equiv 1 \pmod{N}$ when $\gcd(m, N) = 1$. Since $ed \equiv 1 \pmod{\phi(N)}$, we have $ed = k\phi(N) + 1$ for some integer $k$. Therefore:

$$c^d = (m^e)^d = m^{ed} = m^{k\phi(N) + 1} = (m^{\phi(N)})^k \cdot m \equiv 1^k \cdot m = m \pmod{N}$$

An eavesdropper who can factor $N$ into $p$ and $q$ can compute $\phi(N) = (p-1)(q-1)$ and then $d = e^{-1} \bmod \phi(N)$, breaking the scheme. Shor's algorithm does exactly this — in polynomial time on a quantum computer.

Try It Yourself: RSA Key Generation

Choose $p = 61$, $q = 53$. Then $N = 3233$, $\phi(N) = 60 \times 52 = 3120$. Choose $e = 17$ (coprime to 3120). Compute $d = e^{-1} \bmod 3120 = 2753$. Verify: $17 \times 2753 = 46801 = 15 \times 3120 + 1 \equiv 1 \pmod{3120}$. Encrypt $m = 65$: $c = 65^{17} \bmod 3233 = 2790$. Decrypt: $2790^{2753} \bmod 3233 = 65$. ✓

15.1.3 Reduction of Factoring to Period Finding

Shor's key insight (1994) was that factoring reduces to finding the period of the function

$$f(x) = a^x \bmod N$$

where $a$ is a random integer coprime to $N$. The function $f(x)$ is periodic with period $r$ (the order of $a$ modulo $N$):

$$a^{x+r} \equiv a^x \pmod{N}, \quad a^r \equiv 1 \pmod{N}$$

If we can find $r$, and if $r$ is even and $a^{r/2} \not\equiv -1 \pmod{N}$, then:

$$\gcd(a^{r/2} - 1, N) \quad \text{and} \quad \gcd(a^{r/2} + 1, N)$$

are non-trivial factors of $N$.

Why does this work? Since $a^r \equiv 1 \pmod{N}$, we can write $a^r - 1 \equiv 0 \pmod{N}$, which factors as:

$$(a^{r/2} - 1)(a^{r/2} + 1) \equiv 0 \pmod{N}$$

This means $N$ divides $(a^{r/2} - 1)(a^{r/2} + 1)$. If $a^{r/2} \not\equiv \pm 1 \pmod{N}$, then neither factor is $0 \bmod N$, meaning each factor shares a non-trivial common divisor with $N$. The GCD extracts that divisor.

Common Misconception: "Shor's algorithm finds factors directly."

Shor's algorithm does not find factors directly. It finds the order $r$ of an element $a$ modulo $N$ using a quantum subroutine. The factors are then extracted classically via GCD. The quantum part solves period finding; the classical part converts a period into factors.

What fraction of random $a$ yield useful factors? For $N = p^k$ or $N$ with multiple distinct prime factors, at least half of all $a \in [2, N-2]$ coprime to $N$ produce a useful $r$. This means the expected number of trials is at most 2. The proof relies on the structure of $\mathbb{Z}_N^*$ and can be found in standard number theory references.

Worked Example: Factoring $N = 15$

Choose $a = 2$. Compute powers of 2 mod 15: - $2^0 \bmod 15 = 1$ - $2^1 \bmod 15 = 2$ - $2^2 \bmod 15 = 4$ - $2^3 \bmod 15 = 8$ - $2^4 \bmod 15 = 16 \bmod 15 = 1$ ✓

The period is $r = 4$. Since $r$ is even: $a^{r/2} = 2^2 = 4$. Check: $4 \not\equiv -1 \pmod{15}$ (since $-1 \equiv 14 \pmod{15}$). Now: - $\gcd(4 - 1, 15) = \gcd(3, 15) = 3$ ✓ - $\gcd(4 + 1, 15) = \gcd(5, 15) = 5$ ✓

We found $15 = 3 \times 5$.

Worked Example: Factoring $N = 21$

Choose $a = 2$. Compute powers of 2 mod 21: - $2^0 = 1$, $2^1 = 2$, $2^2 = 4$, $2^3 = 8$, $2^4 = 16$, $2^5 = 11$, $2^6 = 1$

The period is $r = 6$. Since $r$ is even: $a^{r/2} = 2^3 = 8$. Check: $8 \not\equiv -1 \pmod{21}$ (since $-1 \equiv 20$). Now: - $\gcd(8 - 1, 21) = \gcd(7, 21) = 7$ ✓ - $\gcd(8 + 1, 21) = \gcd(9, 21) = 3$ ✓

We found $21 = 7 \times 3$.

Worked Example: When the period is odd

Try $N = 35$, $a = 2$. Powers: $2^1 = 2, 2^2 = 4, 2^3 = 8, 2^4 = 16, 2^5 = 32, 2^6 = 29, 2^7 = 23, 2^8 = 11, 2^9 = 22, 2^{10} = 9, 2^{11} = 18, 2^{12} = 1$. Period $r = 12$ (even). $2^6 \bmod 35 = 29$. $\gcd(28, 35) = 7$ and $\gcd(30, 35) = 5$. ✓

But if we had tried $a = 7$: powers are $7^1 = 7, 7^2 = 14, 7^3 = 28, 7^4 = 21, 7^5 = 7, 7^6 = 14, ...$. This repeats with period $r = 4$. However, $7^2 \bmod 35 = 14$. Check: $\gcd(13, 35) = 1$ and $\gcd(15, 35) = 5$. We get one factor but not two — $7$ shares a factor with 35 (since $\gcd(7, 35) = 7 > 1$, we would have caught this at step 4 of the algorithm). This is why step 4 checks $\gcd(a, N)$ first.

Algorithm 15.1: Shor's Algorithm (Classical Reduction)

1. If $N$ is even, return $2$.
2. If $N = p^k$ for prime $p$, use classical root-finding.
3. Choose random $a \in [2, N-2]$.
4. Compute $g = \gcd(a, N)$. If $g > 1$, return $g$ (lucky!).
5. Use the quantum period-finding subroutine to find the order $r$ of $a \bmod N$.
6. If $r$ is odd or $a^{r/2} \equiv -1 \pmod{N}$, go back to step 3.
7. Compute $\gcd(a^{r/2} \pm 1, N)$; these are non-trivial factors.

Recurring Theme: Quantum is Linear Algebra, Not Magic

The quantum part of Shor's algorithm solves a linear algebra problem — finding the period of a function. The QFT is essentially a change of basis (from the computational basis to the Fourier basis) that reveals periodicity. There is nothing magical happening; the quantum computer is performing a specific linear transformation on a high-dimensional vector space that happens to extract the period efficiently. The speedup comes from the fact that an $n$-qubit system naturally lives in a $2^n$-dimensional space, and the QFT transforms this space in $O(n^2)$ operations — exponentially faster than the classical FFT, which requires $O(2^n n)$ operations on a $2^n$-dimensional vector.

15.2 Quantum Period Finding

15.2.1 The Core Idea

We want to find $r$ such that $a^r \equiv 1 \pmod{N}$. The quantum period-finding circuit uses two registers:

• Register 1 (control): $t = 2n$ qubits, where $n = \lceil \log_2 N \rceil$, initialized to $|0\rangle^{\otimes t}$.
• Register 2 (target): $n$ qubits, initialized to $|1\rangle$ (representing $a^0 \bmod N$).

The circuit proceeds in four stages:

1. Superposition: Apply $H^{\otimes t}$ to Register 1.
2. Modular exponentiation: Apply the unitary $U_f$ where $U_f |x\rangle|y\rangle = |x\rangle|y \cdot a^x \bmod N\rangle$.
3. Quantum Fourier transform: Apply $\text{QFT}^\dagger$ to Register 1.
4. Measurement: Measure Register 1 to obtain a phase estimate.

Why $t = 2n$ qubits? The choice of $t$ is crucial. With $n = \lceil \log_2 N \rceil$ qubits, we can represent values up to $2^n$. But we need to resolve the period $r$ to within $\pm 1/2$ with high probability. By the analysis of QPE precision (Chapter 16), we need $t \geq 2n + 1 + \lceil \log_2(2 + 1/\epsilon) \rceil$ qubits for success probability $1 - \epsilon$. The choice $t = 2n$ gives success probability at least $4/\pi^2 \approx 40.5\%$ per trial, which is sufficient since we can repeat.

15.2.2 Step-by-Step State Evolution

Let us trace the quantum state through each stage of the algorithm.

After Step 1 (Superposition):

$$|\Psi_1\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t - 1} |x\rangle |1\rangle$$

After Step 2 (Modular Exponentiation):

$$|\Psi_2\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t - 1} |x\rangle |a^x \bmod N\rangle$$

Since $f(x) = a^x \bmod N$ is periodic with period $r$, we can group terms. Let $m = \lfloor 2^t / r \rfloor$ (the number of complete periods that fit in $[0, 2^t - 1]$). For each residue class $y \in \{0, 1, \ldots, r-1\}$, the values of $x$ giving $a^x \bmod N = a^y$ are $x \in \{y, y+r, y+2r, \ldots, y+(m-1)r\}$ (plus a remainder term if $2^t$ is not divisible by $r$):

$$|\Psi_2\rangle \approx \frac{1}{\sqrt{m r}} \sum_{y=0}^{r-1} \sum_{j=0}^{m-1} |jr + y\rangle |a^y \bmod N\rangle$$

The approximation ignores the incomplete last period; when $2^t \gg r$ (which is ensured by $t = 2n$), this is negligible.

After Step 3 (Inverse QFT on Register 1):

The key insight is that the state $\frac{1}{\sqrt{m}} \sum_{j=0}^{m-1} |jr + y\rangle$ in Register 1 is nearly an eigenstate of the QFT. Applying $\text{QFT}^\dagger_{2^t}$:

$$\text{QFT}^\dagger_{2^t} |jr + y\rangle = \frac{1}{\sqrt{2^t}} \sum_{k=0}^{2^t - 1} e^{-2\pi i (jr+y) k / 2^t} |k\rangle$$

The combined state becomes:

$$|\Psi_3\rangle \approx \frac{1}{\sqrt{r}} \sum_{y=0}^{r-1} \sum_{k=0}^{r-1} e^{2\pi i k y / r} \left| k \cdot \frac{2^t}{r} \right\rangle |a^y \bmod N\rangle$$

Wait — let's be more precise. The measurement outcome on Register 1 will be concentrated near integer multiples of $2^t / r$. Specifically, measuring Register 1 yields some integer $c$ such that:

$$\left| \frac{c}{2^t} - \frac{k}{r} \right| < \frac{1}{2^{t+1}}$$

with probability at least $4/\pi^2 \approx 0.405$ for some $k$ coprime to $r$.

ASCII Diagram: Measurement Probability Distribution for Period Finding
======================================================================

Probability
    |       |       |       |       |       |       |       |
    |  ██   |       |  ██   |       |  ██   |       |  ██   |
    |  ██   |       |  ██   |       |  ██   |       |  ██   |
    |  ██   |  ..   |  ██   |  ..   |  ██   |  ..   |  ██   |
    |  ██   |  ..   |  ██   |  ..   |  ██   |  ..   |  ██   |
    |  ██   |  ..   |  ██   |  ..   |  ██   |  ..   |  ██   |
----+-------+-------+-------+-------+-------+-------+-------+---→ c
    0     2^t/r    2·2^t/r  3·2^t/r          (r-1)·2^t/r  2^t

Peaks at c ≈ k · 2^t/r for k = 0, 1, ..., r-1
The spacing between peaks is 2^t/r, which reveals the period r.

15.2.3 The Modular Exponentiation Circuit

The unitary $U_a$ acts as $U_a |x\rangle|y\rangle = |x\rangle|y \cdot a^x \bmod N\rangle$. Implementing this efficiently requires decomposing $x$ into its binary representation:

$$x = x_0 + 2x_1 + 4x_2 + \cdots + 2^{t-1} x_{t-1}$$

Then:

$$a^x \bmod N = \prod_{j=0}^{t-1} a^{2^j x_j} \bmod N = \prod_{j: x_j = 1} a^{2^j} \bmod N$$

This means we can implement modular exponentiation as a sequence of controlled modular multiplications: for each qubit $j$ in Register 1, if the qubit is $|1\rangle$, multiply Register 2 by $a^{2^j} \bmod N$. The values $a^{2^j} \bmod N$ can be precomputed classically via repeated squaring.

ASCII Circuit Diagram: Shor's Period-Finding Core
=================================================

Register 1 (t qubits)          Register 2 (n qubits)
|0⟩ ──H───●───────────────────●─────●───[QFT†]─── Measure
|0⟩ ──H───│───────────────────●─────│──────────── Measure
|0⟩ ──H───│─────────────────────────●──────────── Measure
  ...      │                    │     │              ...
|0⟩ ──H───│─────────────────────────●──────────── Measure
           │                    │     │
|1⟩ ──────[×a^(2^0)]───[×a^(2^1)]──[×a^(2^{t-1})]─── (ignored)

Each controlled operation multiplies Register 2 by a^(2^j) mod N
when the corresponding control qubit is |1⟩.

Detailed circuit for $N = 15$, $a = 7$:

The constants $a^{2^j} \bmod 15$ are: - $a^{2^0} = 7^1 \bmod 15 = 7$ - $a^{2^1} = 7^2 \bmod 15 = 49 \bmod 15 = 4$ - $a^{2^2} = 7^4 \bmod 15 = 4^2 \bmod 15 = 1$ (since $7^4 \equiv 1 \pmod{15}$, so all higher powers are trivial)

This is very convenient: for $N = 15$, $a = 7$, the order is $r = 4$, so $a^{2^2} \equiv 1 \pmod{15}$ and the controlled operations beyond the first two become trivial.

ASCII Diagram: Modular Multiplication Unitary for N=15, a=7
============================================================

For ×7 mod 15 (4-qubit target):
  |0⟩ → |0⟩      |4⟩ → |13⟩     |8⟩ → |11⟩     |12⟩ → → |9⟩
  |1⟩ → |7⟩      |5⟩ → |5⟩      |9⟩ → |3⟩      |13⟩ → → |1⟩
  |2⟩ → |14⟩     |6⟩ → |12⟩     |10⟩ → |10⟩    |14⟩ → → |14⟩
  |3⟩ → |6⟩      |7⟩ → |4⟩      |11⟩ → |2⟩     |15⟩ → → |15⟩

This permutation can be decomposed into elementary gates.

15.2.4 Full Derivation of the Measurement Probability

Let us derive the probability of obtaining measurement outcome $c$ on Register 1 after the inverse QFT. Starting from:

$$|\Psi_2\rangle = \frac{1}{\sqrt{2^t}} \sum_{x=0}^{2^t - 1} |x\rangle |a^x \bmod N\rangle$$

and discarding Register 2 (since it's entangled with Register 1), the reduced state of Register 1 is:

$$\rho_1 = \frac{1}{2^t} \sum_{x,x'} |x\rangle \langle x'| \cdot \frac{1}{r} \sum_{y=0}^{r-1} e^{i \text{phase}(x-x')}$$

Wait — let's use the simpler route. Since the states $|a^y \bmod N\rangle$ for different $y$ are orthogonal, we can write:

$$|\Psi_2\rangle = \frac{1}{\sqrt{r}} \sum_{y=0}^{r-1} |\psi_y\rangle |a^y \bmod N\rangle$$

where $|\psi_y\rangle = \frac{1}{\sqrt{m}} \sum_{j=0}^{m-1} |jr + y\rangle$ is a "comb" state with teeth spaced by $r$.

After applying $\text{QFT}^\dagger_{2^t}$ to Register 1:

$$\text{QFT}^\dagger |\psi_y\rangle = \frac{1}{\sqrt{2^t \cdot m}} \sum_{c=0}^{2^t - 1} \left( \sum_{j=0}^{m-1} e^{-2\pi i c(jr+y)/2^t} \right) |c\rangle$$

The inner sum is a geometric series:

$$\sum_{j=0}^{m-1} e^{-2\pi i c j r / 2^t} = \frac{1 - e^{-2\pi i c m r / 2^t}}{1 - e^{-2\pi i c r / 2^t}}$$

This is sharply peaked when $c \approx k \cdot 2^t / r$ for integer $k$. The probability of measuring $c$ is:

$$P(c) \approx \frac{1}{r} \cdot \left| \frac{\sin(\pi c r / 2^t)}{r \sin(\pi c / 2^t)} \right|^2$$

normalized over the peak positions. The key point: with high probability, the measured $c$ satisfies $|c/2^t - k/r| < 1/2^{t+1}$ for some $k$.

15.3 The Continued Fraction Step

Given a measurement outcome $c$, we know that:

$$\frac{c}{2^t} \approx \frac{k}{r}$$

for some integer $k$ coprime to $r$. The continued fraction expansion of $c / 2^t$ yields a sequence of convergents $p_i / q_i$, each being the best rational approximation with denominator at most $q_i$. We test each $q_i$ as a candidate for $r$ by checking whether $a^{q_i} \equiv 1 \pmod{N}$.

The continued fraction algorithm computes the continued fraction expansion of a rational number. For $\phi = c/2^t$:

$$\phi = a_0 + \cfrac{1}{a_1 + \cfrac{1}{a_2 + \cfrac{1}{\ddots}}}$$

The convergents $p_i/q_i$ are computed by the recurrence:

$$p_i = a_i p_{i-1} + p_{i-2}, \quad q_i = a_i q_{i-1} + q_{i-2}$$

with $p_{-2} = 0, p_{-1} = 1, q_{-2} = 1, q_{-1} = 0$.

Algorithm 15.2: Continued Fractions

def continued_fraction(numerator, denominator):
    """Return list of convergents (p_i, q_i) for numerator/denominator."""
    convergents = []
    a_terms = []
    a, b = numerator, denominator
    while b != 0:
        a_terms.append(a // b)
        a, b = b, a % b
    p_prev, q_prev = 0, 1
    p_curr, q_curr = 1, 0
    for a_i in a_terms:
        p_next = a_i * p_curr + p_prev
        q_next = a_i * q_curr + q_prev
        convergents.append((p_next, q_next))
        p_prev, q_prev = p_curr, q_curr
        p_curr, q_curr = p_next, q_next
    return convergents

Worked Example: Continued Fractions for $N = 15$, $a = 7$

The order of $7 \bmod 15$ is $r = 4$. With $t = 8$ qubits, we measure $c$ such that $c/2^8 \approx k/4$. Possible measurements: - $c = 0$: $0/256 = 0/4$ → $k = 0$, gives $r = 4$ but $\gcd(0, 4) = 4$, so we'd need $k$ coprime to $r$. If $k = 0$, we get no information. - $c = 64$: $64/256 = 1/4$ → convergent is $1/4$, so $r = 4$. ✓ - $c = 128$: $128/256 = 1/2$ → convergent is $1/2$, so candidate $r = 2$. Check: $7^2 \bmod 15 = 4 \neq 1$. Not the period. Next convergent from $128/256$: reduce to $1/2$, then expand. Actually, we should look at the continued fraction expansion: $128/256 = 1/2$. Convergents: $1/2$. $7^2 = 49 \equiv 4 \pmod{15} \neq 1$. Try next. - $c = 192$: $192/256 = 3/4$ → convergent is $3/4$, so $r = 4$. ✓

Let's trace the continued fraction for $c = 64$: $64/256$: - $a_0 = 256/64 = 4$ remainder $0$ → $64/256 = 1/4$ - Convergents: $(1, 4)$ → test $q = 4$: $7^4 \bmod 15 = 2401 \bmod 15 = 1$. ✓ Found $r = 4$!

Worked Example: Continued Fractions for $N = 21$, $a = 2$

The order of $2 \bmod 21$ is $r = 6$. With $t = 10$ qubits ($2^t = 1024$), expected measurements: $c \approx k \cdot 1024/6$. - $k = 1$: $c \approx 171$, $c/1024 \approx 0.16699...$ - Continued fraction of $171/1024$: $0 + 1/(5 + 1/(117 + 1/3))$. Convergents: $0/1, 1/5, 117/590, 352/1775$... Wait, let me recalculate. - $171/1024$: $a_0 = 0$, $1024/171 = 5$ rem $169$, $171/169 = 1$ rem $2$, $169/2 = 84$ rem $1$, $2/1 = 2$ rem $0$. - Convergents: $0/1, 1/6, 1/7$... Hmm, $171/1024$ is approximately $1/6 = 0.1\overline{6}$. The convergent $1/6$ gives $q = 6$. Check: $2^6 \bmod 21 = 64 \bmod 21 = 1$. ✓ Found $r = 6$.

15.4 Complexity Analysis

15.4.1 Gate Complexity

The overall quantum circuit uses $O((\log N)^2 \log \log N \log \log \log N)$ gates — polynomial in the number of bits. For an $n$-bit number, Shor's algorithm runs in $O(n^3)$ time, compared to the sub-exponential $O(\exp(c n^{1/3}))$ of the best classical algorithm.

Detailed gate count for modular exponentiation: Each controlled modular multiplication by $a^{2^j} \bmod N$ requires $O(n^2)$ gates using the schoolbook multiplication algorithm, or $O(n \log n \log \log n)$ using fast multiplication. With $t = 2n$ such operations, the total is $O(n^3)$ or $O(n^2 \log n \log \log n)$.

Detailed QFT gate count: The QFT on $t = 2n$ qubits requires $O(t^2) = O(n^2)$ gates (each qubit has a Hadamard and $O(t)$ controlled-phase gates). The swaps at the end require $O(t)$ additional SWAP gates.

15.4.2 Qubit Requirements

To factor a 2048-bit RSA modulus, a quantum computer would need roughly $2 \times 2048 = 4096$ logical qubits for the control register plus $2048$ for the target, totaling about 6,000–10,000 logical qubits. With error correction overhead (surface codes at $\sim 10^{-3}$ physical error rate), this translates to millions of physical qubits — beyond current hardware but well within the roadmap for fault-tolerant quantum computers.

More detailed estimate (Gidney & Ekerå, 2021): Factoring a 2048-bit RSA integer requires approximately 20 million noisy qubits running for 8 hours, using: - 6,146 logical qubits (4,096 control + 2,048 target + ancilla) - ~3.5 × 10^6 physical qubits per logical qubit (with surface code at $10^{-3}$ physical error rate) - Total: ~20 million physical qubits - Windowed mode: can be reduced to ~2 million qubits with longer runtime

Recurring Theme: Noise is the Enemy

The 20-million-qubit estimate highlights why error correction is essential. Without it, a single error in billions of gate operations would corrupt the computation. The surface code requires approximately 1,000 physical qubits per logical qubit at current error rates, and the modular exponentiation circuit has depth $O(n^3)$, meaning billions of gate operations for RSA-2048. Each operation must succeed with logical error rate below $10^{-15}$, necessitating multiple levels of concatenation or large code distances.

15.5 The Discrete Logarithm Problem

Shor's original 1994 paper solved both the factoring problem and the discrete logarithm problem. The discrete logarithm problem is: given a cyclic group $G = \langle g \rangle$ of order $p$ and an element $h = g^x$, find $x$.

The quantum algorithm for discrete logarithms follows the same pattern:

1. Encode the problem as period finding: The function $f(a, b) = g^a h^{-b}$ has period $(x, 1)$ in $\mathbb{Z}_p \times \mathbb{Z}_p$, since $f(a+x, b+1) = g^{a+x} h^{-(b+1)} = g^a \cdot g^x \cdot h^{-b} \cdot h^{-1} = g^a h^{-b}$.

2. Use the quantum period-finding subroutine to find $(x, 1)$, which gives the discrete logarithm $x$.

This breaks Diffie-Hellman key exchange, elliptic curve Diffie-Hellman (ECDH), and the Digital Signature Algorithm (DSA) — essentially all widely deployed public-key cryptography except symmetric-key and post-quantum schemes.

Common Misconception: "Shor's algorithm only factors integers."

Shor's algorithm is more general. The period-finding subroutine can be applied to any function that can be efficiently computed by a quantum circuit and has a periodic structure. This includes factoring (period of $a^x \bmod N$), discrete logarithms (period of $g^a h^{-b}$), and Pell's equation. The core insight — that quantum computers can efficiently find periods of functions — is broadly applicable.

15.6 Qiskit Implementation: Factoring $N = 15$

We now implement Shor's algorithm to factor $15 = 3 \times 5$. We choose $a = 7$ (since $\gcd(7, 15) = 1$). The order of $7$ modulo $15$ is $r = 4$, because $7^4 = 2401 \equiv 1 \pmod{15}$.

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from math import gcd, ceil, log2
from fractions import Fraction

def continued_fraction_expansion(phi, max_denom):
    """Return convergents of phi with denominator <= max_denom."""
    convergents = []
    a_terms = []
    x = phi
    for _ in range(20):
        a = int(x)
        a_terms.append(a)
        frac = x - a
        if abs(frac) < 1e-10:
            break
        x = 1.0 / frac
    p0, q0 = 0, 1
    p1, q1 = 1, 0
    for a in a_terms:
        p = a * p1 + p0
        q = a * q1 + q0
        if q > max_denom:
            break
        convergents.append((p, q))
        p0, q0 = p1, q1
        p1, q1 = p, q
    return convergents

def modular_exponentiation_a2j(a, N):
    """Precompute a^(2^j) mod N for j = 0, 1, ..., 2n-1."""
    n = ceil(log2(N))
    t = 2 * n
    values = []
    val = a % N
    for j in range(t):
        values.append(val)
        val = (val * val) % N
    return values

def controlled_modular_multiply(circuit, control_qubit, target_register,
                                 constant, N, n):
    """
    Apply controlled multiplication of target_register by `constant` mod N.
    Uses a simple repeated-addition approach for small N.
    For N=15, we implement explicitly.
    """
    dim = 2**n
    U = np.zeros((dim, dim), dtype=complex)
    for y in range(dim):
        if y < N:
            new_y = (y * constant) % N
            U[new_y][y] = 1.0
        else:
            U[y][y] = 1.0
    from qiskit.circuit.library import UnitaryGate
    ugate = UnitaryGate(U, label=f"x{constant}")
    circuit.append(ugate.control(1), [control_qubit] + target_register[:])

def build_shor_circuit(N, a):
    n = ceil(log2(N))
    t = 2 * n

    ctrl = QuantumRegister(t, 'ctrl')
    targ = QuantumRegister(n, 'targ')
    meas = ClassicalRegister(t, 'meas')
    qc = QuantumCircuit(ctrl, targ, meas)

    for i in range(t):
        qc.h(ctrl[i])

    qc.x(targ[0])

    a_powers = modular_exponentiation_a2j(a, N)
    for j in range(t):
        controlled_modular_multiply(qc, ctrl[j], list(targ),
                                     a_powers[j], N, n)

    for i in range(t // 2):
        qc.swap(ctrl[i], ctrl[t - 1 - i])
    for i in range(t):
        qc.h(ctrl[i])
        for k in range(2, t - i + 1):
            angle = -2 * np.pi / (2 ** k)
            qc.cp(angle, ctrl[i + k - 1], ctrl[i])

    qc.measure(ctrl, meas)
    return qc

N = 15
a = 7
n = ceil(log2(N))
t = 2 * n

qc = build_shor_circuit(N, a)
print(f"Circuit depth: {qc.depth()}")
print(f"Qubit count: {qc.num_qubits}")
print(f"Gate count: {sum(qc.count_ops().values())}")

simulator = AerSimulator()
from qiskit import transpile
qc_compiled = transpile(qc, simulator)
job = simulator.run(qc_compiled, shots=8192)
result = job.result()
counts = result.get_counts()

print("\nTop measurement outcomes:")
for outcome, count in sorted(counts.items(),
                              key=lambda x: -x[1])[:8]:
    phase = int(outcome, 2) / (2**t)
    print(f"  |{outcome}⟩ (phase = {phase:.6f}), count = {count}")

found_factors = set()
for outcome, count in counts.items():
    if count < 50:
        continue
    c = int(outcome, 2)
    phase = c / (2**t)
    convergents = continued_fraction_expansion(phase, N)
    for p, q in convergents:
        if q > 0 and q < N:
            if pow(a, q, N) == 1:
                r = q
                if r % 2 == 0:
                    x = pow(a, r // 2, N)
                    if x != N - 1:
                        f1 = gcd(x + 1, N)
                        f2 = gcd(x - 1, N)
                        if 1 < f1 < N:
                            found_factors.add(f1)
                        if 1 < f2 < N:
                            found_factors.add(f2)

print(f"\nFactors found: {found_factors}")
print(f"Verification: {found_factors} == {{3, 5}}? "
      f"{found_factors == {3, 5}}")

Expected output: The circuit should yield factors {3, 5} with high probability. The measurement histogram will show peaks at phases $k/4$ for $k \in \{0, 1, 2, 3\}$, corresponding to the period $r = 4$.

15.6.1 Qiskit Implementation: Factoring $N = 21$

Let's also implement Shor's algorithm for $N = 21 = 3 \times 7$:

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_aer import AerSimulator
from math import gcd, ceil, log2
import numpy as np

def build_shor_circuit_v2(N, a):
    """Build Shor's period-finding circuit for general N and a."""
    n = ceil(log2(N))
    t = 2 * n

    ctrl = QuantumRegister(t, 'ctrl')
    targ = QuantumRegister(n, 'targ')
    meas = ClassicalRegister(t, 'meas')
    qc = QuantumCircuit(ctrl, targ, meas)

    for i in range(t):
        qc.h(ctrl[i])

    qc.x(targ[0])

    a_powers = modular_exponentiation_a2j(a, N)
    for j in range(t):
        if a_powers[j] != 1:
            controlled_modular_multiply(qc, ctrl[j], list(targ),
                                         a_powers[j], N, n)
        else:
            pass

    for i in range(t // 2):
        qc.swap(ctrl[i], ctrl[t - 1 - i])
    for i in range(t):
        qc.h(ctrl[i])
        for k in range(2, t - i + 1):
            angle = -2 * np.pi / (2 ** k)
            qc.cp(angle, ctrl[i + k - 1], ctrl[i])

    qc.measure(ctrl, meas)
    return qc

N = 21
a = 2
n = ceil(log2(N))
t = 2 * n

print(f"Factoring N = {N} with a = {a}")
print(f"Order of {a} mod {N}: checking...")

order = None
for r in range(1, N):
    if pow(a, r, N) == 1:
        order = r
        break
print(f"Order r = {order}")

if order % 2 == 0:
    x = pow(a, order // 2, N)
    if x != N - 1:
        f1, f2 = gcd(x - 1, N), gcd(x + 1, N)
        print(f"Factors: {f1} × {f2} = {N}")
    else:
        print("a^(r/2) ≡ -1 mod N, try different a")
else:
    print("r is odd, try different a")

qc = build_shor_circuit_v2(N, a)
print(f"\nCircuit: {qc.num_qubits()} qubits, depth {qc.depth()}")

simulator = AerSimulator()
qc_compiled = transpile(qc, simulator)
job = simulator.run(qc_compiled, shots=4096)
counts = job.result().get_counts()

print("\nTop measurement outcomes:")
for outcome, count in sorted(counts.items(), key=lambda x: -x[1])[:8]:
    phase = int(outcome, 2) / (2**t)
    print(f"  |{outcome}⟩ (phase = {phase:.6f}), count = {count}")

found_factors = set()
for outcome, count in counts.items():
    if count < 30:
        continue
    c = int(outcome, 2)
    phase = c / (2**t)
    convergents = continued_fraction_expansion(phase, N)
    for p, q in convergents:
        if 0 < q < N:
            if pow(a, q, N) == 1:
                r = q
                if r % 2 == 0:
                    x = pow(a, r // 2, N)
                    if x != N - 1:
                        f1 = gcd(x + 1, N)
                        f2 = gcd(x - 1, N)
                        if 1 < f1 < N:
                            found_factors.add(f1)
                        if 1 < f2 < N:
                            found_factors.add(f2)
                break

print(f"\nFactors found: {found_factors}")

15.6.2 A Note on Qiskit's Removed Built-in Implementation

Older tutorials (and older editions of the Qiskit textbook) show a built-in Shor class that handled the classical parts automatically:

# ⚠️ DOES NOT WORK on any currently supported Qiskit release.
# from qiskit.algorithms import Shor        # removed in Qiskit 0.25
# from qiskit.utils import QuantumInstance  # removed in Qiskit 1.0

Both qiskit.algorithms.Shor and qiskit.utils.QuantumInstance have been removed from Qiskit — the algorithm library was spun out and the Shor implementation was dropped entirely rather than ported. If you find code like this online, it predates Qiskit 1.0 and will raise ImportError.

Use the manual implementation from §15.6.1 instead. This is the better outcome pedagogically: the built-in class hid exactly the modular-exponentiation circuit construction and continued-fraction post-processing that make Shor's algorithm worth studying.

15.7 Experimental Demonstrations

15.7.1 Factoring 15

In 2001, Vandersypen et al. (IBM Almaden) used a 7-qubit NMR quantum computer to factor $N = 15$ using Shor's algorithm — the first experimental demonstration. They used a simplified circuit requiring only $2n = 8$ qubits for the control register, but clever compilation reduced the requirement to 7 physical qubits.

Key technical details: - Platform: NMR on a custom-synthesized molecule with 7 addressable spins - The 7-qubit molecule was a custom-synthesized perfluorobutadienyl iron complex (five $^{19}$F and two $^{13}$C spins), made for the experiment by IBM chemists Yannoni and Breyta - Pulse sequence: ~300 radiofrequency pulses - Total experiment time: several hours (NMR has long coherence but slow readout) - Result: correctly identified factors 3 and 5 with high confidence

15.7.2 Factoring 21

In 2012, Martín-López et al. (University of Bristol) factored $N = 21$ using an iterative photonic implementation, demonstrating the scalability of the qubit-recycling approach. More recently, superconducting qubit processors have factored 15 and 21 with higher fidelity.

The qubit-recycling technique allows period finding to be performed one qubit at a time, reducing the qubit requirement from $t + n$ to $n + 1$ — at the cost of requiring sequential measurements and feed-forward operations. This is essentially iterative phase estimation (Chapter 16) applied to Shor's algorithm.

15.7.3 The Scaling Challenge

Factoring 15 or 21 is trivial classically. The true test will be factoring a number beyond classical reach — likely requiring thousands of logical qubits. Current estimates suggest that factoring a 2048-bit RSA integer will require $\sim 20$ million physical qubits with surface-code error correction, placing it in the "long-term" category of quantum applications.

ASCII Diagram: The Resource Gap
================================

Physical qubits (log scale)
10^8 ┤ ┌────────────────────────── RSA-2048 (estimated)
    │ │
10^7 ┤ │
    │ │
10^6 ┤ │
    │ │
10^5 ┤ │
    │ │    ┌─────────────────────── FeMoco (150 logical qubits)
10^4 ┤ │    │
    │ │    │
10^3 ┤─│────│──── Current NISQ devices (IBM, Google, etc.)
    │ │    │
10^2 ┤ │    │
    │ │    │
  10 ┤─│────│──── Factored 15 (experimentally)
    │ │    │
   1 ┤─│────│──── Factored 21 (experimentally)
    └─┴────┴──────────────────────────────────────→ Year
    2001  2024  2030?  2040?  2050?

15.8 Implications for Cryptography

Shor's algorithm has profound implications:

1. RSA is broken in a post-quantum world. Any adversary with a large-scale quantum computer can decrypt RSA-ciphertexts and forge RSA signatures.

2. Elliptic curve cryptography (ECC) is also broken by a variant of Shor's algorithm that finds the discrete logarithm in polynomial time.

3. Symmetric cryptography (AES, SHA) is affected only quadratically: Grover's algorithm reduces AES-128 to $2^{64}$ effort, so doubling key sizes to AES-256 restores security.

4. Post-quantum cryptography (PQC) — lattice-based, code-based, hash-based, and multivariate schemes — is being standardized by NIST to replace RSA and ECC. In August 2024, NIST finalized the first PQC standards: ML-KEM (CRYSTALS-Kyber) for key encapsulation and ML-DSA (CRYSTALS-Dilithium) for digital signatures.

Recurring Theme: We're at the Beginning

The transition to post-quantum cryptography is one of the largest cryptographic migrations in history. Every TLS certificate, every SSH key, every digital signature infrastructure needs updating. This process is expected to take 10-15 years. Meanwhile, the "harvest now, decrypt later" threat means that data encrypted today with RSA or ECC could be decrypted by a future quantum computer. This is not a theoretical concern — intelligence agencies are known to be storing encrypted traffic precisely for this purpose.

The "harvest now, decrypt later" threat means that encrypted data intercepted today could be stored and decrypted once quantum computers mature — a concern for long-lived secrets.

NIST Post-Quantum Cryptography Standardization Timeline:

15.9 Shor's Algorithm and the Quantum Fourier Transform: A Deeper Look

15.9.1 Why the QFT Reveals Periodicity

The quantum Fourier transform $\text{QFT}_N$ maps a state with periodicity $r$ to a state concentrated at frequencies that are multiples of $N/r$. This is completely analogous to the classical discrete Fourier transform, which maps a periodic discrete signal to peaks at the harmonic frequencies.

Formally, if $|x\rangle$ appears with amplitude $\alpha_x$ in the computational basis, the QFT transforms:

$$\text{QFT}_N \sum_x \alpha_x |x\rangle = \sum_k \hat{\alpha}_k |k\rangle$$

where $\hat{\alpha}_k = \frac{1}{\sqrt{N}} \sum_x \alpha_x e^{2\pi i kx/N}$. For a periodic state $\alpha_x = \frac{1}{\sqrt{M}} \sum_j \delta_{x, jr + y}$, the Fourier coefficients $\hat{\alpha}_k$ are concentrated at $k$ such that $k$ is a multiple of $N/r$.

15.9.2 Comparison with Classical FFT

The classical FFT operates on a vector of $N = 2^n$ complex numbers and requires $O(N \log N) = O(2^n n)$ operations. The QFT, by contrast, operates on an $n$-qubit state and requires only $O(n^2)$ quantum gates — an exponential reduction in the number of operations. However, this comparison is subtle:

• The classical FFT produces all $N$ Fourier coefficients. The QFT produces a single random Fourier coefficient (upon measurement).
• To extract all Fourier coefficients classically would require $O(N)$ measurements, each collapsing the state.
• Shor's algorithm circumvents this by needing only one good Fourier coefficient to recover the period via continued fractions.

Common Misconception: "The QFT gives an exponential speedup over the classical FFT."

The QFT computes the Fourier transform of a quantum state in $O(n^2)$ gates, versus $O(2^n n)$ for the classical FFT. But these are different tasks: the classical FFT computes all Fourier coefficients, while the QFT produces one random coefficient. The speedup in Shor's algorithm comes from the entire algorithm, not just the QFT step. The quantum advantage is that the periodic state (produced by modular exponentiation) is prepared in superposition, and the QFT extracts the period from this superposition with high probability in a single measurement.

15.10 The Quantum Fourier Transform and Period Finding: A Detailed Analysis

15.10.1 The QFT as a Change of Basis

The quantum Fourier transform $\text{QFT}_N$ on $n = \log_2 N$ qubits is defined by:

$$\text{QFT}_N |j\rangle = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} e^{2\pi i jk/N} |k\rangle$$

This is exactly the discrete Fourier transform, but applied to a quantum state. The key difference from the classical FFT is that the QFT operates on the amplitudes of a quantum state in superposition, producing all $N$ Fourier components simultaneously — but only one can be read out per measurement.

The QFT circuit can be decomposed into $O(n^2)$ elementary gates:

ASCII Diagram: QFT Circuit (n qubits)
========================================

|j_{n-1}⟩ ──H──R₂──R₃──...──R_n──────────────────×─── → |k_{n-1}⟩
                  │    │         │                  │
|j_{n-2}⟩ ──────/─/──H──R₂──...──R_{n-1}──────────×─×─── → |k_{n-2}⟩
                             │         │            │ │
   ...                       /    ...  /            │ │
|j_1⟩     ─────────────────/────/───H──R₂──────────×─×─── → |k_1⟩
                                              │      │ │
|j_0⟩     ───────────────────────────────/──H────×─×─── → |k_0⟩

R_k = controlled phase rotation by angle 2π/2^k
× = SWAP gates (to reverse qubit order)

Gate count: n Hadamard + n(n-1)/2 controlled rotations + n/2 SWAPs = O(n²)

Each controlled-phase rotation $R_k$ has angle $2\pi / 2^k$. The total gate count is $O(n^2) = O((\log N)^2)$, compared to $O(N \log N) = O(2^n n)$ for the classical FFT — an exponential reduction.

Common Misconception: "The QFT gives an exponential speedup over the classical FFT."

This is only true in a specific sense. The QFT computes the Fourier transform of a quantum state in $O(n^2)$ gates, which is exponentially faster than the classical FFT's $O(2^n n)$. But the QFT output is a quantum state — measuring it yields only one Fourier component. To extract all components would require $O(2^n)$ measurements, negating the speedup. Shor's algorithm circumvents this by needing only one good component to recover the period via continued fractions.

15.10.2 Why Period Finding Works: The QFT of a Periodic State

The key insight is that a periodic superposition $\sum_j |jr + y\rangle$ has most of its amplitude concentrated at frequencies $k \cdot N/r$ in the Fourier basis. This is analogous to how a classical periodic signal has its power concentrated at harmonics.

For a state $|\psi\rangle = \frac{1}{\sqrt{m}} \sum_{j=0}^{m-1} |jr + y\rangle$ (a "comb" with teeth spaced by $r$), the QFT gives:

$$\text{QFT}|\psi\rangle = \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} e^{2\pi i ky/r} |k \cdot N/r\rangle$$

Measuring this state yields $k \cdot N/r$ with probability $1/r$ for each $k = 0, 1, \ldots, r-1$. From the measurement $c = k \cdot N/r$, we compute $c/N \approx k/r$ and use continued fractions to recover $r$.

15.10.3 Efficiency of Modular Exponentiation

The most computationally expensive part of Shor's algorithm is the modular exponentiation circuit. For an $n$-bit number $N$, we need to compute $a^x \bmod N$ for $x \in \{0, 1, \ldots, 2^{2n}-1\}$, which requires $2n$ controlled modular multiplications.

Classical precomputation: The values $a^{2^j} \bmod N$ for $j = 0, 1, \ldots, 2n-1$ can be computed classically in $O(n)$ steps via repeated squaring:

$$a^{2^0} \bmod N, \quad a^{2^1} = (a^{2^0})^2 \bmod N, \quad a^{2^2} = (a^{2^1})^2 \bmod N, \quad \ldots$$

Each squaring is a multiplication of two $n$-bit numbers modulo $N$, requiring $O(n^2)$ classical operations.

Quantum implementation: Each controlled modular multiplication $|x\rangle|y\rangle \to |x\rangle|y \cdot a^{2^j} \bmod N\rangle$ must be implemented as a quantum circuit. There are several approaches:

1. Permutation approach (used in our Qiskit code): For small $N$, construct the $2^n \times 2^n$ permutation matrix that maps $|y\rangle \to |y \cdot a^{2^j} \bmod N\rangle$. This requires $O(2^n)$ entries, feasible only for small $N$.

2. Reversible arithmetic: Implement modular multiplication using reversible adders, multipliers, and modular reduction circuits. This requires $O(n^2)$ gates per multiplication and $O(n)$ ancilla qubits.

3. Beauregard's approach: Uses a single multiplication circuit that is reused for all $2n$ controlled multiplications, reducing the total gate count to $O(n^3)$ with $O(n)$ ancilla qubits.

The overall complexity of Shor's algorithm:

Total: $O(n^3)$ gates, $O(n)$ qubits. For a 2048-bit number, this is roughly $10^{10}$ gates on $\sim 6000$ logical qubits.

15.10.4 Shor's Algorithm and the Hidden Subgroup Framework

Shor's algorithm is a special case of the hidden subgroup problem (HSP). Given a group $G$, a subgroup $H \leq G$, and a function $f: G \to S$ that is constant on cosets of $H$ and distinct on different cosets, find $H$.

• Factoring corresponds to the HSP over $G = \mathbb{Z}_N^*$ with $H = \langle a \rangle$, the cyclic subgroup generated by $a$.
• Discrete logarithm corresponds to the HSP over $G = \mathbb{Z}_p^* \times \mathbb{Z}_p^*$.
• Other instances: Pell's equation, principal ideal problem in number fields.

The quantum solution uses the QFT over $G$ (or a sufficiently large abelian subgroup) to extract information about $H$. For abelian groups, this always works efficiently. For non-abelian groups (relevant to graph isomorphism and certain lattice problems), efficient quantum algorithms are not known in general.

Recurring Theme: Quantum is Linear Algebra, Not Magic

The hidden subgroup framework makes clear that Shor's algorithm is not a magical factoring machine but a specific application of the quantum Fourier transform over a group structure. The QFT is a change of basis — from the "spatial" (computational) basis to the "frequency" (Fourier) basis — that reveals the periodic structure of the function $f$. This is exactly analogous to how the classical Fourier transform reveals the frequency content of a signal. The quantum advantage comes from the fact that the QFT operates on exponentially many amplitudes simultaneously.

15.10.5 Complete Probability Analysis of the Order-Finding Measurement

Let us derive the exact probability distribution of the measurement outcome in the order-finding subroutine. This is one of the most important calculations in quantum computing, as it determines the success probability of Shor's algorithm.

After modular exponentiation, the state of the control register (tracing over the target register) is:

$$\rho = \frac{1}{r} \sum_{y=0}^{r-1} |\psi_y\rangle\langle\psi_y|$$

where $|\psi_y\rangle = \frac{1}{\sqrt{m}} \sum_{j=0}^{m-1} |jr + y\rangle$ with $m = \lfloor 2^t / r \rfloor$.

Applying the inverse QFT and measuring gives outcome $c$ with probability:

$$P(c) = \frac{1}{r} \sum_{y=0}^{r-1} |\langle c | \text{QFT}^\dagger | \psi_y \rangle|^2$$

Using the derivation from Section 15.2.2, this evaluates to:

$$P(c) = \frac{1}{2^t r} \left| \frac{\sin(\pi c r / 2^t)}{\sin(\pi c / 2^t)} \right|^2 \cdot \left| \frac{\sin(\pi m r \phi)}{\sin(\pi r \phi)} \right|^2$$

where $\phi = c/2^t - \lfloor c/2^t \cdot r \rfloor / r$ is the deviation from the nearest fraction with denominator $r$.

For the special case where $r | 2^t$ (i.e., $r$ divides $2^t$ exactly), the probability simplifies dramatically:

$$P(c = k \cdot 2^t / r) = \frac{1}{r}$$

for each $k = 0, 1, \ldots, r-1$, and $P(c) = 0$ for all other $c$. This means the measurement always yields a value of the form $k \cdot 2^t / r$, exactly revealing the period.

For the general case where $r \nmid 2^t$, the probability peaks at $c \approx k \cdot 2^t / r$ with width $\sim 2^t / r^2$. The probability of measuring a value close enough to recover $r$ via continued fractions is at least $4/\pi^2$.

Worked Example: Probability Distribution for N=15, a=7, t=8

$N = 15$, $a = 7$, $r = 4$, $t = 8$, $2^t = 256$.

Since $r = 4$ divides $256 = 2^8$, the measurement always yields one of $c \in \{0, 64, 128, 192\}$ (corresponding to $k/4$ for $k = 0, 1, 2, 3$).

• $c = 0$: $c/256 = 0/4 = 0$. Convergent gives $r = 1$... Actually, $k = 0$ gives the trivial phase 0, which doesn't reveal $r$. This happens with probability $1/4 = 25\%$.
• $c = 64$: $c/256 = 1/4$. Continued fractions immediately give $r = 4$. ✓ (probability 25%)
• $c = 128$: $c/256 = 1/2$. Continued fractions give convergent $1/2$, but $7^2 = 4 \not\equiv 1$. Try next convergent: $r = 4$ works since $1/2 = 2/4$ reduces to $1/2$ but $k/r = 2/4$ gives $r = 4$ after checking $7^4 \equiv 1$. ✓ (probability 25%)
• $c = 192$: $c/256 = 3/4$. Convergents give $r = 4$. ✓ (probability 25%)

Total success probability: at least $3/4 = 75\%$ (the $c = 0$ outcome gives no information). If $c = 0$, we simply repeat the algorithm with a different random $a$.

15.11 Shor's Algorithm and Post-Quantum Cryptography: A Deeper Look

15.11.1 The Transition from RSA to Post-Quantum Cryptography

The discovery of Shor's algorithm has triggered a massive effort to develop and standardize post-quantum cryptographic schemes. NIST's Post-Quantum Cryptography Standardization Process began in 2016 and concluded in 2024 with three final standards:

FIPS 203: ML-KEM (Module-Lattice-Based Key Encapsulation Mechanism) - Based on the hardness of the Module-LWE (Learning With Errors) problem - Key sizes: ~1 KB (public key), ~1 KB (ciphertext) - Security level: 128, 192, or 256 bits - Primary replacement for RSA key exchange

FIPS 204: ML-DSA (Module-Lattice-Based Digital Signature Algorithm) - Based on Module-LWE and Module-SIS (Short Integer Solution) - Signature sizes: ~2.5 KB - Primary replacement for RSA/ECDSA signatures

FIPS 205: SLH-DSA (Stateless Hash-Based Digital Signature Algorithm) - Based solely on hash functions (conservative security assumption) - Signature sizes: ~7-50 KB (much larger than ML-DSA) - Backup option with minimal security assumptions

15.11.2 Lattice-Based Cryptography and the LWE Problem

The security of ML-KEM and ML-DSA is based on the Learning With Errors (LWE) problem, which is believed to be hard for both classical and quantum computers:

LWE Problem: Given a matrix $A \in \mathbb{Z}_q^{m \times n}$ and a vector $b = As + e \bmod q$ where $s$ is a secret vector and $e$ is a small error vector, find $s$.

The best known quantum algorithm for LWE runs in time $2^{O(n/\log n)}$, which is sub-exponential but still infeasible for appropriate parameter choices (e.g., $n = 1024$).

Recurring Theme: We're at the Beginning

The transition to post-quantum cryptography is one of the largest cryptographic migrations in history. Every TLS certificate, every SSH key, every digital signature infrastructure needs updating. This process is expected to take 10-15 years. Meanwhile, the "harvest now, decrypt later" threat means that data encrypted today with RSA or ECC could be decrypted by a future quantum computer. This is not a theoretical concern — intelligence agencies are known to be storing encrypted traffic precisely for this purpose.

15.11.3 The "Harvest Now, Decrypt Later" Threat

The "harvest now, decrypt later" attack model works as follows:

1. Today: An adversary records encrypted traffic (e.g., HTTPS sessions, VPN connections, encrypted emails).
2. Years later: When a large-scale quantum computer becomes available, the adversary uses Shor's algorithm to break the RSA/ECC keys and decrypt the stored traffic.

This attack is practical because: - Data with long-term sensitivity (medical records, trade secrets, classified information) remains valuable for decades. - RSA-2048 keys used today will be vulnerable to anyone who records the ciphertexts. - The cost of storage is negligible compared to the value of the decrypted data.

Mitigation: Organizations handling sensitive long-term data should begin migrating to hybrid encryption (combining RSA/ECC with post-quantum algorithms) immediately, even before post-quantum schemes are fully standardized. Hybrid encryption provides security against both classical and quantum adversaries.

15.12 Alternative Approaches to Factoring on Quantum Computers

While Shor's algorithm is the most well-known quantum factoring algorithm, several alternatives have been proposed, each with different resource requirements and tradeoffs.

15.12.1 Variational Quantum Factoring

Variational quantum factoring (VQF) reformulates factoring as an optimization problem. Given $N = pq$, the binary representations of $p$ and $q$ satisfy:

$$N = \sum_{i,j} 2^{i+j} p_i q_j$$

where $p_i, q_j \in \{0, 1\}$ are the binary digits of $p$ and $q$. This is a quadratic unconstrained binary optimization (QUBO) problem that can be encoded as a cost function on a quantum computer and solved using variational methods.

Advantages: VQF uses shallow circuits (depth $O(1)$) compared to Shor's depth $O(n^3)$, making it potentially runnable on NISQ devices.

Disadvantages: VQF is not guaranteed to find the correct factors, and there is no proven quantum advantage. The optimization landscape can have many local minima, and classical heuristics (simulated annealing, branch-and-bound) can solve the same QUBO efficiently for small instances.

15.12.2 Quantum Annealing for Factoring

Quantum annealing (implemented on D-Wave hardware) can also be used for factoring by encoding the multiplication constraint as an Ising Hamiltonian:

$$H_{\text{Ising}} = \sum_i h_i \sigma_i^z + \sum_{i

The ground state of this Hamiltonian encodes the factors $p$ and $q$. However, the number of couplings grows as $O(n^2)$ for an $n$-bit number, and the embedding on D-Wave's sparse connectivity graph requires additional qubits.

Results: D-Wave has factored numbers up to $N = 1{,}099{,}551{,}473{,}989$ (a 40-bit RSA number), but using a hybrid quantum-classical approach where the quantum annealer only solves a small subproblem. The quantum advantage is not established.

15.12.3 Hybrid Classical-Quantum Approaches

Several hybrid approaches combine classical and quantum computation:

1. Quantum-enhanced sieving: Use a quantum algorithm to accelerate the sieving step of the number field sieve. This gives a polynomial speedup over the best classical factoring algorithm, but still sub-exponential.

2. Quantum random walk factoring: Use quantum walks to search for smooth numbers (a key step in the quadratic sieve). This gives a quadratic speedup over the classical search.

3. Grover-accelerated factoring: Apply Grover's search to the trial division method. This gives a square-root speedup ($O(\sqrt{N})$ instead of $O(N)$), which is still exponential and far worse than Shor's polynomial-time algorithm.

Common Misconception: "Grover's algorithm can break RSA."

Grover's algorithm provides a quadratic speedup for brute-force search, reducing the time to break RSA from $O(N)$ to $O(\sqrt{N})$. For a 2048-bit modulus, this is still $O(2^{1024})$ — utterly infeasible. Shor's algorithm provides an exponential speedup ($O((\log N)^3)$), which is why it breaks RSA. Grover's algorithm is relevant for symmetric-key cryptography (reducing AES-128 to $2^{64}$), not for public-key cryptography.

15.13 The Quantum Threat Timeline

When will quantum computers actually be able to break RSA-2048? The answer depends on several factors:

Hardware milestones needed: 1. Error-corrected logical qubits: $\sim 6{,}000$ logical qubits, each requiring $\sim 1{,}000$ physical qubits. 2. Total physical qubits: $\sim 20$ million (Gidney & Ekerå, 2021). 3. Gate fidelity: Logical error rate $< 10^{-12}$ per operation. 4. Coherence time: Long enough to execute $\sim 10^{10}$ logical operations. 5. Runtime: Estimated 8 hours for RSA-2048.

Current state (2024-2025): - Largest superconducting processors: $\sim 1{,}000$ physical qubits - Best two-qubit gate fidelity: $\sim 99.9\%$ - Longest coherence time: $\sim 1$ ms (superconducting), $\sim 10$ s (trapped ions) - No demonstrated error-corrected logical qubit with $< 10^{-6}$ logical error rate

Expert estimates for RSA-2048 factoring: - Optimistic: 2035-2040 (accelerated by new algorithms and hardware breakthroughs) - Moderate: 2040-2050 (assuming steady progress along current roadmaps) - Conservative: 2050+ (accounting for unforeseen engineering challenges)

Key uncertainty: The timeline depends heavily on whether error correction overhead can be reduced below the current $\sim 1{,}000$:1 ratio. Algorithmic improvements (e.g., better factoring circuits, improved surface codes, magic state distillation optimizations) could reduce the overhead by 10-100x, accelerating the timeline.

Recurring Theme: We're at the Beginning

Current quantum computers cannot break any cryptographically relevant code. The gap between current hardware ($\sim 10^3$ physical qubits with $\sim 99.9\%$ gate fidelity) and the hardware needed for RSA-2048 ($\sim 10^7$ physical qubits with $\sim 99.9999\%$ logical gate fidelity) is 4 orders of magnitude in qubit count and 3 orders of magnitude in gate fidelity. Progress has been steady but slow, with qubit counts doubling roughly every 2 years. At this rate, reaching $10^7$ qubits would take $\sim 25$ years — around 2050. However, breakthroughs in error correction, qubit quality, or architecture could accelerate this timeline significantly.