Level 5 · Quantum Computing: Fourier and Phase

Level 4 hid answers in signs, +1 or -1, and used Hadamards to turn them into outcomes. This level does the same for phases of any angle. Its tool is the quantum Fourier transform (QFT), a gate that turns a pattern of turning arrows into a peak you can measure. It does two jobs. It finds the period of a repeating pattern, and in phase estimation it reads the angle by which a gate turns one of its special vectors. Together they make the quantum half of Shor's algorithm, which factors whole numbers, and this level ends with the classical half: continued fractions, powers modulo N and a greatest common divisor.

The first half is the QFT: one basis state, the whole matrix and why it is unitary, numpy's FFT and its opposite sign, the inverse, periodic states and their peaks, the controlled phase gate, and the short circuit that builds the QFT from Hadamards, controlled phases and swaps. The second half is phase estimation, in small steps: eigenphases, controlled gates and phase kickback, the counting register, where its outcomes land, and how to turn an outcome into a period and the period into factors.

The paths meet here. The QFT is the discrete Fourier transform of the Signals & Control path (Level 3) with the arrow turning the other way and a factor 1/\sqrt N, so numpy's np.fft computes it, and the peaks of a periodic state are the bins of a periodic signal. Euclid's algorithm, fast powers and the greatest common divisor come from the Algorithms path, and a continued fraction is nothing but the quotients of Euclid's algorithm. Seeded sampling is the Statistics path's Monte Carlo. numpy is available, and the examples use it for small full matrices.

The QFT of one basis state

A register of n qubits has N = 2^n basis states. The QFT sends each basis state |j\rangle to an even spread over all of them, each amplitude turned by its own angle:

\omega is the N-th root of unity of Level 1, one N-th of a full turn. Amplitude k is an arrow of length 1/\sqrt N turned j\cdot k times by that step. Note the plus sign in the exponent and the factor 1/\sqrt N; this path always uses this convention. For N = 8 and j = 5:

import cmath, math
n, j = 3, 5
N = 2 ** n
amps = [cmath.exp(2j * math.pi * j * k / N) / math.sqrt(N) for k in range(N)]
print("turn of amplitude k, in eighths:", [j * k % N for k in range(N)])
print("probabilities:", [round(abs(a) ** 2, 4) for a in amps])
print("amplitude 1:", complex(round(amps[1].real, 4), round(amps[1].imag, 4)))
plot([a.real for a in amps], [a.imag for a in amps], kind="scatter", title="QFT|5>, N = 8")

turn of amplitude k, in eighths: [0, 5, 2, 7, 4, 1, 6, 3]
probabilities: [0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]
amplitude 1: (-0.25-0.25j)

Press Run to see eight arrows on a circle of radius 1/\sqrt8. Every probability is 1/8: measuring \mathrm{QFT}|j\rangle straight away tells you nothing about j. The information is all in the turns, 5k \bmod 8 eighths for amplitude k. For j = 0 nothing turns and the result is the even spread, the same as a Hadamard on every qubit (Level 4). For one qubit, \omega = e^{i\pi} = -1, so \mathrm{QFT}|0\rangle = |+\rangle and \mathrm{QFT}|1\rangle = |-\rangle: the QFT on one qubit is the Hadamard gate.

Common mistakes: a minus sign in the exponent (that is the inverse, below). Dividing by N instead of \sqrt N, which leaves a state that is not normalised. Expecting the probabilities to show j.

Practice: A Ring of Arrows from a Single Index

The matrix, and why it is unitary

Column j of a gate's matrix is where it sends |j\rangle, so the QFT matrix F_N has the entry \omega^{jk}/\sqrt N in row k, column j. Swapping j and k changes nothing, so F_N is symmetric. It is a legal gate because it is unitary, and the reason is Level 1's fact that roots of unity add up to 0. The inner product of columns j and j' is

For d = 0 every term is 1 and the sum is 1. Otherwise \omega^d \ne 1 is itself a root of unity, multiplying the sum by \omega^d only moves each term to the next, so the sum equals itself times \omega^d, and must be 0. Different columns are orthogonal, and each has length 1.

import numpy as np
N = 8
k = np.arange(N)
F = np.exp(2j * np.pi * np.outer(k, k) / N) / np.sqrt(N)      # row k, column j
print("symmetric:", np.allclose(F, F.T), " unitary:", np.allclose(F.conj().T @ F, np.eye(N)))
w = np.exp(2j * np.pi / N)
print([round(float(abs(sum(w ** (d * kk) for kk in range(N)))), 4) for d in range(N)])
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
F2 = np.exp(2j * np.pi * np.outer([0, 1], [0, 1]) / 2) / np.sqrt(2)
print("F_2 is H:", np.allclose(F2, H))

