How Quantum Algorithms Work

Quantum computing introduces a paradigm shift in how we process information, leveraging the principles of superposition and interference to solve certain problems exponentially faster than classical computers. This article aims to clearly explain how a class of quantum algorithms exploit these principles. Familiarity with qubits, state vectors, the Hadamard gate, tensor products, and Dirac (bra–ket) notation is assumed. Readers looking to build or refresh this foundation may find the interactive essay “Quantum computing for the very curious” especially helpful.

Background

Linear functions

Boolean functions are mappings from $n$ Boolean variables to a Boolean output. A simple example is the $\text{NOT}$ function, which flips a bit: $\text{NOT}(0)=1$ and $\text{NOT}(1)=0$. We can represent this mapping using a truth table.

$$\begin{array}{c|c} x & \text{NOT}(x) \\ \hline 0 & 1 \\ 1 & 0 \end{array}$$

Some 2-variable functions include $\text{OR}$ (inclusive OR) and $\text{XOR}$ (exclusive OR). $\text{OR}$ outputs $1$ if any input is $1$ and $0$ otherwise. $\text{XOR}$ outputs $1$ if inputs differ and $0$ otherwise.

$$\begin{array}{cc|c|c} x_1 & x_2 & \text{OR}(x_1, x_2) & \text{XOR}(x_1, x_2) \\ \hline 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 \\ 1 & 0 & 1 & 1 \\ 1 & 1 & 1 & 0 \end{array}$$

For more than two inputs, exclusive OR is written using the $\oplus$ operator: $$\text{XOR}(x_1, x_2, \dots, x_n) = x_1 \oplus x_2 \oplus \dots \oplus x_n,$$ so $0 \oplus 1 = 1$ but $0 \oplus 1 \oplus 1 = 0$ and $0 \oplus 1 \oplus 1 \oplus 1 = 1$.

Each additional $1$ flips the result, so the order of the inputs does not matter. Thus $x_1 \oplus x_2 \oplus \dots \oplus x_n$ equals $1$ exactly when an odd number of inputs are $1$. This value is called the parity of the input string.

An $n$-variable Boolean function $f$ is linear if there exists some string $\mathbf{s}$ such that $f(\mathbf{x})$ equals $$ \mathbf{x} \boldsymbol{\cdot} \mathbf{s} \text{ mod } 2 = \underbrace{\begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_{n-1}\\ x_n \end{pmatrix}}_{\mathbf{x}} \boldsymbol{\cdot} \underbrace{\begin{pmatrix} s_1\\ s_2\\ \vdots\\ s_{n-1}\\ s_n \end{pmatrix}}_{\mathbf{s}} \enspace \text{mod} \enspace 2 = s_1 x_1 + s_2 x_2 + \dots + s_n x_n \enspace \text{mod} \enspace 2,$$ which can be expressed as $s_1 x_1 \oplus s_2 x_2 \oplus \dots \oplus s_n x_n$ because every linear function is the XOR of a selected subset of its inputs.

Is the inclusive OR function a linear function?

No. A 2-variable linear function must be of the form $f(x_1, x_2) = s_1 x_1 \oplus s_2 x_2$. For the $\text{OR}$ function, we see that $\text{OR}(0,1) = 1$ (implying $s_2 = 1$) and $\text{OR}(1,0) = 1$ (implying $s_1 = 1$). This would mean $\text{OR}(x_1, x_2) = x_1 \oplus x_2$. However, $\text{OR}(1,1) = 1$, while $1 \oplus 1 = 0$. Since the $\text{OR}$ function doesn't match any possible XOR-sum of its inputs, it is not a linear function.

While truth tables are useful for small functions, they obscure structure. Karnaugh maps arrange inputs so adjacent cells differ by one bitThe row and column headers of a Karnaugh map are ordered using a Gray code, where each successive value differs by exactly one bit. Consequently, any horizontal or vertical move between adjacent cells corresponds to flipping exactly one input variable.. In this layout, linear functions appear as simple, repeating geometric patterns.

All eight 3-variable linear functions and their Karnaugh maps, with corresponding product coverings highlighted.
Draw Karnaugh maps for the $4$-variable linear functions corresponding to the hidden strings $0110$ and $1111$.

The Hadamard transform

The Hadamard transform is the application of a series of Hadamard gates to all qubits in a quantum register. Equivalently, it applies $H^{\otimes n}$ on an $n$-qubit system, mapping computational basis states into equal superpositions (with specific relative phases).

General Hadamard

The single-qubit Hadamard transform is simply the Hadamard gate $H$ applied to a qubit.

Single Qubit Hadamard

The Hadamard gate $H$ is represented by the matrix $$ H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}, $$ so it maps $|0\rangle$ to $$\frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix} = \frac{|0\rangle + |1\rangle}{\sqrt{2}} = |+\rangle.$$ and $|1\rangle$ to $$\frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = \begin{pmatrix} 1/\sqrt{2} \\ -1/\sqrt{2} \end{pmatrix} = \frac{|0\rangle - |1\rangle}{\sqrt{2}} = |-\rangle.$$ We can unify these two cases into a single algebraic expression by using the term $(-1)^x$ to control the sign of the $|1\rangle$ component: $$H|x\rangle = \frac{|0\rangle+(-1)^{x}|1\rangle}{\sqrt{2}},$$ which is equivalent to $$\frac{1}{\sqrt{2}} \sum_{y \in \{0, 1\} } (-1)^{x y} |y\rangle. $$

The two-qubit Hadamard transform is equivalent to applying the Hadamard transform to each qubit independently.

Two Qubit Hadamard

We start by splitting the operator across the tensor product: $$ \begin{aligned} |\psi\rangle &= H|x_1\rangle \otimes H|x_2\rangle\\ &= \left( \frac{|0\rangle+(-1)^{x_1}|1\rangle}{\sqrt{2}} \right) \otimes \left( \frac{|0\rangle+(-1)^{x_1}|1\rangle}{\sqrt{2}} \right). \end{aligned} $$ When we distribute the tensor product, the coefficients multiply, resulting in $$ \frac{1}{2} \left(|00\rangle+(-1)^{x_2}|01\rangle+(-1)^{x_1}|10\rangle+(-1)^{x_1 + x_2}|11\rangle\right), $$ which is equivalent to $$\frac{1}{2} \sum_{y_1 \in \{0, 1\} } \sum_{y_2 \in \{0, 1\}} (-1)^{x_1 y_1 + x_2 y_2} |y_1 y_2\rangle. $$

In the general case, we apply the Hadamard gate to n qubits simultaneously. With the input state $|x\rangle = |x_1 x_2 \dots x_n\rangle$, we have $$ \begin{aligned} H^{\otimes n} |x_1 \ldots x_n\rangle &= \bigotimes_{i=1}^n H |x_i\rangle \\ &= \bigotimes_{i=1}^n \left( \frac{1}{\sqrt{2}} \sum_{y_i \in \{0,1\}} (-1)^{x_i y_i} |y_i\rangle \right)\\ &= \frac{1}{\sqrt{2^n}} \bigotimes_{i=1}^n \left( \sum_{y_i \in \{0,1\}} (-1)^{x_i y_i} |y_i\rangle \right). \end{aligned} $$

We now expand the tensor product using linearity, one tensor slot at a time. Applying linearity in the first slot gives: $$ \begin{aligned} |\psi\rangle &= \frac{1}{\sqrt{2^n}} \sum_{y_1 \in \{0,1\}} (-1)^{x_1 y_1} \left( |y_1\rangle \otimes \bigotimes_{i=2}^n \sum_{y_i \in \{0,1\}} (-1)^{x_i y_i} |y_i\rangle \right). \end{aligned} $$

Applying linearity in the second slot gives: $$ \begin{aligned} |\psi\rangle &= \frac{1}{\sqrt{2^n}} \sum_{y_1 \in \{0,1\}} \sum_{y_2 \in \{0,1\}} (-1)^{x_1 y_1} (-1)^{x_2 y_2} \\ &\quad \times \left( |y_1\rangle \otimes |y_2\rangle \otimes \bigotimes_{i=3}^n \sum_{y_i \in \{0,1\}} (-1)^{x_i y_i} |y_i\rangle \right). \end{aligned} $$

Continuing this process for the third slot, fourth slot, and so on, and exhausting all n tensor factors, we obtain: $$ \begin{aligned} |\psi\rangle &= \frac{1}{\sqrt{2^n}} \sum_{y_1,\dots,y_n \in \{0,1\}} \left( \prod_{i=1}^n (-1)^{x_i y_i} \right) \left( |y_1\rangle \otimes \cdots \otimes |y_n\rangle \right)\\ &= \frac{1}{\sqrt{2^n}} \sum_{y \in \{0,1\}^n} \left( \prod_{i=1}^n (-1)^{x_i y_i} \right) |y\rangle\\ &= \frac{1}{\sqrt{2^n}} \sum_{y \in \{0,1\}^n} (-1)^{\sum x_i y_i} |y\rangle. \end{aligned} $$

Finally, we recognize that the sum in the exponent, $\sum x_i y_i$, is exactly the definition of the bitwise dot product modulo 2 (denoted as $x \cdot y$). This gives us the standard compact form of the Hadamard transform: $$ H^{\otimes n}|x\rangle = \frac{1}{\sqrt{2^n}}\sum_{y \in \{0,1\}^n} (-1)^{x \cdot y}|y\rangle. $$ We recognize bitwise dot product modulo 2 as a constraint for a linear function! For any given basis state $\ket{y}$, the sign of the amplitude is positive when $f_x(y)=0$ and negative otherwise.

Phase Oracle

A quantum oracle $U_f$ is a black box that implements the function $f: \{0,1\}^n \rightarrow \{0,1\}$ by acting on the output qubit. The oracle performs a unitary transformation defined by $U_f\ket{x}\ket{y} = \ket{x}\ket{y \oplus f(x)}$The XOR construction ensures the operation is unitary and reversible, a strict requirement for all quantum gates..

An oracle query (or oracle call) refers to the operation of invoking this quantum oracle $U_f$ on a quantum state $\ket{x}\ket{y}$. Specifically, during an oracle query, the input state $\ket{x}$ is processed along with an ancillary qubit $\ket{y}$ such that the output qubit is updated according to the function $f.$

Bernstein-Vazirani algorithm

Given a linear black box Boolean function $f: \{0, 1\}^n \rightarrow \{0, 1\}$, the Bernstein-Vazirani problem is to find the hidden string $\mathbf{s}$ in the fewest evaluations of $f$Ethan Bernstein and Umesh Vazirani formulated this problem in 1993, demonstrating one of the first cases of quantum speedup..

In general, a classical computer requires $n$ oracle calls to find the hidden string of an $n$-variable function. This can be accomplished by evaluating the function at all inputs where a single variable is set to $1$ and all others are set to $0$. For example, a three-variable linear function $f(x_1, x_2, x_3)$ equals $ s_1 x_1 \oplus s_2 x_2 \oplus s_3 x_3.$ Since $f(1,0,0)= s_1$, $f(0,1,0)= s_2$, and $f(0,0,1)= s_3$, we can determine $s_1, s_2, s_3$ with three oracle calls.

All eight 3-variable linear functions and their Karnaugh maps, with one-hot inputs highlighted.

In 1992, Ethan Bernstein and Umesh Vazirani proposed a quantum algorithm to find the hidden string with a single oracle query.[2]

First, the circuit is initialized to the state $\ket{\Psi_0} = \ket{0}^{\oplus n}\ket{1}.$ Second, a Hadamard transform is applied on all qubits, resulting in the state $\ket{\Psi_1} = \ket{+}^{\oplus n}\ket{-}.$

In the next part of the algorithm, an oracle $U_f$ is applied to the state $\ket{\Psi_1}$.

To analyze this, consider the state of an oracle is applied to an output bit initialized to the $\ket{-}$ state.

