<!-- Source: https://web3-lab.annaburd.me/how-quantum-computing-works/quantum-algorithms/ -->

# Quantum algorithms

What a quantum computer does better: query problems from Deutsch to Simon,
the cost of classical arithmetic, phase estimation and Shor's factoring, and Grover's search.

## Query-model algorithms

Two models of computation

Quantum algorithms are often analyzed in the query model, which differs from the ordinary computational model only in how the input is accessed.

Standard model

The algorithm receives the entire input.

inputx

algorithm

output

Query model

The algorithm can only interrogate a black box.

oraclef

query i

answer f(i)

algorithm

output

In the standard model, the complete input $x$ is available from the start. The algorithm may read any part of it whenever it likes, perform arbitrary computations, and eventually produce an output. Cost: the total number of elementary computational steps.

In the query model, the function $f$ is hidden inside a black box called an oracle. The algorithm never sees the function directly. Instead, it repeatedly asks questions of the form “what is $f(i)$?”, receives the answer, performs arbitrary computation, and decides which query to ask next. Cost: the number of oracle queries.

The query model isolates the cost of obtaining information from the cost of computation itself. This makes it possible to compare classical and quantum algorithms in a clean and mathematically precise way.

Examples of query problems

In the query model, the input is not a string—it is an unknown function $f$. The algorithm cannot inspect the function directly. It can only ask questions like “$\text{what is } f(x)\text{?}$” for inputs $x$ of its choice.

Hidden function

Flip the values below to change the hidden function $f$.

$x$$000$$001$$010$$011$$100$$101$$110$$111$$f(x)$

Each problem asks a different question about the same hidden function. The answer updates automatically as you change the function.

Input:

$$
f: \Sigma^n \to \Sigma
$$

Output:

$1$ if there exists a string $x \in \Sigma^n$ for which $f(x) = 1$
$0$ if there is no such string

Does any input satisfy $f(x) = 1$?

Parity

Input:

$$
f: \Sigma^n \to \Sigma
$$

Output:

$0$ if $f(x) = 1$ for an even number of strings $x \in \Sigma^n$
$1$ if $f(x) = 1$ for an odd number of strings $x \in \Sigma^n$

Is the number of inputs with $f(x) = 1$ even or odd?

Minimum

Input:

$$
f: \Sigma^n \to \Sigma^m
$$

Output:

The string $y \in \{ f(x): x \in \Sigma^n \}$ that comes first in the lexicographic ordering of $\Sigma^m$

Which output value comes first in lexicographic order?

$$
0
$$

Unique search

Input:

$$
f: \Sigma^n \to \Sigma
$$

Promise:

Exactly one input $z$ satisfies $f(z) = 1$; all other inputs satisfy $f(x) = 0$

Output:

$$
z
$$

If exactly one input satisfies $f(x) = 1$, which input is it?

$$
101
$$

Query gates

In a circuit model, access to the hidden function is represented by a query gate (or oracle gate). It behaves like an ordinary component, but its behavior is fixed by the unknown $f$: given $x$ on its input wires it outputs $f(x)$. The function is supplied by the problem instance, not the algorithm, and each use counts as one query.

Query gates can be combined with ordinary logic gates just like any other circuit component. The circuit below solves the Parity query problem for a function with two possible inputs, $0$ and $1$. It queries $f(0)$ and $f(1)$, then outputs $1$ exactly when one of the two values is $1$ and the other is $0$ (odd parity).

Why study hidden functions?

Admittedly this model looks a bit weird at first — why lock the input inside a box and count questions instead of simply reading it? But many computational tasks — from searching a database to testing a physical device — can only access information by asking questions. The query model captures exactly this situation by treating the input as an unknown function that can only be queried.

Quantum query gates

Classical query gates output the value $f(x)$ directly. That is convenient for classical circuits, but it cannot be used in quantum circuits.

The reason is that quantum gates must be unitary (and therefore reversible). A gate that simply replaces its input with $f(x)$ is generally not reversible, since many different inputs may produce the same output.

So for the quantum circuit model we choose a different definition that is always unitary. The query gate $U_f$ for any function $f: \Sigma^n \to \Sigma^m$ is defined, for all $x \in \Sigma^n$ and $y \in \Sigma^m$, by its action on basis states:

$$
U_f\bigl(\lvert y\rangle\lvert x\rangle\bigr)=\lvert y\oplus f(x)\rangle\lvert x\rangle
$$

In circuit form, $U_f$ leaves the top register holding $x$ and writes $f(x)$ into the bottom register by XOR. Notice that the function value is added into the second register rather than replacing it — this small change makes the operation reversible for every possible function $f$:

$$
\lvert x\rangle
$$

$$
\lvert y\rangle
$$

$$
\lvert y\oplus f(x)\rangle
$$

Starting the bottom register at $\lvert 0^m\rangle$ makes the gate output $f(x)$ directly, since $0^m \oplus f(x) = f(x)$.

$$
\lvert x\rangle
$$

$$
\lvert 0^m\rangle
$$

$$
\lvert f(x)\rangle
$$

The extra register may seem unnecessary at first, but it is what makes the oracle useful:

- it keeps $U_f$ unitary and reversible for every $f$;
- it lets the oracle act on a superposition of many inputs at once;
- and it preserves the phases that quantum algorithms exploit through interference.

Deutsch’s problem

Deutsch’s problem asks whether a function is constant or balanced. The function takes one bit as input and returns one bit as output, so there are only four possible functions.

Deutsch’s problem

Input:

$$
f: \Sigma \to \Sigma
$$

Output:

$0$ if $f$ is constant, $1$ if $f$ is balanced

Build a function

Flip the two outputs to define a function $f$. There are exactly four possible functions. Try to discover them all.

$a$$0$$1$$f(a)$

Found 0 of 4 possible functions of the form $f: \Sigma \to \Sigma$:

$f_1$Not found

| $a$ | $f_1(a)$ |
| --- | --- |
| $0$ | ? |
| $1$ | ? |

$f_2$Not found

| $a$ | $f_2(a)$ |
| --- | --- |
| $0$ | ? |
| $1$ | ? |

$f_3$Not found

| $a$ | $f_3(a)$ |
| --- | --- |
| $0$ | ? |
| $1$ | ? |

$f_4$Not found

| $a$ | $f_4(a)$ |
| --- | --- |
| $0$ | ? |
| $1$ | ? |

The classical approach

A classical algorithm must evaluate both possible inputs: after seeing only one value, you still cannot distinguish a constant function from a balanced one.

Deterministic classical cost: 2 queries

Deutsch’s algorithm

Deutsch’s algorithm solves the same problem using only one query. It prepares two qubits, performs a single query to the oracle $U_f$, applies one more Hadamard gate, and measures the first qubit — the measurement directly reveals $f(0) \oplus f(1)$, which is:

- $0$ for constant functions
- $1$ for balanced functions