symmetric: True  unitary: True
[8.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
F_2 is H: True

The second line is the sum \sum_k \omega^{dk} for d = 0, \dots, 7: N for d = 0, and 0 for every other d, even when \omega^d goes round the circle more than once (d = 2, 4, 6). np.outer(k, k) is the table of all products j\cdot k at once.

Common mistakes: 1 / N instead of 1 / np.sqrt(N). Raising a rounded \omega to a large power instead of computing each entry with exp; for big registers, reduce j * k % N first.

Practice: The Table of Turns

numpy's FFT and the opposite sign

Signals & Control (Level 3) wrote the discrete Fourier transform with the arrow turning backwards and no factor in front, X_k = \sum_n x_n e^{-2\pi i kn/N}, and its inverse with the arrow turning forwards and a factor 1/N. numpy follows the same convention: np.fft.fft is the first, np.fft.ifft the second. The QFT turns forwards with a factor 1/\sqrt N, so

• the QFT of a list a is np.fft.ifft(a) * np.sqrt(N);
• the inverse QFT of a is np.fft.fft(a) / np.sqrt(N).

On any state, the QFT is linear, so amplitude k of the result is b_k = \frac{1}{\sqrt N}\sum_j a_j\,e^{2\pi i\,jk/N}. Here is that sum next to the numpy routines, for the state 0.6|0\rangle + 0.8i|3\rangle of two qubits:

import numpy as np, cmath, math
def show(v):
    return [complex(round(z.real, 4) + 0.0, round(z.imag, 4) + 0.0) for z in v]

a = np.array([0.6, 0, 0, 0.8j])
N = len(a)
direct = [sum(a[j] * cmath.exp(2j * math.pi * j * k / N) for j in range(N)) / math.sqrt(N) for k in range(N)]
print("the formula:      ", show(direct))
print("ifft(a) * sqrt(N):", show(np.fft.ifft(a) * math.sqrt(N)))
print("fft(a) / sqrt(N): ", show(np.fft.fft(a) / math.sqrt(N)))
print("fft(a):           ", show(np.fft.fft(a)))

the formula:       [(0.3+0.4j), (0.7+0j), (0.3-0.4j), (-0.1+0j)]
ifft(a) * sqrt(N): [(0.3+0.4j), (0.7+0j), (0.3-0.4j), (-0.1+0j)]
fft(a) / sqrt(N):  [(0.3+0.4j), (-0.1+0j), (0.3-0.4j), (0.7+0j)]
fft(a):            [(0.6+0.8j), (-0.2+0j), (0.6-0.8j), (1.4+0j)]

The first two lines agree. fft(a) / sqrt(N) has the right size but the turns go the wrong way: amplitudes 1 and 3 trade places, which is the inverse QFT, a different state. fft(a) alone is not even normalised. The sum costs N^2 operations; the FFT gets the same numbers in about N\log_2 N by divide and conquer. A quantum chip does something else again, with about n^2/2 gates (below), but it never shows you the amplitudes.

Common mistakes: np.fft.fft for the QFT. Forgetting * np.sqrt(N) after ifft. Printing raw numpy arrays: round them and use .tolist().

Practice: Borrowing the Library's Transform

Undoing the QFT

The inverse of a unitary matrix is its conjugate transpose. F_N is symmetric, so only the conjugation is left, and it flips the sign in every exponent: the inverse QFT sends |k\rangle to \frac{1}{\sqrt N}\sum_j e^{-2\pi i\,jk/N}|j\rangle. That is np.fft.fft(a) / np.sqrt(N).

Applying the QFT twice does not undo it. Two forward turns add up instead of cancelling: F_N^2 sends |j\rangle to |{-j} \bmod N\rangle, and only four QFTs in a row give the identity.

import numpy as np
N = 8
k = np.arange(N)
F = np.exp(2j * np.pi * np.outer(k, k) / N) / np.sqrt(N)
Finv = F.conj().T
x = np.array([0.1, 0.5j, -0.3, 0.2, 0, 0.6, 0.4j, -0.3])
x = x / np.linalg.norm(x)
print(np.allclose(Finv @ F, np.eye(N)), np.allclose(Finv @ x, np.fft.fft(x) / np.sqrt(N)))
print("F @ F sends |j> to:", [int(np.argmax(np.abs((F @ F)[:, j]))) for j in range(N)])
print("four QFTs:", np.allclose(np.linalg.matrix_power(F, 4), np.eye(N)))

True True
F @ F sends |j> to: [0, 7, 6, 5, 4, 3, 2, 1]
four QFTs: True

|3\rangle goes to |5\rangle because -3 \bmod 8 = 5; |0\rangle and |4\rangle are their own negatives.

Common mistakes: running the QFT again to undo it. Using ifft for the inverse, which has the plus sign.

Practice: Turning Every Arrow Back

Periodic states: peaks at multiples of N/r

Here is what the QFT is for. A state is periodic with period r (where r divides N) when its amplitudes repeat every r places. The simplest is a comb: equal amplitudes on x = \text{offset}, \text{offset} + r, \text{offset} + 2r, \dots and 0 elsewhere. Take N = 32 and r = 4:

import numpy as np
N, r = 32, 4
a = np.zeros(N, dtype=complex)
a[1::r] = 1                                   # offset 1: x = 1, 5, 9, ..., 29
a = a / np.linalg.norm(a)
p = np.abs(np.fft.ifft(a) * np.sqrt(N)) ** 2  # QFT, then the Born rule
print("outcomes:", [k for k in range(N) if p[k] > 1e-9], " probabilities:", np.round(p[p > 1e-9], 4).tolist())
plot(range(N), p, kind="stem", title="a comb with period 4, after the QFT", xlabel="outcome k", ylabel="probability")

outcomes: [0, 8, 16, 24]  probabilities: [0.25, 0.25, 0.25, 0.25]

Press Run: all the probability sits on the multiples of N/r = 8. The reason is the same cancellation as before. Write x = s + mr, with s the place inside one period and m the period number. Then

The inner sum is over powers of e^{2\pi i\,rk/N}. When k is a multiple of N/r that number is 1 and the sum is N/r; otherwise it is a root of unity other than 1 and its powers add up to 0. So only k = c\cdot N/r survives, and for those k/N = c/r: the period sits in the denominator. A measurement gives a random c, and a later section recovers r from c/r.

The pattern inside one period only decides how the probability is shared among the peaks. And shifting a pattern changes no probability at all: moving every amplitude \Delta places multiplies b_k by e^{2\pi i\,k\Delta/N}, an arrow of length 1.

import numpy as np
N, r = 32, 4
def show(z):
    return complex(round(z.real, 4) + 0.0, round(z.imag, 4) + 0.0)
def periodic(pattern, shift):
    a = np.array([pattern[(x - shift) % r] for x in range(N)], dtype=complex)
    return a / np.linalg.norm(a)

comb = [0, 1, 0, 0]
b1 = np.fft.ifft(periodic(comb, 0)) * np.sqrt(N)
b3 = np.fft.ifft(periodic(comb, 2)) * np.sqrt(N)       # the same comb, two places later
print("amplitude 8:", show(b1[8]), "and", show(b3[8]), " same probabilities:", np.allclose(abs(b1) ** 2, abs(b3) ** 2))
print("shift is a turn:", np.allclose(b3, b1 * np.exp(2j * np.pi * np.arange(N) * 2 / N)))
q = np.abs(np.fft.ifft(periodic([0.2, 0.4j, -0.4, 0.8], 0)) * np.sqrt(N)) ** 2
print("another pattern:", [(k, round(float(q[k]), 4)) for k in range(N) if q[k] > 1e-9])

amplitude 8: 0.5j and -0.5j  same probabilities: True
shift is a turn: True
another pattern: [(0, 0.13), (8, 0.17), (16, 0.29), (24, 0.41)]

Common mistakes: returning amplitudes where probabilities are asked. Keeping outcomes whose probability is rounding noise (1e-33): use a threshold. Using fft instead of ifft, which mirrors the outcomes, k and N - k trading places.

Practice: Where a Repeating Pattern Lands

When the period does not fit

Real periods are rarely powers of 2. When r does not divide N, the ideal peaks sN/r fall between whole numbers, and the probability spreads over the outcomes around each of them. Take a comb with r = 7 in N = 64:

import numpy as np
N, r = 64, 7
a = np.zeros(N, dtype=complex)
a[0::r] = 1                                   # x = 0, 7, 14, ..., 63
a = a / np.linalg.norm(a)
p = np.abs(np.fft.ifft(a) * np.sqrt(N)) ** 2
print("states on the comb:", np.count_nonzero(a), " ideal peaks:", [round(s * N / r, 2) for s in range(r + 1)])
for k in range(24, 32):
    s = (k * r + N // 2) // N                 # the nearest peak, in whole numbers
    print(f"k = {k}  P = {p[k]:.4f}  nearest s = {s}  k*r - s*N = {k * r - s * N:+3d}  within 1: {abs(k * r - s * N) <= r}")
plot(range(N), p, kind="stem", title="r = 7 does not divide N = 64", xlabel="outcome k", ylabel="probability")

states on the comb: 10  ideal peaks: [0.0, 9.14, 18.29, 27.43, 36.57, 45.71, 54.86, 64.0]
k = 24  P = 0.0009  nearest s = 3  k*r - s*N = -24  within 1: False
k = 25  P = 0.0022  nearest s = 3  k*r - s*N = -17  within 1: False
k = 26  P = 0.0068  nearest s = 3  k*r - s*N = -10  within 1: False
k = 27  P = 0.0719  nearest s = 3  k*r - s*N =  -3  within 1: True
k = 28  P = 0.0350  nearest s = 3  k*r - s*N =  +4  within 1: True
k = 29  P = 0.0035  nearest s = 3  k*r - s*N = +11  within 1: False
k = 30  P = 0.0008  nearest s = 3  k*r - s*N = +18  within 1: False
k = 31  P = 0.0001  nearest s = 3  k*r - s*N = +25  within 1: False

Press Run: clusters instead of clean stems, each centred on its ideal peak. The peak at s = r sits at N itself, the same place as 0 on the circle of outcomes, so outcome 63 belongs to it and the cluster at 0 wraps round the end. It is the leakage of the Signals & Control path: a wave that does not fit a whole number of cycles into the recording spreads over the neighbouring bins. The table shows the cluster around 27.43: outcomes 27 and 28 take most of it, and the probability falls off quickly on either side. The comb has 10 states, not 64/7; with another offset it can have 9, so normalise by the number of states on the comb.

To ask "is k within distance 1 of a peak?", stay with whole numbers. |k - sN/r| \le 1 is the same as

which has no division and no rounding question, even at distance exactly 1. In the table, k r - sN for outcome 28 is +4, so it is 4/7 from the peak and counts; for 26 it is -10, more than 1 away. Integer rounding, (k * r + N // 2) // N, finds the nearest s.

Common mistakes: dividing by \sqrt{N/r} instead of the square root of the number of states on the comb. Forgetting the peak at N, which leaves out outcomes just below N.

Practice: When the Period Does Not Fit the Register

The controlled phase gate

The QFT circuit needs one new gate. The controlled phase CP(\lambda) on two qubits multiplies the amplitude of every basis state where both qubits read 1 by e^{i\lambda}; on two qubits its matrix is \mathrm{diag}(1, 1, 1, e^{i\lambda}). It is usually drawn as "apply the phase gate P(\lambda) = \mathrm{diag}(1, e^{i\lambda}) to the target when the control reads 1", but "both read 1" shows that the two qubits play the same role: control and target may be swapped. With full matrices on 3 qubits, using the projectors |0\rangle\langle0| = \mathrm{diag}(1, 0) and |1\rangle\langle1| = \mathrm{diag}(0, 1) to say "if the control reads 0, do nothing; if it reads 1, apply P":

import numpy as np, cmath
lam = 0.9
I, P0, P1 = np.eye(2), np.diag([1, 0]), np.diag([0, 1])
P = np.diag([1, cmath.exp(1j * lam)])
def kron(*ms):
    out = np.eye(1)
    for m in ms:
        out = np.kron(out, m)
    return out

cp_02 = kron(P0, I, I) + kron(P1, I, P)      # control qubit 0, phase gate on qubit 2
cp_20 = kron(I, I, P0) + kron(P, I, P1)      # control qubit 2, phase gate on qubit 0
print("the same gate:", np.allclose(cp_02, cp_20), " diagonal:", np.allclose(cp_02, np.diag(np.diag(cp_02))))
print("turned:", [format(i, "03b") for i in range(8) if abs(cp_02[i, i] - 1) > 1e-12], "by", round(cmath.phase(cp_02[5, 5]), 4))

the same gate: True  diagonal: True
turned: ['101', '111'] by 0.9

Only |101\rangle and |111\rangle turn, the states where qubits 0 and 2 both read 1 (qubit 0 is the leftmost bit, as always). The gate is diagonal: it never moves amplitude, it only turns it. On a simulator it is one pass over the indices with a mask that has both qubits' bits set (Level 3).

Common mistakes: masks 1 << q, which number the qubits from the right. Turning states where only one of the two qubits reads 1. Treating CP as if the order of its qubits mattered.

Practice: A Turn Only When Both Lamps Are Lit

The QFT, one qubit at a time

The QFT of a basis state is a product state: every qubit carries its own turn, and nothing is entangled. Write the outcome k in bits, k = k_0 2^{n-1} + k_1 2^{n-2} + \dots + k_{n-1} (qubit 0 most significant). The exponential of a sum is a product, so

Qubit q is turned by j/2^{q+1} of a turn, and only the fractional part counts. Write j in bits j_0 j_1 \dots j_{n-1} and use binary fractions, 0.b_1b_2b_3 = b_1/2 + b_2/4 + b_3/8. Then qubit 0 gets 0.j_{n-1} (a half turn or none), qubit 1 gets 0.j_{n-2}j_{n-1}, and the last qubit gets 0.j_0j_1\dots j_{n-1} = j/N, which depends on every bit. For j = 6 = 110 in three qubits:

import numpy as np
n, j = 3, 6
N = 2 ** n
formula = np.exp(2j * np.pi * j * np.arange(N) / N) / np.sqrt(N)
turns = [j / 2 ** (q + 1) % 1 for q in range(n)]
print("turn of qubit q:", turns)
state = np.ones(1)
for t in turns:
    state = np.kron(state, np.array([1, np.exp(2j * np.pi * t)]) / np.sqrt(2))
print("product of the three qubits = QFT|6>:", np.allclose(state, formula))

turn of qubit q: [0.0, 0.5, 0.75]
product of the three qubits = QFT|6>: True

Qubit 0 gets 0.0 = 0, qubit 1 gets 0.10 = 0.5, qubit 2 gets 0.110 = 0.75. This product form is what the circuit builds, one qubit at a time.

Common mistakes: forgetting the % 1, which keeps whole turns. Putting j/N on qubit 0: in this bit order it belongs on the last qubit.

The circuit: Hadamards, controlled phases, swaps

How does a circuit put 0.j_0j_1\dots j_{n-1} on a qubit? Start with qubit 0, which holds the bit j_0. A Hadamard gives (|0\rangle + (-1)^{j_0}|1\rangle)/\sqrt2, and (-1)^{j_0} = e^{2\pi i\,0.j_0}: the first binary digit. A CP(\pi/2) between qubit 0 and qubit 1 then adds a quarter turn exactly when both read 1, which is j_1/4 of a turn on the |1\rangle part of qubit 0: the second digit. A CP(\pi/4) with qubit 2 adds j_2/8, the third. In general the controlled phase between qubit q and a later qubit m has angle \pi/2^{m-q}: the further away, the smaller the turn. Here is the stage for qubit 0 with full matrices on 3 qubits:

import numpy as np, cmath, math
n = 3
I, P0, P1 = np.eye(2), np.diag([1, 0]), np.diag([0, 1])
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
def kron(*ms):
    out = np.eye(1)
    for m in ms:
        out = np.kron(out, m)
    return out
def cp(a, b, lam):                              # CP(lam) on qubits a and b of n
    zero, one = [I] * n, [I] * n
    zero[a] = P0
    one[a], one[b] = P1, np.diag([1, cmath.exp(1j * lam)])
    return kron(*zero) + kron(*one)

stage0 = cp(0, 2, math.pi / 4) @ cp(0, 1, math.pi / 2) @ kron(H, I, I)   # H first, then the turns
for j in (6, 3):
    e = np.zeros(2 ** n)
    e[j] = 1
    out = stage0 @ e
    rest = j & 0b011                            # qubits 1 and 2 still hold their input bits
    turn = cmath.phase(out[4 + rest] / out[rest]) / (2 * math.pi) % 1
    print(format(j, "03b"), " turn on qubit 0:", round(turn, 4), " j / 8 =", j / 8)

110  turn on qubit 0: 0.75  j / 8 = 0.75
011  turn on qubit 0: 0.375  j / 8 = 0.375

After its stage, qubit 0 holds 0.j_0j_1j_2 = j/8, the factor that the product form puts on the last qubit. The controlled phases could use j_1 and j_2 because qubits 1 and 2 had not been touched yet. So the circuit works through the qubits in order q = 0, 1, \dots, n-1, each with a Hadamard and then controlled phases with every later qubit. Qubit q ends up holding 0.j_q\dots j_{n-1}, which belongs on qubit n-1-q: everything comes out in reversed order, and a final row of swaps, qubit q with qubit n-1-q, puts it right. In four qubits, for example, qubit 0's stage uses the angles \pi/2, \pi/4, \pi/8 with qubits 1, 2, 3, and the swaps are (0, 3) and (1, 2).

The swaps reverse the bit string of every index: in four qubits, |0001\rangle (index 1) trades places with |1000\rangle (index 8), and |0110\rangle stays. To check a whole circuit, build its matrix column by column, running it on each basis state, and compare with F_N up to a global phase.

The count: n Hadamards, n(n-1)/2 controlled phases and \lfloor n/2\rfloor swaps. For 20 qubits that is 20 + 190 + 10 = 220 gates for a matrix with 2^{40} entries. Most of the controlled phases have tiny angles (\pi/2^{19} is about six millionths of a radian), and the approximate QFT leaves out those below a threshold, at a small, measurable cost in accuracy.

Common mistakes: negative angles, which build the inverse QFT. The angle \pi/2^m instead of \pi/2^{m-q}. Forgetting the swaps, or swapping every pair twice, which undoes them.

Practice: Reading the Register in a Mirror, The Fourier Gate, Wire by Wire, Skipping the Smallest Turns

Eigenphases: how far a gate turns a vector

The second half of the level reads angles off a gate. Some vectors pass through a gate U only multiplied by a number: U|u\rangle = \lambda|u\rangle, an eigenvector with eigenvalue \lambda (Level 1). A unitary never changes lengths, so |\lambda| = 1, a pure turn, and this path writes it in whole turns:

The number \varphi is the eigenphase.

So \varphi = 0.25 is a quarter turn (\lambda = i), and \varphi = 0.5 is \lambda = -1. The rotation R_x(\theta) has the eigenvectors |+\rangle and |-\rangle:

import numpy as np, cmath, math
theta = 1.2
Rx = np.array([[math.cos(theta / 2), -1j * math.sin(theta / 2)],
               [-1j * math.sin(theta / 2), math.cos(theta / 2)]])
for name, u in [("|+>", np.array([1, 1]) / np.sqrt(2)), ("|->", np.array([1, -1]) / np.sqrt(2)), ("|0>", np.array([1, 0]))]:
    w = Rx @ u
    lam = np.vdot(u, w)                          # <u|Uu>, the only candidate for lambda
    phi = cmath.phase(lam) / (2 * math.pi) % 1.0
    print(name, " eigenvector:", np.allclose(w, lam * u), " phi =", round(phi, 4))
print((-1e-16) % 1.0, (-1e-17) % 1.0)

|+>  eigenvector: True  phi = 0.9045
|->  eigenvector: True  phi = 0.0955
|0>  eigenvector: False  phi = 0.0
0.9999999999999999 1.0

R_x(1.2) turns |+\rangle by e^{-0.6i}, which is -0.6/2\pi \approx -0.0955 turns, or 0.9045 once brought into [0, 1) by % 1.0; it turns |-\rangle by +0.0955. |0\rangle is not an eigenvector, so its \varphi means nothing. The last line is a rounding trap: an eigenvalue that should be exactly 1 can come out with a tiny negative angle, and % 1.0 turns it into 0.9999999999999999, or even into 1.0, which is outside [0, 1); neither is 0. Treat a phase within about 1e-9 of 1 as 0.

Common mistakes: giving the angle in radians, or in (-0.5, 0.5]. Dividing w[i] / u[i] where u[i] is 0. Calling a vector an eigenvector because Uu has the same length as u.

Practice: How Far Does the Gate Turn This Vector?

Controlled-U and phase kickback

A controlled-U uses one qubit to decide whether U acts on the others: CU|0\rangle|w\rangle = |0\rangle|w\rangle and CU|1\rangle|w\rangle = |1\rangle U|w\rangle. With the control as qubit 0 and the target as the m qubits after it, the first half of the indices has control 0 and the second half control 1, so the matrix is the block matrix

and the amplitude for control c and target basis state w sits at index c * 2**m + w. Now give the target an eigenvector. The |1\rangle part of the control picks up e^{2\pi i\varphi}, the target does not change, and the eigenphase lands on the control:

That is Level 4's phase kickback, where the eigenvalue was -1; here it can be any turn. With R_x(1.2) and |-\rangle, where \varphi \approx 0.0955:

import numpy as np, cmath, math
theta = 1.2
Rx = np.array([[math.cos(theta / 2), -1j * math.sin(theta / 2)],
               [-1j * math.sin(theta / 2), math.cos(theta / 2)]])
minus = np.array([1, -1]) / np.sqrt(2)
phi = theta / 2 / (2 * math.pi)                 # Rx|-> = e^{2 pi i phi}|->
CU = np.block([[np.eye(2), np.zeros((2, 2))], [np.zeros((2, 2)), Rx]])
c = np.array([0.6, 0.8])
after = CU @ np.kron(c, minus)
kicked = np.array([0.6, 0.8 * cmath.exp(2j * math.pi * phi)])
print("control picks up the turn, target unchanged:", np.allclose(after, np.kron(kicked, minus)))

H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
HI = np.kron(H, np.eye(2))
final = HI @ CU @ HI @ np.kron([1, 0], minus)   # H on the control, CU, H again
p0 = float(sum(abs(final[w]) ** 2 for w in range(2)))   # control reads 0: the first half
print(round(p0, 4), round((1 + math.cos(2 * math.pi * phi)) / 2, 4))

control picks up the turn, target unchanged: True
0.9127 0.9127

The second part is the Hadamard test: Hadamards on the control around CU turn the kicked phase into a probability, p_0 = (1 + \cos 2\pi\varphi)/2. One number cannot tell \varphi from 1 - \varphi (the cosine is the same), and to pin \varphi down to many digits you would need a great many runs. Phase estimation does better.

Common mistakes: applying U to the control-0 half. Expecting the target to change. Forgetting that with the control in |1\rangle alone the turn is a global phase, invisible.

Practice: The Control Qubit Feels the Turn, Two Phases That Look Alike

The counting register: doubling the turn

A power of U turns an eigenvector further: U^2|u\rangle = e^{2\pi i\,2\varphi}|u\rangle, and U^{2^k}|u\rangle = e^{2\pi i\,2^k\varphi}|u\rangle. Phase estimation uses t counting qubits, qubits 0 to t-1, in front of the target:

1. a Hadamard on every counting qubit;
2. counting qubit j controls U^{2^{t-1-j}}: qubit 0, the most significant bit, gets the largest power U^{2^{t-1}}, and the last counting qubit gets U itself.

Each controlled power kicks its turn back onto its own qubit. In the basis state |y\rangle of the counting register, qubit j is worth 2^{t-1-j}, exactly the power it controls, so the turns add up to \varphi\cdot y:

The amplitude for counting value y and target basis state w sits at index y * 2**m + w. With full matrices, t = 2, and the gate U = P(2\pi\cdot0.3), whose eigenvector |1\rangle has \varphi = 0.3:

import numpy as np, cmath, math
phi = 0.3
U = np.diag([1, cmath.exp(2j * math.pi * phi)])
u = np.array([0, 1])
I, P0, P1 = np.eye(2), np.diag([1, 0]), np.diag([0, 1])
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
def kron(*ms):
    out = np.eye(1)
    for m in ms:
        out = np.kron(out, m)
    return out

def counting(power0, power1):
    state = kron(H, H, I) @ np.kron(np.kron([1, 0], [1, 0]), u)   # 2 counting qubits, then the target
    state = (kron(P0, I, I) + kron(P1, I, power0)) @ state        # qubit 0 controls power0
    state = (kron(I, P0, I) + kron(I, P1, power1)) @ state        # qubit 1 controls power1
    return [state[y * 2 + 1] for y in range(4)]                   # |y>|1>: where u lives

reg = counting(U @ U, U)
print("sizes:", [round(float(abs(a)), 4) for a in reg])
print("turns:", [round(cmath.phase(a) / (2 * math.pi) % 1, 4) for a in reg])
print("powers swapped:", [round(cmath.phase(a) / (2 * math.pi) % 1, 4) for a in counting(U, U @ U)])

sizes: [0.5, 0.5, 0.5, 0.5]
turns: [0.0, 0.3, 0.6, 0.9]
powers swapped: [0.0, 0.6, 0.3, 0.9]

The turns are 0, 0.3, 0.6, 0.9: \varphi times y. With the powers the wrong way round, y = 1 and y = 2 trade their turns. Squaring the matrix gives the next power, U, U^2, U^4, \dots, without multiplying U in a long loop.

Now compare with the very first section: \mathrm{QFT}|x\rangle has amplitude e^{2\pi i\,xy/N}/\sqrt N at y. When \varphi = x/2^t exactly, the counting register is \mathrm{QFT}|x\rangle, so the inverse QFT turns it into |x\rangle.

Common mistakes: giving qubit 0 the power U^1. Applying U once per qubit whatever its position. Reading the target's index instead of the counting value.

Practice: Doubling the Turn on Every Wire

Reading the phase: the inverse QFT

The last two steps of phase estimation: the inverse QFT on the counting register, then a measurement of it. Outcome x estimates \varphi \approx x/2^t. When \varphi = 5/8 and t = 3, the outcome is 5 for sure. When \varphi = 0.62, \varphi\cdot 8 = 4.96 is not a whole number, and the probability spreads:

import numpy as np
for phi, t in [(5 / 8, 3), (0.62, 3), (0.62, 4)]:
    N = 2 ** t
    reg = np.exp(2j * np.pi * phi * np.arange(N)) / np.sqrt(N)      # the counting register
    probs = np.abs(np.fft.fft(reg) / np.sqrt(N)) ** 2               # inverse QFT, Born rule
    best = int(np.argmax(probs))
    print(f"phi = {phi}, t = {t}: best outcome {best} (estimate {best / N}) with P = {probs[best]:.4f}, total {probs.sum():.4f}")
plot(range(16), probs, kind="bar", title="phi = 0.62, t = 4", xlabel="outcome x", ylabel="probability")

phi = 0.625, t = 3: best outcome 5 (estimate 0.625) with P = 1.0000, total 1.0000
phi = 0.62, t = 3: best outcome 5 (estimate 0.625) with P = 0.9948, total 1.0000
phi = 0.62, t = 4: best outcome 10 (estimate 0.625) with P = 0.9792, total 1.0000

Press Run to see the last case: 0.62\cdot16 = 9.92, and outcome 10 (estimate 0.625) takes most of the probability, with its neighbours sharing the rest. One more counting qubit halves the step between estimates. The two outcomes on either side of \varphi\cdot 2^t always hold at least 8/\pi^2 \approx 0.81 of the probability together.

If the target is not an eigenvector, write it as a combination \sum_i c_i|u_i\rangle of eigenvectors. Each part runs through the circuit with its own phase, and the outcome distribution becomes the mixture \sum_i |c_i|^2 P_{\varphi_i}(x): the counting register is measured alone, so its probabilities add up over the target's basis states (Level 3's partial measurement).

Common mistakes: the plus sign in the inverse QFT, which moves the peak to 2^t - x. Forgetting one of the two factors 1/\sqrt N, so the list sums to N or 1/N.

Practice: Where the Counting Register Lands, The Whole Phase Reader on a Small Chip

A closed form, and distance on a circle

The probability of each outcome has a formula. Put d = \varphi - x/N. The inverse QFT turns amplitude y back by xy/N, so

a geometric series with ratio e^{2\pi i d}. It adds up to \frac{1}{N}\cdot\frac{1 - e^{2\pi i Nd}}{1 - e^{2\pi i d}}, and with |1 - e^{i\alpha}| = 2|\sin(\alpha/2)|:

One line per outcome instead of a sum of N terms. Checked against the last case of the previous section:

import math
phi, t = 0.62, 4
N = 2 ** t
for x in (9, 10, 11):
    d = phi - x / N
    print(x, round(d, 4), round((math.sin(math.pi * N * d) / (N * math.sin(math.pi * d))) ** 2, 4))

def dist(a, b):                                     # distance between two phases on the circle
    g = abs(a - b) % 1.0
    return min(g, 1 - g)
phi, t = 0.98, 3
N = 2 ** t
for x in (0, 7):
    d = phi - x / N
    print(x / N, round(dist(phi, x / N), 4), round((math.sin(math.pi * N * d) / (N * math.sin(math.pi * d))) ** 2, 4))

9 0.0575 0.0075
10 -0.005 0.9792
11 -0.0675 0.0055
0.0 0.02 0.9198
0.875 0.105 0.0346

Phases live on a circle: \varphi = 0.98 is only 0.02 from the estimate 0.0, closer than to 0.875. The formula already knows, because \sin^2 repeats: d = 0.98 gives the same probability as d = -0.02. A distance between phases must wrap the same way: abs(a - b) % 1, then the shorter way round. More counting qubits do not always help at once: for a fixed tolerance, one more qubit can lower the chance of a good estimate before the next one raises it.

Common mistakes: dividing by \sin(\pi d) when it is 0. Measuring the error as abs(phi - x / N) without wrapping. Forgetting to square.

Practice: How Many Wires for a Sharp Reading?

Continued fractions

Period finding hands back an outcome y with y/2^t \approx s/r, and the question is which simple fraction y/2^t is close to. Continued fractions answer it. Every fraction can be written as

with whole-number terms a_0 \ge 0 and a_1, a_2, \dots \ge 1. Take the whole part, flip what is left, and repeat. With the fraction stored as two integers, "flip what is left" turns the pair (a, b) into (b, a % b): the terms are exactly the quotients of Euclid's algorithm for the greatest common divisor. Take 13653/16384:

y, q = 13653, 16384
a, b = y, q
terms = []
while b:
    terms.append(a // b)
    print(f"{a} = {a // b} * {b} + {a % b}")
    a, b = b, a % b
print("terms:", terms)

x, float_terms = y / q, []
for step in range(8):
    whole = int(x)
    float_terms.append(whole)
    if x == whole:
        break
    x = 1 / (x - whole)
print("with floats:", float_terms)

13653 = 0 * 16384 + 13653
16384 = 1 * 13653 + 2731
13653 = 4 * 2731 + 2729
2731 = 1 * 2729 + 2
2729 = 1364 * 2 + 1
2 = 2 * 1 + 0
terms: [0, 1, 4, 1, 1364, 2]
with floats: [0, 1, 4, 1, 1364, 1, 1, 57823383]

The exact expansion stops after six terms, when the remainder is 0. The float version agrees for five terms and then produces nonsense, although 13653/16384 is stored exactly as a float: each 1 / (x - whole) adds a little rounding error, and the error grows with every flip until it swamps the true value. Use integers (//, %, or fractions.Fraction), which Python keeps exact at any size.

Cutting the expansion after term a_k gives the convergent h_k/k_k. The convergents follow a two-step recurrence, starting from h_{-2} = 0, h_{-1} = 1 and k_{-2} = 1, k_{-1} = 0:

from fractions import Fraction
target = Fraction(13653, 16384)
h_prev, h = 0, 1
k_prev, k = 1, 0
for a in [0, 1, 4, 1, 1364, 2]:
    h_prev, h = h, a * h + h_prev
    k_prev, k = k, a * k + k_prev
    print(f"term {a:4}  convergent {h}/{k}  off by {float(abs(target - Fraction(h, k))):.2e}")

term    0  convergent 0/1  off by 8.33e-01
term    1  convergent 1/1  off by 1.67e-01
term    4  convergent 4/5  off by 3.33e-02
term    1  convergent 5/6  off by 2.03e-05
term 1364  convergent 6824/8189  off by 7.45e-09
term    2  convergent 13653/16384  off by 0.00e+00

The convergents come out in lowest terms, their denominators grow, and each is closer than the one before; the last is the fraction itself. 5/6 is the simple fraction hiding in 13653/16384: it is off by less than 10^{-4}, and the next convergent needs a denominator over 8000 to do better.

Common mistakes: floats. Stopping one step early and losing the last term. Forgetting a_0 = 0 when y < q.

Practice: Peeling a Fraction Down to Whole Numbers

From an outcome to a period

Shor's algorithm factors N through the order of a number a with \gcd(a, N) = 1: the smallest r \ge 1 with a^r \bmod N = 1. Python's pow(a, e, N) computes a^e \bmod N by repeated squaring, fast even when e has hundreds of digits. For N = 91 and a = 3:

N, a = 91, 3
print([pow(a, e, N) for e in range(1, 9)])
r = next(e for e in range(1, N) if pow(a, e, N) == 1)
print("order:", r)
print("ideal outcomes s/r * 2**14:", [round(s / r * 2 ** 14, 2) for s in range(r)])

[3, 9, 27, 81, 61, 1, 3, 9]
order: 6
ideal outcomes s/r * 2**14: [0.0, 2730.67, 5461.33, 8192.0, 10922.67, 13653.33]

The powers repeat with period 6. The quantum part runs phase estimation on the gate |w\rangle \mapsto |a\,w \bmod N\rangle, whose eigenphases are s/r for s = 0, 1, \dots, r - 1. After the controlled powers, each counting value x sits next to the target value a^x \bmod N, which repeats with period r in x: the comb picture of the earlier sections, once for every target value. An outcome y lands near one of the ideal outcomes s\cdot 2^t/r, for an s you do not know. With t = 14 counting qubits, s = 5 gives y = 13653, the fraction of the previous section.

The rule that turns y into a candidate period: take the convergents of y/2^t, keep the last one whose denominator is less than N, and check its denominator r with pow(a, r, N) == 1. The period is less than N, the denominators grow, and the later convergents are the closer ones. A theorem makes it safe: if |y/2^t - s/r| < 1/(2r^2), then s/r (in lowest terms) is one of the convergents, and 2^t \ge N^2 makes the nearest outcomes that close; here 2^{14} = 16384 \ge 91^2 = 8281. Four outcomes:

N, a, t = 91, 3, 14
for y in (13653, 13654, 5461, 0):
    num, den = y, 2 ** t
    dens, k_prev, k = [], 1, 0
    while den:
        term = num // den
        k_prev, k = k, term * k + k_prev
        dens.append(k)
        num, den = den, num % den
    r = [d for d in dens if d < N][-1]
    print(f"y = {y:5}  denominators {dens}  candidate {r}  pow(a, r, N) = {pow(a, r, N)}")

y = 13653  denominators [1, 1, 5, 6, 8189, 16384]  candidate 6  pow(a, r, N) = 1
y = 13654  denominators [1, 1, 6, 4093, 8192]  candidate 6  pow(a, r, N) = 1
y =  5461  denominators [1, 3, 16384]  candidate 3  pow(a, r, N) = 27
y =     0  denominators [1]  candidate 1  pow(a, r, N) = 3

13653 and 13654, the two outcomes next to 5\cdot 16384/6 \approx 13653.33, both give 6. Earlier convergents (denominators 1 and 5) would give wrong candidates, which is why the rule takes the last one below N, not the first. 5461 \approx 16384/3 came from s = 2: 2/6 = 1/3 in lowest terms, the shared factor hides part of the period, and the check fails. y = 0 comes from s = 0 and says nothing. A failed check costs one more run.

Common mistakes: taking the numerator. Taking the last convergent overall (16384 is not a period). Computing a ** r % N, which builds an enormous number first.

Practice: From a Dial Reading to a Repeat Length

Factors from a period, shot by shot

The last step is pure number theory. If r is even, let x = a^{r/2} \bmod N. Then x^2 \bmod N = 1, so N divides x^2 - 1 = (x - 1)(x + 1). If x is neither 1 nor N - 1, N divides neither factor alone, so each factor shares a proper divisor with N, and the greatest common divisor finds it:

import math
N, a, r = 91, 3, 6
x = pow(a, r // 2, N)
print("x =", x, " x*x mod N =", x * x % N, " factors:", math.gcd(x - 1, N), math.gcd(x + 1, N))
for N, a in [(143, 3), (133, 3)]:
    r = next(e for e in range(1, N) if pow(a, e, N) == 1)
    print(N, a, "order", r, " x =", pow(a, r // 2, N) if r % 2 == 0 else "none: r is odd")

x = 27  x*x mod N = 1  factors: 13 7
143 3 order 15  x = none: r is odd
133 3 order 18  x = 132

27^2 = 729 = 8\cdot91 + 1, and \gcd(26, 91) = 13, \gcd(28, 91) = 7: 91 = 7\cdot13. Two things can go wrong. The order can be odd, like 15 for a = 3 modulo 143. Or x can be N - 1, a square root of 1 that says nothing: for a = 3 modulo 133, x = 132, and \gcd(131, 133) = 1, \gcd(133, 133) = 133. Then try another a; when N has two different odd prime factors, a random a works at least half the time.

A whole run is many shots: each shot measures one outcome, and the classical part either succeeds or the next shot is tried. Simulated shots use Level 2's procedure, here with the draw called d, because r is taken by the period: make rng = random.Random(seed) once, and for each shot draw d = rng.random() and pick the first outcome whose running total of probabilities is greater than d (if rounding leaves none, the last outcome with nonzero probability). Every shot draws exactly once, whether or not it then succeeds. Six shots of phase estimation with \varphi = 0.58 and t = 3, where 0.58\cdot 8 = 4.64:

import numpy as np, random
phi, t = 0.58, 3
N = 2 ** t
probs = (np.abs(np.fft.fft(np.exp(2j * np.pi * phi * np.arange(N)) / np.sqrt(N)) / np.sqrt(N)) ** 2).tolist()
rng = random.Random(7)
for shot in range(1, 7):
    d = rng.random()
    total, outcome = 0.0, max(i for i in range(N) if probs[i] > 1e-12)
    for i, p in enumerate(probs):
        total += p
        if total > d:
            outcome = i
            break
    print(shot, round(d, 4), outcome, outcome / N)

1 0.3238 5 0.625
2 0.1508 4 0.5
3 0.6509 5 0.625
4 0.0724 3 0.375
5 0.5359 5 0.625
6 0.3657 5 0.625

The running totals are 0.0136, 0.0267, 0.044, 0.0794, 0.2863, 0.9306, \dots, so d = 0.0724 stops at outcome 3 and d = 0.1508 at outcome 4; most shots read 5, the estimate 0.625 closest to 0.58.

Common mistakes: forgetting the odd case or the x = N - 1 case. A new random.Random(seed) in every shot. Skipping the draw for a shot that fails, which shifts every later shot.

Aside (optional): why factoring matters. The RSA cryptosystem, still common for signatures and key exchange on the web, publishes a number N = pq with two secret primes of about 300 digits each, and its safety rests on nobody being able to factor N. The best classical methods need far more time than the age of the universe for that size; Shor's algorithm needs a number of gates that grows only like a small power of the number of digits. The catch is the hardware. Recent estimates for breaking a 2048-bit RSA key call for around a million physical qubits running for days with error correction, while today's largest chips have on the order of a thousand qubits, too noisy for long circuits, and the numbers honestly factored with Shor's algorithm so far are tiny, like 15 and 21. Because recorded messages could be decrypted later, new "post-quantum" schemes that resist Shor's algorithm are already being deployed.

Practice: A Repeat Length Cracks the Number, Cracking a Small Number, Shot by Shot

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

• Borrowing the Library's Transform
• The Fourier Gate, Wire by Wire
• Where a Repeating Pattern Lands
• Doubling the Turn on Every Wire
• Where the Counting Register Lands
• From a Dial Reading to a Repeat Length

Checklist

You can now:

• Write \mathrm{QFT}|j\rangle = \frac{1}{\sqrt N}\sum_k e^{2\pi i\,jk/N}|k\rangle and its matrix, and explain with roots of unity why it is unitary.
• Compute the QFT as np.fft.ifft(a) * np.sqrt(N) and the inverse as np.fft.fft(a) / np.sqrt(N), and say why the QFT applied twice gives |{-j} \bmod N\rangle.
• Predict the peaks of a periodic state at the multiples of N/r, explain why a shift leaves the probabilities alone, and test closeness to a fractional peak with integers, |kr - sN| \le r.
• Apply the controlled phase CP(\lambda), which turns the states where both qubits read 1.
• Explain the QFT circuit from the product form: a Hadamard and controlled phases \pi/2^{m-q} build 0.j_q\dots j_{n-1} on qubit q, and swaps undo the reversed order.
• Find the eigenphase \varphi\in[0,1) of an eigenvector, and use CU = \mathrm{diag}(I, U) and phase kickback to move it onto a control qubit.
• Build the counting register \frac{1}{\sqrt{2^t}}\sum_y e^{2\pi i\varphi y}|y\rangle with qubit j controlling U^{2^{t-1-j}}, and find where phase estimation's outcomes land, by the inverse QFT or by \sin^2(\pi Nd)/(N^2\sin^2(\pi d)).
• Measure distances between phases around the circle.
• Expand a fraction into a continued fraction with exact integers, compute its convergents, and turn an outcome into a period with the last convergent whose denominator is below N.
• Split N with \gcd(a^{r/2} \pm 1, N), recognise the failing cases, and run seeded shots with one generator.