We see that $\ket{\Psi}$ is $$ \begin{aligned} U_f(\ket{x} \otimes \ket{-}) &= U_f \left(\ket{x} \otimes \left(\frac{\ket{0}-\ket{1}}{\sqrt{2}}\right)\right)\\ &= \frac{1}{\sqrt{2}} U_f \left(\ket{x}\ket{0}-\ket{x}\ket{1}\right)\\ &= \frac{1}{\sqrt{2}} \left(U_f(\ket{x}\ket{0})-U_f(\ket{x}\ket{1})\right)\\ &= \frac{1}{\sqrt{2}} \left(\ket{x}\ket{f(x)}-\ket{x}\ket{1 \oplus f(x)})\right)\\ &= \ket{x}\frac{1}{\sqrt{2}} \left(\ket{f(x)}-\ket{1 \oplus f(x)})\right)\\ &= \ket{x}\begin{cases} \frac{1}{\sqrt{2}} \left(\ket{0}-\ket{1})\right) & \text{if } f(x)=0\\ \frac{1}{\sqrt{2}} \left(\ket{1}-\ket{0})\right) & \text{if } f(x)=1 \end{cases}\\ &= \ket{x}\begin{cases} \ket{-} & \text{if } f(x)=0\\ -\ket{-} & \text{if } f(x)=1 \end{cases}\\ &= (-1)^{f(x)}\ket{x}\ket{-}. \end{aligned} $$

Now, assume a prior Hadamard transform on the upper register.

We see that $\ket{\Psi_2}$ is $$ \begin{aligned} U_f\left(\ket{+}^{\otimes n} \otimes \ket{-}\right) &= U_f\left(\left(\frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} \ket{x})\right) \otimes \ket{-}\right)\\ &= \frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} U_f\left(\ket{x} \otimes \ket{-}\right)\\ &= \frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} (-1)^{f(x)}\ket{x}\ket{-}. \end{aligned} $$

Therefore, the resulting quantum state $\ket{\Psi_2}$ is a uniform superposition of all basis states, but the phase of all basis states for which $f=1$ is negative: $$\ket{\Psi_2}=U_f\left(\ket{\Psi_1}\right)=\frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} (-1)^{f(x)}\ket{x}\ket{-}.$$ Measurement is dependent only on amplitude squared, not phase. Therefore, measuring the state $\ket{\Psi_2}$ is equivalent to $\ket{\Psi_1},$ which would yield a random $\ket{x}.$

Finally, discarding the last qubit and applying yet another Hadamard transform onto the $n$ input qubits yields $$ \begin{aligned} \ket{\Psi_3} &= H^{\otimes n} \left(\frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} (-1)^{f(x)}\ket{x}\right)\\ &= \frac{1}{\sqrt{2^n}}\sum_{y \in \{0,1\}^n}(-1)^{x \cdot y} \left(\frac{1}{\sqrt{2^n}} \sum_{x \in \{0,1\}^n} (-1)^{f(x)}\right)\ket{y}\\ &= \sum_{y \in \{0,1\}^n}\underbrace{\frac{1}{2^n}\sum_{x \in \{0,1\}^n} (-1)^{x \cdot y+f(x)}}_{\text{amplitude of }\ket{y}}\ket{y}. \end{aligned} $$

This expression is a linear combination of all the possible basis states. The amplitude of each basis state is dependent on the similarity between the oracle function and the linear function encoded by the basis state.

It turns out that when the oracle function is linear, the amplitude of all basis states that encode incorrect linear functions is zero. This is because any two different $n$-variable linear functions will have $2^{n-1}$ similar outputs and $2^{n-1}$ different outputs. Conversely, the amplitude of the state encoding the correct linear function is $1$.

Deutsch-Jozsa algorithm

In the Deutsch-Jozsa problemFirst published in 1992, this was the first algorithm to show an exponential speedup over classical deterministic algorithms., we are given a Boolean function $f: \{0, 1\}^n \rightarrow \{0, 1\}.$ The function is known to be either constant or balanced:

The objective of the Deutsch-Jozsa problem is to determine whether the function is constant or balanced using the fewest number of evaluations of the function.

Let us apply the setup of the Bernstein-Vazirani algorithm to the Deutsch-Jozsa problem using the oracle for $f$. Consider the amplitude of the basis state $\ket{0}^{\oplus n}$: $$\frac{1}{2^n}\sum_{x \in \{0,1\}^n} (-1)^{f(x)}.$$ There are two cases:

This is a profound result: we can determine whether a function is constant or balanced using only $1$ evaluation of the function. By contrast, the classical algorithm would require at least $2^{n-1} + 1$ evaluations of the function to determine whether it is constant or balanced. Therefore, the Deutsch-Jozsa algorithm provides an exponential speedup over the best-possible classical algorithm.[3]

Simon's algorithm

Given a Boolean function $f: \{0, 1\}^n \rightarrow \{0, 1\}^n$ and the promise that for all $p, q \in \{0, 1\}^n$: $$f(p) = f(q)$$ if and only if $p \oplus q \in \{0^n, s\}$ for some hidden string $s = \{s_1, s_2, \ldots, s_n\},$ Simon's problem asks us to find $s.$

This promise means the function $f$ is either strictly one-to-one (if $s = 0^n$) or exactly two-to-one (if $s \neq 0^n$). When it is two-to-one, every possible output is generated by exactly two inputs. These paired inputs are separated by a constant bitwise $\text{XOR}$ mask, which is the hidden string $s$. A classical computer requires an exponential number of queries to find two matching outputs. Simon's algorithm solves this with an exponential speedup by leveraging quantum interference to reveal $s$ in $O(n)$ queries.

Simon's problem and its quantum solution was conceived by Dan Simon in 1993.[4]

The algorithm begins by initializing two $n$-qubit registers to the zero state, $\ket{\Psi_0} = \ket{0}^{\oplus n}\ket{0}^{\oplus n}.$ A Hadamard transform is applied to the first register, creating an equal superposition of all possible $2^n$ binary inputs. This prepares the system to evaluate the function globally: $$ \ket{\Psi_1} = \left(\frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} \ket{x}\right) \ket{0^n}. $$ Next, the oracle operation $U_f$ is applied. This operation computes $f(x)$ for every state in the superposition and stores the result in the second register, inextricably entangling the input and output states: $$ \ket{\Psi_2} = \frac{1}{\sqrt{2^n}}\sum_{x \in \{0,1\}^n} \ket{x}\ket{f(x)}. $$

If we measure the second register, it collapses to a single, random output value $z = f(x)$. Because the registers are entangled, this measurement forces the first register to collapse into the specific input states compatible with $z$. Assuming $s \neq 0^n$, the problem's constraints guarantee that exactly two inputs, $p$ and $q$, map to $z$. The first register thus becomes an even superposition of those two inputs: $$ \ket{\Psi_3} = \frac{1}{\sqrt{2}} (\ket{p} + \ket{q}). $$

Measuring the first register at this stage would simply yield either $p$ or $q$ at random, which provides no immediate information about the relationship $p \oplus q = s$. To extract the hidden string $s$, we must interfere the states $\ket{p}$ and $\ket{q}$. Applying a second Hadamard transform to the first register achieves this by rewriting the state in a basis where the phase of each possible measurement outcome $\ket{y}$ depends on its dot product with $p$ and $q$: $$ \frac{1}{\sqrt{2}}(H^{\otimes n}\ket{p} + H^{\otimes n}\ket{q}) = \frac{1}{\sqrt{2^{n+1}}} \sum_{y \in \{0,1\}^n} \left( (-1)^{p \cdot y} + (-1)^{q \cdot y} \right) \ket{y}. $$

The probability amplitude for measuring any specific state $\ket{y}$ is governed by the sum of the two phase terms: $(-1)^{p \cdot y} + (-1)^{q \cdot y}$. If the dot products $p \cdot y$ and $q \cdot y$ have different parities, these terms evaluate to $1$ and $-1$, destructively interfering to zero. The state $\ket{y}$ will only have a nonzero probability of being measured if the terms constructively interfere, which requires the parities to be identical: $$ p \cdot y \bmod 2 = q \cdot y \bmod 2. $$ By substituting the problem's core constraint, $q = p \oplus s$, we can isolate $s$ to see which $y$ values survive this interference: $$ \begin{aligned} p \cdot y \bmod 2 &= (p \oplus s) \cdot y \bmod 2 \\ \cancel{p \cdot y \bmod 2} &= \cancel{(p \cdot y \bmod 2)} \oplus (s \cdot y \bmod 2) \\ s \cdot y \bmod 2 &= 0 \\ y_1 s_1 \oplus y_2 s_2 \oplus \ldots \oplus y_{n-1} s_{n-1} \oplus y_n s_n &= 0. \end{aligned} $$ This equation shows that any measured string $y$ is guaranteed to be orthogonal to the hidden string $s$ modulo 2. While a single measurement of $y$ only provides one equation, running the entire quantum circuit $O(n)$ times will produce a set of $n-1$ linearly independent $y$-values. Classical Gaussian elimination can then be used to solve this system of equations and recover $s$.

As a concrete example, suppose the algorithm is run on a 3-qubit system and yields the $y$-values $001, 110,$ and $111.$ Because $y \cdot s = 0 \bmod 2$ for every measurement, these values correspond to the following system of linear equations: $$ \begin{cases} s_3 = 0\\ s_1 \oplus s_2 = 0\\ s_1 \oplus s_2 \oplus s_3 = 0 \end{cases} $$ This translates into a matrix equation over modulo $2$: $$ \begin{bmatrix} 0 & 0 & 1 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} s_1 \\ s_2 \\ s_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}. $$ To solve for $s$, we perform row operations using $\text{XOR}$ for addition: $$ \begin{aligned} & \left[ \begin{array}{ccc|c} 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \end{array} \right] \\[.5em] & \hspace{1.2em} \Big\downarrow \text{\scriptsize swap 1 and 2} \\[.5em] & \left[ \begin{array}{ccc|c} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 1 & 1 & 0 \end{array} \right] \\[.5em] & \hspace{1.2em} \Big\downarrow \text{\scriptsize subtract 1 from 3} \\[.5em] & \left[ \begin{array}{ccc|c} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 \end{array} \right] \\[.5em] & \hspace{1.2em} \Big\downarrow \text{\scriptsize subtract 2 from 3} \\[.5em] & \left[ \begin{array}{ccc|c} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right] \end{aligned} $$ The reduced matrix yields the simplified equations: $$ \begin{cases} s_1 \oplus s_2 = 0 \\ s_3 = 0 \end{cases} $$ The first equation implies $s_1 = s_2$, and the second explicitly gives $s_3 = 0$. This leaves two valid configurations for $s$: $000$ or $110$. By definition, if the function is two-to-one, $s$ cannot be $000$. Therefore, the hidden string is $110$.

The Fourier transform

Discrete Fourier transform

An $N$th root of unity, where $N$ is a positive integer, is any complex number $z$ that satisfies the equation: $$z^N = 1.$$ Because we are working in the complex plane, there are always exactly $N$ distinct solutions to this equation. Using Euler's formula, these roots can be expressed as $$z_k = \exp\left(\frac{2\pi i k}{N}\right) = \cos\left(\frac{2\pi k}{N}\right) + i\sin\left(\frac{2\pi k}{N}\right)$$ where $k = 0, 1, 2, \dots, N-1$.

Geometrically, if you draw a unit circle on the complex plane, the $N$th roots of unity are $N$ equally spaced points along that circle, always starting at $z = 1$ (when $k = 0$).

Now, consider the following sum: $$\sum_{k=0}^{N-1} z_k^j = \sum_{k=0}^{N-1} (e^{2\pi i k / N})^j.$$ Depending on the value of $j$, the sum evaluates to one of two results.

When $j \neq 0$, you're summing $N$ equally-spaced points on the unit circle, so they cancel perfectly by symmetry. When $j = 0$, every term is $1$, so you get $N$. Therefore, we have $$ \sum_{k=0}^{N-1} e^{2\pi i k j / N} = \begin{cases} N & \text{if } j \equiv 0 \pmod N \\ 0 & \text{otherwise} \end{cases}. $$ This idea is the central to everything that follows.