$$
\textcolor{#6d28d9}{\lvert 0\rangle}
$$

$$
\textcolor{#b45309}{\lvert 1\rangle}
$$

$$
\begin{cases}0 & \text{if } f \text{ is constant}\\[2pt] 1 & \text{if } f \text{ is balanced}\end{cases}
$$

Step through the circuit

1. 1Prepare the two qubits in $\lvert 0\rangle\lvert 1\rangle$.
2. 2Apply a Hadamard to each qubit.
3. 3Apply the query gate $U_f$ — the single query.
4. 4Apply a Hadamard to the top qubit.
5. 5Measure the top qubit to read $f(0) \oplus f(1)$.

The Deutsch–Jozsa circuit

Deutsch’s algorithm works only for the simplest case: a function $f: \Sigma \to \Sigma$, which maps a single input bit to a single output bit. The Deutsch–Jozsa algorithm generalizes this idea to functions of the form $f: \Sigma^n \to \Sigma$ for any $n \geq 1$, allowing the input to consist of any number of bits.

$$
\lvert 0\rangle
$$

$$
\lvert 1\rangle
$$

$$
y \in \Sigma^n
$$

The purpose of the circuit is not to compute $f(x)$, but to extract information about the function $f$ as a whole. After one query, measuring the $n$ query qubits produces a bit string $y \in \Sigma^n$. The meaning of this string depends on the query problem: once we specify what property of $f$ we want to determine, we can interpret $y$ according to an appropriate decision rule.

The Deutsch–Jozsa problem

The Deutsch–Jozsa problem generalizes Deutsch’s problem: for an input function $f: \Sigma^n \to \Sigma$, the task is to output $0$ if $f$ is constant and $1$ if $f$ is balanced.

The Deutsch–Jozsa problem

Input:

$f: \Sigma^n \to \Sigma$ for some $n \geq 1$

Promise:

$f$ is either constant or balanced

Output:

$0$ if $f$ is constant, $1$ if $f$ is balanced

For $n = 1$ these are the only two possibilities, so this is exactly Deutsch’s problem. When $n \geq 2$, however, some functions $f: \Sigma^n \to \Sigma$ are neither constant nor balanced.

Build a function

Flip the four outputs to define any function $f: \Sigma^2 \to \Sigma$ and see which family it falls into.

$x$$00$$01$$10$$11$$f(x)$

$1$ of the four inputs map to $1$ — neither all of them nor half of them — so this function is neither constant nor balanced.

Input functions that are neither constant nor balanced are “don’t care” inputs. The promise excludes them, so on such a function an algorithm may output anything without being considered wrong.

The Hadamard transform

The Hadamard gate acts on the computational basis states like this:

$$
H|0\rangle = \frac{1}{\sqrt{2}}\big(|0\rangle + |1\rangle\big),\qquad H|1\rangle = \frac{1}{\sqrt{2}}\big(|0\rangle - |1\rangle\big).
$$

The only difference between these two equations is the sign of the $|1\rangle$ term. The factor $(-1)^a$ captures this perfectly: it equals $1$ when $a = 0$ and $-1$ when $a = 1$. So we can combine both cases into a single expression:

$$
H|a\rangle = \frac{1}{\sqrt{2}}\big(|0\rangle + (-1)^{a}|1\rangle\big),\qquad a \in \Sigma.
$$

Notice that the only difference between the two terms is their phase. The $|0\rangle$ term always has a positive sign, while the sign of the $|1\rangle$ term depends on the input bit $a$. We can capture both cases with the exponent $ab$, where $b$ labels the basis state in the sum:

- for $b = 0$, we have $ab = 0$, so $(-1)^{ab} = 1$, giving the positive sign of $|0\rangle$;
- for $b = 1$, we have $ab = a$, so $(-1)^{ab} = (-1)^a$, giving the correct sign of $|1\rangle$.

Therefore, both terms can be written as a single summation:

$$
H|a\rangle = \frac{1}{\sqrt{2}}\sum_{b \in \{0,1\}}(-1)^{ab}|b\rangle.
$$

From one Hadamard to many

The one-qubit Hadamard identity extends naturally to a register of $n$ qubits. Take an $n$-bit input string whose bits $x_i$ all lie in $\Sigma = \{0, 1\}$, and write the basis state it labels:

$$
x = x_{n-1}\cdots x_1 x_0, \qquad |x_{n-1}\cdots x_1 x_0\rangle.
$$

Applying a Hadamard gate to every qubit means applying $H$ independently to each bit:

$$
H^{\otimes n}|x_{n-1}\cdots x_1 x_0\rangle = \big(H|x_{n-1}\rangle\big) \otimes \cdots \otimes \big(H|x_0\rangle\big).
$$

Every qubit now becomes a superposition of $|0\rangle$ and $|1\rangle$. In the one-qubit formula the summation index was called $b$, but here we need one such index per qubit, so we rename it to $y_i$ for the $i$-th qubit. So the one-qubit identity, with $a$ renamed to $x_i$ and $b$ renamed to $y_i$, reads:

$$
H|x_i\rangle = \frac{1}{\sqrt{2}}\sum_{y_i \in \Sigma}(-1)^{x_i y_i}|y_i\rangle.
$$

Now substitute the one-qubit identity for each factor $H|x_i\rangle$ in the tensor product, using a separate index $y_i$ for each qubit:

$$
\textcolor{#0f766e}{\big(H|x_{n-1}\rangle\big)} \otimes \cdots \otimes \textcolor{#be185d}{\big(H|x_0\rangle\big)}
$$

$$
= \textcolor{#0f766e}{\left(\frac{1}{\sqrt{2}}\sum_{y_{n-1} \in \Sigma}(-1)^{x_{n-1}y_{n-1}}|y_{n-1}\rangle\right)} \otimes \cdots \otimes \textcolor{#be185d}{\left(\frac{1}{\sqrt{2}}\sum_{y_{0} \in \Sigma}(-1)^{x_{0}y_{0}}|y_{0}\rangle\right)}
$$

$$
= \frac{1}{\sqrt{2^{n}}}\sum_{\textcolor{#0f766e}{y_{n-1}} \in \Sigma}\cdots\sum_{\textcolor{#be185d}{y_{0}} \in \Sigma}\textcolor{#0f766e}{(-1)^{x_{n-1}y_{n-1}}}\cdots\textcolor{#be185d}{(-1)^{x_{0}y_{0}}}\,\big(\textcolor{#0f766e}{|y_{n-1}\rangle}\otimes\cdots\otimes\textcolor{#be185d}{|y_{0}\rangle}\big)
$$

$$
= \frac{1}{\sqrt{2^{n}}}\sum_{\textcolor{#0f766e}{y_{n-1}} \in \Sigma}\cdots\sum_{\textcolor{#be185d}{y_{0}} \in \Sigma}(-1)^{\textcolor{#0f766e}{x_{n-1}y_{n-1}}+\cdots+\textcolor{#be185d}{x_{0}y_{0}}}\,|\textcolor{#0f766e}{y_{n-1}}\cdots\textcolor{#be185d}{y_{0}}\rangle
$$

$$
= \frac{1}{\sqrt{2^{n}}}\sum_{y \in \Sigma^{n}}(-1)^{\textcolor{#0f766e}{x_{n-1}y_{n-1}}+\cdots+\textcolor{#be185d}{x_{0}y_{0}}}\,|\textcolor{#0f766e}{y_{n-1}}\cdots\textcolor{#be185d}{y_{0}}\rangle
$$

$$
= \frac{1}{\sqrt{2^{n}}}\sum_{y \in \Sigma^{n}}(-1)^{x\cdot y}|y\rangle
$$

By the multilinearity of the tensor product, the tensor product distributes over the sums, producing one term for every $n$-bit string $y$. The phase factors multiply together, so their exponents add. Thus $H^{\otimes n}$ maps $|x\rangle$ to an equal superposition of all basis states $|y\rangle$, differing only in their phases.

Walking through the circuit

Let’s follow the state as it passes through each stage of the Deutsch–Jozsa circuit.

$$
\lvert 0\rangle
$$

$$
\lvert 1\rangle
$$

$$
y \in \Sigma^n
$$

Step through the circuit

1. 1Prepare the $n$ query qubits in $|0\rangle$ and the target qubit in $|1\rangle$.
2. 2Apply a Hadamard to every qubit.
3. 3Apply the query gate $U_f$ — the single query.
4. 4Apply a Hadamard to each of the $n$ query qubits.
5. 5Measure the query register to get $y \in \Sigma^n$.

The Bernstein–Vazirani problem

Imagine that someone secretly chooses an $n$-bit string $s$.

You cannot see $s$ directly. Instead, you may query a function $f$. For any input $x$, the function looks only at the positions where the secret string has a $1$. It counts how many of those positions also contain a $1$ in $x$, and returns:

- $1$ if the count is odd,
- $0$ if the count is even.

Bernstein–Vazirani problem

Input:

$$
f: \Sigma^n \to \Sigma
$$

Promise:

there exists a binary string $s = s_{n-1}\cdots s_0$ for which $f(x) = s \cdot x$ for all $x \in \Sigma^n$

Output:

the string $s$

What does $s \cdot x$ mean?

The binary dot product works in two steps.

1. Compare the corresponding bits of $s$ and $x$.
2. Count only the positions where both bits are $1$. If this count is odd, the answer is $1$; if it is even, the answer is $0$.

For example, compare the two strings bit by bit. Only the columns where both bits are $1$ contribute to the dot product. Click any bit to change it.

$$
s
$$

$$
x
$$

$s_i x_i$$1$$0$$0$$1$$0$$1$

$$
\begin{aligned}
f(110101) &= \textcolor{#0369a1}{(1\cdot1)}\oplus\textcolor{#94a3b8}{(0\cdot1)}\oplus\textcolor{#94a3b8}{(1\cdot0)}\oplus\textcolor{#0369a1}{(1\cdot1)}\oplus\textcolor{#94a3b8}{(0\cdot0)}\oplus\textcolor{#0369a1}{(1\cdot1)}\\
&= \textcolor{#0369a1}{1}\oplus\textcolor{#94a3b8}{0}\oplus\textcolor{#94a3b8}{0}\oplus\textcolor{#0369a1}{1}\oplus\textcolor{#94a3b8}{0}\oplus\textcolor{#0369a1}{1}\\
&= 1
\end{aligned}
$$

Mathematically, this is written as follows, where multiplication is ordinary binary multiplication ($1\cdot1 = 1$, otherwise $0$), and $\oplus$ denotes XOR:

$$
s\cdot x = s_{n-1}x_{n-1}\oplus\cdots\oplus s_0x_0.
$$

The quantum algorithm

Unlike Deutsch’s and Deutsch–Jozsa’s problems, where the goal is to learn one property of the function, the Bernstein–Vazirani problem asks for the entire hidden string $s$.

Surprisingly, the quantum algorithm requires no new circuit. It uses exactly the same circuit as Deutsch–Jozsa:

- $n$ query qubits initialized to $|0\rangle$,
- one target qubit initialized to $|1\rangle$,
- Hadamard gates before and after a single query to the oracle $U_f$.

The only difference is the promise on the function. Because $f(x) = s\cdot x$, the measurement no longer reveals whether the function is constant or balanced — it reveals the hidden string $s$ itself.

$$
\lvert 0\rangle
$$

$$
\lvert 1\rangle
$$

$$
s
$$

Step through the circuit

The first three stages are identical to those of the Deutsch–Jozsa algorithm. Since we’ve already derived them, we’ll begin at the state $|\pi_3\rangle$, where the new promise on $f$ finally changes the outcome.

1. 1Prepare the $n$ query qubits in $|0\rangle$ and the target qubit in $|1\rangle$.
2. 2Apply a Hadamard to every qubit.
3. 3Apply the query gate $U_f$ — the single query.
4. 4Apply a Hadamard to each of the $n$ query qubits.
5. 5Measure the query register to read $y = s$.

Simon’s problem

As in Bernstein–Vazirani, someone secretly chooses an $n$-bit string $s$, and the task is to recover it. What changes is how the function hides it.

The function $f$ no longer returns a single bit but a whole string, and no individual value $f(x)$ tells you anything about $s$. Instead, $s$ is written into the pattern of collisions: $f$ gives the same answer on two different inputs exactly when those inputs differ by $s$.

Simon’s problem

Input:

$$
f: \Sigma^n \to \Sigma^m
$$

Promise:

there exists a string $s \in \Sigma^n$ such that

$$
\big[f(x) = f(y)\big]\iff\big[(x = y)\ \text{ or }\ (x\oplus s = y)\big]
$$

for all $x, y \in \Sigma^n$

Output:

the string $s$

The promise says that $x$ and $y$ collide only in the two ways it lists: either they are the same input, or one is the other shifted by $s$. Which of those matters depends on whether $s$ is the all-zero string.

Case 1: $s = 0^n$

Shifting by $s$ changes nothing, since $x \oplus 0^n = x$, so both branches of the promise say the same thing and the condition simplifies to

$$
\big[f(x) = f(y)\big]\iff\big[x = y\big]
$$

This is exactly the definition of one-to-one. On three bits, all eight inputs have different outputs, so a query never repeats a value.

$x$$f(x)$

$000$$101$

$001$$010$

$010$$111$

$011$$001$

$100$$110$

$101$$011$

$110$$100$

$111$$000$

Case 2: $s \neq 0^n$

Now $x \oplus s$ is a genuinely different input from $x$, and the promise forces the two to agree:

$$
f(x) = f(x\oplus s)
$$

Every input is paired with exactly one partner, and the promise also rules out any other coincidence, so different pairs must have different outputs. The function is therefore two-to-one. With $s = 110$, the eight inputs collapse onto four outputs:

$x,\ x \oplus s$$f(x)$

$000,\ 110$$101$

$001,\ 111$$010$

$010,\ 100$$111$

$011,\ 101$$001$

Nothing in a single answer points at $s$. It shows up only once two inputs are found to share an output, and then $s = x \oplus (x \oplus s)$.

The two cases are what makes the problem hard classically. Learning $s$ means finding a collision, and a classical algorithm has no way to force one: it can only keep querying inputs and comparing the answers it has already seen.

Simon’s algorithm

Simon’s algorithm consists of running the following circuit several times, followed by a post-processing step. The circuit is the familiar shape — Hadamards, one query, Hadamards, measurement — with two changes forced by the new function:

- the workspace is now $m$ qubits rather than one, since $f$ returns a string of $m$ bits;
- those qubits start in $|0\rangle$ and carry no gates at all — not even a Hadamard.

$$
\textcolor{#6d28d9}{\lvert 0\rangle}
$$

$$
\textcolor{#b45309}{\lvert 0\rangle}
$$

$$
y \in \Sigma^n
$$

Step through the circuit

1. 1Prepare the $n$ query qubits and the $m$ workspace qubits in $|0\rangle$.
2. 2Apply a Hadamard to each query qubit.
3. 3Apply the query gate $U_f$ — the single query.
4. 4Apply a Hadamard to each query qubit again.
5. 5Measure the query register to read $y \in \Sigma^n$.
6. 6Repeat steps 1–5, then solve the collected equations $y\cdot s=0$ for $s$ — the one step that is classical, not the circuit.

Up to this point the practical value of these algorithms is thin, and the accounting is generous. The speedup is counted in oracle queries while everything around the query is assumed free. Someone still has to build the gate for $f$, which can easily cost more than the queries it saves. The promise has to hold, and rejecting a function that fails it is roughly the problem you started with. Measurements come back noisy and runs have to be repeated. And all of it takes far more thought than the same job written in ordinary 0-1 bits.

Whether that changes further on, we will see. It is too early for disappointment either way. The road so far is a sequence of historical milestones, each adding a piece of the knowledge the later algorithms are built from:

- Deutsch showed quantum computation could outperform classical in principle: one query instead of two.
- Deutsch–Jozsa demonstrated an exponential separation in the query model, under a promise, and only against classical algorithms that must be exactly right every time.
- Bernstein–Vazirani found a separation that randomness cannot close, and applied recursively, a superpolynomial one.
- Simon introduced hidden-period techniques, and gave the first exponential separation against randomized classical algorithms.

All four are statements about the query model, where the only cost counted is the number of calls to the oracle, and where the function comes with a promise attached. Neither assumption holds outside it, and nothing here proves quantum computers are faster on ordinary inputs. The machinery does carry over though. Superposition, phase kickback and interference reappear in algorithms that are handed no black box at all. Whether that finally repays the trouble is an open question at this point.

## The cost of classical algorithms

Measuring cost

In the query model there was exactly one thing to count. Outside it there is no oracle to call, so before any classical and quantum algorithm can be compared on a real problem, we need a yardstick that works for both.

An abstract view of computation

Whatever the computational model, the input and output are binary strings.

inputx

computation

outputy

The middle box could be a Turing machine, a Boolean circuit, a quantum circuit or a Python program. Only the computation changes. Inputs and outputs remain binary strings, and numbers, vectors, matrices, graphs, or molecules all enter the computation through an appropriate binary encoding.

Input length

There is rarely a single standard encoding. We choose one, and the details matter less than they seem: converting between any two reasonable encodings adds only a small overhead. What the choice does determine is the input length: the number of bits in the encoded input. For a nonnegative integer written in binary,

$$
\lg(N)=\begin{cases}1, & N = 0,\\[2pt]1+\lfloor\log_2 N\rfloor, & N \geq 1.\end{cases}
$$

| number | binary encoding | length |
| --- | --- | --- |
| 0 | 0 | 1 |
| 5 | 101 | 3 |
| 12 | 1100 | 4 |
| 1 000 000 | 11110100001001000000 | 20 |
| a 617-digit RSA modulus | 1011…0111 | 2048 |

This is the key idea. The input length grows logarithmically with the number it represents. A 2048-bit input therefore describes a number close to $2^{2048}$. An algorithm that tests every divisor up to $\sqrt{N}$ performs about $2^{n/2}$ operations on an $n$-bit input. It may look efficient when measured against $N$, but it is exponential when measured against the true input size $n$.

Elementary operations

The cost of a circuit is measured by the number of elementary gate applications it performs. Which gates are considered elementary is a modeling choice: we first fix a gate set, and each application of a gate from that set counts as one computational step. The set does not have to be minimal—some gates may themselves be implementable using other gates in the same set.

ANDORNOTFANOUT

We count FANOUT as a gate. It is often treated as free, but making it explicit highlights an important contrast: classical circuits can copy bits freely, whereas quantum circuits cannot.

$X$$Y$$Z$$H$$S$$S^\dagger$$T$$T^\dagger$

CNOTmeasurement

This gate set is universal: any unitary operation can be approximated to arbitrary accuracy using only these gates.

Size and depth

The size of a circuit is the total number of gates in it. Its depth is the largest number of gates on any path from an input wire to an output wire. Size corresponds to sequential running time, while depth corresponds to parallel running time.

Cost as a function of input length

A circuit has a fixed number of input wires, so it accepts inputs of only one length, and its cost is simply its size, $\mathrm{cost}(C)=\mathrm{size}(C)$. An algorithm, however, must work for inputs of arbitrary length. It is therefore represented by a family of circuits $\{C_1, C_2, \ldots\}$, where $C_n$ handles $n$-bit inputs. The cost of the algorithm is then the size of the circuit for each input length:

$$
t(n)=\mathrm{size}(C_n).
$$

For example, a classical factoring algorithm is a family of Boolean circuits, while a quantum factoring algorithm is a family of quantum circuits. Both solve the same problem on $n$-bit inputs, differing only in the gate set they use.

This lets us compare algorithms by how $t(n)$ grows with the input length. An algorithm is considered efficient if $t(n)$ is bounded by a polynomial in $n$.

Cost analysis: integer addition

Now that cost is defined as a function of the input length, we can work through a complete example. The simplest one is integer addition: given two integers $N$ and $M$, compute their sum $N + M$. Both inputs are provided in binary.

The algorithm itself is familiar from elementary school. What changes is the model of computation. Instead of describing the sequence of arithmetic steps, we must build the algorithm as a Boolean circuit from elementary gates and determine its cost. Later, we will construct the same algorithm on a quantum circuit and compare how the resource requirements differ.

The algorithm

Binary addition follows the same schoolbook procedure as decimal addition. Starting with the least significant bit, each column adds the two input bits together with the carry from the previous column, producing a sum bit and a new carry for the next column.

The important observation is that every column performs exactly the same computation. It receives three input bits—the operand bits $x_i$ and $y_i$, and the incoming carry $c_i$—and produces two output bits: the sum $s_i$ and the outgoing carry $c_{i+1}$.

bit 7bit 6bit 5bit 4bit 3bit 2bit 1bit 0carries11111000$N = 156$$+\; M = 107$$N + M = 263$100000111

Building the Boolean circuit

The addition algorithm consists of one operation repeated for every bit position. We therefore start by building a circuit for a single column. Once that building block is complete, the full adder is obtained simply by connecting copies of it together.

Half adder

Every column of the addition except the least significant one—bit $0$, where there is nothing to carry from—must also handle an incoming carry. Starting from a half adder, we add a second half adder to incorporate the carry, then combine the two possible carry outputs with an OR gate. The result is a full adder, implementing the three-input, two-output function performed by every column.

Full adder

Half adder

The complete adder is built by repeating the same full adder circuit, with the carry propagating from one bit to the next.

Half adder

Full adder

Count the gates

An $n$-bit adder is one half adder and $n-1$ full adders, so

$$
t(n)=10+21(n-1)=21n-11.
$$

Whether that constant comes out as 21, or 31, or something else again depends on the gate set and on how the XOR is expanded, and it is not what we are after. What the construction establishes is that there exists a family $\{C_1, C_2, \ldots\}$ of Boolean circuits, where $C_n$ adds two $n$-bit nonnegative integers together, such that $\mathrm{size}(C_n) = O(n)$.

Asymptotic notation

The exact number of gates depends on implementation details such as the gate set or the choice of intermediate operations. These differences affect only constant factors, while the overall growth of the algorithm stays the same. Asymptotic notation describes that growth by ignoring constant multipliers and lower-order terms.

The most commonly used notation is Big O, which gives an upper bound on the growth rate of a function. For two functions $g(n)$ and $h(n)$, we write that $g(n) = O(h(n))$ if there exists a positive real number $c > 0$ and a positive integer $n_0$ such that $g(n) \leq c \cdot h(n)$ for all $n \geq n_0$.

cost (time)input length n

Growth classes

Examples

Polynomial, Subexponential, and Exponential Growth

These three names describe how an algorithm's cost grows with the input size. The chart shades these regions on its log scale.

Polynomial — $O(n^{b})$ for a fixed $b>0$. This is the usual boundary for what we call efficient.

Subexponential — $2^{o(n)}$: the exponent grows more slowly than $n$. A stricter definition requires $O(2^{n^{\varepsilon}})$ for every $\varepsilon>0$. The number field sieve is subexponential under the first definition, but not under this stricter one.

Exponential — $2^{\Theta(n)}$: the exponent grows linearly with $n$. In particular, an algorithm that is not subexponential is not automatically exponential. There is a gap between the two classes.

The [exponential-time hypothesis](https://en.wikipedia.org/wiki/Exponential_time_hypothesis) (ETH) conjectures that NP-complete problems have no subexponential-time algorithms.

Cost analysis: integer multiplication

The next example is one step up from addition: given two integers $N$ and $M$, compute their product $N \cdot M$. Both inputs are again provided in binary. As with addition, the algorithm itself is familiar. The task is to express it as a Boolean circuit and determine how its cost grows with the input length.

The algorithm

Binary long multiplication follows the same procedure as decimal long multiplication. For each bit of $M$, we form a partial product by either copying $N$ or producing a row of zeros, depending on whether that bit is $1$ or $0$. Each partial product is then shifted according to the position of the corresponding bit of $M$. Adding all of these shifted rows gives the final product.

bit 7bit 6bit 5bit 4bit 3bit 2bit 1bit 0$N = 13$$\times\; M = 11$

$M_0\,{=}\,1:\ N \ll 0$1101

$M_1\,{=}\,1:\ N \ll 1$1101

$M_2\,{=}\,0$0000

$M_3\,{=}\,1:\ N \ll 3$1101

$N \cdot M = 143$10001111

The key observation is that every bit of every partial product depends on exactly two input bits: one bit from $N$ and one bit from $M$. This gives us a simple building block for the first stage of the circuit.

Building the Boolean circuit

For a pair of bits $N_i$ and $M_j$, the corresponding partial-product bit is $1$ exactly when both bits are $1$. This is precisely the function computed by an AND gate.

Partial-product bit

We therefore obtain all partial products by arranging these AND gates in a grid. For two $n$-bit inputs, there is one gate for every pair $(i, j)$, giving an $n \times n$ array and therefore $n^2$ AND gates.

Partial-product array

This produces the partial products, but they still have to be added together. Here we can reuse the $n$-bit adder from the previous example. The shifted partial products are added one after another, requiring $n - 1$ such additions.

0 0 0 0 1 1 0 1

0 0 0 1 1 0 1 0

n-bit adder #1

0 0 1 0 0 1 1 1

0 0 0 0 0 0 0 0

n-bit adder #2

0 0 1 0 0 1 1 1

0 1 1 0 1 0 0 0

n-bit adder #3

1 0 0 0 1 1 1 1

Count the gates

Counting the two stages: the array contributes $n^2$ AND gates, and the summation contributes $n - 1$ adders of $O(n)$ gates each. The total is

$$
t(n)=\underbrace{n^2}_{\text{partial products}}+\underbrace{(n-1)\cdot O(n)}_{\text{summation}}=O(n^2).
$$

So there is a family $\{C_1, C_2, \ldots\}$ of Boolean circuits, where $C_n$ multiplies two $n$-bit nonnegative integers, with $\mathrm{size}(C_n) = O(n^2)$. By the standard multiplication algorithm, there are Boolean circuits of size $O(n^2)$ for multiplying $n$-bit integers.

More generally, the same array argument with an $n \times m$ grid gives circuits of size $O(nm)$ for multiplying an $n$-bit integer by an $m$-bit integer.

Faster multiplication: convolution and the Fourier transform

Schoolbook multiplication costs $O(n^2)$, and for a long time that was taken to be optimal. In 1960, Karatsuba showed that it was not, using divide and conquer to reduce multiplication to a smaller number of multiplications.

The same search for structure leads further: the pairwise products of schoolbook multiplication form a convolution, and the Fourier transform provides a way to compute that convolution efficiently.

Multiplying in blocks

An $n$-bit integer can be split into $k$ blocks of $b$ bits, with each block treated as a single digit in base $B = 2^b$. Thus $k = \lceil n/b \rceil$, $N = (a_0, a_1, \ldots, a_{k-1})$ and $M = (c_0, c_1, \ldots, c_{k-1})$.

Block widthb = 2 bits

$8$ bits → $4$ blocks of $2$$B = 2^{2} = 4$$a_i, c_j \in \{0, 1, 2, 3\}$

$$
N
$$

a₀1a₁2a₂1a₃3

$$
N = (3121)_{4} = 217
$$

$$
M
$$

c₀2c₁1c₂3c₃2

$$
M = (2312)_{4} = 182
$$

Schoolbook multiplication of these block digits forms every product $a_i c_j$, giving $k^2$ block products. This is not a saving by itself: larger blocks give fewer products, but each product is a multiplication of wider numbers.

c₀ = 2c₁ = 1c₂ = 3c₃ = 2

a₀ = 1a₀c₀2a₀c₁1a₀c₂3a₀c₃2

a₁ = 2a₁c₀4a₁c₁2a₁c₂6a₁c₃4

a₂ = 1a₂c₀2a₂c₁1a₂c₂3a₂c₃2

a₃ = 3a₃c₀6a₃c₁3a₃c₂9a₃c₃6

The reason for changing to blocks is that they make the structure of the product visible. A cell of the array represents $(a_i B^{\,i})(c_j B^{\,j}) = a_i c_j \, B^{\,i+j}$. Summing over all cells therefore gives $N \cdot M = \big(\textstyle\sum_i a_i B^{\,i}\big)\big(\sum_j c_j B^{\,j}\big) = \sum_i \sum_j a_i c_j \, B^{\,i+j}$. The index sum $i+j$ determines where each product contributes: cells with the same index sum multiply the same power of $B$, so their products can be added together.

Given this structure, the question is therefore how to combine the $k^2$ products more efficiently, rather than compute and handle each one separately.

Convolution

From the product table above, we already know that products with the same index sum $i+j$ belong together. Let $d_l = \sum_{i+j=l} a_i c_j$. Then the product can be written as $N M = \sum_l d_l B^{\,l}$. The sequence $(d_0, d_1, \ldots, d_{2k-2})$ is the convolution of the block sequences $(a_0, a_1, \ldots, a_{k-1})$ and $(c_0, c_1, \ldots, c_{k-1})$.

So the $k^2$ cells of the multiplication array collapse into just $2k-1$ diagonal sums:

d₀ · B⁰ = 2 · 12

d₁ · B¹ = 5 · 420

d₂ · B² = 7 · 16112

d₃ · B³ = 15 · 64960

d₄ · B⁴ = 10 · 2562560

d₅ · B⁵ = 11 · 102411264

d₆ · B⁶ = 6 · 409624576

N × M = 217 × 18239494

Calculating convolution

The same sum can be pictured two ways: as a diagonal in the product table above, or by sliding the reversed $M$ blocks under the $N$ blocks.

a₀1a₁2a₂1a₃3

c₃2c₂3c₁1c₀2

d₀ = a₀c₀ = 2

Nothing has become faster yet. The same $k^2$ pairwise products still appear in the definition of the convolution.

But we have changed what we are trying to compute. Schoolbook multiplication computes every $a_i c_j$ and immediately assigns it to a diagonal. The individual products are discarded after contributing to their diagonal sum. The result we actually need is only $d_0, d_1, \ldots, d_{2k-2}$.

So the problem can now be stated precisely: can we compute all the convolution coefficients $d_l$ without computing all $k^2$ products $a_i c_j$ individually? That is the problem the Fourier transform will solve.

Another way to compute convolution

So far, the coefficients $d_0, d_1, \ldots, d_{2k-2}$ have been a sequence of numbers attached to powers of the base in $NM = \sum_l d_l B^{\,l}$. Instead of fixing the base at $B$, leave it as a variable: the same structure becomes a polynomial, which for the blocks above is $2 + 5x + 7x^{2} + 15x^{3} + 10x^{4} + 11x^{5} + 6x^{6}$.

At the same time, the two multiplicands can be written as polynomials too—$N(x) = \sum_i a_i x^i$ and $M(x) = \sum_j c_j x^j$. Multiplying them gives $D(x) = N(x)\,M(x) = \sum_l d_l x^l$. So the coefficients of $D(x)$ are exactly the convolution coefficients we want: we have simply turned the two input sequences into polynomials and the convolution into their product.

Now comes the useful part: a polynomial can be represented in another way—not by its coefficients, but by its values at enough distinct points. A degree-$d$ polynomial is completely determined by $d+1$ such values. In this representation, multiplication becomes much simpler. At every point $D(x) = N(x)\,M(x)$, so we can evaluate $N$ and $M$, multiply the corresponding values, and obtain the values of $D$. There are no cross terms: just one ordinary multiplication per point.

Since $D$ has degree $2k-2$, $2k-1$ values are enough to recover all its coefficients. The strategy is therefore:

Blocks of Na₀, a₁, …

Blocks of Mc₀, c₁, …

evaluate

pointwise ×

interpolate

Convolutiond₀, d₁, …

at x₀

N(x₀)−0.49

M(x₀)2.05

N(x₀) × M(x₀) = D(x₀)−1

Points → recovered coefficients of D(x)

2, 5, 7, 15, 10, 11, 6

We have recovered all the convolution coefficients—but we have not made the computation faster yet. Evaluating the polynomials and interpolating the result still costs $O(k^2)$ when done naively. The key question is therefore not whether this representation works, but whether we can choose the evaluation points so that the evaluations themselves can be computed efficiently—that is where the Fourier transform enters.

Choosing the evaluation points

We want evaluation points where one evaluation can reuse work from another. A natural pair to try is $x$ and $-x$.

The two points differ only in the sign of $x$. To make that useful, sort the coefficients of $N$ by whether their position is even or odd. Call the two halves $N_{\mathrm{e}}$ and $N_{\mathrm{o}}$; both are polynomials in $x^2$, and $N(x) = N_{\mathrm{e}}(x^2) + x\,N_{\mathrm{o}}(x^2)$—for example, $N(x) = 1 + 2x + x^{2} + 3x^{3} = (1 + x^{2}) + x(2 + 3x^{2})$. Now the advantage is visible: replacing $x$ by $-x$ leaves $x^2$ unchanged, so the even half stays the same while the odd half changes sign, $N(-x) = N_{\mathrm{e}}(x^2) - x\,N_{\mathrm{o}}(x^2)$.

Thus both $N(x)$ and $N(-x)$ can be calculated from the same two quantities, $N_{\mathrm{e}}(x^2)$ and $N_{\mathrm{o}}(x^2)$. Once these are known, the two results require only one multiplication and two additions: form $x \cdot N_{\mathrm{o}}(x^2)$ once, then add it to and subtract it from $N_{\mathrm{e}}(x^2)$. The important part is that $N_{\mathrm{e}}$ and $N_{\mathrm{o}}$ each have half as many coefficients as $N$. The same rule therefore applies to them, and to their halves in turn—each split naming its pieces by the choices that made them, so that $N_{\mathrm{oe}}$ is the even half of $N$’s odd half. Repeating this keeps halving the size of the problem.

The caveat is that we have only used the $(x, -x)$ pairing once. To keep saving work, the points left after that first split must themselves form $(x, -x)$ pairs, so that the same idea can be applied again. Real numbers do not work. Once we square them, all points become nonnegative, so the pairing is lost. We therefore move to the complex plane.

The roots of unity have exactly the structure we need: $\omega_j = e^{2\pi i j/k}$, $j = 0, \ldots, k-1$. They are evenly spaced around the unit circle. Each point has an opposite partner, $\omega_{j+k/2} = -\omega_j$, and squaring sends each pair to the same point. The resulting points are again evenly spaced, so the pairing survives and the process can repeat: $k \rightarrow k/2 \rightarrow k/4 \rightarrow \cdots \rightarrow 1$.

u = x²v = u²q = v²

We now have $8$ evaluation points, paired as $x$ and $-x$. To make each pair share the same work, we separate $N(x)$ into its even- and odd-power terms: $N(x)=N_{\mathrm{e}}(x^2)+x\,N_{\mathrm{o}}(x^2)$.

This rewrite makes $x^2$ the input to both smaller polynomials. For each pair $(x,-x)$, this input is the same because $x^2=(-x)^2$.

Thus the $8$ original points give only $4$ distinct inputs for $N_{\mathrm{e}}$ and $N_{\mathrm{o}}$: $\omega_{0}, \omega_{2}, \omega_{4}, \omega_{6}$.

We can therefore evaluate the two smaller polynomials using just these $4$ points.

From coefficients to values at the roots of unity

We have now split the polynomials down to their individual coefficients. The next step is to reverse the process: combine the pieces back up to obtain the values of $N$ and $M$ at the chosen roots of unity. These values are exactly what we need to multiply the polynomials pointwise.

1. Split down to individual coefficients

We apply the same recursive split to both $N(x)$ and $M(x)$, until every branch contains a single coefficient.

For $N(x)$, the leaves are $\color{#b9785c}{a_{0}=1},\color{#b9785c}{a_{1}=2},\color{#b9785c}{a_{2}=1},\color{#b9785c}{a_{3}=3}$, and $\color{#5f8f7f}{a_{4}},\color{#5f8f7f}{a_{5}},\color{#5f8f7f}{a_{6}},\color{#5f8f7f}{a_{7}}=0$.

For $M(x)$, following the same process as with $N(x)$, we obtain $\color{#b9785c}{c_{0}=2},\color{#b9785c}{c_{1}=1},\color{#b9785c}{c_{2}=3},\color{#b9785c}{c_{3}=2}$, and $\color{#5f8f7f}{c_{4}},\color{#5f8f7f}{c_{5}},\color{#5f8f7f}{c_{6}},\color{#5f8f7f}{c_{7}}=0$.

These $a_k$ and $c_k$ are now the individual coefficients that we recombine upward.

The coefficients themselves have not changed. They are still the original coefficients of the two polynomials. What has changed is how they are organized. At each split, we separate even and odd powers, and the resulting branches record these choices. The tree makes this recursive structure explicit. We can then reverse the same structure to combine the coefficients and evaluate the polynomials at all the roots of unity.

2. Recombine upward

Reverse the same tree. At each level, combine the even and odd pieces until we have evaluated both polynomials at all $8$ chosen roots of unity:

$$
\{a_k\}\xrightarrow{\text{combine upward}}\{N(\omega_j)\},\qquad \{c_k\}\xrightarrow{\text{combine upward}}\{M(\omega_j)\}.
$$

At each point $\omega_j$, multiply the two values:

$$
D(\omega_j)=N(\omega_j)M(\omega_j).
$$

| Point | $N(\omega_j)$ | $M(\omega_j)$ | $D(\omega_j)$ |
| --- | --- | --- | --- |
| $\omega_{0}=1$ | $7$ | $8$ | $56$ |
| $\omega_{1}=0.707+0.707i$ | $0.293+4.536i$ | $1.293+5.121i$ | $-22.849+7.364i$ |
| $\omega_{2}=i$ | $-i$ | $-1-i$ | $-1+i$ |
| $\omega_{3}=-0.707+0.707i$ | $1.707+2.536i$ | $2.707-0.879i$ | $6.849+5.364i$ |
| $\omega_{4}=-1$ | $-3$ | $2$ | $-6$ |
| $\omega_{5}=-0.707-0.707i$ | $1.707-2.536i$ | $2.707+0.879i$ | $6.849-5.364i$ |
| $\omega_{6}=-i$ | $i$ | $-1+i$ | $-1-i$ |
| $\omega_{7}=0.707-0.707i$ | $0.293-4.536i$ | $1.293-5.121i$ | $-22.849-7.364i$ |

3. Interpolate

The $8$ values $D(\omega_j)$ determine the degree-$6$ product uniquely. Interpolating them gives

$$
\boxed{D(x)=2 + 5x + 7x^{2} + 15x^{3} + 10x^{4} + 11x^{5} + 6x^{6}}.
$$

So we have multiplied the two polynomials without forming all $k^2$ coefficient products.

The complete process is:

$$
\{a_k\},\{c_k\}\xrightarrow{\text{evaluate}}\{N(\omega_j)\},\{M(\omega_j)\}\xrightarrow{\text{multiply}}\{D(\omega_j)\}\xrightarrow{\text{interpolate}}D(x).
$$

But there is still one important question: **how did we evaluate the polynomials at all those roots of unity efficiently?**

Naming the operation: DFT and FFT

The operation we have just performed—taking the coefficients of a polynomial and evaluating it at the roots of unity—is the Discrete Fourier Transform (DFT). For $N(x)$, the DFT takes $a_0,a_1,\ldots,a_{k-1}$ and produces $N(\omega_0),N(\omega_1),\ldots,N(\omega_{k-1})$.

The recursive even/odd splitting above is what makes this evaluation fast. Recall that $N(x)=N_{\mathrm{e}}(x^2)+xN_{\mathrm{o}}(x^2)$ while $N(-x)=N_{\mathrm{e}}(x^2)-xN_{\mathrm{o}}(x^2)$. We evaluate the smaller polynomials $N_{\mathrm{e}}$ and $N_{\mathrm{o}}$ once; if their values are $u$ and $v$, the paired results are $u+xv$ and $u-xv$. This combination is one butterfly.

Because the roots of unity are paired as $\omega$ and $-\omega$, squaring them gives the points needed by the smaller transforms. The same split therefore repeats, $8 \rightarrow 4 \rightarrow 2 \rightarrow 1$. The recursion does not merely divide the problem into smaller pieces: each smaller problem has exactly the same structure as the original one.

The Fast Fourier Transform (FFT) is this recursive algorithm for computing the DFT efficiently.

$$
\boxed{\text{DFT}=\text{what we compute}}
$$

$$
\boxed{\text{FFT}=\text{how we compute it quickly}}
$$

At each level, the butterflies combine the results of the two half-size transforms. There are $O(k)$ butterfly operations per level and $\log_2 k$ levels, giving $T(k)=2T(k/2)+O(k)=O(k\log k)$. The same applies to $M(x)$, so we can write the complete multiplication algorithm compactly as

$$
\boxed{\operatorname{FFT}(N),\operatorname{FFT}(M)\;\longrightarrow\;\text{pointwise multiplication}\;\longrightarrow\;\text{inverse FFT}}
$$

The inverse transform takes the values $D(\omega_j)$ back to the coefficients $d_0,d_1,\ldots,d_{2k-2}$, which are exactly the convolution coefficients we wanted. Thus the Fourier transform turns convolution into pointwise multiplication, $\operatorname{DFT}(a*c)=\operatorname{DFT}(a)\odot\operatorname{DFT}(c).$ The $k^2$ pairwise products have been replaced by two fast transforms, $k$ pointwise products, and one inverse transform.

The Schönhage–Strassen algorithm

The Fourier transform gave us a fast way to multiply *polynomials*. But originally we set out to multiply *integers*, and for them there are still open questions:

- Can the transform be made exact? The roots of unity it evaluates at are *complex numbers*, and their coordinates are irrational, so complex arithmetic is only ever approximate—but a product of integers has to come out exactly.
- Can the leftover multiplications be removed? The $k$ pointwise products are still multiplications of integers. The transform has shrunk their operands—to about $n/k$ bits each—but has not made them go away.

And the Schönhage–Strassen algorithm closes both.

**The first fix is to change where the arithmetic happens.** Instead of the complex plane, run the transform inside [modular arithmetic](https://web3-lab.annaburd.me/math-behind-key-pairs/#modular-foundations)—the integers modulo $2^m + 1$. There $2^m \equiv -1$, so $2^{2m} \equiv 1$. The number $2$ therefore behaves as a $2m$-th root of unity: its powers close into a cycle, with the opposite points satisfying $2^{j+m} \equiv -2^{j}$.

exponent

$$
2^{0} \equiv 1 \pmod{2^{4}+1}
$$

The highlighted $\pm$ pair, opposite on the ring:

$$
2^{0} = 1
$$

$$
2^{4} = 16 \equiv -1
$$

Square both. Since $2^{8}\equiv 1$, the extra factor drops and they meet:

$$
(2^{0})^2 = 2^{0} \equiv 1
$$

$$
(2^{4})^2 = 2^{8} \equiv 2^{0} \equiv 1
$$

Both land on the same point, $2^{0}\equiv 1$—just as $\omega$ and $-\omega$ square to $\omega^{2}$. That collapse halves the points, and the transform recurses on what remains.

This change solves two problems at once. The transform is now exact, because nothing is represented by an approximate complex number. And every root of unity is a power of $2$, so multiplying by one is just a shift of the bits, with the part that runs off the top folded back with a minus sign. The transform therefore needs only additions and shifts, and costs $O(n \lg n)$ bit operations.

**The second fix is to recurse.** Each of the $k$ pointwise products is a multiplication of much smaller integers—the very problem we began with, in miniature. So we solve those products using the same algorithm. This recursion is the heart of Schönhage–Strassen.

To analyze the cost, one free parameter remains: **how many blocks should we use?** Cutting an $n$-bit integer into $k$ blocks gives blocks of about $n/k$ bits, and the transform turns the multiplication into $k$ pointwise multiplications of numbers that size. There is a trade-off: fewer blocks mean larger pointwise multiplications, more blocks mean a longer transform. Balancing the two costs gives $k \approx \sqrt{n}$, so each block has about $n/k \approx \sqrt{n}$ bits. At one level of the algorithm we therefore have:

- an $n$-bit integer split into about $\sqrt{n}$ blocks of about $\sqrt{n}$ bits each;
- two forward transforms and one inverse transform, costing $O(n \lg n)$ additions and shifts;
- about $\sqrt{n}$ pointwise multiplications;
- each pointwise multiplication multiplying two $\sqrt{n}$-bit numbers, producing a result of about $2\sqrt{n}$ bits.

That last point is crucial. The recursive problems are multiplications of roughly $2\sqrt{n}$-bit numbers, and there are about $\sqrt{n}$ of them, so the recursive part contains $\sqrt{n} \cdot 2\sqrt{n} = 2n$ bits in total—only a constant factor more than the $n$ input bits. This gives the recurrence

$$
t(n)=O(n\lg n)+\sqrt{n}\,t(2\sqrt{n}),
$$

where the first term is the work done by the transforms and the second is the cost of the $\sqrt{n}$ recursive multiplications. Now look at what happens as we recurse: each level square-roots the operand size, $n \to \sqrt{n} \to \sqrt{\sqrt{n}} \to \cdots$. At first this may look as though the recursive work should become dramatically smaller, but there are more and more subproblems at each level, and the number of bits across all of them grows by the same factor that the logarithm of their size shrinks. So, up to constant factors, each level still costs $O(n \lg n)$:

| Level | Operand size | Bits in total | Transform work |
| --- | --- | --- | --- |
| 0 | $n$ | $n$ | $n\lg n$ |
| 1 | $2\sqrt{n}$ | $2n$ | $2n\cdot\tfrac{1}{2}\lg n=n\lg n$ |
| 2 | $\approx 2\sqrt{\sqrt{n}}$ | $4n$ | $4n\cdot\tfrac{1}{4}\lg n=n\lg n$ |
| r | $\approx 2n^{1/2^{r}}$ | $2^{r}n$ | $2^{r}n\cdot\tfrac{\lg n}{2^{r}}=n\lg n$ |

The only thing left to determine is how many levels there are. Taking a square root halves the exponent, so after $r$ levels the operands are about $n^{1/2^{r}}$ bits wide, and we stop when that reaches constant size. Taking logarithms, the condition reads $\lg n / 2^{r} = O(1)$, or $2^{r} = \Theta(\lg n)$, so $r = O(\lg\lg n)$. Each of those levels costs $O(n \lg n)$, giving

$$
O(n\lg n)\times O(\lg\lg n)=\boxed{t(n)=O(n\lg n\lg\lg n)}
$$

That is where the second logarithm comes from: the transform itself costs only $O(n \lg n)$ bit operations, and the extra $\lg\lg n$ is the price of repeating that work over $O(\lg\lg n)$ levels of recursion.

For the small example above, all of this machinery is obviously overkill. Splitting, padding, transforming, and rebuilding the result carry their own overhead. The advantage appears only for sufficiently large inputs, when replacing $k^2$ pairwise block products with $O(k\lg k)$ transform work saves more than that setup costs.

Beyond Schönhage–Strassen

Schönhage–Strassen remained the asymptotically fastest known integer multiplication algorithm for decades. In 2019, Harvey and van der Hoeven presented [Integer multiplication in time O(n log n)](https://annals.math.princeton.edu/2021/193-2/p04), an algorithm of complexity $O(n \lg(n))$, which is conjectured to be optimal up to constant factors.

Cost analysis: integer division

Given two integers $N$ and $M$, integer division computes a quotient $Q$ and a remainder $R$ such that $N = QM + R$ with $0 \le R < M$. As before, the inputs are given in binary. But division differs from addition and multiplication in an important way: it must make a decision at each step. Given the current partial remainder, it has to determine whether $M$ fits and, depending on the answer, either subtract $M$ or leave the remainder unchanged.

A Boolean circuit cannot branch on this decision. Instead, it has to implement the decision itself using logic gates. This makes the cost of division more interesting to analyze than the straightforward bit-by-bit operations we have seen so far.

The algorithm

Decimal long division is awkward because, at each step, we must determine the next quotient digit from several possibilities. In binary, that choice disappears: the next quotient bit can only be $0$ or $1$. So each step reduces to a single question: **does the divisor fit into the current partial remainder?**

The algorithm—shift and subtract—processes the bits of $N$ from most significant to least significant. Start with $R = 0$. At each step, bring in the next input bit $N_i$ by shifting the current remainder left by one position, $R' = 2R + N_i$, then compare $R'$ with $M$:

- if $R' \ge M$, subtract $M$ and set the quotient bit to $1$;
- if $R' < M$, keep $R'$ unchanged and set the quotient bit to $0$.

The updated value becomes the remainder for the next step.

N = 217M = 11

00010011Q = 19

11<1011→q₇ =0

1111<1011→q₆ =0

110110<1011→q₅ =0

11011101≥1011→q₄ =1

−1011

0010

0010100101<1011→q₃ =0

0101001010<1011→q₂ =0

1010010100≥1011→q₁ =1

−1011

01001

1001110011≥1011→q₀ =1

−1011

R = 801000

After all bits have been processed, the quotient bits form $Q$, and the final value of $R$ is the remainder.

The important point for the cost analysis is that $R'$ never becomes arbitrarily large. Since $R < M$, we have $R' = 2R + N_i < 2M$, so $R'$ needs at most one bit more than $M$. The numbers involved therefore stay within essentially the same width throughout the algorithm.

Building the Boolean circuit

Every step performs the same computation, so we only need to design one circuit and then repeat it once for each bit of $N$. At each step, the circuit must do two things: compute $R' - M$, and decide whether to keep that result or keep $R'$ instead.

The subtraction can reuse the adder from the addition example. Using two’s complement, $R' - M = R' + \overline{M} + 1$, so we invert every bit of $M$ and set the adder’s carry-in to $1$. Its carry-out gives the quotient bit $q_i$: it is $1$ when the subtraction can be kept, and $0$ when we must keep the original remainder.

We still need to implement this choice:

$$
R = \begin{cases} R' - M, & q_i = 1,\\ R', & q_i = 0.\end{cases}
$$

A circuit cannot skip the subtraction when $q_i = 0$; it computes $R' - M$ in every case and then uses $q_i$ to choose which result to keep. For each bit, a small multiplexer selects between the two candidate results, $R_j = \big(q_i \land (R' - M)_j\big) \lor \big(\overline{q_i} \land R'_j\big)$. When $q_i = 1$ the first term passes the subtraction result through; when $q_i = 0$ the second passes the original $R'$ through. The same selection circuit is applied independently to every bit.

One-bit multiplexer

The multiplexer uses a constant number of gates per bit, so for $n$-bit numbers it contributes $O(n)$ gates. Together with the $O(n)$-gate subtractor, one division step therefore still uses only $O(n)$ gates.

0 1 1 0 1

Division step

subtract R′ − M

select on qᵢ

0 0 0 1 0

This one step is the whole algorithm. We repeat the same circuit once per bit of $N$, passing the remainder from one step to the next and collecting the quotient bits as they are produced—just as the adder was built by repeating a full-adder stage.

n-bit divider

step: bit 7

step: bit 6

step: bit 5

step: bit 0

Count the gates

The algorithm performs one division step for each of the $n$ bits of $N$. Each step processes $n$ bits and uses $O(n)$ gates: $O(n)$ for the subtraction and $O(n)$ for the bit-by-bit selection. Repeating this step $n$ times gives

$$
t(n)=n\cdot\big(\underbrace{O(n)}_{\text{subtract}}+\underbrace{O(n)}_{\text{select}}\big)=O(n^2).
$$

So there is a family $\{C_1, C_2, \ldots\}$ of Boolean circuits, where $C_n$ divides one $n$-bit nonnegative integer by another and returns both the quotient and the remainder, with $\mathrm{size}(C_n) = O(n^2)$. Counting the two widths separately, as in the figures above, an $n$-bit dividend and an $m$-bit divisor give $n$ steps of $O(m)$ gates, hence circuits of size $O(nm)$.

A faster algorithm

Schoolbook division has the same $O(n^2)$ cost as schoolbook multiplication. But division does not fundamentally require repeated subtraction: it can be reduced to multiplication. The key is the reciprocal of the divisor. Since $N/M = N \cdot (1/M)$, we can divide by $M$ by first computing $1/M$, then multiplying by $N$; a final correction recovers the exact quotient and remainder.

To compute $1/M$ efficiently, we use [Newton’s method](https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division), $x_{k+1} = x_k(2 - M x_k)$, whose approximation $x_k$ to $1/M$ roughly doubles its number of correct bits each iteration. Since the precision grows as the approximation improves, the resulting costs form a geometric series, $\mathrm{M}(n) + \mathrm{M}(n/2) + \mathrm{M}(n/4) + \cdots = O(\mathrm{M}(n))$, where $\mathrm{M}(n)$ is the cost of multiplying two $n$-bit integers.

Thus division can be performed in $O(\mathrm{M}(n))$ bit operations: **asymptotically, division costs no more than multiplication**. Any fast multiplication algorithm therefore gives a fast division algorithm— [Schönhage–Strassen](https://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm) multiplication brings division to $O(n \lg n \lg\lg n)$, and the more recent [Harvey–van der Hoeven](https://annals.math.princeton.edu/2021/193-2/p04) algorithm improves it to $O(n \lg n)$.

Cost analysis: greatest common divisor

Given two nonnegative integers $a$ and $b$, their greatest common divisor $\gcd(a, b)$ is the largest integer that divides both. The classic method for finding the gcd is the [Euclidean algorithm](https://web3-lab.annaburd.me/math-behind-key-pairs/#euclidean-algorithm). It repeatedly replaces $(a, b) \longrightarrow (b,\, a \bmod b)$ until the remainder becomes $0$. The last nonzero remainder is the gcd.

Two things determine the cost: **the cost of one division** and **the number of divisions**. Each Euclidean step computes a remainder $a \bmod b$, using the division circuit from the previous section, and a division of two $n$-bit numbers costs $O(n^2)$ gates. A simple analysis would say that the algorithm takes $O(n)$ divisions, giving $O(n) \cdot O(n^2) = O(n^3)$. It does indeed take only $O(n)$ steps: every two steps, the current remainder is at most half the value from two steps earlier, so the numbers lose at least one bit every two steps. But $O(n^3)$ is too loose. **Not every division is an $n$-bit division.** As the numbers get smaller, later divisions become cheaper, so we need to account for the size of each quotient.

Suppose the $i$-th division has quotient $q_i$, and let $d_i$ be the number of bits in that quotient. Schoolbook division performs one compare-and-subtract operation for each quotient bit, and each such operation costs $O(n)$ gates, so the $i$-th division costs $O(d_i\, n)$ and the whole run costs $\sum_i O(d_i\, n) = O(n \sum_i d_i)$. The key question is therefore: **how large can the total number of quotient bits $\sum_i d_i$ be?**

Bounding the total quotient size

Let the sequence of values produced by the Euclidean algorithm be $a_0 = a$, $a_1 = b$, $a_2, \ldots, a_k$, where the $i$-th division is $a_{i-1} = q_i a_i + a_{i+1}$ and $a_{i+1}$ is the remainder. Since the remainder is nonnegative, the right-hand side is at least $q_i a_i$ on its own, so $a_{i-1} \ge q_i a_i$: **each step shrinks the current value by at least a factor of $q_i$**. Applying this inequality repeatedly gives

$$
a_0 \ge q_1 a_1 \ge q_1 q_2\, a_2 \ge \cdots \ge q_1 q_2 \cdots q_k\, a_k.
$$

The last nonzero value $a_k$ is the gcd, so $a_k \ge 1$ and the right-hand side is at least the product of the quotients alone: $q_1 q_2 \cdots q_k \le a_0$. This is the crucial bound— **although there may be many divisions, their quotients cannot all be large, because their product is limited by the original input.** Taking logarithms converts that product into a sum:

$$
\sum_i \log_2 q_i = \log_2 (q_1 q_2 \cdots q_k) \le \log_2 a_0 < n,
$$

the last step because $a_0$ has $n$ bits, so $a_0 < 2^n$. Now relate this to the actual number of quotient bits. A number with $d_i$ bits sits between the two neighbouring powers of two, $2^{d_i - 1} \le q_i < 2^{d_i}$, and taking logarithms of the left inequality gives $d_i - 1 \le \log_2 q_i$, which means $d_i \le \log_2 q_i + 1$. The $+1$ is the rounding up to a whole number of bits, and each division pays it once.

Summing over the divisions, $\sum_i d_i \le \sum_i \log_2 q_i + \sum_i 1$. The first sum is less than $n$, and there are only $O(n)$ Euclidean divisions, so the second contributes another $O(n)$, giving $\sum_i d_i = O(n)$. So although the algorithm may perform $O(n)$ divisions, **the total number of quotient bits across all those divisions is only $O(n)$**. The total cost is therefore

$$
\sum_i O(d_i\, n) = O\!\Big(n \sum_i d_i\Big) = \boxed{O(n^2)}.
$$

Thus, with schoolbook division, the Euclidean algorithm costs $O(n^2)$ gates. This is the same asymptotic cost as a single $n$-bit multiplication or division using schoolbook arithmetic. With faster multiplication and division algorithms, a recursive version of the Euclidean algorithm can be implemented in $O(\mathrm{M}(n)\log n)$ bit operations, where $\mathrm{M}(n)$ is the cost of multiplying two $n$-bit integers—a count of bit operations rather than of circuit gates.

Cost analysis: modular exponentiation

Modular exponentiation computes $a^b \bmod N$ for nonnegative integers $a$, $b$ and a modulus $N \ge 2$, each at most $n$ bits long. It is the core operation of [RSA](https://web3-lab.annaburd.me/how-cryptography-works/#asymmetric-cryptography-rsa) and [Diffie–Hellman](https://web3-lab.annaburd.me/how-cryptography-works/#asymmetric-cryptography-diffie-hellman). The question here is not just how to compute it, but how the cost grows with $n$.

Multiplying by $a$ repeatedly is inefficient for two reasons. It takes about $b$ multiplications, and an $n$-bit exponent can be as large as $2^n - 1$, making the number of multiplications exponential in the input length. It also constructs the full integer $a^b$, whose intermediate values can grow exponentially large — even though the final result modulo $N$ is smaller than $N$. An efficient algorithm must avoid both problems.

The efficient method is [square-and-multiply](https://web3-lab.annaburd.me/math-behind-key-pairs/#square-and-multiply). Instead of multiplying by $a$ once for every unit in $b$, we use the binary representation of $b$ to build the required power by repeated squaring. Writing the exponent as $b = \sum_{k=0}^{n-1} b_k 2^k$ with $b_k \in \{0, 1\}$ gives

$$
a^b = a^{\sum_k b_k 2^k} = \prod_{k\,:\,b_k=1} a^{2^k}.
$$

The powers $a,\; a^2,\; a^4,\; a^8,\; \ldots$ are each obtained by squaring the previous one, so generating all $n$ of them takes $n-1$ squarings. We then multiply together only the powers corresponding to the $1$-bits of $b$, requiring at most another $n-1$ multiplications.

Crucially, we reduce modulo $N$ after every multiplication. Every intermediate value therefore stays below $N$, so we never construct the enormous integer $a^b$. The entire computation uses at most $2n-2$ modular multiplications, giving a total of **$O(n)$ modular multiplications**.

Cost of one modular multiplication

A modular multiplication takes two values, multiplies them, and then reduces the product modulo $N$. Because $(x\,y) \bmod N = ((x \bmod N)\,(y \bmod N)) \bmod N$, we can reduce after each multiplication and never need to store values larger than $N - 1$. Since $N$ is represented using at most $n$ bits, each value involved in a modular multiplication has at most $n$ bits.

Multiplying two $n$-bit values produces a product of at most $2n$ bits. From the multiplication circuit above, this costs $O(n^2)$ gates. We then reduce the $2n$-bit product modulo $N$, and the division circuit also costs $O(n^2)$ gates. Therefore one modular multiplication costs $O(n^2) + O(n^2) = O(n^2)$ gates.

Total cost

Square-and-multiply uses $O(n)$ modular multiplications, and each modular multiplication costs $O(n^2)$ gates. Therefore **$t(n) = O(n) \cdot O(n^2) = O(n^3)$**.

Thus there is a family of Boolean circuits $C_1, C_2, \ldots$ such that $C_n$ computes $a^b \bmod N$ for inputs of at most $n$ bits, with $\mathrm{size}(C_n) = O(n^3)$. In other words, modular exponentiation can be implemented by a family of polynomial-size Boolean circuits.

The $O(n^3)$ bound uses the basic $O(n^2)$ multiplication and division circuits described above. Replacing them with faster algorithms improves the bit-operation cost to $O(n\,\mathrm{M}(n))$, where $\mathrm{M}(n)$ is the cost of multiplying two $n$-bit integers.

Cost analysis: integer factorization

Given an integer $N\ge2$, [integer factorization](https://web3-lab.annaburd.me/math-behind-key-pairs/#primes) finds its prime factorization: the unique representation $N=p_1^{e_1}p_2^{e_2}\cdots p_k^{e_k}$, where the $p_i$ are distinct primes and the $e_i$ are positive integers. Here $N$ is given in binary using $n$ bits, and we ask the same question as in the previous blocks: how does the work required to recover the answer grow with the input length $n$?

For the arithmetic problems considered so far, we could construct Boolean circuits whose size grows polynomially with $n$: $O(n)$, $O(n^2)$, or $O(n^3)$, depending on the operation. But for factorization **no polynomial-size circuit family is known** so far—perhaps there is one, but we do not know.

The practical difficulty is illustrated by the [RSA Factoring Challenge](https://en.wikipedia.org/wiki/RSA_Factoring_Challenge). One of its targets, [RSA-1024](https://en.wikipedia.org/wiki/RSA_numbers#RSA-1024), was a 1,024-bit number with a US$100,000 prize. The challenge ended in 2007 with RSA-1024 still unfactored, and the remaining prizes were withdrawn. The largest RSA challenge number that has been factored is [RSA-250](https://en.wikipedia.org/wiki/RSA_numbers#RSA-250), an 829-bit number factored in February 2020 using the general number field sieve. The computation required roughly 2,700 CPU core-years.

Trial division

The simplest approach is trial division: test possible divisors one at a time. If $N=uv$ is composite and both $u$ and $v$ were greater than $\sqrt N$, then their product would be greater than $N$, which is impossible. So every composite $N$ has at least one factor $u\le\sqrt N$. We therefore only need to test prime candidates $d\le\sqrt N$. For each candidate, compute $N\bmod d$: a remainder of $0$ means that $d$ is a factor. Divide it out and repeat the process on the quotient until the remaining factor is prime.

$N =$$\sqrt{899}\approx30.0$

The worst case for trial division is an input with no small factor. Since an $n$-bit input satisfies $N<2^n$, we have $\sqrt N<2^{n/2}$. Trial division may therefore need to test up to $O(2^{n/2})$ candidate divisors. Each test uses the $O(n^2)$-gate division circuit built above, giving

$$
t(n)=\underbrace{O(2^{n/2})}_{\text{candidate divisors}}\cdot\underbrace{O(n^2)}_{\text{division test}}=\boxed{O(n^2 2^{n/2})}.
$$

Testing only prime candidates reduces the number of tests: the number of primes below $2^{n/2}$ is about $\tfrac{2^{n/2}}{(n/2)\ln 2}$. So restricting the search to primes saves roughly a factor of $n$, but the exponential term $2^{n/2}$ remains. Trial division is therefore still exponential in the input length.

A congruence of squares

Trial division looks for a factor directly: try $2$, then $3$, then $5$, and so on. For a hard case—a large integer with no unusually small factor—general-purpose factoring methods take a different approach. Instead of searching for a divisor, they construct a relation from which a divisor can be extracted.

That relation is a congruence of squares. Suppose we find two numbers $x$ and $y$ such that $x^2\equiv y^2\pmod N$. In other words, $x^2$ and $y^2$ leave the same remainder when divided by $N$. Therefore, $N\mid(x^2-y^2)=(x-y)(x+y)$.

So $N$ divides the product of $x-y$ and $x+y$. If $N$ is composite, its factors can be distributed between these two terms, and we can often recover one of them by computing $\gcd(x-y,N)$.

For example, take $N =$

$$
503^2\equiv168^2\equiv218\pmod{737}
$$

different as ordinary integers, equal after reduction

$$
503^2-168^2=(503-168)(503+168)=335\cdot671
$$

so the difference is a multiple of 737

$$
\gcd(335,737)=67
$$

and the gcd with 737 pulls one factor out of the product

$$
737=67\cdot11
$$

a factor, found without dividing by anything

The gcd itself is cheap: the above costs $O(n^2)$ gates. The difficult part is finding the pair $x,y$ in the first place.

The quadratic sieve

We want to find $x^2\equiv y^2\pmod N$. The quadratic sieve approaches this indirectly. Instead of trying to find $y$ directly, it looks for many values of $x^2\bmod N$ that factor completely into small primes. These are called smooth values. Once enough smooth values have been collected, we can combine them so that their product becomes a perfect square. That gives us the second square.

The list of small primes is called the factor base, and every value that factors completely over it is kept as a relation.

The table's parity column records whether the exponent of each factor-base prime is odd or even. With factor base $\{2,3,5,7\}$, for example, $224=2^5\cdot7$ has odd exponents for $2$ and $7$, and even exponents for $3$ and $5$, so its parity vector is $(1,0,0,1)$. Parity is all we keep, because a number is a perfect square exactly when every exponent in its factorization is even — whether an exponent is $2$ or $6$ makes no difference to that question.

Parity rows add the way values multiply. Multiplying two values adds their exponents, so it adds their parity bits mod 2, one prime at a time: $(1,0,0,1)+(1,0,0,1)=(0,0,0,0)$. A prime used an odd number of times in each of the two values is used an even number of times in their product, so the two $1$s cancel. A set of rows adding to all zeros is therefore a set whose values multiply to a perfect square, and finding such a set is the only thing the parity column is for.

Both halves of the congruence come out of that one set. Multiplying the chosen values of $x$ gives the left-hand root; squaring it replaces each one by its residue from the table, so $x^2$ is congruent to the product of those residues — the product just shown to be a square. The right-hand root $y$ is that square's root, and nothing has to search for it: halving every exponent in the factorization writes it down directly, which is possible only because the parities were all even. Reduce $y$ mod $N$ and $x^2\equiv y^2\pmod N$ is in hand, with the gcds left to finish.

$N =$$\sqrt{737}\approx27.1$

factor base$\{2,3,5,7\}$

| $x$ | $x^2\bmod N$ | factors into | parity |
| --- | --- | --- | --- |

| 28 | 47 | $47$ | × |
| --- | --- | --- | --- |
| 29 | 104 | $2^{3}\cdot13$ | × |
| 30 | 163 | $163$ | × |
| 31 | 224 | $2^{5}\cdot7$ | click (1,0,0,1) |
| 32 | 287 | $7\cdot41$ | × |
| 33 | 352 | $2^{5}\cdot11$ | × |
| 34 | 419 | $419$ | × |
| 35 | 488 | $2^{3}\cdot61$ | × |
| 36 | 559 | $559$ | × |
| 37 | 632 | $2^{3}\cdot79$ | × |
| 38 | 707 | $7\cdot101$ | × |
| 39 | 47 | $47$ | × |
| 40 | 126 | $2\cdot3^{2}\cdot7$ | click (1,0,0,1) |
| 41 | 207 | $3^{2}\cdot23$ | × |
| 42 | 290 | $2\cdot5\cdot29$ | × |
| 43 | 375 | $3\cdot5^{3}$ | click (0,1,1,0) |
| 44 | 462 | $2\cdot3\cdot7\cdot11$ | × |
| 45 | 551 | $551$ | × |
| 46 | 642 | $2\cdot3\cdot107$ | × |
| 47 | 735 | $3\cdot5\cdot7^{2}$ | click (0,1,1,0) |
| 48 | 93 | $3\cdot31$ | × |
| 49 | 190 | $2\cdot5\cdot19$ | × |
| 50 | 289 | $289$ | × |
| 51 | 390 | $2\cdot3\cdot5\cdot13$ | × |
| 52 | 493 | $493$ | × |
| 53 | 598 | $2\cdot299$ | × |
| 54 | 705 | $3\cdot5\cdot47$ | × |
| 55 | 77 | $7\cdot11$ | × |
| 56 | 188 | $2^{2}\cdot47$ | × |
| 57 | 301 | $7\cdot43$ | × |
| 58 | 416 | $2^{5}\cdot13$ | × |
| 59 | 533 | $533$ | × |
| 60 | 652 | $2^{2}\cdot163$ | × |
| 61 | 36 | $2^{2}\cdot3^{2}$ | click (0,0,0,0) |
|  |  |  |  |

The demo keeps the numbers small enough to show the bookkeeping: candidate values, smooth relations, parity rows, and the final gcds. A real quadratic sieve uses the same logic at a scale where hand-picking rows is impossible. It locates smooth values with an actual sieve—the operation the algorithm is named for—and then uses linear algebra over $\mathbb F_2$ to find a set of relation rows whose parity sum is zero.

Note the tradeoff: a larger factor base makes smooth values easier to find, since more primes can divide them. But it also makes each parity row wider, increasing the size of the linear system and the number of relations needed before a dependency is guaranteed. The running time therefore depends on choosing a factor-base size that balances the cost of finding relations against the cost of solving the resulting linear system.

The general number field sieve

The [general number field sieve](https://www.ams.org/notices/199612/pomerance.pdf) (GNFS) improves on the quadratic sieve by changing how its smooth relations are constructed. In the quadratic sieve, we look for smooth values among $x^2\bmod N$, which are roughly as large as $N$. As $N$ grows, smooth values become increasingly rare.

GNFS takes a different route: instead of searching for smooth values of roughly size $N$, it constructs smaller values and looks for pairs that are smooth. That makes smooth relations much easier to find, and is the key reason GNFS can handle much larger integers.

The search begins by choosing a polynomial $f$ of degree $d$ with a root $m$ modulo $N$, so that $f(m)\equiv0\pmod N$. Instead of testing single values of $x$, GNFS tests coprime pairs of integers $(a,b)$: it searches over a finite range, sweeping the integer grid within it and keeping the pairs for which $\gcd(a,b)=1$. From each pair, it constructs two integers: $a-bm$ on the ordinary integer side, and $b^{d}f(a/b)$ on the number-field side. These are the two values we test for smoothness. If both factor completely over their respective factor bases, the pair gives a smooth relation and is kept.

Choosing $f$ may look like the difficult part, but constructing a suitable polynomial is surprisingly straightforward. Fix a degree $d$, choose $m=\lfloor N^{1/d}\rfloor$, and write $N$ in base $m$. Use those base-$m$ digits as the coefficients of $f$. Then, by construction, $f(m)=N$, so in particular $f(m)\equiv0\pmod N$. That is exactly the property GNFS needs.

$$
N =
$$

degree3

factor base$\{2,3,5,7\}$

$$
m=\lfloor \sqrt[3]{737}\rfloor=9
$$

The only choices are the degree and the factor base. The degree determines m, and m determines the coefficients.

$$
737=\textcolor{#4338ca}{1}\cdot9^{3}+\textcolor{#0369a1}{8}
$$

$$
f(x)=\textcolor{#4338ca}{1}\cdot x^{3}+\textcolor{#0369a1}{8}
$$

Polynomial form of the line above, with x in place of m.

$$
f(9)=737\equiv0\pmod{737}
$$

This is why we constructed f this way: m = 9 is a root of f modulo N, which links the two sides of the GNFS construction.

shares a factor, repeats a smaller pairnot yet reachedtested, a side left a prime outside the baseboth values factor over the base: a relation

One thing to notice is that GNFS is not necessarily faster than the quadratic sieve on small numbers. It does more work per relation, but that extra cost is offset by its better asymptotic scaling as $N$ grows. Only for sufficiently large $N$ does GNFS become the faster method.

The cost of GNFS

Almost all of the work in GNFS goes into two jobs. The first is collecting relations: sweep the grid of candidate pairs $(a,b)$, test the two values each pair produces, and keep the pairs where both factor completely over the factor bases. The second is the linear algebra: take the relations that survived and find a set of parity rows summing to zero, which is what turns a pile of relations into a congruence of squares. Neither job can be skipped — the first produces the raw material and the second extracts an answer from it — so the running time is the sum of the two.

Both jobs are governed by one number: the smoothness bound $B$, the largest prime allowed in the factor bases. $B$ is the knob we can turn, and it pulls the two jobs in opposite directions. Turn it up and each value has more primes available to factor into, so relations become easier to find; but the factor bases grow with it, and every extra prime is another column in the matrix the second job has to solve.

How many primes is that? There are about $B/\ln B$ primes below $B$. The logarithm moves that count by far less than the choice of $B$ itself does, so from here on we drop it and speak of about $B$ factor-base primes.

Each surviving relation becomes one row of the relation matrix, recording which factor-base primes occur an odd number of times, and each factor-base prime is one column. A set of rows summing to zero is guaranteed once there are more rows than columns, so with about $B$ columns GNFS needs about $B$ relations, plus a small surplus so that the dependency it finds is a usable one. That is where the first job's target comes from: not as many relations as possible, but about $B$ of them.

What one relation costs depends on how often a candidate turns out to be smooth. Write $V$ for the size of the values being tested and set $u=\frac{\ln V}{\ln B}$. This $u$ measures the tested value against the smoothness bound: it is roughly how many factors of size $B$ it takes to build a number of size $V$. When $u$ is small the value is barely larger than the primes allowed to divide it, and smoothness is common. When $u$ is large the value has to be assembled out of many small primes at once, which is rare.

For a random integer of size $V$, the [Dickman function](https://en.wikipedia.org/wiki/Dickman_function) puts the chance of being $B$-smooth at about $u^{-u}$. Turn that probability into work: one success in every $u^{u}$ candidates means about $u^{u}$ candidates tested per relation kept, and about $B\,u^{u}$ candidates tested to collect the $B$ relations we need. This step is heuristic. The values GNFS tests come out of a polynomial and are not random integers, but they behave closely enough to random ones for the estimate to hold up in practice.

The second job works on what the first produced: about $B$ relations, so about $B$ rows, against about $B$ factor-base primes, so about $B$ columns — a matrix roughly $B\times B$. It is a sparse one, since a single relation is divisible by only a handful of factor-base primes and almost every entry in its row is zero. Sparse methods exploit that and cost roughly $B^{2}$, instead of the $B^{3}$ that ordinary elimination would spend on a dense matrix of that size. Real implementations are more delicate; $B^{2}$ is the simplified model we carry through the argument. The two costs together are

$$
t(B)\approx\underbrace{B\,u^{u}}_{\text{finding the relations}}+\underbrace{B^{2}}_{\text{processing them}},\qquad u=\frac{\ln V}{\ln B}.
$$

Everything that follows is an argument about how to choose $B$ in that one equation. Push $B$ down and $B^{2}$ becomes negligible, but $u$ grows and $u^{u}$ grows much faster still: smooth values turn rare and the sieve spends a long time looking for them. Push $B$ up and relations arrive quickly, but the matrix that has to absorb them grows quadratically.

| smoothness bound B | finding relations | linear algebra |
| --- | --- | --- |
| small | expensive — smooth values are rare | cheap — few columns to solve |
| large | cheap — smooth values are common | expensive — many columns to solve |
| balanced | the two costs meet, and the total is as small as it gets |  |

Neither extreme is where the total is smallest. The best $B$ is the one where the two terms are of comparable size, because on either side of that point every saving in one job is paid for by the other.

Locating that point is awkward with the notation we have, because neither term is polynomial in the input length nor exponential in it. Costs in that gap are usually written in L-notation, as $L_N[\alpha,c]$, and only two things about it matter here: the exponent $\alpha$ says where in the gap a cost falls, running from polynomial at $\alpha=0$ to exponential at $\alpha=1$, and the constant $c$ refines the estimate within a scale.

Now measure both costs on that scale. Suppose the values that must be smooth have size $V=L_N[\alpha,\cdot]$, and choose a factor base $B=L_N[\beta,\cdot]$. An $L$-value is an exponential, so taking logarithms leaves $\ln V$ and $\ln B$ as products of powers of $\ln N$ and $\ln\ln N$, and in their ratio the exponents subtract:

$$
u=\frac{\ln V}{\ln B}\approx c\left(\frac{\ln N}{\ln\ln N}\right)^{\alpha-\beta}.
$$

That subtraction drives the rest. Since $\ln u\approx(\alpha-\beta)\ln\ln N$, the search cost $u^{u}=e^{u\ln u}$ has an exponent proportional to $(\ln N)^{\alpha-\beta}(\ln\ln N)^{1-(\alpha-\beta)}$, which is exactly the shape of $L_N[\alpha-\beta,\cdot]$. The leading factor $B$ contributes only $L_N[\beta,\cdot]$, no larger than the search term at the balance point we are heading for, so it leaves the scale alone. On the other side, squaring $B$ doubles a constant but does not touch the exponent, so the linear algebra stays at $L_N[\beta,\cdot]$. The mapping worth remembering is that the size of the numbers being smoothed contributes $\alpha$, the size of the factor base contributes $\beta$, and the search pays the difference between them:

$$
t\approx\underbrace{L_N[\alpha-\beta,\cdot]}_{\text{finding the relations}}+\underbrace{L_N[\beta,\cdot]}_{\text{processing them}}.
$$

A sum of two $L$-terms is set by the larger exponent; the smaller one is swallowed by the slack the notation already carries. So if $\alpha-\beta>\beta$, the cost is relation collection and the matrix was free; if $\beta>\alpha-\beta$, the cost is linear algebra and the relations were free. From either side, moving $\beta$ toward the other case lowers the total, until the two exponents meet:

$$
\alpha-\beta=\beta\quad\Longrightarrow\quad\beta=\frac{\alpha}{2}.
$$

Once $\alpha$ is known, the best factor base is the one that splits that exponent evenly between the two jobs. Nothing in the argument is specific to GNFS: it holds for any method that collects smooth relations and then solves for a dependency among them.

Take the quadratic sieve first. It smooths values of $x^{2}\bmod N$, which are about as large as $N$ itself, and $N=L_N[1,1]$, so $\alpha=1$. Balancing gives $\beta=\frac12$, and a running time of $L_N[\frac12,1]$.

GNFS changes one thing about that calculation: the values it tests are not of size $N$. With the degree chosen as $d\approx\sqrt[3]{\frac{3\ln N}{\ln\ln N}}$, the pair of values $a-bm$ and $b^{d}f(a/b)$ that must both be smooth are heuristically of size $L_N[\frac23,\cdot]$, so $\alpha=\frac23$. The same balancing gives $\beta=\frac13$:

$$
\underbrace{\alpha=1\;\Rightarrow\;\beta=\tfrac12\;\Rightarrow\;L_N[\tfrac12,\cdot]}_{\text{quadratic sieve}}\qquad\underbrace{\alpha=\tfrac23\;\Rightarrow\;\beta=\tfrac13\;\Rightarrow\;L_N[\tfrac13,\cdot]}_{\text{number field sieve}}
$$

The distance between $\frac12$ and $\frac13$ is the whole of the improvement, and it is worth being exact about where it comes from. GNFS does not test smoothness any faster than the quadratic sieve does. It wins because its polynomial construction hands it much smaller numbers to make smooth. That lowers $\alpha$, and every exponent in the analysis is downstream of $\alpha$.

What the balancing argument gives is the shape of the complexity: which power of $\ln N$ appears, and why it is $\frac13$ rather than $\frac12$. The constant in front takes a longer optimization, one that tunes the degree $d$ and the smoothness bound $B$ together rather than one after the other, since $d$ is what sets the size of the values and therefore the rate at which they are smooth. It selects

$$
B=L_N\!\left[\tfrac13,\sqrt[3]{\tfrac89}\right]\quad\Longrightarrow\quad\underbrace{B^{2}}_{\text{linear algebra}}=L_N\!\left[\tfrac13,2\sqrt[3]{\tfrac89}\right]=L_N\!\left[\tfrac13,\sqrt[3]{\tfrac{64}{9}}\right].
$$

Squaring $B$ doubles its constant, and $2\sqrt[3]{8/9}$ is $\sqrt[3]{64/9}$: the familiar constant arrives out of the linear algebra. Relation collection is tuned to cost the same — the balancing argument again, this time with the constants kept — so the total carries that constant too:

$$
t(N)=\underbrace{\exp\!\left(\left(\sqrt[3]{\tfrac{64}{9}}+o(1)\right)(\ln N)^{\tfrac13}(\ln\ln N)^{\tfrac23}\right)}_{\text{heuristic expected time},\ \sqrt[3]{\tfrac{64}{9}}\approx1.923}.
$$

In terms of bits, with $n=\log_2 N$, we have $\ln N=n\ln2$ and $\ln\ln N=\Theta(\log n)$, so the exponent is $\Theta(n^{1/3}(\log n)^{2/3})$ and

$$
t(n)=\boxed{2^{O(n^{\frac13}(\log n)^{\frac23})}}.
$$

Read that exponent both ways. It grows with $n$, so the cost is not polynomial in the input length. But it grows like $n^{1/3}$ rather than like $n$, which leaves it far below the $2^{n/2}$ of trial division. Between the two is what subexponential means, and it is where the best factoring algorithms known today live.

So the fastest known way to factor large integers today is GNFS, with a cost of $2^{O(n^{1/3}(\log n)^{2/3})}$. There might be a faster classical algorithm: no polynomial-time method is known, but factoring has never been shown to be NP-complete either, so as of today this has been neither proved nor disproved. But a faster quantum algorithm is already known—we come to it later.

## Classical circuits as quantum circuits

Classical and quantum computation have been treated separately so far. In practice, quantum algorithms routinely need ordinary classical computation inside a larger quantum circuit: adding two integers, evaluating a function, checking a condition. So can a classical algorithm be run on a quantum computer?

Yes. Any Boolean circuit of size $t$ can be implemented with $O(t)$ quantum gates.

The purpose is compatibility rather than speed. A classical computation run this way is no faster than before, but it now runs coherently: it behaves correctly when its input is part of a superposition, and it leaves its output in a register the rest of the quantum algorithm can use.

Toffoli gates

Classical gates such as AND and OR destroy information. AND takes two input bits and returns one, so its four possible inputs collapse onto two possible outputs: a $0$ at the output could have come from any of $00$, $01$ or $10$, and there is no way to tell which. That is the one obstacle here, because every quantum gate is unitary and therefore reversible. A Boolean circuit has to be rebuilt out of reversible gates before a quantum computer can run it, and the gate that does that work is the Toffoli gate.

Recall that a is a controlled-controlled-NOT: it flips its target qubit only when both control qubits are $1$.

$$
\mathrm{Toffoli}\,\lvert a\rangle\lvert b\rangle\lvert c\rangle=\lvert a\rangle\lvert b\rangle\lvert c\oplus ab\rangle
$$

| in |  | out |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $a$ | $b$ | $c$ | $ab$ | $a$ | $b$ | $c \oplus ab$ |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 1 | 0 | 0 | 0 | 1 |
| 0 | 1 | 0 | 0 | 0 | 1 | 0 |
| 0 | 1 | 1 | 0 | 0 | 1 | 1 |
| 1 | 0 | 0 | 0 | 1 | 0 | 0 |
| 1 | 0 | 1 | 0 | 1 | 0 | 1 |
| 1 | 1 | 0 | 1 | 1 | 1 | 1 |
| 1 | 1 | 1 | 1 | 1 | 1 | 0 |

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert c\rangle
$$

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert c\oplus ab\rangle
$$

Toffoli from elementary gates

The has no three-qubit gate, so one Toffoli gate has to be built from several elementary gates. The circuit below uses 15 elementary gates: two Hadamards, six CNOTs, and seven $T$ or $T^\dagger$ gates. The decomposition is exact — these 15 gates reproduce the Toffoli operation exactly, not approximately.

This illustrates the main cost of translating classical computation into quantum computation. A single classical operation such as AND can correspond to a whole collection of elementary quantum gates. But the number of gates needed is a fixed constant: one Toffoli costs 15 elementary gates, or simply $O(1)$. So a classical circuit with $t$ gates can still be implemented with $O(t)$ elementary quantum gates. The translation introduces overhead, but only a constant-factor overhead.

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert c\rangle
$$

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert c\oplus ab\rangle
$$

Simulating Boolean gates

Now that we can build a Toffoli gate from elementary quantum gates, we can use it to reproduce. The key idea is simple: because quantum gates must be reversible, a Boolean operation is computed into an additional qubit rather than replacing its inputs. With the target initialized to $\lvert 0\rangle$, the quantum circuit can therefore reproduce the same Boolean function on computational-basis states.

$$
\lvert a\rangle
$$

$$
\lvert \neg a\rangle
$$

Nothing to do: NOT is already reversible, and the Pauli-X gate implements exactly the same operation on the two basis states.

FANOUT

$$
\lvert a\rangle
$$

$$
\lvert 0\rangle
$$

$$
\lvert a\rangle
$$

A CNOT with a fresh qubit as its target implements FANOUT on basis states: it copies the input bit to the fresh qubit. This does not violate the no-cloning theorem, because it does not copy an arbitrary quantum state.

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert 0\rangle
$$

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert a\wedge b\rangle
$$

A Toffoli gate with a fresh target qubit computes $ab$ into it: starting from $0$, the target ends as $ab$, which is exactly $a \land b$ for bits $a, b \in \lbrace 0, 1 \rbrace$.

$$
\lvert a\rangle
$$

$$
\lvert b\rangle
$$

$$
\lvert 0\rangle
$$

$$
\lvert \neg a\rangle
$$

$$
\lvert \neg b\rangle
$$

$$
\lvert a\vee b\rangle
$$

By [De Morgan’s law](https://en.wikipedia.org/wiki/De_Morgan%27s_laws) an OR is an AND with everything flipped: flip both inputs, AND them with a Toffoli, then flip the result. The two inputs are left flipped on the way out.

So every Boolean gate can be replaced by $O(1)$ quantum gates, using at most one workspace qubit initialized to $\lvert 0\rangle$. A Boolean circuit with $t$ gates therefore becomes a quantum circuit with $O(t)$ gates and $O(t)$ qubits.

But there is a catch: the quantum version is reversible, so it cannot simply discard the inputs or intermediate values the way a classical circuit does. Instead, they remain in the circuit, leaving behind workspace that must eventually be cleaned up.

Simulating Boolean circuits

Now take a whole circuit rather than one gate. Suppose $C$ is a Boolean circuit of size $t$ computing a function $f: \Sigma^n \to \Sigma^m$:

$$
x
$$

$$
f(x)
$$

Replace each Boolean gate by its quantum simulation, adding a fresh $\lvert 0\rangle$ qubit whenever needed. The resulting quantum circuit $R$ uses $O(t)$ gates and acts on $n + k$ qubits, where the workspace $k = O(t)$. For a basis-state input $x$, the desired $m$-bit output appears in the first $m$ qubits, but the remaining qubits contain leftover intermediate values:

$$
R\bigl(\lvert x\rangle\lvert 0^k\rangle\bigr)=\lvert f(x)\rangle\lvert g(x)\rangle.
$$

Here $g(x)$ is the garbage produced by making the computation reversible.

$$
\lvert x\rangle
$$

$$
\lvert 0^k\rangle
$$

$$
\lvert f(x)\rangle
$$

$$
\lvert g(x)\rangle
$$

Clearing the garbage

The garbage is more than wasted space: if it remains entangled with the result, it can interfere with the quantum algorithm and spoil the interference patterns we rely on. The simple solution is to uncompute it. Because $R$ is made entirely of reversible quantum gates, we can run it backwards using its inverse $R^\dagger$, at the same $O(t)$ cost. This lets us compute the result, use it where needed, and then erase the unwanted intermediate values without erasing the result itself.

The key is that $R$ is deterministic: once we have computed $f(x)$, we can copy that classical result before undoing the computation. We therefore add a fresh $m$-qubit register $\lvert y\rangle$, initially $\lvert 0^m\rangle$, and use $m$ CNOTs to copy the answer into it between $R$ and $R^\dagger$. This is the same FANOUT trick as before: the result wires hold basis-state bits, so copying them does not violate the no-cloning theorem. Then $R^\dagger$ erases the workspace while leaving the copied result untouched.

$$
\lvert x\rangle
$$

$$
\lvert 0^k\rangle
$$

$$
\lvert y\rangle
$$

$$
\lvert x\rangle
$$

$$
\lvert 0^k\rangle
$$

$$
\lvert y\oplus f(x)\rangle
$$

Constructing the query gate

Combine the three circuit segments — the computation of $f(x)$, the XOR of $f(x)$ into the target register, and the uncomputation of the workspace — and call the resulting circuit $Q$. Its cost is

$$
O(t)+m+O(t)=O(t),
$$

since the computation and uncomputation each cost $O(t)$, while the $m$-gate target update is absorbed into $O(t)$.

More importantly, the workspace register is returned to $\lvert 0^k\rangle$ after the uncomputation. Thus the complete circuit acts as

$$
\lvert x\rangle\lvert 0^k\rangle\lvert y\rangle\;\longmapsto\;\lvert x\rangle\lvert 0^k\rangle\lvert y\oplus f(x)\rangle.
$$

$$
\lvert x\rangle
$$

$$
\lvert 0^k\rangle
$$

$$
\lvert y\rangle
$$

$$
\lvert x\rangle
$$

$$
\lvert 0^k\rangle
$$

$$
\lvert y\oplus f(x)\rangle
$$

Because the workspace starts and ends in the fixed state $\lvert 0^k\rangle$, it can be ignored when describing the action of the circuit on the input and target registers. The remaining transformation is exactly the quantum query gate $U_f$ for the function computed by the original Boolean circuit.

In other words, the query-model oracle does not have to be treated as an abstract black box: given a classical circuit for $f$, we can construct its quantum query gate using $O(t)$ gates, only a constant-factor overhead compared with the original circuit.

## Phase estimation and factoring

A quantum state can sometimes pick up a phase when a unitary operation is applied to it. That phase is, but it contains useful information about the operation. Phase estimation is a procedure for extracting that hidden phase.

The spectral theorem

A useful way to understand a matrix is to look for directions that it does not mix with other directions. These are its eigenvectors: if $M\lvert\psi\rangle = \lambda\lvert\psi\rangle$, then applying $M$ to $\lvert\psi\rangle$ does not turn $\lvert\psi\rangle$ into a different direction, it only multiplies it by the number $\lambda$.

$$
\lambda=3
$$

$$
\lambda=2
$$

knocked off its span

$$
\begin{pmatrix}3&1\\0&2\end{pmatrix}
$$

Apply the matrix100%

Try a direction58°

- Eigenvector $(1,0)$ with eigenvalue $\lambda=3$. It keeps its own line.
- Eigenvector $(-1,1)$ with eigenvalue $\lambda=2$. It keeps its own line.
- Any other direction is not an eigenvector: drag the slider and it leaves its dashed line.

For a general matrix, there may not be enough eigenvectors to form a basis. And even when there are enough, they need not be perpendicular to one another. Either way, they are not necessarily convenient as coordinates for the whole space.

The spectral theorem identifies a class of matrices whose eigenvectors can be chosen to form an. This gives us a particularly useful coordinate system: the matrix acts on each direction independently, multiplying it by that direction’s own eigenvalue.

Eigenvalue λ₁1.45

Eigenvalue λ₂0.55

- $\lvert\psi_1\rangle$ and $\lvert\psi_2\rangle$ stay on their own lines. Their eigenvalues only change their lengths.
- Any other vector has components along both eigendirections. Since those components are stretched by different amounts, the vector changes direction as well as length.
- The dashed circle represents all unit vectors. Under $M$, these vectors map to the solid ellipse, whose axes lie along the two eigendirections.

The eigenvalues here are real, so they stretch or shrink the eigenvector directions. A unitary matrix preserves lengths, so its eigenvalues have magnitude 1: in the complex plane, they rotate each direction by a phase instead of changing its length. That phase is what this chapter is after — and it is the one part a real two-dimensional picture cannot show.

The spectral decomposition

A matrix $M$ is normal when it commutes with its:

$$
MM^\dagger=M^\dagger M.
$$

The spectral theorem says that every normal $N \times N$ matrix has an orthonormal basis of eigenvectors $\{\lvert\psi_1\rangle, \ldots, \lvert\psi_N\rangle\}$, together with phases, with corresponding complex eigenvalues $\lambda_1, \ldots, \lambda_N$, such that

$$
M=\sum_{k=1}^{N}\lambda_k\lvert\psi_k\rangle\langle\psi_k\rvert.
$$

Each basis vector satisfies

$$
M\lvert\psi_k\rangle=\lambda_k\lvert\psi_k\rangle.
$$

Writing a matrix in this form is called its spectral decomposition. It says that the entire matrix is determined by an orthonormal set of directions and one complex number for each direction, specifying what $M$ does along it.

Special case: unitary matrices

A unitary matrix satisfies $U^\dagger U = I = UU^\dagger$, so it is normal and the spectral theorem applies. What unitarity adds is a constraint on the eigenvalues. A unitary operation preserves norms, so if $U\lvert\psi_k\rangle = \lambda_k\lvert\psi_k\rangle$, then the output must have the same length as the input. This forces $|\lambda_k| = 1$.

A complex number of modulus one does not change a vector’s length. It only contributes a phase: a rotation in the complex plane. Every such number can be written as $e^{2\pi i\theta}$ for exactly one $\theta \in [0, 1)$.

So suppose $U$ is an $N \times N$ unitary matrix. There exists an orthonormal basis $\{\lvert\psi_1\rangle, \ldots, \lvert\psi_N\rangle\}$, together with phases

$$
\lambda_1=e^{2\pi i\theta_1},\ldots,\lambda_N=e^{2\pi i\theta_N},
$$

such that

$$
U=\sum_{k=1}^{N}\lambda_k\lvert\psi_k\rangle\langle\psi_k\rvert.
$$

Each vector $\lvert\psi_k\rangle$ is an eigenvector of $U$ with eigenvalue $\lambda_k$:

$$
U\lvert\psi_k\rangle=\lambda_k\lvert\psi_k\rangle=e^{2\pi i\theta_k}\lvert\psi_k\rangle.
$$

For a unitary matrix, the spectral decomposition therefore reduces the action of the entire matrix to a collection of phases. Each eigenvector defines an independent direction, and along that direction the matrix does nothing more than multiply by $e^{2\pi i\theta_k}$. The magnitude is fixed at one, so the only information left in each eigenvalue is its phase $\theta_k$.

The phase estimation problem

In the phase estimation problem, we are given two things:

1. A description of a quantum circuit on $n$ qubits implementing a unitary operation $U$.
2. An $n$-qubit quantum state $\lvert\psi\rangle$.

We are promised that $\lvert\psi\rangle$ is an eigenvector of $U$. By the spectral theorem, its eigenvalue has the form $e^{2\pi i\theta}$ for a unique $\theta \in [0, 1)$. The goal is to approximate this phase $\theta$, where

$$
U\lvert\psi\rangle=e^{2\pi i\theta}\lvert\psi\rangle.
$$

The important point is that the eigenvector is given as a quantum state, not as a classical description. We cannot simply read $\theta$ from the circuit, nor can we measure $\lvert\psi\rangle$ to reveal which eigenvector it is. The phase must be extracted by interacting with the state through controlled applications of $U$.

The phase estimate

The phase $\theta$ is a real number, but a quantum measurement can return only finitely many classical bits. We therefore choose a precision $m$: the algorithm will return $m$ bits that specify one of $2^m$ possible approximations to $\theta$.

For example, with $m = 3$, the possible answers are the eight equally spaced points $0, \tfrac{1}{8}, \tfrac{2}{8}, \ldots, \tfrac{7}{8}$.

If the true phase is $\theta = 0.310$, the closest grid point is $\tfrac{2}{8} = 0.250$, so the three-bit answer is $010$, representing the approximation $0.250$. In general, the answer has the form $\theta \approx \tfrac{y}{2^m}$ for $y \in \{0, 1, \ldots, 2^m - 1\}$, and the binary representation of $y$ is the $m$-bit output.

There is one important detail: these points lie on a circle, not on a line. The phases $0$ and $1$ represent the same point, because $e^{2\pi i\cdot 0} = e^{2\pi i\cdot 1} = 1$. So the approximation is understood modulo one. A phase close to $1$ can therefore be approximated by a value close to $0$ when the shortest distance around the circle crosses the boundary.

$$
1
$$

$$
i
$$

$$
-1
$$

$$
-i
$$

$$
2\pi\theta
$$

Phase θ0.310

Precision m3 bits

Angle $2\pi\theta$111.6°

Grid points$2^{3} = 8$

Nearest estimate$\tfrac{2}{8} = 0.250$

Phase error0.060

Angular error21.6°

Phase kickback: making the phase observable

Applying $U$ to $\lvert\psi\rangle$ multiplies the state by $e^{2\pi i\theta}$ and changes nothing else, so measuring the resulting state cannot reveal $\theta$. Phase kickback turns this invisible phase into an observable relative phase: instead of applying $U$ directly, we apply it conditionally on an extra qubit, transferring the phase $e^{2\pi i\theta}$ to the control qubit.

Creating an observable phase

A controlled-$U$ uses an extra qubit to decide whether $U$ is applied: one branch does nothing, while the other applies $U$ to the register. If the control is in a definite state $\lvert 0\rangle$ or $\lvert 1\rangle$ this does not help, because only one branch ever exists and the phase remains global.

The key is to put the control into a superposition. Both branches are then present at once: one where $U$ is applied and one where it is not. Since $\lvert\psi\rangle$ is an eigenvector, it picks up $e^{2\pi i\theta}$ and nothing else, and only in the branch where $U$ acts, so the phase becomes a relative phase between the two branches. A second Hadamard makes those branches interfere, converting the relative phase into measurement probabilities on the control qubit. The register itself is never measured. Everything we learn about $\theta$ comes from the control.

$$
\lvert 0\rangle
$$

$$
\lvert\psi\rangle
$$

Step through the circuit

1. 1Prepare the register in $\lvert\psi\rangle$ and the control qubit in $\lvert 0\rangle$.
2. 2Apply a Hadamard to the control qubit.
3. 3Apply controlled-$U$.
4. 4Apply a Hadamard to the control qubit again.
5. 5Measure the control qubit. The register is never measured.

$$
p_0
$$

$$
p_1
$$

$$
\theta
$$

Phase θ0.310

$p_0=\cos^2(\pi\theta)$0.316

$p_1=\sin^2(\pi\theta)$0.684

What can we learn from one measurement?

The measurement does tell us something about $\theta$: the probabilities change as the phase changes. For example, phases near $0$ tend to produce $\lvert 0\rangle$, while phases near $\tfrac{1}{2}$ tend to produce $\lvert 1\rangle$.

But this is not enough to determine the phase. The same measurement statistics can arise from different phases: $\theta$ and $1-\theta$ are indistinguishable. The probabilities also change very little near $0$ and $\tfrac{1}{2}$, so this measurement gives poor precision there.

So one controlled-$U$ lets us learn something about the phase, but not enough to identify it. To estimate $\theta$ accurately, we need a way to make the measurement more sensitive to different parts of the phase.

Running controlled-U twice

The first experiment was not sensitive enough to distinguish all phases. A natural idea is therefore to apply $U$ more than once. If one application gives the phase $\theta$, then two applications give twice the phase:

$$
U^{2}\lvert\psi\rangle=e^{2\pi i(2\theta)}\lvert\psi\rangle.
$$

So if we put two controlled-$U$ gates on the same control qubit, we get the same experiment as before, but with the phase doubled. This changes how the measurement probabilities respond to $\theta$, giving us information that the single-$U$ experiment could not provide.

$$
\lvert 0\rangle
$$

$$
\lvert\psi\rangle
$$

$$
p_0
$$

$$
p_1
$$

$$
\theta
$$

Phase θ0.310

$p_0=\cos^2(2\pi\theta)$0.136

$p_1=\sin^2(2\pi\theta)$0.864

More sensitivity, more ambiguity

Doubling the phase makes the probabilities change twice as quickly as $\theta$ changes. Phases that were hard to distinguish before can now produce noticeably different probabilities, so the measurement becomes more sensitive to the phase.

But the doubled phase is still read modulo one. In particular, $2\theta$ and $2\theta+1$ represent the same phase, so $\theta$ and $\theta+\tfrac{1}{2}$ produce identical statistics. The original reflection symmetry, $\theta\leftrightarrow 1-\theta$, remains as well. We have therefore gained sensitivity, but also introduced more possible phases that give the same measurement statistics.

This is the central tension in phase estimation: using more applications of $U$ gives finer information about the phase, but also creates more ambiguity about which phase produced it. The solution will be to use several powers of $U$ together, so that the ambiguities from one measurement are resolved by the others.

What do we gain by using both experiments?

We now have two experiments with complementary strengths. One application of $U$ covers the whole range of $\theta$, but resolves it coarsely. Two applications make the probabilities change twice as quickly, but introduce additional ambiguities. It is natural to ask whether the information from the two experiments can be put together to get a better estimate.

The register itself is not the obstacle. Because $\lvert\psi\rangle$ is an eigenvector, each experiment leaves it unchanged and separates it from the control qubit. Measuring the control therefore does not disturb $\lvert\psi\rangle$, so the experiment can be repeated with the same state.

The difficulty is that measurement throws away most of the information available before measurement. Just before measurement, the control qubit has amplitudes whose relative phase depends on $\theta$. Measurement turns those amplitudes into a single classical bit, $0$ or $1$. To learn the corresponding probabilities accurately, we need many repetitions.

So if we run the $U$ and $U^{2}$ experiments separately, we end up with two collections of classical measurement results. We can estimate two probabilities and try to use them together, but each estimate is noisy and each experiment has its own ambiguities.

This raises the next question: can we arrange the experiments so that their phase information is combined before measurement, rather than after?

Two control qubits

Rather than running the two experiments one after another, we can give each of them its own control qubit and run them in a single circuit. The upper control drives one application of $U$, and the lower control drives two.

one U

two U

$$
\textcolor{#6d28d9}{\lvert 0\rangle}
$$

$$
\textcolor{#b45309}{\lvert 0\rangle}
$$

$$
\lvert\psi\rangle
$$

$$
a_0
$$

$$
a_1
$$

Step through the circuit

1. 1Prepare the register in $\lvert\psi\rangle$ and both control qubits in $\lvert 0\rangle$.
2. 2Apply a Hadamard to each control qubit.
3. 3Apply controlled-$U$ once, controlled by $a_0$.
4. 4Apply controlled-$U$ twice, both controlled by $a_1$.

Can we distinguish the phases?

The two controls now carry the control factor of $\lvert\pi_3\rangle$: $\frac{1}{2}\sum\limits_{x=0}^{3}e^{2\pi ix\theta}\lvert x\rangle$.

In general, $\theta$ need not be restricted to a few special values. But to make the problem concrete, let us first pretend that we are promised $\theta=\frac{y}{4}$ for some $y\in\{0,1,2,3\}$. This gives us a smaller problem: can we work out which of these four possible values of $\theta$ we have?

Each possibility gives a different two-qubit state: $\lvert\phi_{y}\rangle=\frac{1}{2}\sum\limits_{x=0}^{3}e^{2\pi i\frac{xy}{4}}\lvert x\rangle$. Explicitly,

$$
\lvert\phi_{0}\rangle=\frac{1}{2}\lvert 0\rangle+\frac{1}{2}\lvert 1\rangle+\frac{1}{2}\lvert 2\rangle+\frac{1}{2}\lvert 3\rangle
$$

$$
\lvert\phi_{1}\rangle=\frac{1}{2}\lvert 0\rangle+\frac{i}{2}\lvert 1\rangle-\frac{1}{2}\lvert 2\rangle-\frac{i}{2}\lvert 3\rangle
$$

$$
\lvert\phi_{2}\rangle=\frac{1}{2}\lvert 0\rangle-\frac{1}{2}\lvert 1\rangle+\frac{1}{2}\lvert 2\rangle-\frac{1}{2}\lvert 3\rangle
$$

$$
\lvert\phi_{3}\rangle=\frac{1}{2}\lvert 0\rangle-\frac{i}{2}\lvert 1\rangle-\frac{1}{2}\lvert 2\rangle+\frac{i}{2}\lvert 3\rangle
$$

Our goal is now clear: determine which of the four states $\lvert\phi_{0}\rangle,\ldots,\lvert\phi_{3}\rangle$ the controls are in. If we can identify the state, we immediately know $y$, and therefore the original phase $\theta=\frac{y}{4}$. And conveniently, notice that all four states are, so they can be distinguished perfectly by a: $\{\lvert\phi_{0}\rangle\langle\phi_{0}\rvert,\ \lvert\phi_{1}\rangle\langle\phi_{1}\rvert,\ \lvert\phi_{2}\rangle\langle\phi_{2}\rvert,\ \lvert\phi_{3}\rangle\langle\phi_{3}\rvert\}$.

Knowing that the four states can be distinguished does not yet give us a way to read out which one we have. We need to change the basis back to the computational basis. Let $V$ be the unitary whose columns are $\lvert\phi_{0}\rangle$, $\lvert\phi_{1}\rangle$, $\lvert\phi_{2}\rangle$, and $\lvert\phi_{3}\rangle$. By construction, $V\lvert y\rangle=\lvert\phi_{y}\rangle$ for every $y\in\{0,1,2,3\}$. In this case,

$$
V=\frac{1}{2}\begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix}
$$

This matrix is the in four dimensions. As a quantum operation, it is called the quantum Fourier transform, or $\mathrm{QFT}_4$.

Now apply the inverse transformation. It takes each of our four states back to the corresponding computational-basis state: $V^\dagger\lvert\phi_{y}\rangle=\lvert y\rangle$.

So instead of building a special measurement for the four $\lvert\phi_{y}\rangle$ states, we can simply apply $V^\dagger$ and then measure the qubits in the computational basis. The measurement gives us $y$, and therefore the phase $\theta=\frac{y}{4}$.

$$
\textcolor{#6d28d9}{\lvert 0\rangle}
$$

$$
\textcolor{#b45309}{\lvert 0\rangle}
$$

$$
\lvert\psi\rangle
$$

At the four promised phases, each curve reaches exactly $1$ at its own quarter and $0$ at the others, so the measurement is certain. Between those phases, the peaks spread out: the outcome is no longer certain, but the nearest quarter remains the most likely.

$$
y=0
$$

$$
y=1
$$

$$
y=2
$$

$$
y=3
$$

$$
\theta
$$

Phase θ0.310

$\Pr[\,y=0\,]$0.043

$\Pr[\,y=1\,]$0.834

$\Pr[\,y=2\,]$0.093

$\Pr[\,y=3\,]$0.030

The quantum Fourier transform

The key idea is to build states whose amplitudes all have the same magnitude but differ in phase.

For example, suppose there are four computational-basis states, labelled $x=0,1,2,3$. The phase can stay constant, or advance by a quarter, half, or three quarters of a full turn each time $x$ increases.

Complex phase

| $x$$y$rows are the frequency y, columns the position x | 0 | 1 | 2 | 3 |
| --- | --- | --- | --- | --- |
| 0 |  |  |  |  |
| 1 |  |  |  |  |
| 2 |  |  |  |  |
| 3 |  |  |  |  |

Frequency $y$1

Position $x$1

$$
e^{2\pi i\cdot\frac{(1)(1)}{4}}=i
$$

As $x$ increases, the phase can advance at different rates. Each rate produces a different pattern, corresponding to a different discrete frequency.

The quantum Fourier transform is the change of basis from the computational-basis states to these frequency patterns. It is the quantum counterpart of the, with the normalization factor $\tfrac{1}{\sqrt{N}}$ that makes the frequency patterns orthonormal and the transformation unitary.

For a positive integer $N$, the quantum Fourier transform $\mathrm{QFT}_N$ is the $N\times N$ unitary defined by

$$
\mathrm{QFT}_N=\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}e^{2\pi i\frac{xy}{N}}\lvert x\rangle\langle y\rvert
$$

Equivalently, its action on a computational-basis state is

$$
\mathrm{QFT}_N\lvert y\rangle=\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}e^{2\pi i\frac{xy}{N}}\lvert x\rangle
$$

The second form is often easier to read. Start with the basis state $\lvert y\rangle$. The transform produces a superposition of all the output basis states $\lvert x\rangle$. Every output basis state has the same amplitude magnitude, $\tfrac{1}{\sqrt{N}}$. What changes with $x$ is the phase $e^{2\pi ixy/N}$.

The phase factor is determined by the product $xy$. For a fixed input $y$, increasing $x$ makes the phase advance in equal steps, and the value of $y$ determines how large those steps are. For example, $y=0$ gives no phase change. $y=1$ advances by one step around the circle — a quarter-turn in the four-state example above. $y=2$ advances twice as far at each step, and so on. Each input basis state $\lvert y\rangle$ is therefore mapped to a different phase pattern.

For an $n$-qubit register, $N=2^{n}$, because that is the number of computational-basis states available. The definition itself does not require $N$ to be a power of two — that restriction comes from applying the transform to a register of whole qubits.

Examples at different sizes

Since $e^{2\pi i\cdot N/N}=1$, only $xy \bmod N$ matters. So, no matter how large $N$ becomes, the entries use only $N$ distinct phases, $e^{2\pi ik/N}$ for $k=0,\ldots,N-1$. Let’s look at a few examples, starting with the smallest transform.

Size $N$1

There is one basis state and one phase: $1$.

$$
\mathrm{QFT}_{1}=\begin{pmatrix}1\end{pmatrix}
$$

Shorthand notation for phase

The same phases keep appearing in every transform. Instead of writing the exponential each time, name the first phase: $\omega_N=e^{2\pi i/N}$. Then every phase is a power of it: $\omega_N^{k}=e^{2\pi ik/N},\qquad\omega_N^{N}=1$.

On the unit circle, $\omega_N$ is one step of $2\pi/N$. Its powers take successive steps around the circle: $1,\,\omega_N,\,\omega_N^{2},\,\ldots,\,\omega_N^{N}=1$. The $N$ distinct powers are the.

A column of the transform follows the same walk. Fixing $y$, its exponents are $0,\,y,\,2y,\,3y,\ldots$, so each row advances by $y$ steps around the circle.

Powers of ω

$$
\mathrm{QFT}_{4}=\frac{1}{2}\begin{pmatrix}1&\textcolor{#0284c7}{1}&1&1\\1&\textcolor{#0284c7}{\omega}&\omega^{2}&\omega^{3}\\1&\textcolor{#0284c7}{\omega^{2}}&1&\omega^{2}\\1&\textcolor{#0284c7}{\omega^{3}}&\omega^{2}&\omega\end{pmatrix}
$$

$$
\omega=\omega_{4}=e^{2\pi i/4}
$$

Size $N$4

Column $y$1

Row $x$1

$$
\omega^{(1)(1)}
$$

Where the arrow lands is a pair of coordinates, written down by Euler’s formula: $\omega_N=e^{2\pi i/N}=\cos\left(\tfrac{2\pi}{N}\right)+i\sin\left(\tfrac{2\pi}{N}\right)$.

So naming $\omega_N$ collapses the definition to a sum of its powers, and the matrix to a table of them:

$$
\begin{aligned}\mathrm{QFT}_N&=\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}\omega_N^{xy}\lvert x\rangle\langle y\rvert\\[6pt]\mathrm{QFT}_N\lvert y\rangle&=\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\omega_N^{xy}\lvert x\rangle\\[12pt]\mathrm{QFT}_N&=\frac{1}{\sqrt{N}}\begin{pmatrix}1&1&1&\cdots&1\\1&\omega_N&\omega_N^{2}&\cdots&\omega_N^{N-1}\\1&\omega_N^{2}&\omega_N^{4}&\cdots&\omega_N^{2(N-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&\omega_N^{N-1}&\omega_N^{2(N-1)}&\cdots&\omega_N^{(N-1)^{2}}\end{pmatrix}\end{aligned}
$$

Turning phase back into a number

Undoing the transform conjugates every phase, so the inverse is the same matrix with the sign of the exponent reversed:

$$
(\mathrm{QFT}_N^\dagger)_{x,y}=\frac{1}{\sqrt{N}}\,\omega_N^{-xy}=\frac{1}{\sqrt{N}}e^{-2\pi ixy/N}
$$

This is the direction used in phase estimation. The controlled-$U$ gates leave the control register in one of the phase patterns above — a Fourier-basis state, not a computational-basis state. That is why measuring the controls directly tells us so little.

$\mathrm{QFT}_N^\dagger$ maps that phase pattern back to the computational basis: $\mathrm{QFT}_N^\dagger\,\mathrm{QFT}_N\lvert y\rangle=\lvert y\rangle$. After that, an ordinary measurement reveals $y$.

Circuits for the QFT

When $N=2^{n}$, the QFT acts on $n$ qubits. Its phase pattern has a simple, repeating structure that we can use: each qubit contributes one level of the pattern, with smaller phase rotations appearing as we move along the qubits. This lets us build the QFT efficiently as a ladder of single-qubit gates and controlled phase rotations, rather than treating every basis state separately.

For a computational-basis input $\lvert y\rangle$, the output can be written as a tensor product of $n$ single-qubit states:

$$
\mathrm{QFT}_{2^{n}}\lvert y\rangle=\bigotimes_{j=1}^{n}\frac{1}{\sqrt{2}}\left(\lvert 0\rangle+e^{2\pi iy/2^{\,n+1-j}}\lvert 1\rangle\right)
$$

The tensor-product symbol $\otimes$ means that we combine these single-qubit states into the full $n$-qubit state. Every factor has the same form, $\tfrac{1}{\sqrt{2}}\left(\lvert 0\rangle+e^{i\varphi}\lvert 1\rangle\right)$, where $\varphi$ is the phase for that qubit. The $\lvert 0\rangle$ term carries no explicit phase because $1=e^{i0}$, so it is the phase reference. The $\lvert 1\rangle$ term carries the relative phase $e^{i\varphi}$.

So each output qubit is an equal superposition of $\lvert 0\rangle$ and $\lvert 1\rangle$, with a phase that depends on $y$ and on which qubit we are looking at. The phases differ by powers of two, giving the QFT its characteristic phase pattern.

Building blocks

The output qubits are not entangled with one another, so we can build the state one qubit at a time.

The Hadamard gate creates the equal superposition $\tfrac{1}{\sqrt{2}}\left(\lvert 0\rangle+\lvert 1\rangle\right)$, which is the basic form of each single-qubit factor above.

A controlled-phase gate adds a phase only to the $\lvert 11\rangle$ state:

$$
\alpha
$$

$$
\mathrm{CP}(\alpha)=\begin{pmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&e^{i\alpha}\end{pmatrix}
$$

The gate is symmetric: it does not matter which qubit is considered the control and which is the target. Both qubits simply need to be $\lvert 1\rangle$ for the phase to be applied. This is why its circuit symbol has two identical dots rather than a separate control and target.

The circuit pattern

The circuit is built from one short pattern repeated across the wires. Each wire gets a Hadamard followed by controlled-phase gates connecting it to the wires below. The phase angles decrease by powers of two: the largest angle, $\pi/2$, connects the wire being worked on to the bottom wire, then $\pi/4$, $\pi/8$, and so on as the connections move upward.

The resulting phase factors appear on the output wires in reverse order. The final swaps reverse the wire order and put them back into the intended positions.

In the picture, the part of the circuit not yet drawn out is folded into a single $\mathrm{QFT}$ box on the left. Unfolding that box reveals another copy of the same pattern.

Qubits $n$5 (N = 32)

Unfoldings1

Cost analysis

Let $s_n$ denote the number of gates we need for $n$ qubits. For $n=1$, a single Hadamard gate is required. For $n\ge 2$, these are the gates required:

- $s_{n-1}$ gates for the QFT on $n-1$ qubits
- $n-1$ controlled-phase gates
- $n-1$ swap gates
- 1 Hadamard gate

$$
s_n=\begin{cases}1 & n=1\\[2pt] s_{n-1}+2n-1 & n\ge 2\end{cases}
$$

This is a recurrence relation with a:

$$
s_n=\sum_{k=1}^{n}(2k-1)=n^{2}
$$

So cost is $n^{2}$ gates for a transform on $N=2^{n}$ amplitudes — quadratic in the number of qubits, for a matrix with $N^{2}$ entries in it.

The swap gates can be reduced. Taken together, they simply reverse the order of the wires, so we need only $\lfloor n/2\rfloor$ swaps if we perform that reversal directly. We can also omit them entirely if we are willing to relabel the wires.

The QFT can also be approximated with fewer gates and lower depth. Its phase angles shrink geometrically: $\pi/2$, $\pi/4$, $\pi/8$, and so on. Once the rotations become small enough, dropping them has little effect while reducing the cost of the circuit.

The inverse QFT

Phase estimation runs this circuit backwards. Reversing the order of the gates and changing every phase angle $\alpha$ to $-\alpha$ gives $\mathrm{QFT}_N^\dagger$ at the same cost.

This is the circuit that turns the phase pattern left behind by the controlled-$U$ gates back into the computational basis, where a measurement can read the encoded number.

Phase estimation with $m$ control qubits

The two-control circuit generalises to $m$ control qubits without changing its basic shape. Each control applies a different power of $U$, so the control register accumulates a phase pattern determined by $\theta$. With $m$ controls, this pattern contains $m$ bits of phase information. It has exactly the form produced by $\mathrm{QFT}_{2^{m}}$ from the corresponding basis state, so we apply $\mathrm{QFT}_{2^{m}}^\dagger$ and measure the controls to recover those bits.

$$
\textcolor{#ffffff}{U}
$$

$$
\textcolor{#ffffff}{U^{2}}
$$

$$
\textcolor{#ffffff}{U^{2^{m-1}}}
$$

$$
\textcolor{#ffffff}{\mathrm{QFT}^{\dagger}_{2^{m}}}
$$

$$
\textcolor{#6d28d9}{\lvert 0^{m}\rangle}
$$

$$
\lvert\psi\rangle
$$

The eigenstate $\lvert\psi\rangle$ is unchanged by every controlled power of $U$, because $\lvert\psi\rangle$ is an eigenvector of $U$ and therefore of every power of it. All the phase information is stored in the control register. Just before measurement, the full state is

$$
\displaystyle\lvert\pi\rangle=\lvert\psi\rangle\otimes\frac{1}{2^{m}}\sum_{y=0}^{2^{m}-1}\sum_{x=0}^{2^{m}-1}e^{2\pi ix(\theta-y/2^{m})}\lvert y\rangle
$$

So the probability of reading $y$ is

$$
\displaystyle p_y=\left\lvert\frac{1}{2^{m}}\sum_{x=0}^{2^{m}-1}e^{2\pi ix(\theta-y/2^{m})}\right\rvert^{2}
$$

Accuracy of a single run

The probability $p_y$ of measuring the control register in $\lvert y\rangle$ depends only on the distance between $\theta$ and the corresponding grid point $y/2^{m}$. If $\theta=y/2^{m}$, every term in the sum is $1$, so $p_y=1$. Otherwise the terms do not line up perfectly, and $p_y$ is smaller.

The possible estimates $y/2^{m}$ are spaced by $2^{-m}$. The nearest grid point is therefore at most half a step from $\theta$, $\lvert\theta-y/2^{m}\rvert\le 2^{-(m+1)}$. For phase differences this small, the probability formula above gives $p_y\ge 4/\pi^{2}\approx 0.405$.

Conversely, if a grid point is at least one full step from $\theta$, $\lvert\theta-y/2^{m}\rvert\ge 2^{-m}$, the same probability formula gives $p_y\le 1/4$.

Thus the nearest grid point has at least a $40.5\%$ chance of appearing in one run, while any grid point at least one full step away has probability at most $25\%$.

$$
\tfrac{4}{\pi^{2}}
$$

$$
\tfrac{1}{4}
$$

$$
\theta
$$

Control qubits $m$3 (2^m = 8)

Phase θ0.310

$\text{nearest }y$2

$\lvert\theta-y/2^{m}\rvert$0.0600

$2^{-(m+1)}$0.0625

$\Pr[\,y\,]$0.443

A single run therefore favors the best approximation but does not guarantee it. Repeating the procedure and taking the mode of the outcomes makes that approximation increasingly likely. The eigenvector $\lvert\psi\rangle$ is unchanged, so it can be reused for every run.

Alternative phase-estimation methods

Standard phase estimation estimates the phase $\theta$ using controlled applications of $U$. Other approaches use different combinations of quantum resources and classical processing:

- Iterative phase estimation extracts the phase bits one at a time, reusing a single control qubit instead of keeping $m$ control qubits at once.
- Kitaev’s phase estimation uses a single control qubit and estimates the phase from interference measurements involving different powers of $U$.
- Maximum-likelihood and Bayesian methods repeat controlled-$U^{k}$ experiments and use classical statistical inference to estimate $\theta$.

These approaches trade off the same basic resources: control qubits, applications of $U$, and classical post-processing.

The order-finding problem: using phase estimation

When working [modulo](https://web3-lab.annaburd.me/math-behind-key-pairs/#modular-foundations) $N$, we only need $N$ possible values, represented by the integers from $0$ to $N-1$. We denote this set by $\mathbb{Z}_N=\{0,1,\ldots,N-1\}$. Thus $\mathbb{Z}_1=\{0\}$, $\mathbb{Z}_2=\{0,1\}$, $\mathbb{Z}_3=\{0,1,2\}$, and so on.

The elements $a \in \mathbb{Z}_N$ that satisfy $\gcd(a, N) = 1$ have an important property: they have a multiplicative inverse modulo $N$. We collect all of them into the set $\mathbb{Z}_N^{*}=\{a\in\mathbb{Z}_N:\gcd(a,N)=1\}$. For $N = 21$, for example, twelve of the twenty-one elements are invertible: $\mathbb{Z}_{21}^{*}=\{1,2,4,5,8,10,11,13,16,17,19,20\}$.

The connection with the greatest common divisor follows from the [Euclidean algorithm](https://web3-lab.annaburd.me/math-behind-key-pairs/#euclidean-algorithm). If $\gcd(a, N) = 1$, it gives integers $x$ and $y$ such that $ax + Ny = 1$. Reducing modulo $N$ gives $ax = 1$, so $x$ is a multiplicative inverse of $a$. Conversely, if $a$ has an inverse modulo $N$, then $\gcd(a, N)$ must be $1$.

Now take any $a \in \mathbb{Z}_N^{*}$ and repeatedly multiply by $a$, producing the powers $a,\ a^{2},\ a^{3},\ldots$ Every one of them is invertible too: if $x$ is the inverse of $a$, then $x^{k}$ is the inverse of $a^{k}$. So all the powers lie in $\mathbb{Z}_N^{*}$, and that set is finite, so they cannot all be different. Two of them must be equal: $a^{i} \equiv a^{j} \pmod{N}$ for some $i < j$. Multiplying both sides by the inverse of $a^{i}$ cancels it and leaves $a^{j-i} \equiv 1 \pmod{N}$, where $j - i$ is positive. So some positive power of $a$ returns to $1$.

So, the smallest positive exponent $r$ for which $a^{r} \equiv 1 \pmod{N}$ is called the order of $a$ in $\mathbb{Z}_N^{*}$.

For elements outside $\mathbb{Z}_N^{*}$ no such exponent exists: if $d = \gcd(a, N) > 1$, then $d$ divides both $a^{r}$ and $N$, so $a^{r} \equiv 1 \pmod{N}$ would force $d$ to divide $1$, which is impossible.

Modulus N21

$a$2

$\gcd(a, N)$1

Size of $\mathbb{Z}_N^{*}$12

Powers of $2$ modulo $21$

2$2^{1}$

4$2^{2}$

8$2^{3}$

16$2^{4}$

11$2^{5}$

1$2^{6}$

back to$2^{1}$

$$
r = 6
$$

The problem

We are given two positive integers $a$ and $N$, with the promise that $\gcd(a, N) = 1$. The task is to find the order of $a$: the smallest positive integer $r$ such that $a^{r} \equiv 1 \pmod{N}$. The two numbers $a$ and $N$ are all we are given. In particular, no factorization of $N$ is provided.

Both numbers are written in binary, so the input length is $n = O(\log N)$ bits. Computing a single power $a^{k} \bmod N$ is efficient: does it using $O(n^3)$ gates. The difficulty is that the order can be almost as large as $N$. Checking the powers one at a time can therefore require $\Omega(N)$ steps, which is exponential in the input length $n$.

The table below runs that scan for $a = 2$. Each modulus is about ten times the one above it, and so is the time.

| $N$ | order $r$ | time |
| --- | --- | --- |
| 9,610,721 | — | — |
| 40,670,489 | — | — |
| 207,335,717 | — | — |
| 4,043,918,803 | — | — |

Scan to measure multiplication speed and estimate the cost at different sizes.

No efficient classical algorithm for order-finding is known. This is significant because order-finding is closely related to integer factorization. In fact, an efficient order-finding algorithm can be used to efficiently, so factorization can be reduced to order-finding.

Multiplication as a unitary operation

We know what we want to find: the length $r$ of the cycle that repeated multiplication by $a$ modulo $N$ runs through. The idea is to turn that repeated multiplication into an operation a quantum computer can apply to a state. For a given element $a \in \mathbb{Z}_N^{*}$, define the operation as $M_{a}\lvert x\rangle=\lvert ax \bmod N\rangle$ for each $x \in \mathbb{Z}_N$.

Because $a$ has a multiplicative inverse modulo $N$, multiplication by $a$ is a bijection on $\mathbb{Z}_N$: every state has exactly one image, and every state has exactly one preimage. In other words, multiplication by $a$ simply permutes the elements of $\mathbb{Z}_N$.

A permutation of the computational basis states is represented by a unitary matrix. This is why $M_{a}$ is a valid quantum operation.

If $d = \gcd(a, N) > 1$, this breaks down. Every product $ax \bmod N$ is divisible by $d$, so the map can reach only a subset of the states. Multiple inputs therefore collide at the same output, while other states are never reached. The map is no longer a permutation, and its matrix is not unitary.

$M_a\text{ on }\mathbb{Z}_8$$a=$

$$
\gcd(a,8)=1
$$

input$\lvert x\rangle$output$\lvert ax\bmod 8\rangle$

$$
M_{3}\lvert 1\rangle=\lvert 3\rangle
$$

input

output

0123456701·······1···1····2······1·3·1······4····1···5·······16··1·····7·····1··

A permutation can be decomposed into cycles: starting from any state, repeatedly applying $M_{a}$ eventually returns to that state. For multiplication by $a$, these cycles are determined by the repeated powers of $a$ modulo $N$.

For example, take $N = 8$ and $a = 3$. The state $\lvert 0\rangle$ remains fixed, while starting from $\lvert 1\rangle$, repeated application of $M_{3}$ gives $\lvert 1\rangle\to\lvert 3\rangle\to\lvert 1\rangle$. The cycle therefore has length $2$. Equivalently, $3^{2} \equiv 1 \pmod{8}$, and no smaller positive power gives $1$, so the order of $3$ modulo $8$ is $r = 2$.

The remaining states form cycles of their own: $\lvert 2\rangle\to\lvert 6\rangle\to\lvert 2\rangle$ and $\lvert 5\rangle\to\lvert 7\rangle\to\lvert 5\rangle$, while $\lvert 4\rangle$ is fixed. Together, these cycles make up the full permutation implemented by $M_{3}$.

This is the key connection: the order we want is encoded as the length of a cycle in the permutation $M_{a}$. The remaining challenge is to extract that cycle length from the unitary using quantum phase estimation.

From the cycle to eigenvalues

At this point the order $r$ is hidden as the number of positions in a cycle. Phase estimation does not measure that cycle length directly. It measures an eigenphase, so the goal is to encode the cycle length $r$ into an eigenphase of the form $j/r$.

The cycle containing $\lvert 1\rangle$ consists of the states $\lvert 1\rangle, \lvert a\rangle, \ldots, \lvert a^{r-1}\rangle$. On this part of the state space, $M_{a}$ has one simple action: move everything one position forward, wrapping the last position back to the first:

$$
\lvert 1\rangle\xrightarrow{M_{a}}\lvert a\rangle\xrightarrow{M_{a}}\lvert a^{2}\rangle\xrightarrow{M_{a}}\cdots\xrightarrow{M_{a}}\lvert a^{r-1}\rangle\xrightarrow{M_{a}}\lvert 1\rangle
$$

A basis state does not have the property we need. For example, $M_{a}\lvert 1\rangle=\lvert a\rangle$, so applying $M_{a}$ changes it into a different basis state. Instead, consider a superposition of the states in the cycle. With the right pattern of phases, the shift preserves this superposition and changes only its overall phase. Such a state is an eigenvector, and the corresponding phase change is its eigenvalue.

Begin with $\lvert\psi_0\rangle$, the equal superposition of all positions in the cycle, with every amplitude having the same phase:

$$
\lvert\psi_{0}\rangle=\frac{1}{\sqrt{r}}\left(\lvert 1\rangle+\lvert a\rangle+\cdots+\lvert a^{r-1}\rangle\right)
$$

Applying $M_{a}$ moves every term one position forward. The last state wraps back to $\lvert 1\rangle$, so the same $r$ terms appear again, only in a different order. The state is therefore unchanged: its eigenvalue is $1$, corresponding to eigenphase $\theta_0=0$. This is a valid eigenvector, but its phase contains no information about $r$.

$$
\begin{aligned}
M_{a}\lvert\psi_{0}\rangle&=\frac{1}{\sqrt{r}}\left(\lvert a\rangle+\lvert a^{2}\rangle+\cdots+\lvert a^{r}\rangle\right)\\[4pt]
&=\frac{1}{\sqrt{r}}\left(\lvert a\rangle+\cdots+\lvert a^{r-1}\rangle+\lvert 1\rangle\right)=\lvert\psi_{0}\rangle
\end{aligned}
$$

We need eigenvectors with nonzero eigenphases. The simplest way to get one is to let the amplitudes acquire a phase difference from one position to the next. Because the cycle contains $r$ positions, this phase difference must fit consistently when the cycle closes: after $r$ steps, the phase must return to its starting value. A natural choice is therefore $1/r$ of a full turn per step. Writing this phase step as $\omega_{r}=e^{2\pi i/r}$, we have $\omega_r^{r}=1$.

Now look at $\lvert\psi_1\rangle$. We assign successive positions phases that differ by $1/r$ of a turn, so position $k$ carries the factor $\omega_r^{-k}$ (the minus sign is a convention):

$$
\lvert\psi_{1}\rangle=\frac{1}{\sqrt{r}}\left(\lvert 1\rangle+\omega_{r}^{-1}\lvert a\rangle+\cdots+\omega_{r}^{-(r-1)}\lvert a^{r-1}\rangle\right)
$$

Applying $M_{a}$ shifts every position forward by one step. The phase pattern shifts with the states, and when the last term wraps back to $\lvert 1\rangle$, its phase factor becomes $\omega_r^{-(r-1)}=\omega_r$. Rearranging the terms shows that every amplitude has acquired the same extra factor $\omega_r$. The phase pattern is therefore unchanged, while the whole state gains the eigenphase $\theta_1=1/r$.

$$
\begin{aligned}
M_{a}\lvert\psi_{1}\rangle&=\frac{1}{\sqrt{r}}\left(\lvert a\rangle+\omega_{r}^{-1}\lvert a^{2}\rangle+\cdots+\omega_{r}^{-(r-1)}\lvert a^{r}\rangle\right)\\[4pt]
&=\frac{1}{\sqrt{r}}\left(\omega_{r}\lvert 1\rangle+\lvert a\rangle+\omega_{r}^{-1}\lvert a^{2}\rangle+\cdots+\omega_{r}^{-(r-2)}\lvert a^{r-1}\rangle\right)\\[4pt]
&=\omega_{r}\cdot\frac{1}{\sqrt{r}}\left(\lvert 1\rangle+\omega_{r}^{-1}\lvert a\rangle+\omega_{r}^{-2}\lvert a^{2}\rangle+\cdots+\omega_{r}^{-(r-1)}\lvert a^{r-1}\rangle\right)\\[4pt]
&=\omega_{r}\lvert\psi_{1}\rangle
\end{aligned}
$$

By the same logic, we can choose different phase steps to obtain a whole family of eigenvectors. The state $\lvert\psi_j\rangle$ is an equal superposition of all $r$ basis states in the cycle through $\lvert 1\rangle$, with only their phases differing. Each component has magnitude $1/\sqrt r$. The label $j$ determines the phase difference between neighbouring positions: the phase advances by $j/r$ of a turn from one position to the next. Thus, the component on $\lvert a^k\rangle$ carries the phase factor $\omega_r^{-jk}=e^{-2\pi i jk/r}$:

$$
\lvert\psi_{j}\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\omega_{r}^{-jk}\lvert a^{k}\rangle\qquad\text{for }j\in\{0,1,\ldots,r-1\}
$$

Every state in this family is an eigenvector of $M_{a}$, with eigenvalue $\omega_r^{\,j}=e^{2\pi i j/r}$.

$$
M_a\lvert\psi_j\rangle=\omega_r^{\,j}\lvert\psi_j\rangle
$$

Note that there are different ways to choose the phase pattern and construct eigenvectors. For this problem, however, these particular eigenvectors are useful because their eigenphases are $\theta_j=j/r$, so the unknown cycle length $r$ appears directly in the denominator.

The phase pattern of an eigenstate

cycle length $r =$

phase step $j =$

$$
\lvert 1\rangle
$$

$$
\lvert a\rangle
$$

$$
\lvert a^{2}\rangle
$$

$$
\lvert a^{3}\rangle
$$

$$
\lvert a^{4}\rangle
$$

0°

288°

216°

144°

72°

$$
\lvert\psi_{1}\rangle=\frac{1}{\sqrt{5}}\left(\lvert 1\rangle+\omega_{5}^{-1}\lvert a\rangle+\omega_{5}^{-2}\lvert a^{2}\rangle+\omega_{5}^{-3}\lvert a^{3}\rangle+\omega_{5}^{-4}\lvert a^{4}\rangle\right)
$$

$$
\theta_{1}=\frac{j}{r}=\frac{1}{5}\text{ turn}
$$

The phase pattern is what makes these states useful for phase estimation. Under every controlled power of $M_{a}$, an eigenvector remains the same target state while its phase accumulates in the control register.

From eigenphase to order

Among the eigenvectors $\lvert\psi_{j}\rangle$ above, start with $j = 0$. Its eigenphase is $0$, which carries no information about the unknown order $r$. The next choice, $j = 1$, is exactly what we need: its eigenphase is $1/r$, putting the unknown order directly in the denominator:

$$
M_{a}\lvert\psi_{1}\rangle=\omega_{r}\lvert\psi_{1}\rangle=e^{2\pi i\frac{1}{r}}\lvert\psi_{1}\rangle
$$

This gives us a direct route from phase estimation to the order. If we can prepare $\lvert\psi_{1}\rangle$, phase estimation gives an estimate of its eigenphase, which in this case is $1/r$. We can then invert that estimate to obtain $r$.

1. Perform phase estimation on $\lvert\psi_{1}\rangle$ using a quantum circuit implementing $M_{a}$, with $m$ control qubits. The controlled powers of $M_{a}$ accumulate the phase $e^{2\pi i\frac{k}{r}}$ in the control register. The inverse QFT converts this accumulated phase into an estimate of the eigenphase. Measuring the control register gives an integer $y$. Dividing by $2^{m}$ turns that $m$-bit readout into a phase estimate $y/2^{m}$ in $[0,1)$. And since the eigenphase is $1/r$, $y/2^{m}\approx 1/r$.
2. Recover the order by inverting the phase estimate and rounding it to the nearest integer: $r\approx\operatorname{round}\!\left(\frac{2^{m}}{y}\right)=\left\lfloor\frac{2^{m}}{y}+\frac{1}{2}\right\rfloor$.

How accurate does the phase estimate need to be?

The estimate of $1/r$ must be accurate enough to distinguish it from the phase corresponding to any other possible order $r'$. Since both $r$ and $r'$ are smaller than $N$, the smallest possible separation between two such phases is $\lvert 1/r-1/r'\rvert=\lvert r'-r\rvert/(rr')>1/N^{2}$.

Therefore, if the phase estimate is within half of this minimum separation from the true phase, it cannot be mistaken for the phase of a different possible order. In other words, it is enough to have $\lvert y/2^{m}-1/r\rvert\le 1/(2N^{2})$.

With $m$ control qubits, the phase-estimation grid has spacing $1/2^{m}$, so the nearest grid point is at most $1/2^{m+1}$ away from the true phase. Choosing $m = 2\lceil\lg N\rceil + 1$ makes this error at most $1/(4N^{2})$, comfortably within the required precision. Thus $O(\log N)$ control qubits are enough.

A single run produces the nearest grid point with probability at least $4/\pi^{2}$, about 40%. Repeating the procedure independently increases the probability of obtaining the correct phase. After $k$ runs, the probability that at least one run produces the nearest grid point is at least $1-(1-4/\pi^{2})^{k}$, so a constant number of repetitions gives any fixed desired success probability, while $O(\log(1/\varepsilon))$ repetitions give failure probability at most $\varepsilon$.

So, we choose enough qubits so that the useful region around the true phase is narrow enough to identify $r$. And adding more qubits makes the grid finer and the phase estimate more precise, while repetitions can further boost the probability of obtaining a sufficiently accurate estimate.

When the eigenphase is a random fraction

The previous procedure assumed that we could start with $\lvert\psi_{1}\rangle$, whose eigenphase is $1/r$. But there is nothing special about $j = 1$: suppose instead that we are given $\lvert\psi_{j}\rangle$ for a random choice of $j \in \{0,\ldots,r-1\}$. Its eigenphase is $j/r$, so phase estimation now returns that fraction rather than $1/r$:

$$
M_{a}\lvert\psi_{j}\rangle=\omega_{r}^{\,j}\lvert\psi_{j}\rangle=e^{2\pi i\frac{j}{r}}\lvert\psi_{j}\rangle
$$

We can estimate $j/r$ as follows:

1. Perform phase estimation on the state $\lvert\psi_{j}\rangle$ using a quantum circuit implementing $M_{a}$, with $m$ control qubits. The outcome is an integer $y$ such that $y/2^{m}$ approximates $j/r$.
2. Find the fraction $u/v$ in lowest terms, with $u, v \in \{0,\ldots,N-1\}$ and $v \neq 0$, that is closest to $y/2^{m}$. The continued fraction algorithm finds this fraction efficiently.

The same precision bound is enough. Two distinct fractions with denominators below $N$ are more than $1/N^{2}$ apart, so an estimate within half that gap identifies $j/r$ uniquely:

$$
\left\lvert\frac{y}{2^{m}}-\frac{j}{r}\right\rvert\le\frac{1}{2N^{2}}\quad\Longrightarrow\quad\frac{u}{v}=\frac{j}{r}
$$

Thus the same choice $m = 2\lceil\lg N\rceil + 1$ makes the correct fraction likely to be recovered.

There is one complication: continued fractions return the fraction in lowest terms. Suppose, for example, that the true eigenphase is $j/r = 2/6$. The algorithm sees only the value $1/3$, so it returns $u/v = 1/3$ rather than $2/6$. In general, if $j$ and $r$ share a common factor, the denominator returned is only $v = r/\gcd(j,r)$, a proper divisor of $r$. A single run therefore may not reveal the order.

Repeating the procedure solves this problem. Each run gives a denominator $r/\gcd(j,r)$ for an independently chosen $j$. Taking the least common multiple of the denominators observed across several runs recovers $r$ with high probability.

The continued fraction algorithm

At this point, phase estimation has given us an integer measurement outcome $y\in\{0,1,\ldots,2^{m}-1\}$. We turn it into a number in $[0,1)$ by dividing by $2^{m}$: $x=y/2^{m}$. The value $x$ is our estimate of the eigenphase $j/r$. The denominator $2^{m}$ is determined entirely by the number of control qubits, so it tells us nothing about the unknown order $r$.

What we want is a fraction $u/v$ that is close to $x$, with a denominator $v < N$. This bound comes from the problem itself: the order satisfies $r < N$. If the phase estimate satisfies $\lvert y/2^{m}-j/r\rvert\le 1/(2N^{2})$, then $x$ is close enough to the true eigenphase that this reduced fraction is uniquely determined among fractions with denominators below $N$.

This is where [continued fractions](https://mathshistory.st-andrews.ac.uk/Biographies/Cataldi/) enter. Starting from $x$, the continued fraction algorithm repeatedly divides with remainder. These divisions produce a sequence of integers $a_{0}, a_{1}, a_{2},\ldots$ called the continued-fraction terms. From these terms we construct fractions $p_{0}/q_{0},\ p_{1}/q_{1},\ p_{2}/q_{2},\ldots$ called the convergents. Each convergent is a rational approximation to $x$. As we move through the sequence, the approximations become better while their denominators grow.

Outcome $y$341

Precision $m$11

Bound $N$21

| Euclidean division | term | $p_i=a_i\,p_{i-1}+p_{i-2}$ | $q_i=a_i\,q_{i-1}+q_{i-2}$ | test |
| --- | --- | --- | --- | --- |
| $341=0\cdot2048+341$ | $a_{0}=0$ | $0$ | $1$ | $q_i<N$ |
| $2048=6\cdot341+2$ | $a_{1}=6$ | $1$ | $6$ | $q_i<N$kept |
| $341=170\cdot2+1$ | $a_{2}=170$ | $170$ | $1021$ | $q_i\ge N$stop |

phase estimate

$x=\frac{y}{2^{m}}=\frac{341}{2048}$≈ 0.16650

recovered fraction

$\frac{u}{v}=\frac{1}{6}$≈ 0.16667

distance $\le 1/(2N^{2})$

$\left\lvert x-\frac{u}{v}\right\rvert$≈ 0.00016 ≤ 0.00113

The state we can prepare

So far, the procedure was described as if we had to start with a particular eigenvector such as $\lvert\psi_{1}\rangle$. But preparing $\lvert\psi_{1}\rangle$, or any other $\lvert\psi_{j}\rangle$, would require knowing the order $r$ in advance — exactly what we are trying to find.

Fortunately, we can start with a state we already know how to prepare: $\lvert 1\rangle$. On the cycle containing $\lvert 1\rangle$, this basis state is an equal superposition of all $r$ eigenvectors: $\lvert 1\rangle=(1/\sqrt r)\sum_{j=0}^{r-1}\lvert\psi_j\rangle$.

To see this, substitute the definition of $\lvert\psi_j\rangle$: $(1/\sqrt r)\sum_{j=0}^{r-1}\lvert\psi_j\rangle=(1/r)\sum_{j=0}^{r-1}\sum_{k=0}^{r-1}\omega_r^{-jk}\lvert a^{k}\rangle$.

For $k=0$, every phase factor is $1$, so the sum over $j$ gives $r$. For every $k\neq 0$, the factors $1,\omega_r^{-k},\omega_r^{-2k},\ldots,\omega_r^{-(r-1)k}$ run through all $r$th roots of unity and sum to zero. All terms with $k\neq 0$ therefore cancel, leaving only the $k=0$ term: $(1/r)\,r\lvert 1\rangle=\lvert 1\rangle$.

This is exactly what we need. Starting with $\lvert 1\rangle$ means that phase estimation runs simultaneously on all the eigenvectors $\lvert\psi_{j}\rangle$. Each one contributes its own eigenphase $\theta_j=j/r$, and the measurement selects one of these phases.

The important point is that we do not need to know which $j$ was selected. Whatever phase we obtain has the form $j/r$, so the continued fraction step can recover its reduced denominator. Repeating the procedure with fresh copies of $\lvert 1\rangle$ gives several such denominators, whose least common multiple reveals the unknown order $r$ with high probability.

Implementation

Every piece is now in place. We know a state we can prepare, $\lvert 1\rangle$, and an operation $M_{a}$ whose eigenphases are the fractions $j/r$. We also know how to read one of those fractions off the control register. Putting them together gives the circuit that finds the order of $a \in \mathbb{Z}_N^{*}$.

$$
\textcolor{#ffffff}{M_{a}}
$$

$$
\textcolor{#ffffff}{M_{a}^{2}}
$$

$$
\textcolor{#ffffff}{M_{a}^{2^{m-1}}}
$$

$$
\textcolor{#ffffff}{\mathrm{QFT}^{\dagger}_{2^{m}}}
$$

$$
\textcolor{#6d28d9}{\lvert 0^{m}\rangle}
$$

$$
\lvert\psi\rangle
$$

What does one run cost? Write $n$ for the number of bits of $N$. The control register holds $m = 2\lceil\lg N\rceil + 1 = O(n)$ qubits, so the circuit opens with $O(n)$ Hadamard gates and closes with an inverse Fourier transform over $2^{m}$, which costs $O(n^{2})$ gates.

The controlled unitaries are the expensive part, and they are cheaper than they look. Nothing forces us to apply $M_{a}$ repeatedly: both $a$ and $N$ are known in advance, so each power $b = a^{k} \bmod N$ for $k = 1, 2, 4, 8, \ldots, 2^{m-1}$ can be worked out classically by before the circuit is built. What the circuit runs is then a single multiplication $M_{b} = M_{a}^{k}$, at cost $O(n^{2})$. With $O(n)$ of them, the controlled unitaries cost $O(n^{3})$, and that dominates the total: the whole circuit runs in $O(n^{3})$ gates.

This is the payoff of the whole construction. Searching for the order classically means walking through the powers of $a$ one at a time, and the order can be almost as large as $N$, so the walk can run to $\Omega(N) = \Omega(2^{n})$ steps — exponential in the input length. The circuit above answers the same question with a number of gates that grows like $n^{3}$.

Factoring through order-finding

Order-finding may seem far removed from factoring: it tells us about the exponents that make powers of $a$ repeat, not about the divisors of $N$. The key idea is that this periodicity contains exactly the information we need. From the order $r$ of a suitable $a$ modulo $N$, a few lines of classical arithmetic can reveal a non-trivial factor of $N$.

The reduction works under a few conditions. We take $N$ to be odd and composite, and choose $a$ so that $\gcd(a,N)=1$. For such an $N$, the [standard analysis](https://doi.org/10.1017/CBO9780511976667) guarantees that a randomly chosen $a$ produces useful factors with probability at least $1/2$. If an attempt fails, we simply choose another $a$ and repeat.

So, before running order-finding, we deal with the easy cases classically. If $N$ is **even**, we immediately have the factor $2$. If $N$ is **prime**, there is nothing to factor. Classical [Miller–Rabin](https://doi.org/10.1145/800116.803773) and [AKS](https://annals.math.princeton.edu/2004/160-2/p12) primality tests can detect this efficiently. And if $N$ is a **perfect power** $N=p^{k}$, taking successive roots hands us $p$ directly. What is left is an **odd composite that is not a prime power**, and that is the only case the quantum procedure is needed for.

How a repeating power reveals a factor

Suppose the order $r$ we get back is even. Then $r/2$ is a whole number, so we may halve the exponent and set $x=a^{r/2}\bmod N$. Squaring $x$ puts the exponent back to $r$, and $a^{r}\equiv1$ by definition of the order. So the halved power is a square root of $1$ modulo $N$:

$$
x=a^{r/2}\bmod N,\qquad x^{2}\equiv1\pmod N,\qquad N\mid(x-1)(x+1)
$$

The last step is where the factor comes from. Saying that $x$ squares to $1$ is the same as saying that $N$ divides $x^2-1$, and a difference of squares splits that quantity into the two brackets $x-1$ and $x+1$. So $N$ divides a product of two numbers that differ by only $2$ — and that is a sharp constraint on where the prime factors of $N$ can be hiding.

Because $N$ is odd, no prime of $N$ can divide two numbers that differ by $2$, so each prime power making up $N$ has to sit entirely in one bracket or the other. If they all sit in the same bracket, that bracket is a multiple of $N$ and we learn nothing. But if they are shared between the two, then $\gcd(x-1,N)$ picks up precisely the parts on the left and $\gcd(x+1,N)$ precisely the parts on the right. Both are proper factors of $N$, and Euclid’s algorithm produces them in a moment.

Modulus N21

$a$2

$r$6

$x=a^{r/2}$8

$x=2^{3}=8$$x^2\equiv1\pmod{21}$

$$
\gcd(x-1,21)=7
$$

$$
\gcd(x+1,21)=3
$$

$$
21=7\times3
$$

Putting it all together

We now have all the pieces of [Shor’s algorithm](https://arxiv.org/abs/quant-ph/9508027). The full run makes one thing especially clear: **almost all of the work is classical**. Choosing $a$, checking the conditions,, and are all classical operations. The only quantum step is finding the order $r$ — the crucial part that makes the whole approach useful.

Modulus N21

1. 1
   
   Draw $a$ at random from $\{2,\ldots,N-1\}$
   
   a = 4
2. 2
   
   Take $d=\gcd(a,N)$
3. quantum step 3
   
   Find the order $r$: the least $r>0$ with $a^{r}\equiv1\pmod N$
4. 4
   
   Check the parity of $r$
5. 5
   
   Halve the exponent: $x=a^{r/2}\bmod N$
6. 6
   
   Take $\gcd(x-1,N)$ and $\gcd(x+1,N)$

Candidates for a

splits $N$shares a factorredraw

## Grover's algorithm

Unstructured search

Let $\Sigma = \{0, 1\}$ denote the binary alphabet. Suppose we are given a function we can compute efficiently, $f: \Sigma^n \to \Sigma$, and our goal is to find a solution: a binary string $x \in \Sigma^n$ for which $f(x) = 1$.

Search

Input:

$$
f: \Sigma^n \to \Sigma
$$

Output:

A string $x \in \Sigma^n$ satisfying $f(x) = 1$, or “no solution” if no such string exists

This is unstructured search because $f$ is arbitrary. There is no promise attached to it, so there is no structure to exploit—no ordering, periodicity, or gradient to follow. Learning that $f(x) = 0$ for one string rules out that string but tells us nothing about any other.

A PIN that opens a lock is easy to check, but a failed attempt gives no clue about the next one. The same pattern appears whenever we have a cheap way to test a candidate but no useful information about where to look next. In every case, $f$ is the cheap checker, and Search asks us to find something it accepts.

Unique search

Input:

$$
f: \Sigma^n \to \Sigma
$$

Promise:

There is exactly one string $z \in \Sigma^n$ for which $f(z) = 1$

Output:

$$
z
$$

Hereafter, let $N = 2^n$ denote the number of strings in $\Sigma^n$, so that we can express costs in terms of the size of the search space rather than the number of bits. It is also useful to name the two sets into which the strings are divided: $A_1 = \{x \in \Sigma^n: f(x) = 1\}$ and $A_0 = \{x \in \Sigma^n: f(x) = 0\}$. Let $s = \lvert A_1\rvert$ be the number of solutions. Search asks us to produce an element of $A_1$, while Unique search is the special case $s = 1$.

By iterating through all $x \in \Sigma^n$ and evaluating $f$ on each one, we can solve Search with $N$ queries, and no deterministic algorithm can guarantee a solution with fewer. Probabilistic algorithms can do slightly better on average by stopping as soon as a solution is found, but they still require a number of queries that is linear in $N$.

Search by hand

Number of solutions

Queries used0Classical average33Grover6

Grover’s algorithm is a quantum algorithm for Search requiring $O(\sqrt{N})$ queries. Compared with Shor’s exponential speedup, a quadratic saving sounds modest, and it is. Whether it offers a practical advantage is a separate question. But Grover is still important: it applies to completely unstructured search, with no promise or hidden structure, and its quadratic speedup is the largest possible in the query model.

Phase query gates

So far, queries have been made through the, which writes the answer into a workspace qubit: $U_f\bigl(\lvert a\rangle\lvert x\rangle\bigr)=\lvert a\oplus f(x)\rangle\lvert x\rangle$.

Grover’s algorithm is easier to describe using a second form of query, which records the answer as a phase rather than in a qubit. For a function $f: \Sigma^n \to \Sigma$, the phase query gate $Z_f$ is the $n$-qubit operation defined by $Z_f\lvert x\rangle=(-1)^{f(x)}\lvert x\rangle$ for every $x \in \Sigma^n$.

This is the first of the two phase gates Grover’s algorithm needs: it marks the solutions by reversing their sign, and leaves every non-solution exactly as it was. A measurement cannot see that mark directly—a sign is not a probability—which is why the rest of the algorithm is needed. A single query marks every solution at once, and Grover’s algorithm is the machinery that turns those marks into amplitude a measurement can find.

Each gate from the other

$$
\lvert x\rangle
$$

$$
(-1)^{f(x)}\lvert x\rangle
$$

$$
\lvert -\rangle
$$

$$
Z_f
$$

The phase query is not a different oracle. It is the same $U_f$ query used with the workspace qubit prepared in $\lvert -\rangle$. In that case, the workspace qubit returns to $\lvert -\rangle$, while the value of $f(x)$ appears as a phase on $\lvert x\rangle$ (phase kickback): $U_f\bigl(\lvert -\rangle\lvert x\rangle\bigr)=(-1)^{f(x)}\lvert -\rangle\lvert x\rangle$.

The construction runs in the other direction too, so an algorithm counted in $Z_f$ queries and one counted in $U_f$ queries are counted on the same scale.

The second phase gate is not a query to the problem’s function $f$. It is a fixed operation that we can build directly into the circuit, based on the simple, known function $\mathrm{OR}: \Sigma^n \to \Sigma$ defined by

$$
\mathrm{OR}(x)=\begin{cases}0,&x=0^n\\[2pt]1,&x\neq 0^n.\end{cases}
$$

Its phase query gate is

$$
Z_{\mathrm{OR}}\lvert x\rangle=\begin{cases}\lvert x\rangle,&x=0^n\\[2pt]-\lvert x\rangle,&x\neq 0^n.\end{cases}
$$

Unlike $Z_f$, which queries the unknown $f$, $Z_{\mathrm{OR}}$ is a fixed, known operation that does not depend on the search problem. Its circuit can be built once and for all as an $(n-1)$-fold controlled-$Z$ gate with $X$ gates before and after it. It provides the fixed reflection used in each Grover iteration, adding gates but no queries. Thus, only $Z_f$ counts toward the query count.

$$
\lvert x_1\rangle
$$

$$
\lvert x_2\rangle
$$

$$
\lvert x_n\rangle
$$

$$
Z_{\mathrm{OR}}
$$

Grover’s algorithm

Grover’s algorithm repeatedly applies one fixed Grover operation:

$$
G=H^{\otimes n}\,Z_{\mathrm{OR}}\,H^{\otimes n}\,Z_f
$$

Each application consists of a phase query $Z_f$, followed by a fixed sequence of gates that amplifies the amplitudes of the solutions. Only $Z_f$ depends on the unknown function $f$ and therefore counts as a query, so with the other three gates fixed, $t$ iterations use exactly $t$ queries.

The operation is applied to the uniform superposition, where $N = 2^n$:

$$
\lvert u\rangle=H^{\otimes n}\lvert 0^n\rangle=\frac{1}{\sqrt N}\sum_{x\in\Sigma^n}\lvert x\rangle
$$

Before any amplification, every string has the same probability of being measured. If there are $s$ solutions, a measurement finds one with probability $s/N$—the quantum equivalent of one blind guess.

Grover’s algorithm uses the repeated application of $G$ to move probability amplitude from non-solutions to solutions. After the right number of iterations, no more and no less, measuring the register is much more likely to produce a solution. The walkthrough below runs the circuit one stage at a time, with the state after each stage and the picture that goes with it.

$$
\lvert 0^n\rangle
$$

$$
x \in \Sigma^n
$$

Step through the circuit

1. 1Prepare. Start the $n$ qubits in $\lvert 0^n\rangle$ and apply a Hadamard to each, producing the uniform superposition $\lvert u\rangle$.
2. 2Iterate. Apply the Grover operation $G$ exactly $t$ times.
3. 3Measure. Measure all $n$ qubits in the standard basis and return the resulting string. One classical evaluation of $f$ checks whether it is a solution.

Cost analysis

Two quantities determine the cost of Grover’s algorithm: the number of iterations and the success probability after those iterations. For $s$ solutions among $N$ possible strings, the optimal iteration count is approximately $t\approx\tfrac{\pi}{4\theta}-\tfrac12$, with success probability close to $1$ when the solutions are sparse.

The iteration count is where the speedup appears. Since $\sin\theta=\sqrt{s/N}$, for small $\theta$ we have $\theta\approx\sqrt{s/N}$. Therefore, $t\approx\tfrac{\pi}{4}\sqrt{N/s}=O\bigl(\sqrt{N/s}\bigr)$. Each Grover iteration uses one query to $f$, so the same expression gives the query complexity.

For Unique search, $s=1$, so Grover’s algorithm needs approximately $\tfrac{\pi}{4}\sqrt{N}$ queries. The corresponding success probability is very high: the failure probability is about $1/N$.

A classical search needs $N/2$ queries on average, so Grover’s algorithm reduces the number of queries from $O(N)$ to $O(\sqrt{N})$.

| Search space $N$ | Classical, on average | Grover iterations |
| --- | --- | --- |
| $2^{10}$ | $512$ | $25$ |
| $2^{20}$ | $524,288$ | $804$ |
| $2^{40}$ | $5.5 \times 10^{11}$ | $823,550$ |
| $2^{80}$ | $6.0 \times 10^{23}$ | $8.6 \times 10^{11}$ |

The rotation picture also explains an important limitation. The success probability is periodic in the iteration count: stopping too early leaves probability on the non-solution side, while continuing past the optimum rotates the state away from the solution direction again. Unlike classical search, where checking more candidates cannot reduce your chances, running more Grover iterations can make the result worse.

The number of solutions also changes the rotation speed. More solutions make $\theta$ larger, so each iteration rotates farther and fewer iterations are needed. The $\sqrt{N/s}$ dependence captures this directly: increasing the number of solutions makes the search easier.

There is one extreme case to handle separately. If more than half of the strings are solutions, then $\theta>\pi/4$. The usual iteration formula gives zero iterations, which is reasonable: if most strings are solutions, simply choosing a string at random already succeeds with probability greater than one half.

The useful regime for Grover’s speedup is therefore the sparse-search regime, where $s\ll N$. There, the algorithm reduces the query complexity from classical $O(N/s)$ to $O(\sqrt{N/s})$.

A natural question is whether a cleverer quantum algorithm could do better than Grover’s quadratic speedup. For unstructured search, the answer is no. Bennett, Bernstein, Brassard and Vazirani proved that every quantum algorithm solving Search with a black-box $f$ requires $\Omega(\sqrt N)$ queries.

This matches Grover’s $O(\sqrt N)$ query complexity up to a constant factor, so Grover’s algorithm is optimal. The constant $\pi/4$ in the unique-search case is optimal as well.

Unknown number of solutions

The optimal iteration count depends on $s$, the number of solutions, but $s$ is not always given to the algorithm. Fortunately, a measured candidate is easy to verify: evaluate $f(x)$ classically. A failed attempt therefore does not produce a wrong answer—it only means we need to try again.

Instead of choosing one precise iteration count, we can choose the count at random and repeat. With an appropriate randomized schedule, the algorithm finds a solution in $O\bigl(\sqrt{N/s}\bigr)$ queries when solutions exist, and $O(\sqrt{N})$ queries when there are none.

A poorly chosen iteration count can waste a single attempt because the state may have rotated past the solution direction. Randomizing the count prevents the algorithm from repeatedly getting stuck at the wrong point in the rotation. The lack of knowledge about $s$ therefore costs only a constant factor, not the quadratic speedup.

Qubits nN = 1,024

Solutions s (hidden)1

1. 1
   
   Draw $t$ at random from $\{1,\ldots,\lfloor\pi\sqrt{N}/4\rfloor\}$
2. 2
   
   Run $t$ Grover iterations on $\lvert u\rangle$
3. 3
   
   Measure all $n$ qubits, giving a string $x$
4. 4
   
   Check $f(x)$ classically: accept it or draw again

$$
\lvert A_0\rangle
$$

$$
\lvert A_1\rangle
$$

Iteration counts tried

Start searching to draw an iteration count and try it. The ceiling here is $\lfloor\pi\sqrt{N}/4\rfloor$ = 25.

Queries this run—

Knowing $s$, it would take25

Classical, on average513

The demo draws from a fixed ceiling, $\lfloor\pi\sqrt{N}/4\rfloor$, which is generous whenever solutions turn out to be plentiful. A more sophisticated approach grows the ceiling instead: set $T=1$, draw $t$ uniformly from $\{1,\ldots,T\}$, and on a failure raise $T$ and try again—stopping when the classical check accepts a string, or reporting “no solution” once $T$ has climbed past $\sqrt{N}$.

The rate of increase has to be carefully balanced. Raise $T$ too slowly and the run piles up long shots that were never likely to land, so the queries mount. Raise it too quickly and each attempt overshoots the count it was looking for, and the success probability drops. Growing by a fifth at a time, $T\leftarrow\lceil\tfrac54 T\rceil$, works.

From query complexity to real cost

The $O(\sqrt N)$ bound counts only oracle queries. In practice, each query is a full reversible implementation of $f$, and the Grover iterations must run sequentially for a deep, coherent computation. Classical search, by contrast, is easy to distribute across many machines.

So the quadratic speedup is real, given that someone eventually builds a large and stable enough quantum processor. Its cryptographic consequence is simple: Grover effectively halves the security exponent of symmetric key search and hash preimage search: $\text{AES-128}\to 2^{64},\quad\text{AES-256}\to 2^{128},\quad\text{SHA-256 preimage}\to 2^{128}.$

But a halved exponent is answered by a doubled key. Moving symmetric keys and hashes to the larger sizes restores the original margin exactly, and that is [already the standing advice](https://web3-lab.annaburd.me/how-cryptography-works/#asymmetric-cryptography-pqc), so the practical significance of the speedup keeps shrinking.

Shor’s algorithm is the sharper threat. It exploits the structure underlying RSA and elliptic-curve cryptography and breaks them outright, where no key size helps. Even there, though, standardised replacements already exist—lattice-based ML-KEM and ML-DSA, hash-based SLH-DSA—and TLS 1.3 already ships hybrid key exchange. The [risk table on the cryptography page](https://web3-lab.annaburd.me/how-cryptography-works/#asymmetric-cryptography-pqc) sorts the primitives along exactly this line—broken by Shor, weakened by Grover, or believed resistant to both.

So for now, quantum algorithms are a great deal more interesting as research than as a practical threat. The theory is settled well ahead of the hardware, and the cryptography that would be affected mostly knows what to do about it already.

These notes are based on John Watrous's [Understanding Quantum Information and Computation](https://www.youtube.com/playlist?list=PLOFEBzvs-VvqKKMXX4vbi4EB1uaErFMSO), made with IBM Quantum, with additional explanations, derivations and interactive demos.