Lecture 15
Order Finding and Shor's Algorithm I

Week 8 | factoring reduces to finding the order of a number mod N, and order finding is phase estimation; Lab 4 launches | reading: Watrous Ch. 7 (order finding and factoring); de Wolf Ch. 5

the powers of 7 mod 15 cycle 1, 7, 4, 13, 1, ...: the length of that cycle is enough to factor 15

In this lecture

Why factoring matters
The number theory we need
From the order to the factors
Order finding is phase estimation
Reading s/r: continued fractions
Shor's algorithm, end to end
Lab 4 launch
Try it in code
Summary and exercises

Lecture 14 ended with a promise: phase estimation applied to multiplication mod $N$ finds the order of a number, and the order factors $N$. This lecture keeps it. By the end you should be able to explain why factoring breaks RSA, compute orders and use them to factor by hand, prove the reduction from factoring to order finding, explain why phase estimation on the easy state $|1\rangle$ returns a random fraction $s/r$, recover $r$ from a measured $y$ with continued fractions, and carry out Shor's algorithm end to end for $N=15$ and $N=35$. Lecture 16 builds the modular-arithmetic circuits and extends the method to discrete logarithms.

1 — Why factoring matters

crypto connection — RSA in one paragraph

Choose two large primes $p,q$ and publish $N=pq$ and an exponent $e$ coprime to $(p-1)(q-1)$. The private key is $d$ with $ed\equiv1\pmod{(p-1)(q-1)}$. Encryption is $m\mapsto m^e\bmod N$ and decryption $c\mapsto c^d\bmod N$; signatures run the same maps the other way. Anyone who factors $N$ knows $(p-1)(q-1)$, computes $d$ with Euclid's algorithm, and can read every message and forge every signature. The best classical algorithm, the general number field sieve, takes time $\exp\big((1.92+o(1))(\ln N)^{1/3}(\ln\ln N)^{2/3}\big)$: subexponential, but still far out of reach for 2048-bit $N$. The largest RSA challenge number factored so far, RSA-250 (829 bits, in 2020), took about 2,700 CPU-core-years. Shor's algorithm factors in time polynomial in $\log N$.

Shor's algorithm is mostly classical. The quantum computer does one job, finding the order of a number mod $N$, and everything before and after it is number theory that runs on a laptop.

2 — The number theory we need

Greatest common divisors. Euclid's algorithm computes $\gcd(x,N)$ by repeated division with remainder: $\gcd(x,N)=\gcd(N\bmod x,\,x)$. It needs $O(\log N)$ steps, so it is fast even for 2048-bit numbers. A random $x$ shares a factor with $N=pq$ only with probability about $1/p+1/q$, which is negligible; finding such an $x$ is exactly as hard as factoring.