The discrete Fourier transform (DFT) maps a time-domain sequence $x = (x_0, x_1, \ldots, x_{N-1})$ to its corresponding frequency-domain sequence $X = (X_0, X_1, \ldots, X_{N-1})$ $$X_k = \sum_{n=0}^{N-1} x_n \, \omega_{N}^{-nk},$$ where $\omega_{N}$ is the primitive $N$th root of unity: ${\displaystyle \omega_{N}=e^{\frac{i 2\pi}{N}}}$ and $(\omega_{N})^N=1$.

Suppose you have a discrete time-domain sequence $x = (x_0, x_1, \ldots, x_{N-1})$ and you're told that it is a "pure signal," meaning that $$ x_k = c \cdot e^{2\pi i u k / N}$$ for some complex constant $c$ and integer $u$ with $0 \le u < N$. The value $u$ is called the frequency of the signal. Visually, a pure discrete signal is nothing more than a single point jumping around a circular orbit with a relentlessly steady rhythm.

Consider the sum over all samples, tested against frequency $j$: $$ X_j = \sum_{k=0}^{N-1} x_k \, e^{-2\pi i j k / N} = c \sum_{k=0}^{N-1} e^{2\pi i k (u - j) / N}. $$ By the root-of-unity fact, $X_j = cN$ when $j = u$ and $X_j = 0$ for every other frequency. A pure signal at frequency $u$ produces a single spike in the frequency domain.

A general periodic signal with period $r$ is one where the samples repeat every $r$ steps: $x_k = x_{k+r}$ for all $k$. Since there are $N$ total samples and the pattern tiles exactly, we need $r$ to divide $N$.

How does periodicity relate to pure signals? A pure signal $e^{2\pi i u k / N}$ completes one full cycle when $k$ advances by $N/u$ steps, so its period is $N/u$. A pure signal fits inside a period-$r$ pattern precisely when $N/u$ divides $r$, i.e. when $u$ is a multiple of $N/r$. The allowed frequencies are then $$u \in \left\{0,\; \frac{N}{r},\; \frac{2N}{r},\; \dots,\; \frac{(r-1)N}{r}\right\},$$ giving exactly $r$ distinct values. This count is no coincidence: a period-$r$ signal is completely determined by its $r$ values in one period, so it has exactly $r$ degrees of freedom. The $r$ allowed pure signals match this exactly, and they are mutually orthogonal: as we saw above, testing one pure signal against a different frequency causes the samples to spread evenly around the unit circle and cancel to zero. Each pure signal is therefore invisible to all the others, making them independent. With $r$ independent signals covering $r$ degrees of freedom, they span the entire space of period-$r$ signals, and any such signal can be expressed as a unique combination of them.

By linearity of the DFT, and since each pure signal at frequency $u$ contributes a single spike at $X_u$ and zero everywhere else, the DFT of any period-$r$ signal is non-zero only at those $r$ multiples of $N/r$. The remaining $N - r$ coefficients are exactly zero.

Quantum Fourier transform

The quantum Fourier transform (QFT) applies the discrete Fourier transform to a quantum state, which usually has dimension $N=2^n$ where $n$ is the total number of qubits.[1] It maps the computational basis state $\ket{j}=\ket{j_1 j_2 \ldots j_n}$ to its corresponding Fourier basis state $$F_N \ket{j} = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \omega_{N}^{jk} \ket{k},$$ where $\omega_{N}$ is an $N$th root of unity: ${\displaystyle \omega_{N}=e^{\frac{i 2\pi}{N}}}$ and $(\omega_{N})^N=1$.

What is the unitary matrix that defines the single-qubit quantum Fourier transform?

Since $n=1,$ we have $N=2^n=2^1=2.$ Therefore, we get $$F_2 \ket{0} = \frac{1}{\sqrt{2}} \sum_{k=0}^{1} \omega_{2}^{0(k)} \ket{k} = \frac{1}{\sqrt{2}}\left(e^{2\pi i{\frac {0 \cdot 0}{2}}}\ket{0}+e^{2\pi i{\frac {0 \cdot 1}{2}}}\ket{1}\right) = \frac{1}{\sqrt{2}}\left(\ket{0}+\ket{1}\right)$$ and $$F_2 \ket{1} = \frac{1}{\sqrt{2}} \sum_{k=0}^{1} \omega_{2}^{1(k)} \ket{k} = \frac{1}{\sqrt{2}}\left(e^{2\pi i{\frac {1 \cdot 0}{2}}}\ket{0}+e^{2\pi i{\frac {1 \cdot 1}{2}}}\ket{1}\right) = \frac{1}{\sqrt{2}}\left(\ket{0}-\ket{1}\right).$$ Thus, the single-qubit QFT performs the transformation $$\begin{pmatrix} \alpha_0\\ \alpha_1\end{pmatrix} \to \alpha_0\frac{1}{\sqrt{2}} \begin{pmatrix} 1\\ 1\end{pmatrix} + \alpha_1\frac{1}{\sqrt{2}} \begin{pmatrix} 1\\ -1\end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} \alpha_0\\ \alpha_1\end{pmatrix},$$ which is the matrix representation of the single-qubit Hadamard gate!

Theorem. Every abelian group is a direct product of cyclic groups.

For a finite abelian group $G = \mathbb{Z}_{N_1} \times \mathbb{Z}_{N_2} \times \cdots \times \mathbb{Z}_{N_m}$, its QFT is the tensor product ($\otimes$) of the QFTs on each cyclic component: $$F_G = F_{N_1} \otimes F_{N_2} \otimes \cdots \otimes F_{N_m}.$$

Describe the unitary transformation that defines the quantum Fourier transform over $G=\mathbb{Z}_{2}^n.$

For a finite abelian group $G=\mathbb{Z}_{2}^n = \overbrace{\mathbb{Z}_{2} \times \mathbb{Z}_{2} \times \cdots \times \mathbb{Z}_{2}}^{\text{$n$ times}},$ we have $F_G = \overbrace{F_{2} \otimes F_{2} \otimes \cdots \otimes F_{2}}^{\text{$n$ times}}.$ We know that $F_2$ is equivalent to a single-qubit Hadamard transformation $H$. Therefore, the quantum Fourier transform over the group $G=\mathbb{Z}_{2}^n$ is an $n$-qubit Hadamard transformation!

What is the unitary matrix that defines the two-qubit quantum Fourier transform?

Since $n=2$, we have $ N=2^n=2^2=4.$ Therefore, we get

$F_4 \ket{j}$ $\displaystyle \frac{1}{\sqrt{4}} \sum_{k=0}^{3} \omega_N^{jk} \ket{k}$ $\frac{1}{2}\left(\omega_N^{j \cdot 0}\ket{00}+\omega_N^{j \cdot 1}\ket{01}+\omega_N^{j \cdot 2}\ket{10}+\omega_N^{j \cdot 3}\ket{11}\right)$
$F_4 \ket{00}$ $\displaystyle \frac{1}{\sqrt{4}} \sum_{k=0}^{3} e^{i 2\pi \frac{0 \cdot k}{4}} \ket{k}$ $\frac{1}{2}\left(e^{2\pi i{\frac {0 \cdot 0}{4}}}\ket{00}+e^{2\pi i{\frac {0 \cdot 1}{4}}}\ket{01}+e^{2\pi i{\frac {0 \cdot 2}{4}}}\ket{10}+e^{2\pi i{\frac {0 \cdot 3}{4}}}\ket{11}\right)$
$F_4 \ket{01}$ $\displaystyle \frac{1}{\sqrt{4}} \sum_{k=0}^{3} e^{i 2\pi \frac{1 \cdot k}{4}} \ket{k}$ $\frac{1}{2}\left(e^{2\pi i{\frac {1 \cdot 0}{4}}}\ket{00}+e^{2\pi i{\frac {1 \cdot 1}{4}}}\ket{01}+e^{2\pi i{\frac {1 \cdot 2}{4}}}\ket{10}+e^{2\pi i{\frac {1 \cdot 3}{4}}}\ket{11}\right)$
$F_4 \ket{10}$ $\displaystyle \frac{1}{\sqrt{4}} \sum_{k=0}^{3} e^{i 2\pi \frac{2 \cdot k}{4}} \ket{k}$ $\frac{1}{2}\left(e^{2\pi i{\frac {2 \cdot 0}{4}}}\ket{00}+e^{2\pi i{\frac {2 \cdot 1}{4}}}\ket{01}+e^{2\pi i{\frac {2 \cdot 2}{4}}}\ket{10}+e^{2\pi i{\frac {2 \cdot 3}{4}}}\ket{11}\right)$
$F_4 \ket{11}$ $\displaystyle \frac{1}{\sqrt{4}} \sum_{k=0}^{3} e^{i 2\pi \frac{3 \cdot k}{4}} \ket{k}$ $\frac{1}{2}\left(e^{2\pi i{\frac {3 \cdot 0}{4}}}\ket{00}+e^{2\pi i{\frac {3 \cdot 1}{4}}}\ket{01}+e^{2\pi i{\frac {3 \cdot 2}{4}}}\ket{10}+e^{2\pi i{\frac {3 \cdot 3}{4}}}\ket{11}\right)$

where ${\displaystyle \omega_{N}=e^{\frac{i 2\pi}{N}}}.$ Thus, the two-qubit QFT performs the transformation $$\begin{pmatrix} \alpha_0\\ \alpha_1\\ \alpha_2\\ \alpha_3\end{pmatrix} \to \alpha_0\frac{1}{2} \begin{pmatrix} 1\\ 1\\ 1\\ 1\end{pmatrix} + \alpha_1\frac{1}{2} \begin{pmatrix} 1\\ i\\ -1\\ -i\end{pmatrix} + \alpha_2\frac{1}{2} \begin{pmatrix} 1\\ -1\\ 1\\ -1\end{pmatrix} + \alpha_3\frac{1}{2} \begin{pmatrix} 1\\ -i\\ -1\\ i\end{pmatrix},$$ which means $$ F_4 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{pmatrix} \begin{pmatrix} \alpha_0\\ \alpha_1\\ \alpha_2\\ \alpha_3\end{pmatrix}.$$

What is the unitary matrix that defines the quantum Fourier transform over the cyclic group $\mathbb{Z}_N$? To describe the matrix, we can find a formula for the $(l,j)$-th entry of the matrix $F_N.$ The matrix element $F_N(l,j)=\braket{l}{F_N | j}$ because column $j$ corresponds to vector $F_N \ket{j}$. By the definition of the QFT, we have $$F_N \ket{j} = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \omega_N^{jk} \ket{k} \implies \braket{l}{F_N | j}= \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \omega_N^{jl} \braket{l}{k}.$$ Therefore, the QFT matrix $F_N$ has elements given by $$F_N(l, j) = \frac{1}{\sqrt{N}} \omega_N^{lj}$$ and thus, we have $$ F_{N}={\frac {1}{\sqrt {N}}}{\begin{bmatrix}1&1&1&1&\cdots &1\\1&\omega &\omega ^{2}&\omega ^{3}&\cdots &\omega ^{N-1}\\1&\omega ^{2}&\omega ^{4}&\omega ^{6}&\cdots &\omega ^{2(N-1)}\\1&\omega ^{3}&\omega ^{6}&\omega ^{9}&\cdots &\omega ^{3(N-1)}\\\vdots &\vdots &\vdots &\vdots & \ddots &\vdots \\[4pt]1&\omega ^{N-1}&\omega ^{2(N-1)}&\omega ^{3(N-1)}&\cdots &\omega ^{(N-1)(N-1)}\end{bmatrix}} $$ where $\omega = \omega_N=e^{\frac{2\pi i}{N}}$

Circuit implementation

In order to realize any unitary operation on most quantum computers, we must decompose the matrix operation into a circuit composed of single-qubit gates and two-qubit gates. Now, how might we construct such a circuit that performs the quantum Fourier transform on $n$-qubits?[2]

Express the $2$-qubit state $F_4(\ket{j})$ as a tensor product of single-qubit states. That is, find qubit states $\ket{\phi_1}$ and $\ket{\phi_2}$ such that $F_4(\ket{j}) = \ket{\phi_1} \otimes \ket{\phi_2}.$

