The Fourier Family

The classical discrete Fourier transform takes a vector of $N$ complex numbers $x_0,\dots ,x_{N-1}$ and returns another, $y_k=\frac{1}{\sqrt{N}}{\sum }_j{x}_j{e}^{2\pi ijk/N}$. It is the single most useful transform in all of computing — it diagonalizes convolution, exposes periodicity, and underlies signal processing wholesale. The fast Fourier transform computes it in $O(N\log N)$ operations. The Quantum Fourier Transform is the same linear map, but acting on the amplitudes of an $n$-qubit state with $N=2^n$:

extended linearly to superpositions. Here is the first miracle. A classical FFT on $N=2^n$ numbers costs $O(N\log N)=O(n{2}^n)$ — exponential in the number of qubits. The QFT circuit costs $O(n^2)$ gates, and the reason deserves a derivation, not an assertion: the dense-looking transform matrix secretly factors into a tensor product. Write the input index in binary, $j=j_0{j}_1\cdots j_{n-1}$ with $j_0$ the most significant bit (the book's convention), and expand the output index the same way, $k={\sum }_{l=0}^{n-1}{k}_l{2}^{n-1-l}$. Substituted into the phase of (7.1), the exponential splits over the bits of $k$,

so the single sum over all $2^n$ values of $k$ in (7.1) factors into $n$ independent two-term sums, one per output qubit. The sum over $k_l\in \{0,1\}$ is just $(|0⟩+e^{2\pi ij/2^{l+1}}|1⟩)/\sqrt{2}$; and because a phase only matters modulo 1, the number $j/2^{l+1}$ enters only through the last $l+1$ bits of $j$ — the binary fraction $0.j_{n-1-l}\cdots j_{n-1}$ (where $0.b_1{b}_2\cdots$ means $b_1/2+b_2/4+\cdots$). Reading the factors off directly — the $l$-th tensor slot is output qubit $l$, which the phase $e^{2\pi ij{k}_l/2^{l+1}}$ pairs with the fraction $j/2^{l+1}$ — the QFT of a basis state is an unentangled product:

qubit by qubit: the shortest fraction $0.j_{n-1}$ belongs to output qubit 0, and the longest $0.j_0{j}_1\cdots j_{n-1}$ to output qubit $n-1$. Note the reversal, because it is about to cost us a layer of gates: the fraction on output qubit $l$ starts at input bit $j_{n-1-l}$, the mirror-image wire. Each factor is cheap to build in place: a Hadamard on the wire holding $j_l$ produces $(|0⟩+e^{2\pi i0.j_l}|1⟩)/\sqrt{2}$ (since $e^{2\pi i0.j_l}=(-1{)}^{j_l}$), and each further bit of the fraction is stamped on by a two-qubit controlled phase rotation

controlled on the wire that holds the bit in question. One Hadamard per wire and a rotation for each pair: $n+(n-1)+\cdots +1=O(n^2)$ gates. And the derivation hands us the final swaps for free: the circuit builds the $0.j_l\cdots j_{n-1}$ factor on wire $l$ (where the controls it needs live), but (7.2) says that factor belongs on qubit $n-1-l$ — every factor lands on the mirror-image wire. The product comes out in reversed wire order, and the closing layer of swaps exists purely to reverse it back.

But now the second fact, the one that keeps this honest: you cannot read the transform out. The QFT rearranges $2^n$ amplitudes, but a measurement still returns just $n$ bits, sampled by the Born rule (Chapter 1's Holevo ceiling has not gone anywhere). The QFT is not a subroutine that hands you a Fourier spectrum to inspect — it is a change of basis, useful only when the thing you care about becomes a measurable label in the new basis. The rest of the chapter is two spectacular instances of exactly that: a hidden phase, and a hidden period, each engineered to land on a basis state the QFT can expose.

Here is the problem the QFT was born to solve. You are given a unitary $U$ and one of its eigenstates $|u⟩$, so that

and you may apply $U$ (and its powers) as a controlled operation. The task is to estimate the phase $\phi$. The trick is phase kickback from Chapter 6: put a register of $t$ counting qubits into the uniform superposition with Hadamards, and apply $U^{2^{t-1-w}}$ controlled on counting qubit $w$. Because $|u⟩$ is an eigenstate, each controlled application does nothing to the target but stamps a phase onto the control, leaving

Look at the counting register: whenever $2^t\phi$ is an integer $k_0$, this is exactly the state $\mathrm{QFT}|k_0⟩$, so applying the inverse QFT and measuring returns $k_0=2^t\phi$ with probability 1 — $\phi$ read off in binary. For a generic $\phi$ the register holds the same geometric sum with a non-integer “frequency” — not the QFT of any basis state — and the inverse QFT concentrates its weight near $2^t\phi$ instead of on a single value. The concentration is quantitative: writing ${\delta }_k=\phi -k/2^t$ for the offset of outcome $k$ from the true phase, a two-line geometric sum gives the outcome distribution

the Fejér kernel — the squared, normalized Dirichlet kernel, and the one that is everywhere non-negative, as a probability must be. A sharp central lobe of width one bin around $2^t\phi$, with sidelobes falling off like $1/{\delta }_k^2$. From (7.6) one can prove the best outcome $\hat{\phi }=k/2^t$ satisfies

Each extra counting qubit doubles the precision; a few repetitions beat the constant success probability down to nothing. Drive it yourself — move the true phase and watch the register concentrate, spiking to a single bar exactly when $\phi$ becomes a $t$-bit fraction:

What you are watching is the Fejér kernel (7.6) — the same ${\sin }^2(N\theta )/{\sin }^2\theta$ that governs a diffraction grating. Constructive interference piles amplitude onto the bin nearest $\phi$; destructive interference scours the rest. Phase estimation is the engine under the hood of an astonishing amount of quantum computing — chemistry energies, linear-systems solvers, and, next, period-finding.

Fix a composite $N$ and a number $a$ coprime to it. The powers $a,a^2,a^3,\dots (\mathrm{mod}N)$ eventually cycle back to 1; the smallest exponent that does so is the order $r$,

Finding $r$ is believed to be classically hard, and it is exactly a period-finding problem — perfect for phase estimation. The unitary is modular multiplication, $U_a|y⟩=|ay\mathrm{mod}N⟩$. Where are its eigenstates? Follow the orbit of 1: $U_a$ walks the length-$r$ cycle $|1⟩\to |a⟩\to |a^2⟩\to \cdots \to |a^{r-1}⟩\to |1⟩$, and the eigenvectors of a cyclic shift are its Fourier modes: for each $s=0,\dots ,r-1$ define

One line verifies the claim: $U_a$ shifts $|a^k⟩\to |a^{k+1}⟩$ inside the sum, and re-indexing $k+1\to k$ pulls out a single factor, $U_a|u_s⟩=e^{2\pi is/r}|u_s⟩$ — an eigenstate with eigenphase exactly $s/r$. (To be precise, these $r$ states span only the orbit of 1; on the rest of the space $U_a$ has other eigenvector blocks with other phases, but a computation started from $|1⟩$ never leaves the orbit, so those blocks never enter.) And the state we can actually prepare is their equal superposition: sum (7.9) over $s$ and the phases cancel except at $k=0$, leaving $|1⟩=\frac{1}{\sqrt{r}}{\sum }_{s=0}^{r-1}|u_s⟩$. So phase estimation on $U_a$ started from $|1⟩$ returns a measured value $k\approx 2^{t}s/r$ for an $s$ drawn uniformly at random.

Turning a $t$-bit estimate of $s/r$ into $r$ itself is a classical theorem about continued fractions (Legendre): a fraction that approximates an unknown rational closely enough is forced to show up among its convergents,

The counting register exists at the size it does precisely to trigger this criterion. What (7.10) needs is $2^{-(t+1)}<1/(2r^2)$, that is $2^{t+1}>2r^2$; and since $r<N$ it suffices to demand $2^{t+1}>2N^2$, i.e. $t>2{\log }_{2}N$. For $N=15$ that is $t\ge 8$ and for $N=21$ it is $t\ge 9$ — exactly the registers the lab uses. (You will often see the rule quoted as $t\ge 2{\log }_{2}N+1$, which for these two $N$ asks for 9 and 10 instead. That form rounds up and carries a spare bit of slack, bought to lift the success probability rather than to satisfy Legendre; both register sizes are valid, and the lab takes the tighter one so the state vectors stay small.) Checking each convergent denominator $q$ against $a^q\stackrel{?}{\equiv }1(\mathrm{mod}N)$ then certifies $r$. That test is exact, but not for the reason one might guess: the convergent denominators are not bounded by $r$. In the lab's own $N=21$, $r=6$, $t=9$ run, the expansion $85/512=[0;6,42,2]$ has denominators 1, 6, 253, 512 — two of them larger than 6. What makes the filter exact is that (7.10) puts $s/r$ in lowest terms among the convergents, with denominator $r/\gcd (s,r)$, a divisor of $r$; and $a^q\equiv 1$ holds exactly when $r$ divides $q$. So the test throws out every proper divisor of $r$ and accepts only multiples of it, which is $r$ itself whenever $\gcd (s,r)=1$. The explicit $q\le N$ guard, not any bound on the denominators, is what keeps the search finite. And when $s$ shares a factor with $r$ so that no $q$ passes, another measurement supplies one that does.

The last step is pure classical number theory, and it is the reason factoring falls. Suppose $r$ is even and $a^{r/2}\not\equiv -1(\mathrm{mod}N)$. Then $a^r-1=(a^{r/2}-1)(a^{r/2}+1)\equiv 0$, so $N$ divides the product but (by the two conditions) neither factor alone — which means

are nontrivial factors of $N$, computable in a blink by Euclid's algorithm. If the two conditions fail — odd $r$, or $a^{r/2}\equiv -1$ — you simply pick another $a$ and retry; one can show a random $a$ works with probability at least $1/2$. That is the whole of Shor's algorithm: a classical shell (pick $a$, take gcds) around a single quantum subroutine (find the order $r$). Worked concretely: for $N=15$, $a=7$ the order is $r=4$, so $a^{r/2}=49\equiv 4$, and $\gcd (3,15)=3$, $\gcd (5,15)=5$. There is 15.

Be precise about the claim. Shor's algorithm is a genuine super-polynomial speedup over the best known classical factoring algorithm (the number field sieve, sub-exponential but super-polynomial). That is not a proven separation — no one has shown factoring is classically hard, only that decades of effort have not cracked it — but it is a real and enormous advantage against the state of the art, and it is structured: it works because factoring hides a periodicity that the Fourier transform is built to find. Not every problem has such structure, as the very next chapter's unstructured search will show.

And the honest caveat. The circuits here run on a laptop for $N=15$ and $N=21$ because those need only a dozen or so ideal qubits. Factoring a 2048-bit RSA modulus needs thousands of logical qubits and billions of gates run coherently — which, on hardware with today's error rates, means millions of physical qubits marshalled by the error correction of Part III. Present-day NISQ devices — noisy, a few hundred physical qubits, no fault tolerance — cannot run it at any cryptographically interesting size. The demonstrations you see factoring 15 are real, and they are also toys. The gap between the toy and the threat is precisely the subject of Part III, and it is why the field's most important engineering problem is not a better algorithm but a longer-lived qubit.

The QFT is one short function — a Hadamard on each wire, a controlled phase rotation for each later wire, and a reversal by swaps. This is the circuit of Figure 7.1, verbatim, including the sign convention that pins (7.1):

1/// The Quantum Fourier Transform on a contiguous list of wires, given
2    /// MOST-significant first. Hadamard on each wire, then controlled phase
3    /// rotations R_m = diag(1, e^{2*pi*i/2^m}) from every less-significant
4    /// wire, then a reversal by swaps. O(n^2) gates. Sign convention:
5    /// |j> -> (1/sqrt(N)) sum_k e^{+2*pi*i*j*k/N} |k>.
6    fn qft(&mut self, w: &[usize]) {
7        let n = w.len();
8        for i in 0..n {
9            self.apply_1q(&g_h(), w[i]);
10            for j in (i + 1)..n {
11                let m = (j - i + 1) as u32;
12                let angle = 2.0 * PI / (1u64 << m) as f64;
13                self.cphase(w[j], w[i], angle);
14            }
15        }
16        for i in 0..n / 2 {
17            self.swap_wires(w[i], w[n - 1 - i]);
18        }
19    }

The strongest referee on it is a two-path cross-check in the spirit of Chapter 2: build the explicit $N\times N$ DFT matrix $\frac{1}{\sqrt{N}}[e^{+2\pi ijk/N}]$, apply it to random states, and demand the circuit agree with it. The inverse circuit gets the same treatment against the conjugate matrix — and that reference is the load-bearing part, because the obvious check is worthless. The lab's iqft is qft reversed with negated angles, so a round-trip $\mathrm{iQFT}(\mathrm{QFT}|\psi ⟩)=|\psi ⟩$ would undo any bug the two circuits share and still print machine dust. Delete the swap layer from both circuits and that is exactly what happens: the round-trip stays down at machine epsilon while the conjugate-matrix referee jumps to order 1. Both achieved values are in the panel below. Phase estimation is the QFT run backwards on a kicked-back phase, in full:

1/// Return the probability distribution over the t-bit counting register after
2/// phase estimation of a phase gate with eigenphase phi. Layout: t counting
3/// wires (0..t) then 1 target wire (t), target prepared in |1> (its own
4/// eigenstate). Since |1> is an eigenstate the counting register ends pure, so
5/// P(k) = |amp[2k+1]|^2 (target bit is the single low bit).
6fn phase_estimate_dist(t: usize, phi: f64) -> Vec<f64> {
7    let n = t + 1;
8    let target = t; // the last wire
9    let mut s = State::basis(n, 1); // index 1 = target (LSB) is |1>, rest |0>
10    for w in 0..t {
11        s.apply_1q(&g_h(), w);
12    }
13    // Controlled-U^{2^{t-1-w}}: since the target holds |1>, U^p contributes a
14    // phase e^{2*pi*i*phi*p}. Apply it as a controlled phase on (w, target).
15    for w in 0..t {
16        let p = 1u64 << (t - 1 - w);
17        let angle = 2.0 * PI * phi * p as f64;
18        s.cphase(w, target, angle);
19    }
20    let counting: Vec<usize> = (0..t).collect();
21    s.iqft(&counting);
22    (0..(1usize << t)).map(|k| s.amps[2 * k + 1].norm_sqr()).collect()
23}

And Shor is a thin classical wrapper — measure, run continued fractions on what came back, then take two gcds — around that quantum subroutine:

1/// Recover the order r from a measured phase k/2^t: among the CF convergents
2/// with denominator <= N, return the smallest positive q with a^q = 1 mod N.
3///
4/// Why that filter is exact. Legendre's criterion guarantees that s/r IN LOWEST
5/// TERMS is one of the convergents of k/2^t; its denominator is r/gcd(s,r), a
6/// divisor of r. And a^q = 1 mod N holds exactly when r divides q. So the a^q
7/// test rejects every convergent whose denominator is a PROPER divisor of r and
8/// accepts only multiples of r — which means that when gcd(s,r) = 1 the first
9/// accepted denominator IS r. (What it does NOT rest on is any bound q <= r:
10/// convergent denominators routinely exceed r. In this lab's own N=21, r=6,
11/// t=9 case, 85/512 = [0;6,42,2] has denominators 1, 6, 253, 512 — two of them
12/// larger than 6. The search is bounded by the explicit q <= N guard below.)
13fn order_from_measurement(k: u64, two_t: u64, a: u64, n: u64) -> Option<u64> {
14    for (_p, q) in cf_convergents(k, two_t) {
15        if q > 0 && q <= n && modpow(a, q, n) == 1 {
16            return Some(q);
17        }
18    }
19    None
20}
21
22// …
23
24/// Order-finding as a real machine would run it: draw an outcome k from the
25/// quantum distribution, hand k/2^t to the continued fractions, and repeat
26/// until CONFIRM runs have certified a candidate or the shot budget is spent.
27/// Returns the smallest certified candidate and the number of runs consumed.
28///
29/// Nothing here inspects the distribution as a whole. Every bit of information
30/// arrives through sampled measurements, so wrong amplitudes cost real shots —
31/// which is the whole point. The previous version of this routine scanned every
32/// outcome whose probability exceeded 1e-6; for N=21, t=9 that is all 512 of
33/// them, and the continued-fraction filter alone then certifies r=6 even from a
34/// uniform distribution. It was refereeing classical number theory, not the
35/// quantum circuit. Sampling narrows, but does not close, that hole: fed a
36/// deliberately uniform distribution this routine still certifies r most of the
37/// time — it just burns more shots. The referees that fail DECISIVELY on flat
38/// amplitudes are the peak-offset and peak-mass checks below; they, not this
39/// routine, carry the amplitude sensitivity of the suite.
40fn recover_order_sampled(
41    dist: &[f64],
42    t: usize,
43    a: u64,
44    n_mod: u64,
45    rng: &mut StdRng,
46) -> Option<(u64, usize)> {
47    let two_t = 1u64 << t;
48    let mut best: Option<u64> = None;
49    let mut confirmed = 0usize;
50    for shot in 1..=MAX_SHOTS {
51        let k = sample_outcome(dist, rng.gen::<f64>());
52        if let Some(q) = order_from_measurement(k as u64, two_t, a, n_mod) {
53            best = Some(match best {
54                Some(b) => b.min(q),
55                None => q,
56            });
57            confirmed += 1;
58            if confirmed >= CONFIRM {
59                return Some((best.unwrap(), shot));
60            }
61        }
62    }
63    best.map(|r| (r, MAX_SHOTS))
64}
65
66// …
67
68/// The full Shor post-processing given a witness a and the order r found by the
69/// quantum subroutine. Returns Ok((f1,f2)) on success, or Err(reason) when the
70/// classical guards reject this witness (odd order, or a^{r/2} = -1 mod N).
71fn shor_factors(a: u64, n: u64, r: u64) -> Result<(u64, u64), String> {
72    if r % 2 != 0 {
73        return Err(format!("order r={r} is odd — retry with another a"));
74    }
75    let x = modpow(a, r / 2, n);
76    if x == n - 1 {
77        return Err(format!("a^(r/2) = -1 mod N — retry with another a"));
78    }
79    let f1 = gcd(x + 1, n);
80    let f2 = gcd(if x >= 1 { x - 1 } else { n - 1 }, n);
81    // pick two nontrivial factors
82    for (p, q) in [(f1, n / f1), (f2, n / f2)] {
83        if p > 1 && p < n && p * q == n {
84            return Ok((p.min(q), p.max(q)));
85        }
86    }
87    Err(format!("gcds gave only trivial factors ({f1}, {f2})"))
88}

The middle function is worth dwelling on, because an earlier version of this lab got it wrong in a way that is easy to miss and is worth recognizing elsewhere. It did not sample at all: it scanned every outcome whose probability exceeded a fixed threshold and handed each one to the continued fractions. For $N=21$, $t=9$ every one of the 512 outcomes clears that threshold, so the scan covered the whole register — and the number-theoretic filter alone then certifies $r=6$ from a uniform distribution. The referee built on it could not see the amplitudes at all; it was grading Euclid's algorithm, not the quantum circuit. Drawing Born-distributed outcomes one at a time and retrying — which is what a real machine forces on you anyway — puts the interference back inside the loop.

The lab's referees are written to JSON and read live by the panel below; nothing here is hardcoded. The QFT and its inverse are graded against the explicit DFT matrices. Phase estimation pins exact $t$-bit phases to machine precision, and for generic phases it is graded on the margin it leaves against the $2^{-(t+1)}$ bound and the $4/{\pi }^2$ floor — a hinged “violation” would report exactly zero whether the bound was cleared by a hair or by a mile, which measures nothing. Two further referees check where the order-finding distribution actually puts its weight: the $r$ most likely outcomes must sit exactly on the nearest $t$-bit fractions of the $s/r$, and the mass on them must clear the same $4/{\pi }^2$ floor. Those two are the ones that fail loudly — hundreds of bins off, and a negative margin — if you strip the inverse QFT out of the order-finder. Finally, Shor's factors for 15 and 21 are compared not against $N$ but against an independent trial division: comparing them to $N$ could never fail, since shor_factors only ever returns a pair whose product is $N$ already. A bad witness ($a=14$, where $a\equiv -1$) is shown being rejected by the guards and retried.

Run it yourself with cargo run --release in Rust-QML/ch07-fourier. You have now built, from amplitudes upward, the algorithm that made the world take quantum computing seriously — and seen exactly where its power comes from (periodicity plus interference) and exactly what it still lacks (the coherence of a fault-tolerant machine).