Units. The numbers $a\in\{1,\dots,N-1\}$ with $\gcd(a,N)=1$ form the group $\mathbb Z_N^*$ under multiplication mod $N$: products of units are units, and every unit has an inverse (from Euclid's algorithm). For $N=15$, $\mathbb Z_{15}^*=\{1,2,4,7,8,11,13,14\}$, eight elements.

definition 15.1 — order

The order of $a\in\mathbb Z_N^*$ is the least $r\ge1$ with $a^r\equiv1\pmod N$. The powers $1,a,a^2,\dots$ then repeat with period exactly $r$: $a^x\equiv a^{x'}$ if and only if $x\equiv x'\pmod r$.

What is the order of 4 modulo 51?

List the powers $4,4^2,4^3,\dots$ reducing mod 51 each time, and stop at the first one that is 1.

$4,\ 16,\ 64\equiv13,\ 13\cdot4=52\equiv1$. So $r=4$.

example 15.1 — the order of 7 mod 15

Find the order of 7 modulo 15.

$7^1=7$, $7^2=49\equiv4$, $7^3\equiv4\cdot7=28\equiv13$, $7^4\equiv13\cdot7=91\equiv1$. So $r=4$: the powers cycle $1,7,4,13,1,7,\dots$, the loop in the cartoon.

The order always divides $|\mathbb Z_N^*|$ (Lagrange's theorem), which is $(p-1)(q-1)$ for $N=pq$. But for a 2048-bit $N$ the order is typically a number of about 2048 bits too, and no known classical method finds it much faster than factoring $N$. Computing a single power $a^x\bmod N$ is easy, though: repeated squaring computes $a,a^2,a^4,a^8,\dots$ by squaring mod $N$ and multiplies together the ones selected by the bits of $x$, using $O(\log x)$ multiplications of $n$-bit numbers. That is $O(n^3)$ bit operations for an $n$-bit exponent with schoolbook multiplication, and the quantum circuit in Lecture 16 does the same thing reversibly.

3 — From the order to the factors

theorem 15.1 — an even order gives a factor

Let $a\in\mathbb Z_N^*$ have even order $r$, and suppose $a^{r/2}\not\equiv-1\pmod N$. Then $\gcd(a^{r/2}-1,\,N)$ and $\gcd(a^{r/2}+1,\,N)$ are nontrivial factors of $N$.

Let $h=a^{r/2}\bmod N$. Then $h^2=a^r\equiv1$, so $N$ divides $h^2-1=(h-1)(h+1)$. Now $N\nmid h-1$, because $h\equiv1$ would mean $a^{r/2}\equiv1$ with $r/2\lt r$, contradicting minimality of $r$; and $N\nmid h+1$ by assumption. So $N$ divides the product $(h-1)(h+1)$ but neither factor. If $\gcd(h-1,N)$ were 1, then $N\mid(h-1)(h+1)$ would force $N\mid h+1$, which is false; and it is not $N$, since $N\nmid h-1$. So $\gcd(h-1,N)$ is a proper factor, and likewise $\gcd(h+1,N)$.

$a\in\mathbb Z_N^*$, $N$ odd, has even order $r$, and $h=a^{r/2}\bmod N$. Which of these can never happen?

$h\equiv-1\pmod N$This does happen: $a=14$ mod 15 has $r=2$ and $h=14\equiv-1$ (Example 15.2). It is one of the two ways the method fails.
$h\equiv1\pmod N$
$\gcd(h-1,N)=1$This happens exactly when $h\equiv-1$: then $h-1\equiv-2$, and $\gcd(N-2,N)=1$ for odd $N$.
$h$ is neither $1$ nor $-1$ mod $N$That is the good case of Theorem 15.1, for example $h=4$ for $a=7$, $N=15$.

Look at the first line of the proof: what does $r$ being the least exponent with $a^r\equiv1$ say about $a^{r/2}$?

$h\equiv1$ would mean $a^{r/2}\equiv1$ with $r/2\lt r$, contradicting the minimality of $r$. The only real failure for even $r$ is $h\equiv-1$.

In words: $h$ is a square root of 1 mod $N$ other than $\pm1$. For a prime modulus only $\pm1$ exist, but for $N=pq$ there are four square roots of 1, and the "extra" ones split $N$.

For $N=91$, some unit $a$ gives $h=a^{r/2}\equiv64$. Check that $64^2\equiv1\pmod{91}$, then compute the two gcds of Theorem 15.1.

$64^2=4096=45\cdot91+1$. Now factor $h-1=63$ and $h+1=65$ and see what each shares with $91$.

$\gcd(63,91)=7$ and $\gcd(65,91)=13$, and indeed $91=7\times13$. The square root 64 of 1 is neither $1$ nor $-1\equiv90$, so it splits $N$.

example 15.2 — every unit mod 15

For each $a\in\mathbb Z_{15}^*$, find its order and decide whether Theorem 15.1 factors 15.

$a$	1	2	4	7	8	11	13	14
order $r$	1	4	2	4	4	2	4	2
$a^{r/2}\bmod15$	—	4	4	4	4	11	4	14
result	$r$ odd	3, 5	3, 5	3, 5	3, 5	5, 3	3, 5	$-1$: fails

For example $a=7$: $r=4$, $h=7^2\bmod15=4$, and $\gcd(3,15)=3$, $\gcd(5,15)=5$. For $a=11$: $h=11$, $\gcd(10,15)=5$, $\gcd(12,15)=3$. Six of the eight units work; $a=1$ has odd order and $a=14\equiv-1$ gives $h\equiv-1$.

example 15.3 — N = 35, a = 2

Use the order of 2 modulo 35 to factor 35.

The powers of 2 mod 35 are $2,4,8,16,32,29,23,11,22,9,18,1$, so $r=12$. Then $h=2^6=64\equiv29\not\equiv-1$, and $\gcd(28,35)=7$, $\gcd(30,35)=5$. Notice that we needed twelve powers to find $r$ by brute force; for an $n$-bit $N$ the order can be exponentially large in $n$, and that is the step the quantum computer replaces.

theorem 15.2 — a random a works at least half the time

Let $N$ be odd with at least two distinct prime factors. If $a$ is chosen uniformly from $\mathbb Z_N^*$, then $r$ is even and $a^{r/2}\not\equiv-1$ with probability at least $\tfrac12$. (More precisely, at least $1-2^{1-k}$ when $N$ has $k$ distinct prime factors.)

Proof sketch. By the Chinese remainder theorem, choosing $a$ mod $N=p^\alpha q^\beta\cdots$ is the same as choosing independent random units mod each prime power. Let $2^{d_i}$ be the largest power of 2 dividing the order of $a$ mod the $i$-th prime power. The method fails exactly when all the $d_i$ are equal (which includes "$r$ odd", all $d_i=0$), and each $d_i$ takes any particular value with probability at most $\tfrac12$, because the group mod an odd prime power is cyclic of even order. See Nielsen–Chuang, Appendix 4, for the details.

The bound is tight: for $N=33$ and $N=77$ exactly half of the units work. For other small $N$ it is better (computed in verify_L15.py; the code box in Section 8 checks $N=35$):

$N$	15	33	35	39	51	55	77	91
units that work	6 / 8	10 / 20	18 / 24	18 / 24	30 / 32	30 / 40	30 / 60	54 / 72

So a factoring algorithm needs only an order-finding subroutine: check that $N$ is odd and not a prime power (both easy classically), pick random $a$, take $\gcd(a,N)$ in case we are lucky, find the order, and apply Theorem 15.1; expect to repeat about twice.

$N$ is odd and has exactly three distinct prime factors. What lower bound does Theorem 15.2 give for the probability that a random unit $a$ leads to a factor?

Use the precise form of the bound, $1-2^{1-k}$, with $k$ the number of distinct primes.

$1-2^{1-3}=1-\tfrac14=\tfrac34$. More prime factors make the bad event (all the $d_i$ equal) less likely.

4 — Order finding is phase estimation

Let $n=\lceil\log_2N\rceil$, so the numbers $0,\dots,N-1$ fit in an $n$-qubit register.

definition 15.2 — modular multiplication

For $a\in\mathbb Z_N^*$, $U_a$ is the linear map on $n$ qubits with

$$U_a|y\rangle=\begin{cases}|ay\bmod N\rangle,&0\le y\lt N,\\[2pt] |y\rangle,&N\le y\lt2^n.\end{cases}$$

proposition 15.3 — U_a is unitary, and its powers are cheap

$U_a$ permutes the basis states, so it is unitary. Moreover $U_a^k=U_{a^k\bmod N}$; in particular $U_a^{2^j}=U_c$ with $c=a^{2^j}\bmod N$, a single multiplication by a number computed classically.

Multiplication by a unit is a bijection of $\{0,\dots,N-1\}$: its inverse is multiplication by $a^{-1}\bmod N$. Extended by the identity on $y\ge N$, it permutes all $2^n$ basis states, and a permutation matrix is unitary. Applying $U_a$ $k$ times multiplies by $a$ $k$ times, which is multiplying once by $a^k\bmod N$.

The orbit of 1 under $U_a$ is the cycle $1,a,a^2,\dots,a^{r-1}$ of length $r$: $U_a$ acts on it as a cyclic shift. Lecture 13 showed that shifts are diagonalised by Fourier vectors, and that is where the eigenvectors come from.

theorem 15.4 — the spectrum of U_a on the orbit of 1

Let $a$ have order $r$. For $s=0,1,\dots,r-1$ the states

$$|u_s\rangle=\frac{1}{\sqrt r}\sum_{k=0}^{r-1}e^{-2\pi i\,sk/r}\,|a^k\bmod N\rangle$$

are orthonormal eigenvectors of $U_a$ with $U_a|u_s\rangle=e^{2\pi i\,s/r}|u_s\rangle$, and

$$|1\rangle=\frac{1}{\sqrt r}\sum_{s=0}^{r-1}|u_s\rangle .$$

Proof. HW8. (verify_L15.py checks both statements numerically for $N=15$, $a=7$ and $N=35$, $a=2$.)

Take $N=51$, $a=4$, whose order you found above. How many eigenvectors $|u_s\rangle$ does Theorem 15.4 give, and what is the eigenphase $\theta\in[0,1)$ of $|u_3\rangle$?

There is one $|u_s\rangle$ for each $s=0,\dots,r-1$, and $U_a|u_s\rangle=e^{2\pi i\,s/r}|u_s\rangle$.

$r=4$, so there are four, $|u_0\rangle,\dots,|u_3\rangle$, and $|u_3\rangle$ has eigenphase $s/r=\tfrac34$, eigenvalue $e^{3\pi i/2}=-i$.

idea — we never prepare an eigenvector

Preparing $|u_s\rangle$ would require knowing $r$. But the state $|1\rangle$, the number one in binary, is an equal superposition of all $r$ eigenvectors, and Proposition 14.6 says phase estimation on a superposition simply samples one of the eigenphases. Each eigenphase is a fraction $s/r$ with the unknown $r$ as denominator.

corollary 15.5 — what phase estimation returns

Run phase estimation of $U_a$ with $m$ control qubits on the work register $|1\rangle$. The outcome $y$ is distributed as follows: $s$ is uniform on $\{0,\dots,r-1\}$, and given $s$, $y$ has the distribution of Theorem 14.3 for $\theta=s/r$. In particular $y/2^m$ is the nearest $m$-bit fraction to $s/r$ with probability at least $4/\pi^2$.

Each run of phase estimation of $U_a$ on $|1\rangle$ returns $y$ with $y/2^m\approx s/r$, where $s$ is .

Theorem 15.4 writes $|1\rangle$ as $\tfrac1{\sqrt r}\sum_s|u_s\rangle$. What does Proposition 14.6 do with a superposition of eigenvectors?

Every $|u_s\rangle$ has weight $|1/\sqrt r|^2=1/r$, so phase estimation samples $s$ uniformly. That is why we only learn $r$ as the denominator of a random fraction.

Order finding. The controlled powers are controlled multiplications by the classically precomputed numbers $a^{2^j}\bmod N$ (Proposition 15.3). The quantum part ends at the measurement; everything to the right is classical post-processing (Section 5).

There is an equivalent picture that connects to Lecture 13. The controlled multiplications compute $\sum_x|x\rangle|a^x\bmod N\rangle$. If we measured the work register and saw $a^{x_0}$, the control register would collapse to the periodic state $\sum_j|x_0+jr\rangle$ (Definition 15.1), and the inverse QFT of a periodic state peaks near multiples of $2^m/r$ (Proposition 13.5). Measuring the work register is unnecessary, but it explains why the output reveals the period.

5 — Reading s/r: continued fractions

We observe $y$, and $y/2^m$ is very close to $s/r$ with $r\lt N$. Many fractions are close to $y/2^m$; we want the one with a small denominator. Continued fractions find it.

definition 15.3 — continued fractions and convergents

Every rational $x\ge0$ can be written as

$$x=c_0+\cfrac{1}{c_1+\cfrac{1}{c_2+\cfrac{1}{\ddots+\cfrac{1}{c_k}}}}=[c_0;c_1,\dots,c_k]$$

with integers $c_0\ge0$ and $c_1,\dots,c_k\ge1$, found by Euclid's algorithm: $c_0=\lfloor x\rfloor$, then repeat on $1/(x-c_0)$. Truncating after $c_j$ gives the convergent $p_j/q_j=[c_0;c_1,\dots,c_j]$, computed by $p_j=c_jp_{j-1}+p_{j-2}$, $q_j=c_jq_{j-1}+q_{j-2}$, starting from $p_{-1}=1,q_{-1}=0,p_{-2}=0,q_{-2}=1$.

Expand $\tfrac{7}{16}=[0;c_1,c_2,c_3]$ as a continued fraction.

Run Euclid's algorithm: $16=c_1\cdot7+\text{remainder}$, then divide 7 by that remainder, and so on.

$16=2\cdot7+2$, $7=3\cdot2+1$, $2=2\cdot1$. So $\tfrac7{16}=[0;2,3,2]$, with convergents $\tfrac12,\ \tfrac37,\ \tfrac7{16}$.

theorem 15.6 — good approximations are convergents (Legendre)

If $s/r$ is in lowest terms and $\big|x-\tfrac sr\big|\lt\tfrac{1}{2r^2}$, then $s/r$ is one of the convergents of $x$.

Proof. See Hardy and Wright, An Introduction to the Theory of Numbers, Theorem 184.

This dictates the number of control qubits. If $2^m\ge N^2$ and $y$ is the nearest $m$-bit fraction to $s/r$, then $\big|\tfrac{y}{2^m}-\tfrac sr\big|\le\tfrac{1}{2^{m+1}}\le\tfrac{1}{2N^2}\lt\tfrac{1}{2r^2}$, so $s/r$, in lowest terms, is a convergent of $y/2^m$. So the standard choice is $m\approx2n$ control qubits: $m=8$ for $N=15$ ($256\ge225$) and $m=11$ for $N=35$ ($2048\ge1225$). Among the convergents, take the last one whose denominator is below $N$, and check the candidate $q$ by computing $a^q\bmod N$ classically. (No later convergent can also have a denominator below $N$: two different fractions with denominators below $N$ differ by more than $1/N^2$, but both would lie within $1/(2N^2)$ of $y/2^m$.)

What is the smallest number of control qubits $m$ with $2^m\ge N^2$ for $N=91$?

$91^2=8281$. Compare with nearby powers of 2: $2^{13}=8192$.

$2^{13}=8192$ falls just short of 8281, so $m=14$ ($2^{14}=16384$). Here $n=7$, and $m=2n$ exactly.

example 15.4 — N = 15, a = 7, m = 8

Find the output distribution, and say what each outcome tells us.

Here $r=4$ divides $2^8=256$, so each $s/r=s/4$ is an exact 8-bit fraction, and by Theorem 14.2 and Corollary 15.5 the outcome is $y=64s$ for $s=0,1,2,3$, each with probability $\tfrac14$. Reading them: $64/256=1/4$ gives the candidate 4, and $7^4\equiv1$: correct. $192/256=3/4$ also gives 4. $128/256=1/2$ gives the candidate 2, but $7^2\equiv4\not\equiv1$, so it is rejected: $s=2$ shares a factor with $r=4$, and the fraction reduced. $y=0$ says nothing. So each run succeeds with probability exactly $\tfrac12$.

Example 15.4, computed exactly. Four equally likely outcomes; two of them reveal $r=4$.

example 15.5 — N = 35, a = 2, m = 11, outcome y = 853

A run of order finding for $N=35$, $a=2$ returns $y=853$. Recover the order.

Run Euclid's algorithm on $853/2048$:

$$\begin{aligned}2048&=2\cdot853+342, &853&=2\cdot342+169,\\[4pt] 342&=2\cdot169+4, &169&=42\cdot4+1,\\[4pt] 4&=4\cdot1. &&\end{aligned}$$

so $853/2048=[0;2,2,2,42,4]$, with convergents

$$\frac01,\ \frac12,\ \frac25,\ \frac5{12},\ \frac{212}{509},\ \frac{853}{2048}.$$

The last denominator below 35 is 12, and $2^{12}\bmod35=1$, so $r=12$, as in Example 15.3. The outcome came from $s=5$: $5/12\cdot2048=853.33$.

Order finding for $N=35$, $a=2$, $r=12$, $m=11$, computed from Corollary 15.5. Since 12 does not divide 2048, most peaks are split between two neighbouring outcomes. The four green peaks ($s=1,5,7,11$) give $r=12$ directly; the others give a proper divisor of 12, such as 6 from $s=2$ or 4 from $s=3$.

What can go wrong. Three things, all detectable classically by checking $a^q\equiv1$:

$s=0$ gives no information.
$\gcd(s,r)\gt1$: the convergent is $s/r$ in lowest terms, whose denominator is a proper divisor of $r$. This happens with probability $1-\varphi(r)/r$, where $\varphi(r)$ counts the $s\lt r$ coprime to $r$. Remedy: two runs giving denominators $q_1,q_2$ usually have $\operatorname{lcm}(q_1,q_2)=r$.
The outcome is not the nearest one (probability at most about $1-4/\pi^2$ for each $s$, usually much less).

Suppose $r=8$ and the run happened to sample $s=6$, with the nearest outcome observed. Continued fractions return the candidate order .

Continued fractions find $s/r$ in lowest terms. Reduce $6/8$.

$6/8=3/4$, so the method sees the denominator 4, a proper divisor of 8. The check $a^4\equiv1$ fails and we run again (or take an lcm with another candidate).

proposition 15.7 — one run succeeds with reasonable probability

With $2^m\ge N^2$, one run of order finding followed by continued fractions returns $r$ itself with probability at least $\tfrac{4}{\pi^2}\cdot\tfrac{\varphi(r)}{r}$. Since $\varphi(r)/r$ is at least of order $1/\log\log r$, $O(\log\log N)$ runs suffice on average; in practice a handful.

Proof sketch. $s$ is uniform, so $\gcd(s,r)=1$ with probability $\varphi(r)/r$. Given such an $s$, the nearest outcome has probability at least $4/\pi^2$ (Theorem 14.4), and then Theorem 15.6 makes $s/r$ a convergent with denominator exactly $r\lt N$. For the bound on $\varphi(r)/r$, see Hardy–Wright, Theorem 328.

For $N=35$, $a=2$ the exact probability that one run returns $r=12$ is 0.321 (verify_L15.py), slightly below $\varphi(12)/12=\tfrac13$ because a little probability leaks to outcomes beyond the two nearest.

6 — Shor's algorithm, end to end

algorithm — Shor's factoring algorithm (Shor, 1994)

Input: an odd $N$ with at least two distinct prime factors (check classically that $N$ is not a prime power).

Pick $a$ uniformly from $\{2,\dots,N-1\}$. If $g=\gcd(a,N)\gt1$, output $g$ (lucky).
Quantum: run order finding for $a$ with $m\approx2\log_2N$ control qubits on the work register $|1\rangle$; measure $y$.
Find the last convergent $p/q$ of $y/2^m$ with $q\lt N$. If $a^q\not\equiv1\pmod N$, go to step 2 (or combine with an earlier candidate by an lcm).
Set $r=q$. If $r$ is odd or $a^{r/2}\equiv-1\pmod N$, go to step 1.
Output $\gcd(a^{r/2}-1,\,N)$.

The expected number of repetitions is a small constant times $\log\log N$ (Theorem 15.2 and Proposition 15.7). The cost of each quantum run is dominated by the controlled modular multiplications: $m\approx2n$ of them, each costing $O(n^2)$ gates with schoolbook arithmetic, for $O(n^3)$ gates overall. Lecture 16 builds those circuits, and Lecture 17 turns gate counts into physical resource estimates for RSA-2048.

idea — where the quantum speedup really is

Everything in this section except step 2 is classical and efficient. The quantum computer does not try all factors, and it does not "try all exponents in parallel" in any useful sense: it prepares an interference pattern whose peaks sit at multiples of $1/r$, and a single sample of that pattern carries enough information to recover $r$. The structure that makes this possible, a hidden period in a group, is exactly what Lecture 16 generalises to the hidden subgroup problem, which also breaks discrete logarithms and hence elliptic-curve cryptography.

7 — Lab 4 launch: QFT, phase estimation and Shor for N = 15

lab 4 — QFT, phase estimation, and Shor for N = 15 (due in 10 days)

Goal. Implement the Fourier machinery and complete Shor's algorithm end to end on a simulator.

QFT. Build the QFT from Hadamards and controlled phases and verify it against the DFT matrix in NumPy.
Phase estimation. Estimate $\theta=1/3$ for $U=P(2\pi\theta)$ with $t=3,\dots,8$ control qubits. Plot the distributions and the probability of the best estimate against $t$.
Order finding. Using the provided modular-multiplication gates for $N=15$ and $a\in\{2,7,8,11,13\}$, run order finding, apply continued fractions, and factor 15.
Statistics. Record the success rate over 50 runs and compare it with the theoretical probability that a random $s/r$ reveals $r$.
Discussion. Why do published "hardware Shor" demonstrations for 15 and 21 usually compile the modular exponentiation away, and what would an honest demonstration need?

Hardware budget: none; this lab runs on simulators. (Section 8's first code box shows why $N=15$ is a misleading test case: most of its controlled multiplications are the identity.)

8 — Try it in code

tool — order finding for N = 15, a = 7 (tested with qiskit 2.5)

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library import QFTGate, UnitaryGate
from qiskit.quantum_info import Statevector

def mult_gate(c, N, n):
    """|y> -> |c*y mod N> for y < N (identity for y >= N), on n qubits: a permutation matrix"""
    P = np.zeros((2 ** n, 2 ** n))
    for y in range(2 ** n):
        P[(c * y) % N if y < N else y, y] = 1
    return UnitaryGate(P, label=f"x{c}")

N, a, n, m = 15, 7, 4, 8                          # work register: n = 4 qubits; m = 8 control qubits
qc = QuantumCircuit(m + n)
work = list(range(m, m + n))
qc.x(work[0])                                     # work register = |1>
qc.h(range(m))
for j in range(m):
    c = pow(a, 2 ** j, N)                         # a^(2^j) mod N, by repeated squaring (classical)
    qc.append(mult_gate(c, N, n).control(1), [j] + work)
qc.append(QFTGate(m).inverse(), range(m))
p = Statevector.from_instruction(qc).probabilities(qargs=list(range(m)))
print({y: round(float(p[y]), 4) for y in range(2 ** m) if p[y] > 1e-6})
print("multipliers a^(2^j) mod 15:", [pow(a, 2 ** j, N) for j in range(m)])

{0: 0.25, 64: 0.25, 128: 0.25, 192: 0.25}
multipliers a^(2^j) mod 15: [7, 4, 1, 1, 1, 1, 1, 1]

The simulation reproduces Example 15.4. The second line exposes a weakness of $N=15$ as a test case: since $7^4\equiv1$, every multiplier from $j=2$ on is 1, so six of the eight controlled multiplications are the identity and can be deleted. Here the permutation matrices are handed to Qiskit as UnitaryGates, which hides their cost; Lecture 16 builds them from Toffoli gates.

tool — the classical half: which a factor 35?

from math import gcd

def order(a, N):                                  # brute force: fine classically only for tiny N
    r, x = 1, a % N
    while x != 1:
        x, r = x * a % N, r + 1
    return r

N = 35
good = []
for a in range(1, N):
    if gcd(a, N) != 1:
        continue                                  # (a lucky a would already share a factor with N)
    r = order(a, N)
    if r % 2 == 0 and pow(a, r // 2, N) != N - 1:
        h = pow(a, r // 2, N)
        good.append((a, r, gcd(h - 1, N), gcd(h + 1, N)))
units = sum(1 for a in range(1, N) if gcd(a, N) == 1)
print(f"{len(good)} of {units} units a mod {N} lead to a factorisation")
print("first few (a, r, gcd(a^(r/2)-1, N), gcd(a^(r/2)+1, N)):", good[:5])

18 of 24 units a mod 35 lead to a factorisation
first few (a, r, gcd(a^(r/2)-1, N), gcd(a^(r/2)+1, N)): [(2, 12, 7, 5), (3, 12, 7, 5), (4, 6, 7, 5), (6, 2, 5, 7), (8, 4, 7, 5)]

Theorem 15.2 for $N=35$: three quarters of the units work. Python's built-in pow(a, e, N) uses repeated squaring, and math.gcd uses Euclid's algorithm, so both are fast for 2048-bit numbers; only order is hopeless at that size.

9 — Summary and exercises

summary

RSA's security rests on factoring $N=pq$; the best classical algorithm is subexponential, Shor's is polynomial.
If $a$ has even order $r$ and $a^{r/2}\not\equiv-1$, then $\gcd(a^{r/2}\pm1,N)$ are factors. A random unit works with probability at least $\tfrac12$.
$U_a|y\rangle=|ay\bmod N\rangle$ has eigenphases $s/r$, its powers are single multiplications by $a^{2^j}\bmod N$, and $|1\rangle$ is an equal superposition of its eigenvectors on the orbit of 1.
Phase estimation on $|1\rangle$ with $2^m\ge N^2$ returns $y$ with $y/2^m\approx s/r$ for a random $s$; continued fractions recover $s/r$ in lowest terms, and $a^q\equiv1$ checks the candidate.
Each run succeeds with probability about $\varphi(r)/r$; a few runs suffice. The quantum cost is dominated by modular exponentiation, $O(n^3)$ gates.

Next time: Shor II: the modular-arithmetic circuits, discrete logarithms and the hidden subgroup problem. Lab 4 is due in ten days; HW8 covers Lectures 15–16.

exercises 15

Find the order of 4 modulo 15 and use it to factor 15.
Find the order of 5 modulo 33, and use Theorem 15.1 to factor 33.
Expand $1365/2048$ as a continued fraction and list its convergents. If this outcome came from order finding for $N=35$, $a=2$, which candidate order does it give, and why is it not 12?
Why must $\gcd(a,N)=1$ for $U_a$ to be unitary? What does $y\mapsto ay\bmod N$ do when $a=5$, $N=15$?
For $N=35$, $a=2$, what is the probability that a uniformly random $s\in\{0,\dots,11\}$ is coprime to 12? Why does the widget's success rate for this $a$ come out a little lower?
Suppose order finding for $N=35$, $a=2$ were run with only $m=6$ control qubits. The outcome nearest to $64\cdot5/12$ is $y=27$. Expand $27/64$ as a continued fraction. Is $5/12$ a convergent? Which condition of Section 5 fails?

$4^2=16\equiv1$, so $r=2$ and $h=4\not\equiv-1$. $\gcd(3,15)=3$ and $\gcd(5,15)=5$.
Powers of 5 mod 33: $5,25,26,31,23,16,14,4,20,1$, so $r=10$. Then $h=5^5\equiv23\not\equiv-1\equiv32$, and $\gcd(22,33)=11$, $\gcd(24,33)=3$.
$2048=1\cdot1365+683$, $1365=1\cdot683+682$, $683=1\cdot682+1$, $682=682\cdot1$. So $1365/2048=[0;1,1,1,682]$, with convergents $0/1,\,1/1,\,1/2,\,2/3,\,1365/2048$. The last denominator below 35 is 3, so the candidate is 3, and $2^3=8\not\equiv1$: rejected. The outcome came from $s=8$, and $8/12=2/3$ reduces, so continued fractions can only see the denominator 3.
If $g=\gcd(a,N)\gt1$, multiplication by $a$ maps every $y$ to a multiple of $g$, so it is not onto and not invertible. For $a=5$, $N=15$, the images of $0,\dots,14$ are only $0,5,10$: many basis states collide, and the matrix is not unitary. (Such an $a$ factors $N$ by a gcd anyway.)
$\varphi(12)/12=4/12=\tfrac13$ ($s=1,5,7,11$). The exact per-run success probability is 0.321, a little lower, because for these $s$ a small amount of probability falls on outcomes other than the two nearest to $2048s/12$, and those can give wrong convergents. Over 50 runs the widget's rate scatters around 0.32 by several percentage points.
$64=2\cdot27+10$, $27=2\cdot10+7$, $10=1\cdot7+3$, $7=2\cdot3+1$, $3=3\cdot1$, so $27/64=[0;2,2,1,2,3]$ with convergents $0/1,\,1/2,\,2/5,\,3/7,\,8/19,\,27/64$. $5/12$ is not among them, and the method returns the candidate 19, which fails the check. Legendre's condition fails: $|27/64-5/12|=1/192$, which is larger than $1/(2\cdot12^2)=1/288$, because $2^6=64$ is far below $N^2=1225$.

Reading and sources

J. Watrous, Understanding Quantum Information and Computation, Ch. 7, the parts on order finding and factoring (arXiv:2507.11536; IBM course Fundamentals of Quantum Algorithms, Lesson 3).
R. de Wolf, Quantum Computing: Lecture Notes, Ch. 5, "Shor's factoring algorithm" (arXiv:1907.09415).
P. W. Shor, "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer", SIAM J. Comput. 26 (1997), 1484–1509 (arXiv:quant-ph/9508027); conference version, FOCS 1994.
M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, §5.3 and Appendix 4 (number theory).
G. H. Hardy and E. M. Wright, An Introduction to the Theory of Numbers, Ch. X (continued fractions).
F. Boudot, P. Gaudry, A. Guillevic, N. Heninger, E. Thomé and P. Zimmermann, "Comparing the difficulty of factorization and discrete logarithm: a 240-digit experiment", CRYPTO 2020, and the RSA-250 announcement (February 2020).

all numbers and figures in these notes are computed in python (verify_L15.py) and follow the page theme