Using definition 6.1, we have $$F_4(\ket{j}) = \frac{1}{2} \sum_{k=0}^{3} \omega_{4}^{jk} \ket{k}$$ where the state $\ket{k}=\ket{k_1 k_2}$ and the number $k = 2 k_1 + k_2,$ so $$ \begin{aligned} F_4(\ket{j}) &= \frac{1}{2} \sum_{k\in \{0, 1\}^2} \left(\omega_{N}^{j\left(2 k_1 + k_2\right)} \right)\ket{k_1 k_2}\\ &= \frac{1}{2} \left(e^{\frac{i 2\pi}{4} j(0)} \ket{00} + e^{\frac{i 2\pi}{4} j(1)} \ket{01} + e^{\frac{i 2\pi}{4} j(2)} \ket{10} + e^{\frac{i 2\pi}{4} j(3)} \ket{11}\right)\\ &= \frac{1}{2} \left(e^{\frac{i 2\pi}{4} j(0)} \ket{00} + e^{\frac{i 2\pi}{4} j(2)} \ket{10}\right) + \frac{1}{2} \left( e^{\frac{i 2\pi}{4} j(1)} \ket{01} + e^{\frac{i 2\pi}{4} j(3)} \ket{11}\right)\\ &= \frac{1}{2} \left(\left(e^{\frac{i 2\pi}{4}}\right)^{j(0)} \ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(2)} \ket{1}\right) \otimes \ket{0} + \frac{1}{2} \left(\left(e^{\frac{i 2\pi}{4}}\right)^{j(1)} \ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(3)} \ket{1}\right) \otimes \ket{1}\\ &= \frac{1}{2} \left[\left(\left(e^{\frac{i 2\pi}{4}}\right)^{j(0)} \ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(2)} \ket{1}\right) \otimes \ket{0} + \left(\left(e^{\frac{i 2\pi}{4}}\right)^{j(0)} \ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(2)} \ket{1}\right) \otimes e^{j\left(\frac{i 2\pi}{4}\right)} \ket{1}\right]\\ &= \frac{1}{2} \left[\left(\left(e^{\frac{i 2\pi}{4}}\right)^{j(0)} \ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(2)} \ket{1}\right) \otimes \left(\ket{0} + e^{j\left(\frac{i 2\pi}{4}\right)} \ket{1}\right)\right]\\ &= \frac{1}{\sqrt{2}} \left(\ket{0} + \left(e^{\frac{i 2\pi}{4}}\right)^{j(2)} \ket{1}\right) \otimes \frac{1}{\sqrt{2}}\left(\ket{0} + e^{j\left(\frac{i 2\pi}{4}\right)} \ket{1}\right)\\ &= \underbrace{\frac{1}{\sqrt{2}} \left(\ket{0} + e^{j \frac{i 4\pi}{4}}= \ket{1}\right)}_{\ket{\phi_1}} \otimes \underbrace{\frac{1}{\sqrt{2}}\left(\ket{0} + e^{j\left(\frac{i 2\pi}{4}\right)} \ket{1}\right)}_{\ket{\phi_2}} \end{aligned} $$

Now, let us try to express the $n$-qubit state $F_N(\ket{j})$ as a tensor product of single-qubit states. Using definition 6.1, we have $$F_N(\ket{j}) = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} \omega _{N}^{jk} \ket{k}$$ where the state $\ket{k}=\ket{k_1 k_2 \ldots k_n}$ and the number $k=2^{n-1} k_1 + 2^{n-2} k_2 + \cdots + 2^0 k_n,$ so $$ \begin{aligned} F_N(\ket{j}) &= \frac{1}{\sqrt{N}} \sum_{k\in \{0, 1\}^n} \left(\omega_{N}^{jk} \right)\ket{k_1 k_2 \cdots k_n}\\ &= \frac{1}{\sqrt{N}} \sum_{k_1=0}^{1} \sum_{k_2=0}^{1} \cdots \sum_{k_n=0}^{1} \left(\omega_{N}^{j\left(2^{n-1} k_1 + 2^{n-2} k_2 + \cdots + 2^0 k_n\right)} \right)\ket{k_1 k_2 \cdots k_n}\\ &= \frac{1}{\sqrt{N}} \sum_{k_1=0}^{1} \sum_{k_2=0}^{1} \cdots \sum_{k_{n-1}=0}^{1} \left(\omega_{N}^{j\left(2^{n-1} k_1 + 2^{n-2} k_2 + \cdots + 2^0 (0)\right)} \right)\ket{k_1 k_2 \cdots k_{n-1}} \otimes \ket{0}\\ & \hspace{105.75pt} +\left(\omega_{N}^{j\left(2^{n-1} k_1 + 2^{n-2} k_2 + \cdots + 2^0 (1)\right)} \right)\ket{k_1 k_2 \cdots k_{n-1}} \otimes \ket{1}\\ &= \frac{1}{\sqrt{N/2}} \sum_{k_1=0}^{1} \sum_{k_2=0}^{1} \cdots \sum_{k_{n-1}=0}^{1} \left(\omega_{N/2}^{j\left(2^{n-1} k_1 + 2^{n-2} k_2 + \cdots + 2^1 k_{n-1}\right)/2} \right)\ket{k_1 k_2 \cdots k_{n-1}} \otimes \left(\frac{\ket{0}+\omega_{N}^{j}\ket{1}}{\sqrt{2}}\right)\\ &= \frac{1}{\sqrt{N/2}} \sum_{k_1=0}^{1} \sum_{k_2=0}^{1} \cdots \sum_{k_{n-1}=0}^{1} \left(\omega_{N/2}^{j\left(2^{n-2} k_1 + 2^{n-3} k_2 + \cdots + 2^0 k_{n-1}\right)} \right)\ket{k_1 k_2 \cdots k_{n-1}} \otimes \left(\frac{\ket{0}+e^{\frac{i 2 \pi j}{N}}\ket{1}}{\sqrt{2}}\right) \end{aligned} $$ Let the number $k'=2^{n-2} k_1 + 2^{n-3} k_2 + \cdots + 2^0 k_{n-1}$ and the state $\ket{k'}=\ket{k_1 k_2 \ldots k_{n-1}}.$ Since $$j=2^{n-1} j_1 + 2^{n-2} j_2 + \cdots + 2^0 j_n,$$ we have $$ \begin{aligned} \omega_{N/2}^{jk'} = \omega_{N/2}^{\left(2^{n-1} j_1 + 2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'} &= \left(\omega_{N/2}^{\left(2^{n-1} j_1\right)k'} \right) \left(\omega_{N/2}^{\left(2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'}\right)\\ &= \left(e^{\frac{i2\pi}{N/2}2^{n-1}j_1 k'}\right) \left( \omega_{N/2}^{\left(2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'}\right)\\ &= \underbrace{\left(e^{i2\pi j_1 k'}\right)}_{1} \left( \omega_{N/2}^{\left(2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'}\right)\\ &= \omega_{N/2}^{\left(2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'}. \end{aligned} $$ Additionally, notice that $\displaystyle \frac{j}{N}=\frac{2^{n-1} j_1 + 2^{n-2} j_2 + \cdots + 2^0 j_n}{2^n}=[0.j_{1}j_{2}\ldots j_{n}].$ Therefore, we have $$ \begin{aligned} F_N(\ket{j}) &= \underbrace{\frac{1}{\sqrt{N/2}} \sum_{k'\in \{0, 1\}^{n-1}} \left(\omega_{N/2}^{\left(2^{n-2} j_2 + \cdots + 2^0 j_n\right)k'} \right)\ket{k'}}_{\displaystyle F_{N/2}(\ket{j_2 j_3 \cdots j_n})} \otimes \left(\frac{\ket{0}+e^{\frac{i 2 \pi [0.j_{1}j_{2}\ldots j_{n}]}{N}}\ket{1}}{\sqrt{2}}\right). \end{aligned} $$ Repeating once, we get $$F_{N/2}(\ket{j_2 j_3 \cdots j_n})=F_{N/4}(\ket{j_3 j_4 \cdots j_n}) \otimes \left(\frac{\ket{0}+e^{i 2 \pi [0.j_{2}j_3\ldots j_{n}]}\ket{1}}{\sqrt{2}}\right).$$ Completing this process $n$ times, we get $$ \begin{aligned} F_N(\ket{j})&=\left(\frac{\ket{0}+e^{i 2 \pi [0.j_{n}]}\ket{1}}{\sqrt{2}}\right) \otimes \cdots \otimes \left(\frac{\ket{0}+e^{i 2 \pi [0.j_{2}j_3\ldots j_{n}]}\ket{1}}{\sqrt{2}}\right) \otimes \left(\frac{\ket{0}+e^{i 2 \pi [0.j_{1}j_{2}\ldots j_{n}]}\ket{1}}{\sqrt{2}}\right)\\ &=\frac{1}{\sqrt{N}}\left(\ket{0}+e^{i 2 \pi [0.j_{n}]}\ket{1}\right) \otimes \cdots \otimes \left(\ket{0}+e^{i 2 \pi [0.j_{2}j_3\ldots j_{n}]}\ket{1}\right) \otimes \left(\ket{0}+e^{i 2 \pi [0.j_{1}j_{2}\ldots j_{n}]}\ket{1}\right) \end{aligned} $$ The resulting quantum state after an $n$-qubit quantum Fourier transform bears a striking resemblance to the Hadamard transform, save for the complex phase factors.

The quantum Fourier transform circuit uses two key ingredients:

  • the Hadamard gate $H,$ such that $\displaystyle H\ket{x}=\frac{1}{\sqrt{2}}\left(\ket{0}+e^{\frac{i2\pi}{2}x}\ket{1}\right)$
  • the phase gate $P_k=\begin{pmatrix} 1 & 0 \\ 0 & e^{i 2 \pi /2^k} \end{pmatrix},$ so that $P_k\ket{0}=\ket{0}$ and $P_k\ket{1}=e^{i 2 \pi /2^k}\ket{1}$

Assume you are able to construct a quantum circuit realizing an $(n-1)$-qubit QFT. Extend the circuit to realize an $n$-qubit QFT by adding sub-circuit $E_n$.

Circuit for Problem 1.6.5: extending an (n-1)-qubit QFT with sub-circuit E_n

A good starting point for the circuit $E$ would be the Hadamard gate. The Hadamard gate alone produces the state $\frac{1}{\sqrt{2}} \left(\ket{0}+e^{i 2 \pi [0.x_{1}]}\ket{1}\right)$ in the first qubit, but the correct state would be $\frac{1}{\sqrt{2}} \left(\ket{0}+e^{i 2 \pi [0.x_{1}\ldots x_{n}]}\ket{1}\right).$ The missing relative phase is $$e^{i 2 \pi [0.0 x_{2}\ldots x_{n}]}=e^{i 2 \pi \left({2^{n-2} x_2 + 2^{n-3} x_3 + \cdots + 2^0 x_n}\right)/{2^n}}=e^{i 2 \pi \left({x_2}/{2^2}\right)} \times e^{i 2 \pi \left({x_3}/{2^3}\right)}\times \cdots \times e^{i 2 \pi \left({x_n}/{2^n}\right)}.$$ Remember we have a phase gate $P_k$ such that $P_k\ket{0}=\ket{0}$ and $P_k\ket{1}=e^{i 2 \pi /2^k}\ket{1}.$ Therefore, we can add each component phase $$e^{i 2 \pi \left({x_i}/{2^i}\right)}=\begin{cases} 1, & \text{if $x_i=0$}\\ e^{i 2 \pi /{2^i}}, & \text{if $x_i=1$} \end{cases}$$ separately using a sequence of $P_{i}$ gates with each gate $P_i$ controlled by $x_i.$

Solution for Problem 1.6.5: E_n sub-circuit with Hadamard and controlled phase gates

Now, we have enough information to construct the quantum circuit implementing the general $n$-qubit QFT. We know that the single-qubit QFT circuit is the Hadamard gate. Adding extension circuit $E_2$ to the single-qubit QFT circuit gives the two-qubit QFT circuit. We can systematically build the $n$-qubit QFT by iteratively adding extensions circuits $E_2, E_3, \ldots, E_{n-2}, E_{n-1}, E_{n}$.

Figure 1.6.1: Quantum circuit for the n-qubit quantum Fourier transform

Unfortunately, the output states of the qubits are in the reverse order of Eq. \ref{eq:qft}. We can easily invert the order of the output states by applying a sequence of SWAP gates.

Quantum circuit for the n-qubit QFT with SWAP gates

Multi-valued Bernstein-Vazirani algorithm

While alternative positional number systems, such as ternary (base 3), have been considered in classical computing hardware, these systems have not seen widespread adoption. Several practical obstacles prevent their adoption such as the increased complexity of building non-binary hardware and the entrenched dominance of binary-based technology.

Quantum computing, however, offers a unique opportunity to revisit the concept of multi-valued logic. Quantum computing is still in its infancy and unlike most classical hardware technologies, several quantum hardware technologies naturally accommodate more than two states. Companies and research institutions are increasingly exploring multi-valued quantum systems, which have the potential to transcend the binary paradigm.

At the core of multi-valued quantum computing is the concept of the qudit (quantum digit), a generalization of the qubit. While a qubit is restricted to two states, a qudit can exist in $d$ states, where $d>2$. For example, a ternary qudit, or qutrit, has three basis states: $\ket{0}, \ket{1},$ and $\ket{2}.$ This expanded state space allows quantum algorithms to encode, process, and retrieve significantly more information per quantum unit. The benefit of using qudits grows exponentially; while an $n$-qubit system has a state space of $2^n,$ an $n$-qudit system has a state space of $d^n$. Consequently, multi-valued quantum computing offers a far greater advantage compared to multi-valued classical computing, which is only a linear improvement over its binary counterpart.

The Hadamard gate is a single-qubit quantum Fourier transform, i.e., a quantum Fourier transform over the cyclic group $\mathbb{Z}_2.$ The ternary equivalent of the Hadamard gate is called the Chrestenson gate, which is a single-qutrit quantum Fourier transform over the cyclic group $\mathbb{Z}_3$. Determine the matrix for the Chrestenson gate.

Since $N=3,$ we have $\displaystyle \omega_N=\omega_3=e^{\frac{2\pi i}{3}}.$ Substituting into Eq. \ref{eq:qftmatrix} gives $$ F_{3}={\frac {1}{\sqrt {3}}}{\begin{bmatrix}1&1&1\\ 1&\omega_3 &\omega_3^{2}\\ 1&\omega_3^{2}&\omega_3^{4}\end{bmatrix}}={\frac {1}{\sqrt {3}}}{\begin{bmatrix}1&1&1\\ 1&\omega_3 &\omega_3^{2}\\ 1&\omega_3^{2}&\omega_3\end{bmatrix}}={\frac {1}{\sqrt {3}}}{\begin{bmatrix}1&1&1\\ 1&e^{\frac{2\pi i}{3}} & e^{\frac{4\pi i}{3}}\\ 1&e^{\frac{4\pi i}{3}}& e^{\frac{2\pi i}{3}}\end{bmatrix}}. $$

Let $C$ be the matrix for the single-qutrit Chrestenson gate defined above.

A qutrit is initialized to an input state $\ket{x},$ which is either $\ket{0},$ $\ket{1},$ or $\ket{2}.$ A Chrestenson gate is applied, mapping the state $\ket{x}$ to $\ket{\psi}.$ Define an expression for $\ket{\psi}$ in terms of $\ket{x}.$

$$ \begin{aligned} \ket{\psi} = C\ket{x} &= \begin{cases} \frac{1}{\sqrt{3}}\left(\ket{0}+\ket{1}+\ket{2}\right) & \text{if } \ket{x}=\ket{0}\\ \frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3\ket{1}+\omega_3^{2}\ket{2}\right) & \text{if } \ket{x}=\ket{1}\\ \frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3^{2}\ket{1}+\omega_3\ket{2}\right) & \text{if } \ket{x}=\ket{2} \end{cases} \\ &= \frac{\ket{0}+(\omega_3)^{x}\ket{1}+(\omega_3^2)^{x}\ket{1}}{\sqrt{3}}\\ &= \frac{1}{\sqrt{3}}\sum_{y \in \{0,1,2\}} (\omega_3)^{xy}\ket{y} \end{aligned} $$

The general $n$-qutrit Chrestenson transform maps an arbitrary computational basis state $\ket{x} = \ket{x_1 \ldots x_n}$ to $\ket{\psi}.$ Find $\ket{\psi}.$

$$ \begin{aligned} \ket{\psi } = C^{\otimes n}\ket{x_1 \ldots x_n} = \bigotimes_{i=1}^n C\ket{x_i} &= \bigotimes_{i=1}^n \left(\frac{\ket{0}+(\omega_3)^{x_i}\ket{1}+(\omega_3^2)^{x_i}\ket{1}}{\sqrt{3}}\right)\\ &= \prod_{i=1}^n \left(\frac{1}{\sqrt{3}}\sum_{y_i \in \{0,1,2\}} (\omega_3)^{x_i y_i}\ket{y_i}\right)\\ &= \frac{1}{\sqrt{3^n}}\sum_{y_i \in \{0,1,2\}} \left(\prod_{i=1}^n (\omega_3)^{x_i y_i}\right)\ket{y_i}\\ &= \frac{1}{\sqrt{3^n}}\sum_{y \in \{0,1,2\}^n} (\omega_3)^{\sum_{i=1}^n x_i y_i}\ket{y} \\ &= \frac{1}{\sqrt{3^n}}\sum_{y \in \{0,1,2\}^n} (\omega_3)^{x \boldsymbol{\cdot} y}\ket{y} \end{aligned} $$

Given a black box function $f: \{0, 1, 2\}^n \rightarrow \{0, 1, 2\}$ such that $$ f(\mathbf{x}) = \mathbf{x} \boldsymbol{\cdot} \mathbf{s} \text{ mod } 3 = \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_n \end{pmatrix} \boldsymbol{\cdot} \begin{pmatrix} s_1\\ s_2\\ \vdots\\ s_n \end{pmatrix} \text{mod } 3 = (s_1 x_1) \oplus_3 (s_2 x_2) \oplus_3 \ldots \oplus_3 (s_n x_n) $$ the ternary Bernstein-Vazirani problem is to find the hidden string $\mathbf{s} \in \{0, 1, 2\}^n.$

The quantum circuit associated with the ternary Bernstein-Vazirani algorithm is identical to the binary Bernstein-Vazirani algorithm, but replacing all Hadamard gates with Chrestenson gates.

Evaluate the progression of the quantum state from $\ket{\Psi_0}$ to $\ket{\Psi_3}.$ Then, analyze the result of the measurement of $\ket{\Psi_3}.$

Circuit for Problem 1.7.4: ternary Bernstein-Vazirani algorithm
  1. The circuit is initialized to the state $\ket{\Psi_0} = \ket{0}^{\oplus n}\ket{1}.$
  2. The first part of the ternary Bernstein-Vazirani algorithm is a Hadamard transform on all qubits, resulting in the state $$\ket{\Psi_1} = \left(\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} \ket{x}\right) \otimes \left(\frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3\ket{1}+\omega_3^{2}\ket{2}\right)\right)$$
  3. In the next part of the algorithm, an oracle $U_f$ is applied to the state $\ket{\Psi_1}$. Therefore, the resulting quantum state $\ket{\Psi_2}$ is a uniform superposition of all basis states, but phase is induced on all basis states for which $f\neq0$: $$ \begin{aligned} \ket{\Psi_2} = U_f\left(\ket{\Psi_1}\right) &=\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} U_f\left(\ket{x}\otimes \frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3\ket{1}+\omega_3^{2}\ket{2}\right)\right)\\ &=\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} \ket{x}\otimes \frac{1}{\sqrt{3}}\left(\ket{0 \oplus_3 f(x)}+\omega_3\ket{1\oplus_3 f(x)}+\omega_3^{2}\ket{2\oplus_3 f(x)}\right)\\ &=\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} \ket{x}\otimes \begin{dcases} \frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3\ket{1}+\omega_3^{2}\ket{2}\right) & \text{if } \ket{f(x)}=\ket{0}\\ \frac{1}{\sqrt{3}}\left(\ket{1}+\omega_3\ket{2}+\omega_3^{2}\ket{0}\right) & \text{if } \ket{f(x)}=\ket{1}\\ \frac{1}{\sqrt{3}}\left(\ket{2}+\omega_3\ket{0}+\omega_3^{2}\ket{1}\right) & \text{if } \ket{f(x)}=\ket{2} \end{dcases}\\ &=\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} (\omega_3^2)^{f(x)} \ket{x}\otimes \left(\frac{1}{\sqrt{3}}\left(\ket{0}+\omega_3\ket{1}+\omega_3^{2}\ket{2}\right)\right)\\ \end{aligned} $$ Measurement is dependent only on amplitude squared, not phase. Therefore, measuring the state $\ket{\Psi_2}$ is equivalent to $\ket{\Psi_1},$ which would yield a random $\ket{x}.$
  4. Finally, discarding the last qutrit and applying yet another Hadamard transform onto the $n$ input qutrits yields $$ \begin{aligned} \ket{\Psi_3} &= C^{\otimes n} \left(\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} (\omega_3^2)^{f(x)} \ket{x}\right)\\ &= \frac{1}{\sqrt{3^n}}\sum_{y \in \{0,1,2\}^n}(\omega_3)^{x \cdot y} \left(\frac{1}{\sqrt{3^n}}\sum_{x \in \{0,1,2\}^n} (\omega_3^{-1})^{f(x)} \right)\ket{y}\\ &= \sum_{y \in \{0,1,2\}^n}\underbrace{\frac{1}{3^n}\sum_{x \in \{0,1,2\}^n} (\omega_3)^{x \cdot y-f(x)}}_{\text{amplitude of }\ket{y}}\ket{y}. \end{aligned} $$

Since $f(x_1 \ldots x_n) = s \cdot x \text{ mod } 3,$ the amplitude of $\ket{y}$ is $$ \frac{1}{3^n} \sum_{x \in \{0,1\}^n} (\omega_3^2)^{x \cdot y-s\cdot x}=\frac{1}{3^n}\sum_{x \in \{0,1,2\}^n} (\omega_3^2)^{x(y-s)}, $$ which equals $1$ if $y$ is the correct linear form because $$y = s \implies y - s = 0 \implies \frac{1}{3^n}\sum_{x \in \{0,1\}^n} \left(1\right) = \frac{3^n}{3^n} = 1,$$ meaning that measuring the output register yields $s.$

We can generalize these results from ternary to arbitrary radix $r.$

The general $n$-qudit quantum Fourier transform over the group $\mathbb{Z}_r^n$ maps an arbitrary $r$-valued computational basis state $\ket{x} = \ket{x_1 \ldots x_n}$ to $\ket{\psi}.$ Find $\ket{\psi}.$

$$ \begin{aligned} \ket{\psi } = F_r^{\otimes n}\ket{x_1 \ldots x_n} = \bigotimes_{i=1}^n F_r\ket{x_i} &= \bigotimes_{i=1}^n \left( \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \omega_{r}^{x_i k} \ket{k}\right)\\ &= \prod_{i=1}^n \left(\frac{1}{\sqrt{r}} \sum_{k_i=0}^{r-1} \omega_{r}^{x_i k_i} \ket{k_i}\right)\\ &= \frac{1}{\sqrt{r^n}}\sum_{k_i=0}^{r-1} \left(\prod_{i=1}^n (\omega_r)^{x_i k_i}\right)\ket{k_i}\\ &= \frac{1}{\sqrt{r^n}}\sum_{k \in \{0,\ldots,r-1\}^n} (\omega_r)^{\sum_{i=1}^n x_i k_i}\ket{k} \\ &= \frac{1}{\sqrt{r^n}}\sum_{k \in \{0,\ldots,r-1\}^n} (\omega_r)^{x \boldsymbol{\cdot} k}\ket{k} \end{aligned} $$

The quantum circuit associated with the $r$-valued Bernstein-Vazirani algorithm is shown in Figure 1.7.1.

  1. The circuit is initialized to the state $\ket{\Psi_0} = \ket{0}^{\oplus n}\ket{1}.$
  2. The first part of the $r$-ary Bernstein-Vazirani algorithm is a single-qudit quantum Fourier transform applied on each qubit, resulting in the collective state $$\ket{\Psi_1} = \left(\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} \ket{x}\right) \otimes \left( \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \omega_{r}^{k} \ket{k}\right)$$
  3. In the next part of the algorithm, an oracle $U_f$ is applied to the state $\ket{\Psi_1}$. Therefore, the resulting quantum state $\ket{\Psi_2}$ is a uniform superposition of all basis states, but phase is induced on all basis states for which $f\neq0$: $$ \begin{aligned} \ket{\Psi_2} &= U_f\left(\ket{\Psi_1}\right)\\ &=\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} U_f\left(\ket{x}\otimes \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \omega_{r}^{k} \ket{k}\right)\\ &=\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} \ket{x}\otimes \frac{1}{\sqrt{r}} \left(\ket{0 \oplus_r f(x)} + \omega_{r} \ket{1 \oplus_r f(x)} + \omega_{r}^{2} \ket{2 \oplus_r f(x)} + \ldots + \omega_{r}^{r-1} \ket{(r-1) \oplus_r f(x)}\right)\\ &=\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} \ket{x}\otimes \begin{dcases} \frac{1}{\sqrt{r}} \left(\ket{0} + \omega_{r} \ket{1} + \omega_{r}^{2} \ket{2} + \ldots + \omega_{r}^{r-1} \ket{r-1}\right) & \text{if } \ket{f(x)}=\ket{0}\\ \frac{1}{\sqrt{r}} \left(\ket{1} + \omega_{r} \ket{2} + \omega_{r}^{2} \ket{3} + \ldots + \omega_{r}^{r-1} \ket{0}\right) & \text{if } \ket{f(x)}=\ket{1}\\ \frac{1}{\sqrt{r}} \left(\ket{2} + \omega_{r} \ket{3} + \omega_{r}^{2} \ket{4} + \ldots + \omega_{r}^{r-1} \ket{1}\right) & \text{if } \ket{f(x)}=\ket{2} \\ \hspace{90pt}\ldots & \text{if } \hspace{15pt}\ldots \end{dcases}\\ &=\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r^{r-1})^{f(x)} \ket{x}\otimes \frac{1}{\sqrt{r}} \sum_{k=0}^{r-1} \omega_{r}^{k} \ket{k}\\ \end{aligned} $$ Measurement is dependent only on amplitude squared, not phase. Therefore, measuring the state $\ket{\Psi_2}$ is equivalent to $\ket{\Psi_1},$ which would yield a random $\ket{x}.$
  4. Finally, discarding the last qudit and again applying a single-qubit quantum Fourier transform onto each of the $n$ input qudits yields $$ \begin{aligned} \ket{\Psi_3} &= F_r^{\otimes n} \left(\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r^{r-1})^{f(x)} \ket{x}\right)\\ &= \frac{1}{\sqrt{r^n}}\sum_{k \in \{0,\ldots,r-1\}^n}(\omega_r)^{x \cdot k} \left(\frac{1}{\sqrt{r^n}}\sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r^{-1})^{f(x)} \right)\ket{k}\\ &= \sum_{k \in \{0,\ldots,r-1\}^n}\underbrace{\frac{1}{r^n}\sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r)^{x \cdot k+f(x)}}_{\text{amplitude of }\ket{k}}\ket{k}. \end{aligned} $$

Since $f(x_1 \ldots x_n) = s \cdot x \text{ mod } 3,$ the amplitude of $\ket{k}$ is $$ \frac{1}{r^n} \sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r)^{x \cdot k-s\cdot x}\ket{k}=\frac{1}{r^n}\sum_{x \in \{0,\ldots,r-1\}^n} (\omega_r)^{x(k - s)}\ket{k}, $$ which equals $1$ if $k$ is the correct linear form because $$k = s \implies k - s = 0 \implies \frac{1}{r^n}\sum_{x \in \{0,\ldots,r-1\}^n} \left(1\right) = \frac{r^n}{r^n} = 1,$$ meaning that measuring the output register yields $s.$

Finding the eigenvalues of matrices is important across a variety of domains. Consider a unitary matrix $U$ and a complex-valued eigenvalue $\lambda$ such that $U \mathbf{v} = \lambda \mathbf{v}$ for some non-zero eigenvector $\mathbf{v}$. What are the possible values for $\lambda$?

Because $U$ preserves the magnitude of vectors, we have $$ \|U \mathbf{v}\|^2 = (U \mathbf{v})^\dagger U \mathbf{v} = \mathbf{v}^\dagger U^\dagger U \mathbf{v} = \mathbf{v}^\dagger \mathbf{v} = \|\mathbf{v}\|^2.$$ and therefore, $$\|\mathbf{v}\|^2 = \|U \mathbf{v}\|^2 = \|\lambda \mathbf{v}\|^2 = |\lambda|^2 \|\mathbf{v}\|^2. $$ Since $ \|\mathbf{v}\| \neq 0 $ (because $\mathbf{v}$ is an eigenvector), it follows that $ |\lambda|^2 = 1 \implies |\lambda| = 1.$ Thus, the eigenvalues of a unitary matrix are complex numbers with magnitude 1.

Phase estimation

Given a unitary operator $U$ with eigenstate $\ket{\psi}$ and eigenvalue $e^{2 \pi i \theta}$ such that $$U\ket{\psi} = e^{2 \pi i \theta}\ket{\psi},$$ the quantum phase estimation (QPE) problem is to estimate the phase $\theta.$[1]

How might we go about designing a circuit to find $\theta$?

Let's start by evaluating the progression of the quantum state in the following quantum circuit:

The initial quantum state $\ket{\Psi_0}$ is $\ket{0}\ket{\psi}.$ After the Hadamard gate is applied on the first qubit, the state becomes $$\ket{\Psi_1} = \ket{+}\ket{\psi} = \frac{1}{\sqrt{2}}\Bigl(\ket{0}+\ket{1}\Bigr)\ket{\psi} = \frac{1}{\sqrt{2}}\Bigl(\ket{0}\ket{\psi}+\ket{1}\ket{\psi}\Bigr).$$ Next, the controlled-$U$ gate is applied, inducing a relative phase onto the quantum state $$\ket{\Psi_2} = \frac{1}{\sqrt{2}}\Bigl(\ket{0}\ket{\psi}+\ket{1}U\ket{\psi}\Bigr) = \frac{1}{\sqrt{2}}\Bigl(\ket{0}\ket{\psi}+\ket{1}e^{2 \pi i \theta}\ket{\psi}\Bigr),$$ yet $\ket{\Psi_2}$ is still indistinguishable from $\ket{\Psi_1}$ with measurement. Only when the second Hadamard gate is applied on the first qubit, creating an interference pattern, does the phase factor $\theta$ influence measurement probabilities: $$ \begin{aligned} \ket{\Psi_3} &= \frac{1}{\sqrt{2}}\Bigl(H\ket{0}\ket{\psi}+H\ket{1}U\ket{\psi}\Bigr)\\ &= \frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{2}}\Bigl(\ket{0}+\ket{1}\Bigr)\ket{\psi}+\frac{1}{\sqrt{2}}\Bigl(\ket{0}-\ket{1}\Bigr)e^{2 \pi i \theta}\ket{\psi}\right)\\ &= \frac{1}{2}\left(\ket{0}\Bigl(\ket{\psi}+e^{2 \pi i \theta}\ket{\psi}\Bigr)+\ket{1}\Bigl(\ket{\psi}-e^{2 \pi i \theta}\ket{\psi}\Bigr)\right)\\ &= \left(\frac{1+e^{2 \pi i \theta}}{2}\right)\ket{0}\ket{\psi}+\left(\frac{1-e^{2 \pi i \theta}}{2}\right)\ket{1}\ket{\psi}. \end{aligned} $$

In the problem above, the probability of measuring $0$ is $$ \begin{aligned} \overbrace{\left|\frac{1+e^{2 \pi i \theta}}{2}\right|^2 = \frac{1}{4}\left(1+e^{2 \pi i \theta}\right)\left(1+e^{-2 \pi i \theta}\right)}^{\left|z\right|^2 z\overline{z} \quad \forall z \in \mathbb{Z}} &= \frac{1+e^{2 \pi i \theta}+e^{-2 \pi i \theta}+1}{4}\\ &= \underbrace{\frac{1+\cos(2\pi \theta)}{2}=\textcolor{red}{\cos^2(\pi \theta)}}_{\text{cosine half-angle identity}} \end{aligned} $$ and the probability of measuring $1$ is $$ \begin{aligned} \left|\frac{1-e^{2 \pi i \theta}}{2}\right|^2 = \frac{1}{4}\left(1-e^{2 \pi i \theta}\right)\left(1-e^{-2 \pi i \theta}\right) &= \frac{1-e^{2 \pi i \theta}-e^{-2 \pi i \theta}+1}{4}\\ &= \underbrace{\frac{1-\cos(2\pi \theta)}{2} = \textcolor{blue}{\sin^2(\pi \theta)}}_{\text{sine half-angle identity}} \end{aligned}$$ If we plot a graph of the probability of measurement of $0$ and $1$ against $\theta,$ we get

Interestingly, we always measure $0$ when $\theta=0$ and $1$ when $\theta=1/2.$ If $\theta$ is somehow restricted to either $0$ or $0.5,$ the circuit shown in the figure above solves the quantum phase estimation problem! In other words, this simple circuit approximates $\theta$ to $1$-bit of accuracy: $$\theta = 0.x =\begin{cases} 0, & x=0\\ 1/2, & x=1 \end{cases}$$ where $x$ is the measurement outcome. Now, let's try to increase the accuracy.

Consider modifying the circuit such that there are two controlled-$U$ operations in sequence, effectively doubling the relative phase that is applied. How would the measurement probabilities be affected? Would the accuracy of determining $\theta$ increase?

Phase estimation

Two controlled-$U$ operations, each applying a phase factor $\theta$, is equivalent to a single controlled-$U^2$ operation applying a phase factor $2\theta.$ Therefore, the probability of measuring $0$ would be $\cos^2(2\pi \theta)$ and the probability of measuring $1$ is $\sin^2(2\pi \theta).$ If we plot a graph of the probability of measurement of $0$ and $1$ against $\theta,$ we get

We still have $1$-bit accuracy. Instead of measuring $x_1,$ we now measure $x_2.$ We leave as an exercise to the reader to prove the general case: that using $2^n$ controlled-$U$ gates in sequence allows us to measure the $n$th bit or $x_n$.

By measuring the outcomes from our first circuit and its modified version with double the phase, we can achieve $2$-bit accuracy. It's not immediately obvious, however, how we can achieve this result in a single quantum algorithm.

Let's see what happens when we combine them without adding the final Hadamard transform

The initial quantum state $\ket{\Psi_0}$ is $\ket{00}\ket{\psi}.$ After the Hadamard gate is applied on the first qubit, the state becomes $$\ket{\Psi_1} = \ket{++}\ket{\psi} = \frac{1}{2}\Bigl(\ket{00}+\ket{01}+\ket{10}+\ket{11}\Bigr)\ket{\psi} = \frac{1}{2}\Bigl(\ket{00}\ket{\psi}+\ket{01}\ket{\psi}+\ket{10}\ket{\psi}+\ket{11}\ket{\psi}\Bigr).$$ Next, the controlled-$U$ gates are applied, inducing relative phase onto the quantum state $$ \begin{aligned} \ket{\Psi_2} &= \frac{1}{2}\Bigl(\ket{00}\ket{\psi}+\ket{01}U\ket{\psi}+\ket{10}U\ket{\psi}+\ket{11}U^2\ket{\psi}\Bigr)\\ &= \frac{1}{2}\Bigl(\ket{00}\ket{\psi}+\ket{01}e^{2 \pi i \theta}\ket{\psi}+\ket{10}e^{2 \pi i (2\theta)}\ket{\psi}+\ket{11}e^{2 \pi i (3\theta)}\ket{\psi}\Bigr). \end{aligned} $$

  1. Assume the phase factor $\theta$ of unitary operator in Problem 1.7.4 is restricted to $\{0, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}\}.$ Evaluate the state of the first two qubits unique to every possible $\theta$
  2. The set of these states forms an orthonormal basis. What is this basis? What transformation maps this basis back to the computational basis?
  1. $$ \begin{aligned} \ket{0.00}&=\frac{1}{2}\ket{00}+\frac{1}{2}\ket{01}+\frac{1}{2}\ket{10}+\frac{1}{2}\ket{11}\\ \ket{0.01}&=\frac{1}{2}\ket{00}+\frac{i}{2}\ket{01}-\frac{1}{2}\ket{10}-\frac{i}{2}\ket{11}\\ \ket{0.10}&=\frac{1}{2}\ket{00}-\frac{1}{2}\ket{01}+\frac{1}{2}\ket{10}-\frac{1}{2}\ket{11}\\ \ket{0.11}&=\frac{1}{2}\ket{00}-\frac{i}{2}\ket{01}-\frac{1}{2}\ket{10}+\frac{i}{2}\ket{11}\\ \end{aligned} $$
  2. The basis $\{\ket{0.00}, \ket{0.01}, \ket{0.10}, \ket{0.11}\}$ is the Fourier basis! The quantum Fourier transform $\mathcal{QFT}_{2}$ maps the $2$-qubit computational basis to the Fourier basis. In order to realize an inverse mapping, we use the inverse operation: the inverse quantum Fourier transform $\mathcal{QFT}_{2}^{-1}.$
Single-qubit phase estimation

The quantum Fourier transform (QFT) is efficient for identifying the periodicity of quantum amplitudes, which is essential for phase estimation. If you are unfamiliar with QFT, I recommend reading my guide on the quantum Fourier transform first.

The circuit for the general quantum phase estimation algorithm is depicted in Fig. 1.8.3.

Circuit for Problem 1.8.3

For some input to the controls $x=x_1 x_2 \ldots x_n,$ the action of $U_{\text{total}}$ on the target register can be described as $$U^{x_n} \left(U^2\right)^{x_{n-1}} \ldots \left(U^{2^n}\right)^{x_1} = U^{\left(2^n x_1 + 2^{n-1} x_2 + \ldots + x_1 \right)}=U^x,$$ which allows us to reduce the quantum circuit shown in Fig. 1.8.3 to the schematic version in Fig. 1.8.4.

  1. The circuit is initialized to the state $$\ket{\Psi_0} = \ket{0}^{\oplus m}\ket{\psi}.$$
  2. First, a Hadamard transform is applied on the first register, resulting in an even superposition of the first register: $$\ket{\Psi_1} = \left(H^{\otimes m}\ket{0}^{\oplus m}\right)\ket{\psi}=\left(\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}\ket{x}\right)\ket{\psi}.$$
  3. Next, the sequence of controlled unitary operations $U_{\text{total}}$ is applied to the state $\ket{\Psi_1}$. The resulting quantum state $$ \begin{aligned} \ket{\Psi_2}=U_{\text{total}}\ket{\Psi_1}&=U_{\text{total}}\left(\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}\ket{x}\right)\ket{\psi}\\ &=\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}\ket{x}U^x\ket{\psi}\\ &=\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}\ket{x}e^{2\pi i (x \theta)}\ket{\psi} \end{aligned} $$ From the point of view of measurement of the first register, the state $\ket{\Psi_2}$ is equivalent to $\ket{\Psi_1},$ producing a random $\ket{x}.$
  4. Finally, we apply an inverse quantum Fourier transform on the $m$-qubit first register: $$ \begin{aligned} \ket{\Psi_3} &= \mathcal{QFT}_{2^m}^{-1} \left(\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}e^{2\pi i (x \theta)}\ket{x}\right)\\ &= \sum_{k = 0}^{2^m-1} \frac{1}{\sqrt{2^m}} \left(\frac{1}{\sqrt{2^m}}\sum_{x=0}^{2^m-1}e^{2\pi i (x \theta)}\right)e^{-\frac{2\pi i}{2^m}xk}\ket{k}\\ &= \sum_{k = 0}^{2^m-1} \left(\frac{1}{2^m} \sum_{x=0}^{2^m-1}e^{2\pi i x\left( \theta -\frac{k}{2^m}\right)}\right)\ket{k} \end{aligned} $$

In the phase estimation problem, it is only given that $0 \leq \theta < 1.$ Making a further assumption that $2^m\theta$ is an integer, determine the probability of measuring $\ket{2^m\theta}.$

If $k=2^m\theta \implies \theta -\frac{k}{2^m} = 0,$ then the probability of measuring $\ket{k}$ is $$ \left|\frac{1}{2^m} \sum_{x=0}^{2^m-1}e^{2\pi i x\left(0\right)}\right|^2 = 1. $$

Let $a$ be the nearest integer to $2^m \theta.$ Determine its probability of measurement. Do not assume that $2^m\theta$ is an integer.

If the margin of error $\displaystyle \Delta = \theta -\frac{a}{2^m},$ then $$\left|2^m \theta - a \right| \leq \frac{1}{2} \implies \left|2^m \left(\Delta +\frac{a}{2^m}\right) - a \right| \leq \frac{1}{2} \implies \left|\Delta \right| \leq \frac{1}{2^{m+1}}.$$ The probability of measuring $a$ is $$ \left|\frac{1}{2^m} \sum_{x=0}^{2^m-1}e^{2\pi i x \Delta}\right|^2 =\frac{1}{2^{2m}}\left|\left(1+e^{2\pi i \Delta}+e^{4\pi i \Delta}+\ldots\right)\right|^2=\frac{1}{2^{2m}}\left|\frac{1-e^{2\pi i (2^m \Delta)}}{1-e^{2\pi i \Delta}}\right|^2 $$ because of the well-known geometric series formula $\displaystyle 1+r+r^2+\ldots=\frac{1-r^n}{1-r}.$ In general, the expression $\displaystyle \left|1-e^{2ix}\right|^2 = (1-e^{2ix})(1-e^{-2ix}) = 1-e^{2ix}-e^{-2ix}+1.$ Converting to general form, we get $\displaystyle 1-(\cos(2x)+i\sin(2x))-(\cos(2x)- i\sin(2x))+1=2-2\cos(2x),$ which by the cosine double-angle identity equals $\displaystyle 2-2(1-2\sin^2(x))=4\sin^2(x).$ Therefore, $$\displaystyle \frac{1}{2^{2m}}\left|\frac{1-e^{2\pi i (2^m \Delta)}}{1-e^{2\pi i \Delta}}\right|^2 = \frac{1}{2^{2m}}\left|\frac{4\sin(\pi 2^m \Delta)}{4\sin(\pi \Delta)}\right|^2 = \frac{1}{2^{2m}}\left|\frac{\sin(\pi 2^m \Delta)}{\sin(\pi \Delta)}\right|^2.$$ Now, we can plot the probability against a normalized $\Delta.$

Notice that there is a lower bound to the probability of measurement of the correct state. Because $|\sin(\pi \Delta)| \leq |\pi \Delta|$ and $|2 (2^m \pi \Delta)| \leq |\sin(\pi 2^m \Delta)|$ when $|\Delta| \leq 1/2,$ $$\frac{1}{2^{2m}}\left|\frac{\sin(\pi 2^m \Delta)}{\sin(\pi \Delta)}\right|^2 \geq \frac{1}{2^{2m}}\frac{|2\cdot 2^m \Delta|^2}{|\pi \Delta|^2}=\boxed{\frac{4}{\pi^2} \approx 0.405}.$$ Extra qubits can bring this probability arbitrarily closer to $1.$

Shor's algorithm

The prime factorization problem is to decompose a number $N$ into its prime factors. To solve this, we can solve a simpler problem: find any two integers $p$ and $q$ such that $p \cdot q = N.$ If $p$ or $q$ is composite, we can further factor them into primes. This process continues recursively until all factors are prime.

The problem of finding two integers $p$ and $q$ such that $p \cdot q = N$ is very difficult for classical computers. In fact, it was considered so difficult that modern internet security was based on it. RSA (Rivest–Shamir–Adleman) is one of the most widely used public-key cryptosystems. The process begins when the receiver selects two very large prime numbers, $p$ and $q,$ and computes their product $N = p\cdot q,$ which is then shared publicly. To encode a message, you only need the public number $N$. However, to decode the message, you need the original primes $p$ and $q$. This allows two parties to communicate securely over a public channel without needing to share any secret information beforehand.

Classical reduction to order-finding

The greatest common divisor (GCD) of two numbers is the largest number that divides both of them without leaving a remainder. We denote the greatest common divisor of $a$ and $b$ as $\gcd(a,b).$ When $\gcd(a,b)=1,$ we say $a$ and $b$ are coprime.

Euclid's algorithm is a classic procedure for finding $\gcd(a,b)$ where $a > b$:

Two-qubit phase estimation

Find the greatest common divisor of $48$ and $18$ using Euclid's algorithm.

$$ \begin{aligned} 48 &\bmod 18 = 12\\ 18 &\bmod 12 = 6\\ 12 &\bmod 6 = 0 \end{aligned} $$ Therefore, $\gcd(48, 18)=6.$ Notice that $48 = 6 \cdot 8,$ $18 = 6 \cdot 3,$ and $\gcd(8,3)=1.$

Find the greatest common divisor of $77$ and $29$ using Euclid's algorithm. Do the numbers $77$ and $29$ share prime factors?

$$ \begin{aligned} 77 &\bmod 29 = 19\\ 29 &\bmod 19 = 10\\ 19 &\bmod 10 = 9\\ 10 &\bmod 9 = 1\\ 9 &\bmod 1 = 0 \end{aligned} $$ Since $\gcd(77, 29)=1,$ the numbers $77$ and $29$ are coprime, which means they share no prime factors.

The first quantum algorithm to get significant attention was Shor's algorithm for factoring, developed by Peter Shor in 1994.[1] This algorithm demonstrated that quantum computers could efficiently factor large integers, a task that is not only extremely challenging for classical computers but also has profound implications for cryptography and data security.

Shor's algorithm relies on a remarkable property of prime numbers: for any product of primes $N=p\cdot q$ and any positive integer $a$, there exists an integer $r$ such that $a^r \equiv 1 \pmod{N}.$ The "quantum part" of Shor's algorithm is responsible for finding the smallest such $r.$

Given a random integer $a$ where $1 < a < N$ and $\gcd(a,N)=1,$ the order finding problem is to find the smallest positive integer $r$ (called the "order" of $a$ modulo $N$) such that $\displaystyle a^r \equiv 1 \pmod{N}.$

Example. Let $N = 5\cdot 7 = 35.$ If we pick $a=11,$ then $r=3$:

$x$ $11^x$ $11^x \bmod 35$
1 11 11
2 121 16
3 1331 1
4 14641 11
5 161051 16
6 1771561 1
7 19487171 11
8 214358881 16
9 2357947691 1
... ... ...

Given an arbitrary number $a$ where $2\leq a

So how does solving order finding relate to factorization? Because $a^r=mN+1$ for some integer $m,$ our number $N$ is a divisor of $$mN=a^r-1=(a^{r/2}-1)(a^{r/2}+1),$$ which is a product of whole numbers, assuming $r/2$ is even. Clearly, $N$ cannot be a divisor of $a^{r/2}-1$ because that would imply that $r/2$ is the actual order. There are then two cases, depending on whether or not $N$ is a divisor of $a^{r/2}+1$:

  1. If neither $a^{r/2} - 1$ nor $a^{r/2} + 1$ is divisible by $N$, yet their product still equals $mN,$ then $N$ divides the product, but not either factor individually. Thus, at least one of these factors must share a nontrivial divisor with $N$, meaning: $$ \gcd(a^{r/2} - 1, N) \quad \text{or} \quad \gcd(a^{r/2} + 1, N) $$ must yield a factor of $N$ that is neither 1 nor $N$.
  2. If $a^{r/2} + 1$ is divisible by $N$, then $\gcd(a^{r/2} + 1, N)=N$ and $\gcd(a^{r/2} - 1, N)=1.$ Since this does not produce any nontrivial factors of $N,$ we must rerun Shor's algorithm with a different $a.$
Quantum phase estimation

Try to factor $N=5\cdot7=35$ using a single iteration of Shor's algorithm. Use $10, 11, 34, 6,$ and $8$ for each guess $a.$

  1. First, we compute $\gcd(10, 35)=5\neq1,$ so $\displaystyle N=5 \cdot \frac{35}{5}=5\cdot 7.$
  2. First, we compute $\gcd(11, 35)=1.$ Next, we find $r=3,$ which is an odd number, so we must repeat the process with a different $a.$
  3. First, we compute $\gcd(34, 35)=1.$ Next, we find $r=2,$ which is an even number. We then compute $\gcd(34^{2/2}-1, 35)=\gcd(33, 35)=1.$ This iteration did not produce any nontrivial factors, so we must repeat the process with a different $a$.
  4. First, we compute $\gcd(6, 35)=1.$ Next, we find $r=2,$ which is an even number. We then compute $\gcd(6^{2/2}-1, 35)=\gcd(5, 35)=5.$ Therefore, $\displaystyle N=5 \cdot \frac{35}{5}=5\cdot 7.$
  5. First, we compute $\gcd(4, 35)=1.$ Next, we find $r=6,$ which is an even number. We then compute $\gcd(4^{6/2}-1, 35)=\gcd(63, 35)=7.$ Therefore, $\displaystyle N=7 \cdot \frac{35}{7}=7\cdot 5.$
  6. First, we compute $\gcd(8, 35)=1.$ Next, we find $r=2,$ which is an even number. We then compute $\gcd(8^{2/2}-1, 35)=\gcd(7, 35)=7.$ Therefore, $\displaystyle N=7 \cdot \frac{35}{7}=7\cdot 5.$

Now, how exactly does Shor's algorithm seem to find $r$ almost magically efficiently?

Quantum order-finding subroutine

Quantum phase estimation: simplified

The quantum subroutine in Shor's algorithm is a form of quantum phase estimation where the unitary in question is $$U\ket{y}=\begin{cases} \ket{a \cdot y \bmod N}, & 0\leq y < N\\ \ket{y}, & N \leq y < 2^n \end{cases}$$ where $n=\lceil \log_2(N) \rceil.$

What are the eigenvalues of $U$?

We know that $$U^r\ket{y}=\begin{cases} \ket{a^r \cdot y \bmod N}, & 0\leq y < N\\ \ket{y}, & N \leq y < 2^n \end{cases}=\ket{y}$$ because $\displaystyle a^r \equiv 1 \pmod{N}.$ Since $U^r=I,$ we have $$ \begin{aligned} U^r \ket{\psi} &=\lambda^r \ket{\psi}\\ I\ket{\psi} &=\lambda^r \ket{\psi}\\ 1 &=\lambda^r, \end{aligned} $$ which implies that the eigenvalues $\lambda$ of $U$ are $r$-th roots of unity: $$\lambda_k=e^{2\pi i k/r}$$ for $k=0,1,2,\ldots,r-1.$

Verify that the states $$\ket{\psi_k}=\frac{1}{\sqrt{r}}\sum_{j=0}^{r-1} e^{-2\pi ikj/r}\ket{a^j \bmod N}$$ for $k = 0,1,2,\ldots,r-1$ are eigenstates of $U.$

Since $$ \begin{aligned} U\ket{\psi_k} &= \frac{1}{\sqrt{r}}\sum_{j=0}^{r-1} e^{-2\pi ikj/r}\ket{a^{j+1} \bmod N}\\ &= \frac{1}{\sqrt{r}}\left(\ket{a \bmod N}+e^{-2\pi ik/r}\ket{a^{2} \bmod N}+\ldots+e^{-2\pi ik(r-1)/r}\ket{a^{r} \bmod N}\right)\\ &= \frac{1}{\sqrt{r}}\left(\ket{a \bmod N}+e^{-2\pi ik/r}\ket{a^{2} \bmod N}+\ldots+e^{-2\pi ik(r-1)/r}\ket{1}\right)\\ &= \frac{1}{\sqrt{r}}\left(e^{-2\pi ik(r-1)/r}\ket{1}+\ket{a \bmod N}+\ldots+e^{-2\pi ik(r-2)/r}\ket{a^{r-1} \bmod N}\right)\\ &= e^{2\pi ik/r}\frac{1}{\sqrt{r}}\left(e^{-2\pi ik(r)/r}\ket{1}+e^{-2\pi ik/r}\ket{a \bmod N}+\ldots+e^{-2\pi ik(r-1)/r}\ket{a^{r-1} \bmod N}\right)\\ &= e^{2\pi ik/r}\frac{1}{\sqrt{r}}\left(\ket{1}+e^{-2\pi ik/r}\ket{a \bmod N}+\ldots+e^{-2\pi ik(r-1)/r}\ket{a^{r-1} \bmod N}\right)\\ &= e^{2\pi ik/r}\left(\frac{1}{\sqrt{r}}\sum_{j=0}^{r-1} e^{-2\pi ikj/r}\ket{a^j \bmod N}\right)\\ &= e^{2\pi ik/r}\ket{\psi_k}, \end{aligned} $$ we know that $\ket{\psi_k}$ are eigenstates of $U$ with corresponding eigenvalues of $e^{2\pi ik/r}.$

Compute the normalized sum of the eigenstates of $U$.

$$ \begin{aligned} \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\ket{\psi_k}&=\frac{1}{\sqrt{r}}\sum_{k=0}^{r- 1}\left(\frac{1}{\sqrt{r}}\sum_{j=0}^{r-1} e^{-2\pi ikj/r}\ket{a^j \bmod N}\right)\\ &=\frac{1}{r}\sum_{k=0}^{r-1}\left(\sum_{j=0}^{r-1} e^{-2\pi ikj/r}\right)\ket{a^j \bmod N}\\ &=\ket{1}+\frac{1}{r}\sum_{k=1}^{r-1}\left(1+e^{-2\pi ik/r}+e^{-4\pi ik/r}+\ldots+e^{-2\pi ik(r-1)/r}\right)\ket{a^j \bmod N}\\ &=\ket{1}+\frac{1}{r}\sum_{k=1}^{r-1}\left(1\frac{1-(e^{-2\pi i k/r})^r}{1-(e^{-2\pi i k/r})}\right)\ket{a^j \bmod N}\\ &=\ket{1} \end{aligned} $$

Let $y_k$ be the nearest integer to $2^m k/r.$ The quantum phase estimation algorithm maps the input state $\ket{0}^{\otimes m}\ket{\psi_k}$ to $\ket{y_k}\ket{\psi_k}$ with probability greater than $4/\pi^2,$ which can be made arbitrarily close to $1$ by adding extra qubits and truncating the output. However, initializing the input register to an eigenstate $\ket{\psi_k}$ paradoxically requires knowing $r,$ the value we seek to discover, in advance!

Luckily, there's a workaround. Applying quantum phase estimation to the input state $\ket{0^m}\ket{1}$ results in the evolution $$\ket{0^m}\ket{1}=\ket{0^m} \left(\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\ket{\psi_k}\right)=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\ket{0^m} \ket{\psi_k} \to \frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\ket{y_k} \ket{\psi_k}.$$

Our measured $y$ has a $1/r$ chance of being any given $y_k.$ The value $y/2^m$ serves as a decimal approximation for $k/r.$ Now, how is $r$ extracted from the value $k/r$ when $k$ is unknown?

Classical post-processing

Theorem. For any given real number $\alpha \in (0,1)$ and integer $N \geq 1,$ there exists at most one pair of coprime integers $p$ and $q$ with $0 \leq p < q < N$ such that $$\left|\alpha - \frac{p}{q}\right|<\frac{1}{Nq}.$$ This theorem is called Dirichlet's theorem on Diophantine approximation.

Given $\alpha \in (0,1)$ and $N \geq 1,$ the continued fraction algorithm finds a pair of coprime integers $p$ and $q$ with $0 \leq p < q < N$ that approximates $\alpha$ such that $$\left|\alpha - \frac{p}{q}\right|<\frac{1}{Nq}.$$ if such a pair exists.

Example. Let $\alpha = 0.328125$ and $N=35.$ First, we compute $$0.328125=\frac{328125}{1000000}.$$ Next, we write it as a mixed fraction and invert the fractional part to get $$0+\frac{1}{\frac{1000000}{328125}}$$ Then, we express $1000000/328125$ as a mixed fraction, resulting in $$0+\frac{1}{3+\frac{15625}{328125}}\approx \frac{1}{3},$$ which is less than $\displaystyle \frac{1}{35 \cdot 3} \approx 0.01$ away from $\alpha,$ so $p=1$ and $q=3.$

Example. Let $\alpha = 0.8359375$ and $N=35.$ First, we compute $$0.8359375=\frac{8359375}{10000000}.$$ Next, we write it as a mixed fraction and invert the fractional part to get $$0+\frac{1}{\frac{10000000}{8359375}}$$ Then, we express $10000000/8359375$ as a mixed fraction, resulting in $$0+\frac{1}{1+\frac{1640625}{8359375}}\approx 1,$$ which is unacceptably far from $\alpha = 0.8359375,$ so we repeat the process: $$\frac{1}{1+\frac{1640625}{8359375}} \to \frac{1}{1+\frac{1}{\frac{8359375}{1640625}}} \to \frac{1}{1+\frac{1}{5+\frac{156250}{1640625}}} \approx \frac{1}{1+\frac{1}{5}}=\frac{5}{6},$$ which is less than $\displaystyle \frac{1}{35 \cdot 6} \approx 0.005$ away from $\alpha,$ so $p=1$ and $q=6.$ If $5/6$ was too far away, we would keep repeating the process until some other fraction $p/q$ is close enough.

Running the continued fraction algorithm yields $$\frac{p}{q}=\frac{k}{r},$$ but that doesn't necessarily imply that $p=k$ or $q=r.$ Because $p$ and $q$ must be coprime, $p$ and $q$ have no nontrivial factors whereas $k$ is free to have common factors with $r$. This means that $p/q$ may be a reduced form of $k/r.$

The least common multiple (LCM) of two numbers is the smallest positive integer that is divisible by both of the numbers. In other words, it's the smallest number that both numbers divide into without leaving a remainder.

If we run the quantum subroutine $s$ many times, we get a list of fraction approximations $$\frac{p_1}{q_1}, \frac{p_2}{q_2}, \ldots, \frac{p_s}{q_s}.$$

Each $q_i$ will be different depending on the common factors of $k_i$ and $r.$ To calculate $r,$ we can compute the least common multiple (LCM)of a ll the $c_i$ values: $$\text{lcm}(c_1, c_2, \ldots, c_s).$$ Note that $\text{lcm}(a,b,c)=\text{lcm}(\text{lcm}(a,b), c)$ and the relationship between LCM and GCD is given by the formula: $$\text{lcm}(a, b) = \frac{|a \cdot b|}{\gcd(a, b)}$$ where $|a \cdot b|$ denotes the absolute value of the product of $a$ and $b$.

Example. If $N=5 \cdot 7 =35$ and we pick $a=4,$ then $r=6.$ After running the quantum subroutine in Shor's algorithm with a $7$-qubit counting register (extra qubits may be added with that output truncated for nearly perfect accuracy), we have the final quantum state $$\frac{1}{\sqrt{6}}\left(\ket{0}+\ket{21}+\ket{43}+\ket{64}+\ket{85}+\ket{107}\right).$$

We compute $$\text{lcm}(3,6,2)=\text{lcm}(\text{lcm}(3,6),2)=\text{lcm}(6,2)=6,$$ revealing that $r=6.$