Chapter 2: Quantum Algorithms

This section goes over the basics of Quantum Algorithms, starting with simple demonstrations of quantum effects.

An accompanying Jupyter notebook with the Qiskit code used to the figures in this chapter is available here.

This chapter and the accompanying Jupyter notebook are currently a work in progress.

So far, the first section introduces simple demonstrations of how quantum algorithms can utilise quantum effects, and the second section covers the central concept of the Fourier Transform and how it applies to quantum algorithms.

Simple Demonstrations of Quantum effects

This section will go through demonstrations of quantum effects in quantum circuits. While these are of limited practical use, they are worth going through to develop a foundational grasp of quantum algorithms.

The first and second subsection cover the Quantum Teleportation and Superdense Coding quantum communication protocols which both make use of Bell states of entangled qubits in different ways. The last subsection describes the Deutsch-Jozsa algorithm which leverages the nature of quantum systems to form superpositions of states.

Quantum Teleportation

Quantum Teleportation is a protocol for sharing quantum information over a distance by using the quantum property of entanglement.

If two parties, Alice and Bob, share a pair of entangled qubits they can use this to 'teleport' the state of another qubit from one party to the other, as shown in the sketch below.

To summarise what's occurring in the quantum teleportation protocol:

1. A pair of qubits $A$ and $B$ are prepared in the Bell state $\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ and shared among two parties, referred to as 'Alice' and 'Bob'.
2. Alice interacts a qubit $Q$, in the state $\alpha\ket{0}+\beta\ket{1}$, with her qubit $A$. This is done by performing a controlled-not operation on $A$ by using qubit $Q$ as the control qubit, followed by performing the Hadamard gate on qubit $Q$.
3. Alice measures her two qubits, then sends the results to Bob.
4. Depending on the result of the measurement, Bob performs an X and/or Z gate on his qubit $B$, which modifies its state to recover qubit state $\alpha\ket{0}+\beta\ket{1}$. Specifically, if Alice measures $A$ as $\ket{1}$ then the $X$ gate is applied on qubit $B$, and if Alice measures $Q$ as $\ket{1}$ then the $Z$ gate is applied on qubit $B$.

To get a better grasp of how the Quantum Teleportation protocol functions, it is worth looking at the annotated Qiskit circuit below and go through the states at different points along the circuit.

At the beginning of the circuit, before the Bell state has been initialised, the qubit can be represented by the composite state $\ket{\psi_0}$ $$\psi_0=(\alpha\ket{0}+\beta\ket{1})\otimes\ket{0}\otimes\ket{0}$$

Bell state then is generated on $A\otimes B=\ket{00}$.

First, the Hadamard gate is applied on qubit $A$ $$\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})(\ket{0})$$ Which can be reformulated for ease of calculation as $$\frac{1}{\sqrt{2}}(\ket{00}+\ket{10})$$

Then, a $C_{NOT}$ gate is applied to qubit $B$ using qubit $A$ as the control qubit, resulting in qubits $A$ and $B$ forming the following Bell state $$\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$$

This results in the following prepared state for the Quantum Teleportation protocol $$\psi_1 = \frac{1}{\sqrt{2}}(\alpha\ket{0}+\beta\ket{1})(\ket{00}+\ket{11})$$ Which can be expanded as $$\psi_1 = \frac{1}{\sqrt{2}}(\alpha\ket{000} + \alpha\ket{011} + \beta\ket{100} + \beta\ket{111})$$

As the first part of the protocol, Alice applies the $C_{NOT}$ gate on the second qubit $A$, using the first qubit $Q$ as the control qubit. $$\frac{1}{\sqrt{2}}(\alpha\ket{000} + \alpha\ket{011} + \beta\ket{110} + \beta\ket{101})$$ Which again can be rearranged as $$\frac{1}{\sqrt{2}}(\alpha\ket{0}(\ket{00} + \ket{11}) + \beta\ket{1}(\ket{10} + \ket{01}))$$

Then the $H$ gate is applied on the first qubit leading to $\ket{\psi_2}$ $$\ket{\psi_2} = \frac{1}{2}\alpha(\ket{0}+\ket{1})(\ket{00} + \ket{11}) + \beta(\ket{0}-\ket{1})(\ket{10} + \ket{01}))$$

$\ket{\psi_2}$ can be rearranged by grouping together the first and second qubits $Q$ and $A$, both of which Alice possesses.

The first two qubits are measured by Alice, also breaking this Bell pair, and this classical data is shared with Bob, as noted by the measurement gates connected to the classical channels $a$ and $b$. This results in the composite qubit state $\psi_3$ being one of the four following states, each with a probability $(\frac{1}{2})^{2}=\frac{1}{4}$

An interesting thing can be seen from these states; the third qubit $B$ can be changed to the initial state $\alpha\ket{0} + \beta\ket{1}$ through the conditional application of $X$ and $Z$ gates. By applying a phase-flip on the third qubit $B$ via the $Z$ gate if the first qubit $Q$ is measured as $\ket{1}$, and a bit-flip on $B$ via the $X$ gate if the second qubit $A$ is measured as $\ket{1}$, we can recover $\alpha\ket{0} + \beta\ket{1}$ on qubit $B$. With this, we have achieved a teleportation of states from Alice to Bob.

Now of course, since this uses a classical channel to share information between Alice and Bob, the speed of transfer is limited by the speed of this classical communication, with an upper limit of the speed of light.

Superdense Coding

Another quantum communication protocol which makes use of Bell states is Superdense Coding.

Superdense Coding involves two parties, Alice and Bob, initially sharing one qubit each of the Bell state $\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$. To encode two bits of classical information Alice performs two gate operations on her qubit, before sending her qubit to Bob, who after simple gate operation can measure the two bit Alice encoded. This protocol is called 'Superdense' because while only one qubit is sent, two bits of information are transferred from Alice to Bob, contrary to what we would expect classically. It is visualised in the diagram below.

To briefly explain what happens in Superdense Coding:

1. A pair of qubits $A$ and $B$ are prepared in the Bell state $\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ and shared among two parties, referred to as 'Alice' and 'Bob'.
2. Alice encodes information on $A$ by applying $X$ and $Z$ gate conditionally based on the classical information she wishes to send. If she wishes to encode $1$ for her first classical bit $a$ she applies an $X$ gate (and does not apply it if she wishes to encode $0$), and similarly if she wishes to encode $1$ for her second classical bit $b$ she applies an $Z$ gate (and does not apply it is she wishes to encode $0$).
3. Alice sends her qubit $A$ to Bob.
4. Bob transforms the qubits $A$ and $B$ by applying a $C_{NOT}$ gate and $H$ gate; first he applies a $C_{NOT}$ gate with $A$ as the target qubit and $B$ as the control qubit, then he applies $H$ on the second qubit $B$.
5. Bob measures qubits $A$ and $B$ in the computational basis to recover classical bits $a$ and $b$ respectively.

To get a better grasp of how the Superdense Coding protocol functions, it is worth looking at the annotated Qiskit circuit below and understand the evolution of the Bell state.

Similarly to the Quantum Teleportation protocol, a Bell state $\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ is prepared beforehand with a $H$ gate applied on Alice's qubit $A$, followed by a $C_{NOT}$ gate applied on Bob's qubit $B$ with $A$ as the control qubit.

The conditional application of $X$ and $Z$ gates changes the Bell state $\beta_{00}$ to other Bell states.

If Alice decides to encode classical bits $a=0$ and $b=0$, then she needn't do anything, and we are left with $$\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$$

If she instead encodes $a=1$ and $b=0$, then the $X$ gate is applies on $A$ to cause a bit-flip, leading to the $\beta_{10}$ state $$\beta_{10}=\frac{1}{\sqrt{2}}(\ket{10}+\ket{01})$$

If Alice encodes $a=0$ and $b=1$, then the $Z$ gate is applies on $A$ to cause a phase-flip, leading to the $\beta_{01}$ state $$\beta_{01}=\frac{1}{\sqrt{2}}(\ket{00}-\ket{11})$$

Finally if $a=1$ and $b=1$ is encoded then both $X$ and $Z$ gates are applied on qubit $A$, leading to the $\beta_{11}$ state $$\beta_{11}=\frac{1}{\sqrt{2}}(\ket{01}-\ket{10})$$

After this, Alice sends qubit $A$ to Bob, who now has full access to both qubits in the Bell state. Since the Bell states form an orthonormal basis, Bob can distinguish each state. To measure in the computational basis though, he first has to transform them with unitary gates: a $C_{NOT}$ applied on the first qubit $A$ with the second qubit $B$ as the control qubit, followed by a Hadamard gate applied on $B$.

For $\beta_{00}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ the $C_{NOT}$ leads to $$\frac{1}{\sqrt{2}}(\ket{00}+\ket{01})=\ket{0}\otimes\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})=\ket{0}\ket{+}$$ Which with a Hadamard gate applied to $B$ leads to $$\ket{0}\ket{0}$$ Which can be be measured in the computational basis to return $a=0$ and $b=0$ for $A$ and $B$ respectively.

For $\beta_{10}=\frac{1}{\sqrt{2}}(\ket{10}+\ket{01})$ the $C_{NOT}$ applied to $A$ leads to $$\frac{1}{\sqrt{2}}(\ket{10}+\ket{11})=\ket{1}\otimes\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})=\ket{1}\ket{+}$$ Which with $H$ applied to $B$ leads to $$\ket{1}\ket{0}$$ Which upon measurement leads to $a=1$ and $b=0$.

For $\beta_{01}=\frac{1}{\sqrt{2}}(\ket{00}-\ket{11})$ the $C_{NOT}$ applied to $A$ results in $$\frac{1}{\sqrt{2}}(\ket{00}-\ket{01})=\ket{0}\otimes\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})=\ket{0}\ket{-}$$ Which with $H$ applied to $B$ leads to $$\ket{0}\ket{1}$$ Which upon measurement leads to $a=0$ and $b=1$.

Finally for $\beta_{11}=\frac{1}{\sqrt{2}}(\ket{01}-\ket{10})$ the $C_{NOT}$ applied to $A$ results in $$\frac{1}{\sqrt{2}}(\ket{11}-\ket{10})=-\ket{1}\otimes{\sqrt{2}}(\ket{0}-\ket{1})=-\ket{1}\ket{-}$$ Which with $H$ applied to $B$ leads to $$-\ket{1}\ket{1}$$ Which upon measurement leads to $a=1$ and $b=1$.

Deutsch-Jozsa Algorithm

To conclude this section on simple demonstrations of quantum effects, we'll cover the Deutsch-Jozsa algorithm. The Deutsch-Jozsa algorithm demonstrates the advantage quantum circuits to place qubits to be in a superposition of states, which can be operated on in parallel and also interfere with each other.

Though first, we'll build up to the algorithm from circuits to develop an understanding of it. The algorithm deals with a unitary 'black boxes' operation acting on multiple qubits. Such a black box is also known as an 'oracle' and can be formulated as seen in the sketch below.

The unitary black box $U_f$ acts on the target qubit $\ket{y}$, with a function dependent on the data qubit $\ket{x}$.

Note, that the $\oplus$ is modulo two addition operation, which only returns '1' if exactly one of the two bits it is acting on is '1' (and the other is '0'). This also means that the following property holds regardless of whether $k$ is $0$ or $1$ $$0 \oplus k = k$$

In the case of the initial state $\ket{x}=\ket{+}$ and $\ket{y}=\ket{0}$ the resulting state is

Such a state simultaneously holds information on the results $f(0)$ and $f(1)$, as if $f(x)$ was evaluated for two possible values in parallel. Such a feature is known as quantum parallelism.

Though there is a caveat to keep in mind here: upon measurement this superposition collapses, into $\ket{0}\ket{f(0)}$ or $\ket{0}\ket{f(1)}$.

The Hadamard gate is crucial for setting up state in a superposition to exploit quantum parallelism, and a general operation on multiple qubits is is known as the Hadamard Transform (or Walsh-Hadamard Transform). The Hadamard transform involves simply applying a $H$ gate on each qubit to place them in a superposition as seen below with three qubits initially in a $\ket{0}$ state

The resulting state is

Generalising this for $n$ qubits, we can represent the operation as $H^{\otimes n}$ acting on a state $\ket{x_1, x_2, ... x_n}$ resulting in a superposition of states $\ket{z_1, z_2, ... z_n}$. Such a superposition can be represented by the following sum $$H^{\otimes n} \ket{x_1, ... x_n} = \frac{1}{\sqrt{2^n}} \sum_{z} (-1)^{x_1 z_1 + ... + x_n z_n} \ket{z_1, ... z_n}$$ Which can be represented more succinctly by dropping the subscripts and using a dot product between wavevector state values $x \cdot z$ $$H^{\otimes n} \ket{x} = \frac{1}{\sqrt{2^n}} \sum_{z} (-1)^{x \cdot z} \ket{z}$$

The coefficient of $(-1)^{x \cdot z}$ in the sum above can be explained by the following: when $z_i=0$ then the coefficient is $+1$ such as $\ket{0}$ within $\ket{+}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ and $\ket{-}=\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})$ (for $n=1$), negative factors only appear when $x_i=1$ such as $\frac{-\ket{1}}{\sqrt{2}}$ within $H\ket{1}=\ket{-}=\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})$, and the coefficient is only negative when the dot product $x \cdot z$ is odd since an even number of $-1$ factors cancel out to $+1$.

Now, let's move onto the Deutsch Algorithm which builds upon quantum parallelism by including interference between superpositions. The Deutsch Algorithm in particular determines whether an unknown function $f(x)$ with possible inputs $x={0,1}$, and two possible outputs of $0$ and $1$, is constant or balanced. If a function is constant, then this results in the same output for all inputs, so in this case $f(0)=f(1)$. Otherwise if a function is balanced, then there is an equal likelihood of all outputs, which in this case with two possible inputs leads to $f(0)\neq f(1)$.

The advantage of the Deutsch algorithm is that we only have to evaluate $f(x)$ once, as opposed to twice for each possible input $x={0,1}$ to determine this global property of this function, as shown in the circuit diagram shown below.

The input state is $$\ket{\psi_0} = \ket{01}$$ Upon applying the two Hadamard gates it becomes $$\ket{\psi_1} = \frac{\ket{0}+\ket{1}}{\sqrt{2}} \frac{\ket{0}-\ket{1}}{\sqrt{2}}$$

The unitary gate $U_f$ then changes the target qubit $\ket{y}$ into $\ket{y\oplus f(x)}$, applying a function $f(x)$ dependent on the data qubit state $\ket{x}$. $$\ket{\psi_2} = \ket{x}\ket{y\oplus f(x)}$$ In the case of this algorithm $\ket{x}=\frac{1}{\sqrt{2}}(\ket{0}+\ket{1})$ and $\ket{y}=\frac{1}{\sqrt{2}}(\ket{0}-\ket{1})$, so $$\ket{\psi_2}=\frac{1}{\sqrt{2}} (\ket{0} \frac{\ket{0\oplus f(0)} - \ket{1\oplus f(0)}}{\sqrt{2}} + \ket{1} \frac{\ket{0\oplus f(1)} - \ket{1\oplus f(1)}}{\sqrt{2}} )$$

To simplify this we can use the following relation $$\frac{1}{\sqrt{2}} (\ket{0\oplus k} - \ket{1\oplus k} ) = (-1)^{k} \ket{-}$$ Which holds for $k=0$

And for $k=1$

Substituting the relation into $\ket{\psi_2}$ results in

So

Upon applying the Hadamard gate to the data qubit, which is the first qubit, this leads to $$\ket{\psi_3} = \begin{cases} \pm \ket{0} \ket{-} & \text{if } f(0) = f(1)\\ \pm \ket{1} \ket{-} & \text{if } f(0) \neq f(1) \end{cases}$$

Upon measurement of the data qubit $\ket{x'}$ in $\ket{\psi_4}=\ket{x'}\ket{-}$, we can determine whether the function $f(x)$ is constant ($f(0)=f(1)$) or balanced ($f(0)\neq f(1)$) $$x' = \begin{cases} 0 & \text{if } f(0) = f(1) \text{, hence constant} \\ 1 & \text{if } f(0) \neq f(1) \text{, hence balanced} \end{cases}$$

Thanks to the interference between the two possible function outcomes $f(0)$ and $f(1)$, in one measurement we can determine the global property of whether the function is constant or balanced. This interference can be seen in the $(-1)^{f(1)-f(0)}$ term within $\ket{\psi_2}$. This means that after transformation to $\ket{\psi_3}$ and measurement of the data qubit, if the output is $0$ then $f(x)$ is constant, and if the output is $1$ then $f(x)$ is balanced.

Now, let's generalise this from merely one data qubit to $n$ data qubits in the Deutsch-Jozsa algorithm, as depicted below which now depends on $x$ being a value which varies from $0$ to $2^{n}-1$.

The input state is $$\ket{\psi_0}=\ket{0}^{\otimes n} \ket{1}$$

After the Hadamard transform it becomes $$\ket{\psi_1}= H^{\otimes n}\ket{0}^{\otimes n} \otimes \frac{\ket{0}-\ket{1}}{\sqrt{2}}$$

We can substitute in the sum of superpositions we explained earlier $$H^{\otimes n} \ket{x} = \frac{1}{\sqrt{2^n}} \sum_{z} (-1)^{x \cdot z} \ket{z}$$ Which for $n$ qubits in the $\ket{0}$ state leads to an equal superposition of states $$H^{\otimes n} \ket{0}^{\otimes n} = \frac{1}{\sqrt{2^n}} \sum_{z} \ket{z}$$ So by substituting this in, and relabeling $z$ as $x$, we get the state $$\ket{\psi_1}= \frac{1}{\sqrt{2^n}} \sum_{x} \ket{x} \frac{\ket{0}-\ket{1}}{\sqrt{2}}$$

Upon application of the unitary gate $U_f$, the function $f(x)$ is evaluated so $\ket{x, y}$ evolves into $\ket{x, y\oplus f(x)}$ as seen in the Deutsch algorithm, this time leading to $$\ket{\psi_2} = \frac{1}{\sqrt{2^n}}\sum_{x} \ket{x} \frac{\ket{0\oplus f(x)}-\ket{1 \oplus f(x)}}{\sqrt{2}} $$ Similarly we can again use the following relation $$\frac{1}{\sqrt{2}} (\ket{0\oplus k} - \ket{1\oplus k} ) = (-1)^{k} \ket{-}$$ So the state $\ket{\psi_2}$ is $$\ket{\psi_2} = \frac{1}{\sqrt{2^n}}\sum_{x} (-1)^{f(x)} \ket{x} \ket{-}$$

We can again use the following relation $$H^{\otimes n} \ket{x} = \frac{1}{\sqrt{2^n}} \sum_{z} (-1)^{x \cdot z} \ket{z}$$ So that after the application of $n$ Hadamard gates, $H^{\otimes n}$, the state becomes $$\ket{\psi_3} = \frac{1}{\sqrt{2^n}}\sum_{x} (-1)^{x\cdot z f(x)} \ket{x} \ket{-}$$

Then $n$ measurements are carried out to in the computational basis. If we get a bitstring of all $0$s then the function $f(x)$ is constant, otherwise if there is even one $1$ then the function is balanced.

Fourier Transform

Classical Fourier Transform

One way of solving mathematical problems is to transform it into a problem which is easier to solve, such as the Fourier transform.

The Fourier transform is a centrally important and widely used mathematical tool. In the case of signal processing of an audio or electromagnetic wave, it can convert a signal which varies in time (also called the 'time domain') into a visualisation of the component frequencies (called the 'frequency domain').

So, for a time varying equation $f(t)$ can be transformed into a frequency varying $F(\xi)$ through the Fourier transform, where $\xi$ is frequency and $t$ is time, as shown below

$$F(\xi) = \int_{-\infty}^{\infty} f(t) e^{2\pi i\xi t} dt$$

The inverse Fourier transform inverts this, yielding the time domain $f(t)$ from the Frequency domain $F(\xi)$

$$f(t) = \int_{-\infty}^{\infty} F(\xi) e^{-2\pi i \xi t} dt$$

To better understand this, it's worth demonstrating the Fourier transform of a superposition of waveforms as seen below.

Before performing Fourier Transforms on the above waveforms, it's worth noting a rather useful integral identity

$$\delta (a - b) = \int_{-\infty}^{\infty} e^{2 \pi i (a-b) x } dx $$

Where $\delta (a-b)$ is the Dirac delta, which is only equal to 1 if its argument, $a-b$, is equal to $0$, and is otherwise zero.

Let's start with the Fourier Transform of the first waveform which is a sine wave with a frequency of $\xi_1=5$ Hz.

$$F_{1}(\omega) = \int_{-\infty}^{\infty} \text{sin}(2\pi \xi_1 t) e^{2 \pi i \xi t} dt$$

The sine equation can also be re-expressed as $$\text{sin}(\theta) = \frac{1}{2 i} (e^{i\theta}-e^{-i\theta})$$ Which when substituted into $F_{1}(\xi)$ leads to

\begin{align} F_{1}(\xi) &= \frac{1}{2 i} \int_{-\infty}^{\infty} (e^{2\pi i\xi_{1}t}-e^{-2\pi i\xi_{1} t}) e^{2\pi i\xi t} dt \\ &= \frac{1}{2 i}\int_{-\infty}^{\infty} e^{2\pi i(\xi+\xi_{1})t} dt - \frac{1}{2 i}\int_{-\infty}^{\infty} e^{2\pi i(\xi-\xi_{1})t} dt \\ \end{align}

Then, we can substitute the aforementioned integral identity leading to the Dirac delta $$F_{1}(\omega) = \frac{1}{2i} \delta(\xi+\xi_{1}) - \frac{1}{2i}\delta(\xi-\xi_{1})$$ Which when using $\frac{1}{i}=-i$ turns into $$F_{1}(\omega) = -\frac{i}{2}\delta(\xi+\xi_{1}) + \frac{i}{2}\delta(\xi-\xi_{1})$$

When this Fourier transform $F_{1}$ is plotted, it leads to

A similar conclusion can be drawn from the third waveform $w_{3}=1.2 \text{sin}(60\pi t)$ since it is also a sine wave, except with an amplitude factor of $1.2$ which leads to peaks of magnitude $1.2$ at $|\xi|=\xi_{3}=30$.

For cosine waveforms, such as $w_{2}=\text{cos}(30\pi t)$ we have to instead substitute in the following identity $$\text{cos}(\theta) = \frac{1}{2} (e^{i\theta}+e^{-i\theta})$$ Following a similar derivation as above, and using $\omega_{2}=30\pi$ leads to

\begin{align} F_{2}(\xi) &= \int_{-\infty}^{\infty} \text{cos}(2\pi\xi_{2} t) e^{2\pi i\xi t} dt\\ &= \frac{1}{2} \int_{-\infty}^{\infty} (e^{2\pi i\xi_{2}t}+e^{-2\pi i\xi_{2} t}) e^{2\pi i\xi t} dt \\ &= \frac{1}{2} \int_{-\infty}^{\infty} e^{2\pi i(\xi+\xi_{2})t} dt + \frac{1}{2} \int_{-\infty}^{\infty} e^{2\pi i(\xi-\xi_{2})t} dt\\ &= \frac{1}{2} \delta(\xi+\xi_{2}) + \frac{1}{2} \delta(\xi-\xi_{2}) \end{align}

With the linear nature of Fourier transform, we can sum the Fourier transforms of the component waveforms to arrive at the Fourier transform of the total waveform $F_{t}(\omega)$

\begin{align} w &= w_1 + w_2 + w_3 \\ F(w) &= F(w_1) + F(w_2) + F(w_3) \\ F_{t}(\omega) &= F_{1}(\omega) + F_{2}(\omega) + F_{3}(\omega) \\ F_{t}(\omega) &= -0.5i\delta(\xi-5) + 0.5i\delta(\xi+5) \\ &+ 0.5 \delta(\xi-15) + 0.5 \delta(\xi+15) \\ &-0.6 i\delta(\xi-30) + 0.6 i\delta(\xi+30) \end{align}

Due to the complex nature of this equation, containing both imaginary and real terms, the following visualisation is of the magnitude $|F_{t}(\omega)|$

While we were working with continuous integrals with infinite limits to grasp the basics of Fourier Transforms, in reality when applying the Fourier Transform to real data we will work with a discrete and finite dataset with $N$ elements.

For this, we need to discretise the continuous variables $t$ and $\xi$.

For $t$, this is $$t = j \Delta t \text{ for } j=0, 1, ..., N-1$$

As for the frequency term, we can define its discrete spacing as $\Delta \xi=\frac{1}{N\Delta t}$

Then, for $\xi$ it is discretised as $$\xi = \frac{k}{N\Delta t} \text{ for } k=0, 1, ..., N-1$$

Finally, we need to change from a continuous integral into a discrete sum, so for the Fourier Transform we start with $$F_{k} = \sum_{j=0}^{N-1} f_{j} e^{- 2\pi i \Delta \xi \Delta t}$$

Substituting in $\xi$ and $t$ results in \begin{align} F_{k} &= \sum_{j=0}^{N-1} f_{j} e^{2\pi i \frac{k}{N\Delta t} j \Delta t} \Delta t\\ F_{k} &= \Delta t \sum_{j=0}^{N-1} f_{j} e^{2\pi i \frac{k j}{N}} \end{align}

Repeating this with the inverse Fourier Transform results in \begin{align} f_{j} &= \sum_{k=0}^{N-1} F_{k} e^{-2\pi i \frac{k}{N\Delta t} j \Delta t} \Delta \xi \\ &= \frac{1}{N \Delta t} \sum_{k=0}^{N-1} F_{k} e^{-2\pi i \frac{k j}{N}} \end{align}

To make this Fourier and inverse Fourier transform pair unitary, we will modify these so that both sums share the same factor of $\frac{1}{\sqrt{N}}$ \begin{align} F_{k} &= \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1} f_{j} e^{2\pi i \frac{k j}{N}}\\ f_{j} &= \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} F_{k} e^{-2\pi i \frac{k j}{N}} \end{align}

Quantum Fourier Transform

In this subsection we will go over the quantum Fourier transform (QFT) circuit, building atop the introduction to the Fourier transform in the preceding subsection. The QFT circuit is an efficient algorithm for applying the Fourier transform to the amplitudes of qubits, and is pivotal to further applications as explored in the following subsections.

Within this subsection, we will first derive QFT, demonstrate how a 3 qubit state evolves under QFT, then derive the inverse QFT.

Beyond this, a demonstration of the QFT and inverse QFT in Qiskit is available for further viewing in the following Jupyter notebook.

Derivation of QFT Circuit

QFT effectively applies the forward discrete Fourier transform by acting as an operator $$\ket{j} \rightarrow \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1} e^{2\pi i \frac{k j}{N}} \ket{k} $$

When applied to a more general state, QFT generalises to $$\sum_{j=0}^{N-1} x_{j}\ket{j} \rightarrow \sum_{k=0}^{N-1} y_{k}\ket{k}$$ Where $$y_{k} = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1} x_{j} e^{2\pi i \frac{k j}{N}} $$

For QFT performed on $n$ qubits, $N=2^{n}$ and the transformation occurs with regards to an orthonormal basis state ranging from $\ket{0}$ to $\ket{2^{n}-1}$.

Represented in operator form, QFT is then $$\hat{F}_{QFT}= \frac{1}{\sqrt{2^n}}\sum^{2^{n}-1}_{k=0} \sum^{2^{n}-1}_{j=0} e^{2\pi i \frac{jk}{2^{n}}} \ket{k}\bra{j}$$

The number $j$ represented by a state $\ket{j}$ can be presented in binary form, as $j = j_1 j_2 \ldots j_n$ which can be expanded out into $j = j_1 2^{n-1} + j_2 2^{n-2} + \ldots + j_n 2^{0}=\sum^{n}_{l=1} j_l 2^{n-l}$. Binary representation also extends past the decimal point, with $0.j_{l} j_{l+1} \ldots j_{m} = \frac{j_{l}}{2^1} + \frac{j_{l+1}}{2^2} + \ldots + \frac{j_m}{2^{m-l+1}} = \sum_{p=l}^{m} j_{p} 2^{-(p-l+1)}$.

Using this binary representation for $k=k_1 k_2 \ldots k_n=\sum^{n-1}_{l=0} k_l 2^{n-l} $ , one can rearrange the state resulting from application of QFT to $\ket{j}$ as follows: \begin{align} \hat{F}_{QFT} \ket{j} &= \frac{1}{\sqrt{2^n}} \sum^{2^{n}-1}_{k=0} \exp \left( 2\pi i \frac{k j}{2^{n}}\right) \ket{k} \\ &= \frac{1}{\sqrt{2^n}} \sum^{2^{n}-1}_{k=0} \exp \left(2\pi i \frac{j}{2^{n}} \sum^{n}_{l=1} k_l 2^{n-l} \right) \ket{k} \\ &= \frac{1}{\sqrt{2^n}} \sum^{2^{n}-1}_{k=0} \exp \left(2\pi i j\sum^{n}_{l=1} k_l 2^{-l} \right) \ket{k} \\ &= \frac{1}{\sqrt{2^n}} \sum^{2^{n}-1}_{k=0} \prod^{n}_{l=1} \left[ \exp \left(2\pi i j k_l 2^{-l} \right) \right]\ket{k} \end{align} Above, the exponential term was rearranged and the product symbol $\prod$ was used, which acts as $\prod_{l=1}^{m}x_{l} = x_{1} \times x_{2} \times \ldots \times x_{m}$. Continuing on with the rearrangement, this time targeting the sum symbol and ket vector to replace $k$ with its binary representation: \begin{align} \hat{F}_{QFT} \ket{j} &= \frac{1}{\sqrt{2^n}} \sum^{1}_{k_{1}=0} \sum^{1}_{k_{2}=0} \ldots \sum^{1}_{k_{n}=0} \prod^{n}_{l=1} \left[ e^{2\pi i j k_l 2^{-l}} \right]\ket{k_{1} k_{2}\ldots k_{n}}\\ &= \frac{1}{\sqrt{2^n}} \sum^{1}_{k_{1}=0} \sum^{1}_{k_{2}=0} \ldots \sum^{1}_{k_{n}=0} \bigotimes^{n}_{l=1} \left[ e^{2\pi i j k_l 2^{-l}} \ket{k_{l}}\right] \end{align} Above $\bigotimes$ is a tensor product symbol, which multiplies together the scalar terms ($e^{2\pi i j k_l 2^{-l}}$) and performs a tensor of the component vectors $\ket{k_l}$ (which reassemble $\ket{k}$). As a reminder, tensor products were explained in the previous chapter. The sums over $k_l$ can then be moved into the tensor product and expanded as follows: \begin{align} \hat{F}_{QFT} \ket{j} &= \frac{1}{\sqrt{2^n}} \bigotimes^{n}_{l=1} \left[ \sum^{1}_{k_{l}=0}e^{2\pi i j k_l 2^{-l}} \ket{k_{l}}\right]\\ &= \frac{1}{\sqrt{2^n}} \bigotimes^{n}_{l=1} \left[ \ket{0} + e^{2\pi i j 2^{-l}} \ket{1}\right] \end{align}

The exponential coefficient of $\ket{1}$ can be expanded by substituting in $j=\sum^{n}_{p=1} j_p 2^{n-p}$ , as \begin{align} \exp \left( 2\pi i j 2^{-l} \right) &= \exp \left( 2\pi i \sum^{n}_{p=1} j_p 2^{n-p} 2^{-l} \right) \\ &= \exp \left( 2\pi i \sum^{n}_{p=1} j_p 2^{(n-l)-p} \right) \\ &= \exp \left( 2\pi i \left[ \sum^{n-l}_{p=1} j_p 2^{(n-l)-p} + \sum^{n}_{p=n-l+1} j_p 2^{(n-l)-p} \right] \right) \end{align} We can observe that $e^{2\pi i x} = 1$ for integer values of $x$, and the first sum only contains integer values (ranging from $j_1 2^{n-l-1}$ to $j_{n-l} 2^{0}$) so the whole sum effectively cancels out as a factor of unity $1$ and can be removed. \begin{align} \exp \left( 2\pi i j 2^{-l} \right) &= \exp \left( 2\pi i j \left[ \sum^{n}_{p=n-l+1} j_p 2^{(n-l)-p} \right] \right) \\ &= \exp \left( 2\pi i j \left[ \frac{j_{n-l+1}}{2^{1}} + \frac{j_{n-l+2}}{2^{2}} + \ldots + \frac{j_{n}}{2^{l}} \right] \right) \end{align} Now, using the binary fraction representation we are left with different exponential terms depending on the value of $l$. Keep in mind that $j_p$ values can only go up to $j_n$, and that further values such as $j_{n+1}$ are dropped. \begin{align} \text{if } l=1 \text{: } & \exp \left( 2\pi i j \left[\frac{j_{n}}{2} \right] \right)\\ &= e^{(2\pi i) 0.j_{n}}\\ \text{if } l=2 \text{: } & \exp \left( 2\pi i j \left[\frac{j_{n-1}}{2} + \frac{j_{n}}{2^{2}} \right] \right)\\ &= e^{(2\pi i) 0.j_{n-1}j_{n}}\\ \text{if } l=n \text{: } & \exp \left( 2\pi i j \left[\frac{j_{1}}{2} + \frac{j_{2}}{2^{2}} + \ldots \frac{j_{n}}{2^{n}} \right] \right)\\ &= e^{(2\pi i) 0.j_{1}j_{2}\ldots j_{n}} \end{align} Substituting this exponential term back into $\hat{F}_{QFT}\ket{j}$, it can be expanded out as follows \begin{align} \hat{F}_{QFT}\ket{j} &= \frac{1}{\sqrt{2^n}} \left( \ket{0} + e^{(2\pi i) 0.j_{n} } \ket{1}\right) \left( \ket{0} + e^{(2\pi i) 0.j_{n-1}j_{n} } \ket{1}\right) \ldots \left( \ket{0} + e^{(2\pi i) 0.j_{1}j_{2}\ldots j_{n} } \ket{1}\right) \end{align}

Now that we have an equation demonstrating how QFT evolves an $n$ qubit state $\ket{j}$, the next step is implement this evolution with gates.

Before starting with this, let's define a control-phase gate which applies the following operation $R_{k}$ on a target qubit if the control qubit is $1$: \begin{align} R_k = \begin{bmatrix} 1 & 0 \\ 0 & e^{2 \pi i /2^{k}} \end{bmatrix} \end{align}

Starting with the state $$\ket{j} = \ket{j_{1} j_{2} \ldots j_{n}} = \ket{j_1} \otimes \ket{j_2 \ldots j_n}$$

Let's first apply a Hadamard gate to the first qubit, resulting in $$H\ket{j_1} = \frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_1}\ket{1}\right) \otimes \ket{j_2 \ldots j_n}$$

The above operation can be seen to hold for two possible values of $j_1$. If $j_1=1$ then $e^{2 \pi i \frac{j_1}{2}}=e^{2 \pi i \frac{1}{2}} =e^{\pi i} = -1$ so $H\ket{1} = \frac{1}{2}(\ket{0}-\ket{1})$ as expected of the Hadamard gate. If instead $j_1=0$ then $e^{2 \pi i \frac{j_1}{2}}=e^{2 \pi i \frac{0}{2}} =e^{0} = 1$, leading to $H\ket{0} = \frac{1}{2}(\ket{0}+\ket{1})$ as also expected of the Hadamard gate.

Now let's apply $R_2$ (that is with $k=2$) onto the 1st qubit with the second qubit state $\ket{j_2}$ as the control qubit. To keep things simple, states $\ket{j_2 \ldots j_n}$ won't be shown here: \begin{align} R_{2} \frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_1}\ket{1} \right) &= \begin{bmatrix} 1 & 0 \\ 0 & e^{2 \pi i /2^{2}} \end{bmatrix} \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ e^{2 \pi i 0.j_1} \end{bmatrix} \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ e^{2 \pi i 0.j_1 } e^{2 \pi i 0.0j_{2}} \end{bmatrix} \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ e^{2 \pi i (0.j_1 + 0.0j_{2})} \end{bmatrix} \\ &= \frac{1}{\sqrt{2}} \left( \ket{0} + \ket{1} e^{2 \pi i 0.j_1 j_2} \right) \\ \end{align}

This controlled phase operation repeats for the 3rd qubit, except with a controlled $R_3$ gate. This repeats as $p$ increments to $n$, for control qubit $p$ and a controlled $R_p$ gate, leading to the following state for the first qubit $$\frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_1 j_2\ldots j_n} \ket{1} \right)$$

The same course of action occurs for the 2nd qubit, starting with a Hadamard gate and then applying a controlled $R_2$ gate with the 3rd qubit state $\ket{j_3}$ as the control state. Similarly this continues as $p$ increments all the way to $p=n$ for a controlled $R_{p-1}$ operation controlled by state $\ket{j_{p}}$, leading to $$\frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_2 \ldots j_n}\ket{1} \right)$$

As this repeats for following qubits, the number of controlled phase operations decreases, with the $n-1$ qubit state of $\ket{j_{n-1}}$ only having a single $R_2$ operation controlled by $\ket{j_n}$ following the application of a Hadamard, leading to $$\frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_{n-1} j_n} \ket{1} \right)$$

And for the last qubit only a Hadamard gate is need to transform it $$H\ket{j_n} = \frac{1}{\sqrt{2}} \left( \ket{0} + e^{2 \pi i 0.j_n} \ket{1} \right) $$

This results in the following composite state, including all $n$ qubits $$\frac{1}{\sqrt{2^{n}}} \left( \ket{0} + e^{2 \pi i 0.j_1 \ldots j_n} \ket{1} \right) \left( \ket{0} + e^{2 \pi i 0.j_2 \ldots j_n} \ket{1} \right) \ldots \left( \ket{0} + e^{2 \pi i 0.j_n} \ket{1} \right) $$

Note, the order of the states is inverted compared to what $\hat{F}_{QFT}\ket{j}$ results in, so after this swap gates are used to invert the order of the qubit states, resulting in the expected state: $$\hat{F}_{QFT}\ket{j} = \frac{1}{\sqrt{2^n}} \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{n} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{n-1}j_{n} } \right) \ldots \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{1}j_{2}\ldots j_{n} } \right)$$

The gate operations so far, excluding the swap operations (to keep the size of the figure down), are shown below:

Demonstration of QFT acting on a 3 qubit state

To better understand QFT let's go through how QFT evolves a 3 qubit state, $\ket{101}$.

Fundamentally, QFT as a Fourier Transform changes the way data is represented. Mathematically for a general state $\ket{j}$ this is $$\hat{F}_{QFT}\ket{j} = \frac{1}{\sqrt{2^n}} \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{n} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{n-1}j_{n} } \right) \ldots \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{1}j_{2}\ldots j_{n} } \right)$$

This can be visualised as a change from the computational basis states {$\ket{0}$,$\ket{1}$} as seen below

After the QFT this information is transformed to be encoded into states on the 'x-y' plane (like {$\ket{-}$,$\ket{+}$}) which rotate along the 'z' axis of the Blochsphere, as seen below

Specifically, information is encoded as rotations from the $\ket{+}$ state. As a brief reminder, such rotations can be achieved with the $R_Z$ gate, which for rotations of $\frac{\pi}{2}$ and $\frac{\pi}{4}$ rotations (with $2\pi$ being a full rotation) are the S and T gates which are demonstrated below

For $\ket{101}$ the state after QFT is \begin{align} \hat{F}_{QFT}\ket{j_{1} j_{2} j_{3}} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{3} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{2}j_{3} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.j_{1}j_{2}j_{3} } \right)\\ \hat{F}_{QFT}\ket{1 0 1} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{(2\pi i) 0.1 } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.01 } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.101 } \right)\\ \hat{F}_{QFT}\ket{1 0 1} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{1}{2} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{1}{4} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{5}{8} } \right) \end{align}

As can be seen in the phases of different qubits, what effectively happens is that as you move from the last qubit to earlier qubits the phase doubles for each qubit, but whole numbers multiples of $2 \pi$ are dropped since they are equivalent to staying in place.

For the 3rd qubit, $\frac{5}{8}$ rotation (of $2 \pi$) occur. Then for the 2nd qubit this doubles to $\frac{5}{4}=1+\frac{1}{4}$ rotations. Finally for the 1st qubit $\frac{5}{2}=2 + \frac{1}{2}$ rotations occur.

To further develop our intuition, let's increment from $\ket{101}$, where the value of the binary representation within the ket vector is $j= 1\times 2^{2} + 0\times 2^{1} + 1\times 2^{0} = 5$, to $\ket{110}$ which has a corresponding value of $j=6$ and perform QFT again. \begin{align} \hat{F}_{QFT}\ket{1 1 0} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{(2\pi i) 0.0 } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.10 } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) 0.110 } \right)\\ \hat{F}_{QFT}\ket{1 1 0} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{0} \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{2}{4} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{6}{8} } \right)\\ \hat{F}_{QFT}\ket{1 1 0} &= \frac{1}{\sqrt{2^3}} \left( \ket{0} + \ket{1}e^{0} \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{1}{2} } \right) \left( \ket{0} + \ket{1}e^{(2\pi i) \frac{3}{4} } \right) \end{align}

So for qubit 3 a further rotation of $\frac{1}{2^3} 2 \pi = \frac{1}{8} 2 \pi $ occurs, for qubit 2 a further rotation of $\frac{1}{2^2} 2 \pi = \frac{1}{4} 2 \pi $ occurs, and qubit 2 a further rotation of $\frac{1}{2^1} 2 \pi = \frac{1}{2} 2 \pi $ occurs.

This results in the following Bloch states:

Inverse QFT circuit

Now that we are familiar with the circuit formulation of the QFT circuit, let's look at the circuit for the inverse QFT. By definition, the inverse QFT reverts the changes applied by QFT circuit through its component gates. This would mean placing inverted (Hermitian conjugate) gates in the opposite order which the gates were applied in QFT, starting from the end with the swap gates and ending with the Hadamard applied on the 1st qubit.

Before we start going through the reverse course of operations, it's worth recognising what these inverse gates are. The swap and Hadamard gates are both Hermitian, that is they are equal to their inverse, as described in the section on quantum gates in the previous chapter.

Conversely, the inverse of a controlled-phase gate which applies $R_k$ onto a target qubit $\ket{b}$ with a control qubit $\ket{a}$ is a similar operation except with $R_{\tilde{k}}$ which is defined as \begin{align} R_{\tilde{k}} &= \begin{bmatrix} 1 & 0 \\ 0 & e^{-2 \pi i /2^{k}} \end{bmatrix} \end{align} This can be seen by applying both gates and returning to the initial state as follows:

First the initial state is $$\ket{a}\ket{b}$$ Then $R_k$ is applied onto $\ket{b}$ with a control qubit $\ket{a}$ $$\ket{a} (R_{k})^{a} \ket{b}$$ Then $R_{\tilde{k}}$ is applied onto the second qubit with the first as the control qubit $$\ket{a} (R_{\tilde{k}})^{a} (R_{k})^{a} \ket{b}$$ $$\ket{a} (R_{\tilde{k}}R_{k})^{a} \ket{b}$$ Above $R_{\tilde{k}} R_{k}=I$, as shown by the following \begin{align} R_{\tilde{k}} R_{k} &= \begin{bmatrix} 1 & 0 \\ 0 & e^{-2 \pi i /2^{k}} \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & e^{2 \pi i /2^{k}} \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 \\ 0 & e^{-2 \pi i /2^{k}} e^{2 \pi i /2^{k}} \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\\ &= I \end{align} Substituting this back into the composite state leads to $$\ket{a} (I)^{a} \ket{b}$$ Which is simply equivalent to returning to the starting state $$\ket{a}\ket{b}$$

Now that we understand inverses, let's start deriving the inverse QFT circuit, beginning with the end result of QFT $$\frac{1}{\sqrt{2^n}} \left( \ket{0} + e^{(2\pi i) 0.j_{n} } \ket{1}\right) \left( \ket{0} + e^{(2\pi i) 0.j_{n-1}j_{n} } \ket{1}\right) \ldots \left( \ket{0} + e^{(2\pi i) 0.j_{1}j_{2}\ldots j_{n} } \ket{1}\right)$$

The last operation carried out during QFT were global swaps, so the inverse of that will be the first operation we carry out now. The inverse of a swap is just another swap, so this similarly reverses the order of qubit states $$\frac{1}{\sqrt{2^{n}}} \left( \ket{0} + e^{2 \pi i 0.j_1 \ldots j_n} \ket{1} \right) \left( \ket{0} + e^{2 \pi i 0.j_2 \ldots j_n} \ket{1} \right) \ldots \left( \ket{0} + e^{2 \pi i 0.j_n} \ket{1} \right) $$

Now let's focus on the $n$th qubit state and apply the inverse of a Hadamard, $H^{\dagger}=H$, another Hadamard \begin{align} H \frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.j_n} \ket{1} \right) &= \frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.j_n} H \ket{1} \right)\\ &= \frac{1}{\sqrt{2}}\left( \frac{1}{\sqrt{2}}(\ket{0}+\ket{1}) + e^{2 \pi i 0.j_n} \frac{1}{\sqrt{2}}(\ket{0}-\ket{1})\right)\\ &= \frac{1}{2}\left( (1+e^{2 \pi i 0.j_n})\ket{0} + (1-e^{2 \pi i 0.j_n})\ket{1} \right) \end{align} We can recover $\ket{j_{n}}$ by recognising that it holds true for either possible value of $j_n$.

If $j_n=0$, so $\ket{j_n}=\ket{0}$, then \begin{align} H \frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.0} \ket{1} \right) &= \frac{1}{2}\left( (1+e^{0})\ket{0} + (1-e^{0})\ket{1} \right)\\ =& \frac{1}{2}\left( (1+1)\ket{0} + (1-1)\ket{1} \right)\\ =& \frac{1}{2}\left( 2\ket{0} \right)\\ =& \ket{0} \end{align} If instead $j_n=1$, so $\ket{j_n}=\ket{1}$, then \begin{align} H \frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.1} \ket{1} \right) &= \frac{1}{2}\left( (1+e^{\pi i})\ket{0} + (1-e^{\pi i})\ket{1} \right)\\ =& \frac{1}{2}\left( (1-1)\ket{0} + (1-(-1))\ket{1} \right)\\ =& \frac{1}{2}\left( 2\ket{1} \right)\\ =& \ket{1} \end{align} Hence the $n$th qubit is $\ket{j_n}$

Now let's include the $(n-1)$th qubit together with the $n$th qubit, which together are $$\frac{1}{\sqrt{2}} \left( \ket{0} + e^{(2\pi i) 0.j_{n-1}j_{n}}\ket{1}\right) \ket{j_n}$$

Applying $R_{\tilde{2}}$ with the $n$th qubit ($\ket{j_n}$) as the control qubit and the $(n-1)$th qubit as the target qubit results in $$\frac{1}{\sqrt{2}} \left( \ket{0} \ket{j_n}+ e^{(2\pi i) 0.j_{n-1}j_{n}} R_{\tilde{2}}\ket{1} \ket{j_n} \right) $$ This is a result of recognising that $R_{\tilde{p}}\ket{0}=\ket{0}$, where because the target qubit is $\ket{0}$ the application of $R_{\tilde{p}}$ does not accumulate phase.

Just like with the Hadamard gate we can go through what happens for the two possible values of $\ket{j_n}$. If $\ket{j_n}=\ket{0}$, then \begin{align} & \frac{1}{\sqrt{2}} \left( \ket{0} \ket{0} + e^{(2\pi i) 0.j_{n-1}0} \ket{1} \ket{0} \right)\\ &= \frac{1}{\sqrt{2}} \left( \ket{0} + e^{(2\pi i) 0.j_{n-1}} \ket{1} \right) \ket{0} \end{align} If instead $\ket{j_n}=\ket{1}$ \begin{align} & \frac{1}{\sqrt{2}} \left( \ket{0} \ket{1} + e^{(2\pi i) 0.j_{n-1} 1 } e^{-(2\pi i) 0.01 }\ket{1} \ket{1} \right)\\ &= \frac{1}{\sqrt{2}} \left( \ket{0} + e^{(2\pi i) ( 0.j_{n-1}1 - 0.01) } \ket{1} \right) \ket{1}\\ &= \frac{1}{\sqrt{2}} \left( \ket{0} + e^{(2\pi i) ( 0.j_{n-1}) } \ket{1} \right) \ket{1}\\ \end{align} Which together leaves us with $$\frac{1}{\sqrt{2}} \left( \ket{0} + e^{(2\pi i) ( 0.j_{n-1}) } \ket{1} \right) \ket{j_n}$$

Applying the Hadamard gate to $(n-1)$th qubit results in same operation as for the $n$th qubit, leading to $$\ket{j_{n-1}}\ket{j_n}$$

Now let's repeat the same series of actions with the $(n-2)th$ qubit to get a feel for how this scales up as the number of controlled-phase gates increases. We start with $$\frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.j_{n-2}j_{n-1}j_n } \ket{1} \right) \ket{j_{n-1}}\ket{j_n} $$

Then we apply $R_{\tilde{3}}$ with the $n$th qubit ($\ket{j_n}$) as the control qubit and the $(n-2)$th qubit as the target qubit, which similar to before effectively removes a factor of $e^{2 \pi i \frac{j_n}{2^3}}=e^{2 \pi i 0.00j_n}$, leading to $$\frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.j_{n-2}j_{n-1} } \ket{1} \right) \ket{j_{n-1}}\ket{j_n} $$

Then we apply $R_{\tilde{2}}$ with the $n$th qubit ($\ket{j_{n-1}}$) as the control qubit and the $(n-2)$th qubit as the target qubit, which removes a factor of $e^{2 \pi i \frac{j_{n-1}}{2^2}}=e^{2 \pi i 0.0j_{n-1}}$, resulting in $$\frac{1}{\sqrt{2}}\left( \ket{0} + e^{2 \pi i 0.j_{n-2} } \ket{1} \right) \ket{j_{n-1}}\ket{j_n} $$

Again like before, the Hadamard is applied to the $(n-2)$th qubit transforming it into $$\ket{j_{n-2}}\ket{j_{n-1}}\ket{j_n}$$

These steps continue all the way until the 1st qubit, decrementing down the number of qubits. To summarise, for the $p$th qubit controlled-phase operations ranging from $R_{\tilde{n+1-p}}$ to $R_{\tilde{2}}$ are applied, with a control state $\ket{j_l}$ for applying a gate $R_{\tilde{l-p}}$ to qubit $p$, followed by a Hadamard gate which transforms into $\ket{j_p}$.

The resulting inverse QFT circuit, excluding the global swaps at the start, is shown below

A demonstration of the QFT and inverse QFT in Qiskit code is available in the following Jupyter notebook.

Quantum Fourier Transform Applications

This section is currently a work in progress.

QFT is a key tool that can be leveraged by quantum algorithms to potentially speed up calculations compared to classical algorithms, including the famed Shor's algorithm for factoring large numbers. This section will go through how exactly this achieved in multiple steps, split among subsections. First, the use of phase estimation to extract the exponential power within the coefficient of $\ket{1}$ will be worked through. Then phase estimation will be utilised in the quantum solution to the problem of order-finding, which itself will be extended towards the factoring of large numbers. With all of these pieces together, Shor's Algorithm will be demonstrated, and the risk to the widely used RSA encryption standard will be elaborated on. Finally, the general applications of QFT will be explained.

Phase Estimation

Let's look at an application for QFT, specifically inverse QFT, as a process used in larger quantum algorithms: 'phase estimation'.

Phase estimation approximates an exponential factor, or 'phase', $\varphi$ within the eigenvalue $e^{2\pi i \varphi}$ of a given unitary operator $U$ when it acts on the eigenvector $\ket{u}$

As part of conditions for performing phase estimation, we assume that we have 'black boxes'/'oracles' which are capable of preparing the $\ket{u}$ state and performing controlled-$U^{2^{j}}$ operations.

The purpose of phase estimation then, is to discover the phase $\phi$ hidden inside the 'black box' of the unitary operator $U$.

For the controlled-$U^{2^{j}}$ gate, when the control qubit state is $\ket{1}$ and the target qubits are $\ket{u}$, the control qubit state is multiplied by the eigenvalue of $e^{2\pi i 2^{j}}$, as demonstrated below with a control qubit in the state $\ket{+}$,

Similar to the case of the two qubit Deutsch's algorithm, the state which ends up being modified by a controlled-unitary operation is the control qubit state. This mechanism is called 'phase-kickback', and introduces a factor of $e^{2\pi i 2^{j}}$ into the control qubit in this case.

We can make use of the phase-kickback mechanism of controlled-$U^{2^{j}}$ gates for different values of $j$ to effectively extract phase data from the eigenstate of $\ket{u}$, as seen in the circuit below where the $n$th control qubit out of $t$ controls for gate $U^{2^{t-n}}$:

This results in the target qubits preserving the $\ket{u}$ state, and the control qubits becoming $$\frac{1}{\sqrt{2^t}} (\ket{0}+ e^{2\pi i 2^{t-1}\varphi}) (\ket{0}+ e^{2\pi i 2^{t-2}\varphi}) \ldots (\ket{0}+ e^{2\pi i 2^{0}\varphi}) $$

We can estimate $\varphi$ with $t$ qubits as $\tilde{\varphi} =0.\varphi_1 \varphi_2 \ldots \varphi_t$, and substitute it into the above state, which results in $$\frac{1}{\sqrt{2^t}} (\ket{0}+ e^{2\pi i 0.\varphi_t}) (\ket{0}+ e^{2\pi i 0.\varphi_{t-1}\varphi_t }) \ldots (\ket{0}+ e^{2\pi i 0.\varphi_1 \varphi_2 \ldots \varphi_t }) $$

Note, that like with the previous subsection on QFT, we have removed factors of $e^{2\pi i}$ since that just results in a multiplication by $1$.

This state is effectively just $\hat{F}_{QFT}\ket{\tilde{\varphi}}$ (as described at the end of the QFT subsubsection) so by performing the inverse QFT $\hat{F^{\dagger}}_{QFT}$ operation we can achieve $\tilde{\varphi}$. A simplified complete phase estimation circuit, including the measurement of $\tilde{\varphi}$ is shown below

$\tilde{\varphi}$ is an estimate so we need to decide how many $n$ qubits of precision to estimate to. Measurement of quantum states is probabilistic, which adds to the number of qubits required by $\text{log}_{2} \left( 2 + \frac{1}{\epsilon} \right)$ to measure the correct phase with a probability of at least $1-\epsilon$, resulting in $$t=n + \text{log}_{2} \left( 2 + \frac{1}{\epsilon} \right)$$

A detailed derivation of the bounding of probabilities which leads to this value can be found in subsection 5.2.1 of 'Quantum Computation and Quantum Information'.

A demonstration of phase estimation, along with how increasing measurement qubits increases precision and accuracy is shown in the following Jupyter notebook.

Order Finding

A key application of phase estimation is order finding, which is the problem of finding the smallest positive integer power $r$ needed to raise a positive integer $x$ until it equals 1 modulo $N$, where $x$ and $N$ share no common factors greater than 1 (otherwise know as being co-prime). Another were of expressing this is $$x^{r}=1 \ (\mathrm{mod}\ N)= 1+kN$$ Where $k$ is either 0 or an integer. As can be seen above, order finding makes use of modular arithmetic.

Modular arithmetic deals with remainders, so for example $$50 \ (\mathrm{mod}\ 7) = 1 + 49 \ (\mathrm{mod}\ 7) = 1 + 7\times 7 \ (\mathrm{mod}\ 7) = 1 $$

Another way of expressing modular arithmetic is as follows $$a + kN \ (\mathrm{mod}\ N) = a$$

An example of order finding for $x=3$ and $N=10$ is \begin{align} 3^{1} \ (\mathrm{mod}\ 10) &= 3\\ 3^{2} \ (\mathrm{mod}\ 10) &= 9\\ 3^{3} \ (\mathrm{mod}\ 10) &= 27 \ (\mathrm{mod}\ 10) = 7\\ 3^{4} \ (\mathrm{mod}\ 10) &= 81 \ (\mathrm{mod}\ 10) = 1\\ \end{align} So here for $x=3$ and $N=10$ the order is $r=4$.

An interesting property for modular exponential is that it is cyclic with a period of the order $r$ $$x^{r+s}=r^{s}\ (\mathrm{mod}\ N)$$ So for example $$3^7\ (\mathrm{mod}\ 10) = 2187 \ (\mathrm{mod}\ 10) = 7 \ (\mathrm{mod}\ 10) = 3^{3} \ (\mathrm{mod}\ 10) $$ This periodic behaviour is also represented in the following plot for $3^{k} \ (\mathrm{mod}\ 10)$

As explored further in following subsections, order-finding is a key problem demonstrating an advantage for quantum algorithm through the application of QFT, which can be applied to the problem of factoring and thus Shor's algorithm, which is a risk to RSA encryption. The periodic nature we see here, also hints at a class of general applications of QFT which can be accelerated by quantum algorithms utilising QFT.

Now that we've ascertained the importance of Quantum order-finding, let's explore how it works, starting with the following oracle $U$ and its eigenvector $\ket{u_s}$ \begin{align} U\ket{y} &= \ket{ xy \ (\mathrm{mod}\ N) }\\ \ket{u_s} &= \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{sk}{r}} \ket{x^{k} \ (\mathrm{mod}\ N) } \end{align}

The eigenvalue can derived by applying $U$ to $\ket{u_s}$ as follows \begin{align} U \ket{u_s} &= \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{sk}{r}} U \ket{x^{k} \ (\mathrm{mod}\ N) } \\ &= \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{2\pi i \frac{s}{r}} e^{-2\pi i \frac{s(k+1)}{r}} \ket{x \times x^{k} \ (\mathrm{mod}\ N) } \\ &= e^{2\pi i \frac{s}{r}} \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{s(k+1)}{r}} \ket{x^{k+1} \ (\mathrm{mod}\ N) } \end{align} Which upon substitution of $k^\prime=k+1$ becomes $$ U \ket{u_s} = e^{2\pi i \frac{s}{r}} \frac{1}{\sqrt{r}} \sum^{r}_{k=1} e^{-2\pi i \frac{sk^{\prime}}{r}} \ket{x^{k^{\prime}} \ (\mathrm{mod}\ N) } $$

The one thing stopping us from directly extracting an eigenvalue is that the limits of the sum have changed, though this can be resolved by rearranging the last term of the sum for $k^{\prime}=r$, $e^{-2\pi i \frac{sr}{r}} \ket{x^{k^{r}} \ (\mathrm{mod}\ N) }$ to correspond to $k^{\prime}=0$. For the exponential factor, we can recognise that $e^{-2\pi i \frac{sr}{r}} = e^{-2\pi i s} = e^{0}$ since $s$ is an integer, so it can be rearranged as $e^{-2\pi i \frac{s \times 0 }{r}}$ to correspond to $k^{\prime}=0$. For the ket vector argument, we instead utilise the periodic nature of the modular power to get $\ket{x^{k^{r}} \ (\mathrm{mod}\ N) } = \ket{x^{k^{0+r}} \ (\mathrm{mod}\ N) } = \ket{x^{k^{0}} \ (\mathrm{mod}\ N) }$. Together, this means that the limits of the sum can be changed so \begin{align} U \ket{u_s} &= e^{2\pi i \frac{s}{r}} \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{sk^{\prime}}{r}} \ket{x^{k^{\prime}} \ (\mathrm{mod}\ N) } \\ U \ket{u_s} &= e^{2\pi i \frac{s}{r}} \ket{u_s} \end{align} So the eigenvalue is $e^{2\pi i \frac{s}{r}}$, which contains the order $r$ in the power. Phase estimation can be then be used to extract $\frac{s}{r}$. This is where quantum order-finding utilises QFT.

We can then extract order $r$ from the denominator of the fraction through classical post-processing it into a continued fraction.

The idea of the continued fractions algorithm is to express a fraction in terms of a finite list of positive integers $[a_0, a_1 , \ldots, a_N]$, with the exception of $a_0$ which can also be 0, as $$[a_0, a_1 , \ldots, a_N] = a_0 + \frac{1}{a_1 + \frac{1}{a_2 + \frac{1}{\ldots + \frac{1}{a_N}}}} $$

And then to utilise this sequence of numbers to iteratively determine a fraction $\frac{p_{N}}{q_{N}}$ with the following initial values \begin{align} p_0 &= a_0, &\quad q_0 &= 1 \\ p_1 &= 1 + a_{0} a_{1}, &\quad q_1 &= a_1 \end{align} and the following iterative rules for $2 \leq n \leq N $ \begin{align} p_n &= a_{n} p_{n-1} + p_{n-2}\\ q_1 &= a_{n} q_{n-1} + q_{n-2} \end{align}

To grasp how the continued fractions algorithm functions, let's go through step-by-step with $\frac{31}{13}$ (which is $2.384615384615384615\ldots = 2.\overline{384615}$ ). Here we will utilise $b_{n}$ to define the remainder at each step for ease of understanding, while also noting and updating the value of $p_n$ and $q_n$.

Let's start with the first step for $n=0$ \begin{align} \frac{31}{13} &= 2 + \frac{5}{13} \\ a_{0} &= 2 \\ b_{0} &= \frac{5}{13} \\ p_0 &= a_0=2\\ q_0 &= 1 \end{align} After this, we will need to invert the remainder $b_{0}$ so the numerator becomes $1$ as follows $$b_{0} = \frac{1}{ \frac{13}{5}}$$ We can then split the denominator of $b_{0}$ $$b_{0} = \frac{1}{ 2+ \frac{3}{5}}$$ $b_0$ can then be substituted back in $$ \frac{31}{13} = 2 + \frac{1}{ 2+ \frac{3}{5}} $$ The values of $p_0$ and $q_0$ are \begin{align} p_0 &= a_0= 2\\ q_0 &= 1 \end{align}

For $n=1$ we will repeat this whole process, except this time using $p_1= 1+a_{0} a_{1}$ \begin{align} \frac{31}{13} &= 2 + \frac{1}{2+\frac{3}{5}}\\ a_{1} &= 2 \\ b_{1} &= \frac{3}{5} \\ &= \frac{1}{ \frac{5}{3}} \\ &= \frac{1}{ 1 + \frac{2}{3}}\\ p_1 &= 1+ a_{0} a_{1} = 1 + 2 \times 2 = 5\\ q_1 &= a_{1} = 2 \end{align}

For $n=2$, we then have to remember to use the iterative rules of $p_n$ and $q_n$ for $2 \leq n \leq N $ \begin{align} \frac{31}{13} &= 2 + \frac{1}{2+\frac{1}{ 1 + \frac{2}{3}}}\\ a_{2} &= 1 \\ b_{2} &= \frac{2}{3} \\ &= \frac{1}{ \frac{3}{2}} \\ &= \frac{1}{ 1 + \frac{1}{2}}\\ p_2 &= a_{2} p_{1} + p_{0} = 1 \times 5 + 2 = 7\\ q_2 &= a_{2} q_{1} + q_{0} = 1 \times 2 + 1 = 3 \end{align}

For $n=3$, \begin{align} \frac{31}{13} &= 2 + \frac{1}{2+\frac{1}{ 1 + \frac{1}{ 1 + \frac{1}{2}}}}\\ a_{3} &= 1 \\ b_{3} &= \frac{1}{2} \\ p_3 &= a_{3} p_{2} + p_{1} = 1 \times 7 + 5 = 12\\ q_3 &= a_{3} q_{2} + q_{1} = 1 \times 3 + 2 = 5 \end{align} We find that there is no need to invert $b_{3}$ to get a numerator of 1, so there is also no need to split the denominator in the next step of $n=4$, making $n=4$ the final step.

For $n=4$, the final step: \begin{align} \frac{31}{13} &= 2 + \frac{1}{2+\frac{1}{ 1 + \frac{1}{ 1 + \frac{1}{2}}}}\\ a_{4} &= 2 \\ b_{4} &= 0 \\ p_4 &= a_{4} p_{3} + p_{2} = 2 \times 12 + 7 = 31\\ q_4 &= a_{4} q_{3} + q_{2} = 2 \times 5 + 3 = 13 \end{align} So finally we have the sequence $[2, 2, 1, 1, 2]$ and the determined fraction $\frac{p_{N}}{q_{N}}= \frac{p_{4}}{q_{4}} = \frac{31}{13}$.

Let's get back to looking at the operator $U$ and eigenvector $\ket{u_s}$.

$U$ gates can be used to evolve a state $\ket{y \ (\mathrm{mod}\ N)}$ to $\ket{x^{z} y \ (\mathrm{mod}\ N) }$ for an arbitrary $z=z_{n} 2^{n-1} + z_{t-1} 2^{n-2} + \ldots z_{1} 2^{0}$. This can be achieved by applying $n$ gates: $$U^{z} = U^{z_{n} 2^{n-1}} U^{ z_{t-1} 2^{n-2}} \ldots U^{z_{1} 2^{0}} $$

We can't be assumed to know $\ket{u_s}$ since we need $r$ to define it $$\ket{u_s} = \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{sk}{r}} \ket{x^{k} \ (\mathrm{mod}\ N) } $$ Thankfully we can use $\ket{1}$ instead, because $\ket{u_s}$ can be summed to equal $\ket{1}$ \begin{align} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} \ket{u_s} &= \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} \frac{1}{\sqrt{r}} \sum^{r-1}_{k=0} e^{-2\pi i \frac{sk}{r}} \ket{x^{k} \ (\mathrm{mod}\ N) }\\ &= \frac{1}{r} \sum^{r-1}_{k=0} \left(\sum^{r-1}_{s=0} e^{-2\pi i \frac{sk}{r}} \right) \ket{x^{k} \ (\mathrm{mod}\ N) } \end{align}

Above, $\sum^{r-1}_{s=0} e^{-2\pi i \frac{sk}{r}}=r \delta_{k,0}$, where $\delta_{k,0}$ is the Kronecker delta which is $1$ only if $k=0$ (and $0$ otherwise). This conclusion can be arrived at by recognising that if $k=0$ then the exponential is $e^{-2\pi i \frac{s\times 0}{r}}=e^{0}=1$ so the sum is $\sum^{r-1}_{s=0} 1 = r$, but if $k\neq0$ then the sum is $0$. The latter is a little more complicated to derive, but makes use of the sum of geometric series, as described briefly before.

We can define a variable $c^{s}=e^{-2\pi i \frac{sk}{r}} = (e^{-2\pi i \frac{k}{r}})^{s}$, and the sum is now $\sum^{r-1}_{s=0} c^{s}$, a sum of a geometric series. By subtracting the product of $c$ and $\sum^{r-1}_{s=0} c^{s}$ from $\sum^{r-1}_{s=0} c^{s}$, there are only two terms which don't cancel out, namely $-c^{r-1} \times c = -c^{r}$ from the first sum and $c^{0}=1$ from the second sum, so we have \begin{align} \sum^{r-1}_{s=0} c^{s} - c \sum^{r-1}_{s=0} c^{s} &= 1 - c^{r} \\ (1-c) \sum^{r-1}_{s=0} c^{s} &= 1 - c^{r} \\ \sum^{r-1}_{s=0} c^{s} &= \frac{1-c^{r}}{1-c} \end{align} Here, $c^{r}=e^{-2\pi i \frac{r k}{r}} = e^{-2\pi i k} = e^{0} = 1$ (because $k$ is an integer) so the numerator $ 1 - c^{r} = 0$, and therefore the sum itself is zero. On the other hand, the denominator does not go to zero, because $k\neq 0$ so $ c = e^{-2\pi i \frac{k}{r}} \neq e^{0}=1$ and hence $1-c \neq 0$. This all means that for $k\neq 0$ the sum is $0$.

Going back to the sum of $\ket{u_s}$ by substituting in the Kronecker delta, we have \begin{align} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} \ket{u_s} &= \frac{1}{r} \sum^{r-1}_{k=0} r \delta_{k,0} \ket{x^{k} \ (\mathrm{mod}\ N)} \\ &= \sum^{r-1}_{k=0} \delta_{k,0} \ket{x^{k} \ (\mathrm{mod}\ N)} \\ &= \ket{x^{0} \ (\mathrm{mod}\ N)}\\ &= \ket{1 \ (\mathrm{mod}\ N)}\\ &= \ket{1} \end{align} Note here, that $\ket{1 \ (\mathrm{mod}\ N)}$ is a multi-qubit state with all qubits initialised to $1$. As we will see soon, we can rearrange quantum states to substitute in our sum of $\ket{u_s}$.

Let's remember that we're utilising the phase estimation algorithm for quantum order-finding, so there will be some similar steps.

We begin with a $n$ qubit list control qubits initialised to $0$, and a $L$ qubit list of target qubits initialised to $1$. After applying Hadamard gates to each of the control qubits we arrive at $$\frac{1}{\sqrt{2^{n}}} \left( \sum^{2^{n}-1}_{k=0} \ket{k} \right) \ket{1} $$ The application of $U^{k}$ gates controlled by $k$ then leads to $$\frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \left(\ket{k} U^{k} \ket{1} \right) = \frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \ket{k} \ket{x^{k} \ (\mathrm{mod}\ N) } $$

We can then rearrange the left side of the above equation, and substitute in our sum of $\ket{u_s}$ to achieve \begin{align} \frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \ket{k} U^{k} \left( \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} \ket{u_s} \right) \\ = \frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \ket{k} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} U^{k} \ket{u_s} \end{align}

Applying $U^{k}$ to the eigenvector $\ket{u_s}$ results in the eigenvalue of $U$ with $\ket{u_s}$, $e^{2\pi i \frac{s}{r}}$ , being applied $k$ times, leading to a multiplication by $e^{2\pi i \frac{s}{r}\times k } = e^{2\pi i k \frac{s}{r}}$. This means that the quantum state is \begin{align} \frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \ket{k} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} e^{2\pi i k \frac{s}{r}} \ket{u_s} \\ = \frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} e^{2\pi i k \frac{s}{r}} \ket{k} \ket{u_s} \end{align}

Like with phase estimation algorithm, we can perform an inverse QFT to extract an estimation of the phase factor $\frac{s}{r}$, for which we can apply the continued fraction algorithm to estimate the order $r$.

The precision, and thus number of bits used, is defined by the number of bits taken up by $N$ (which is co-prime to $x$), $L=\log_{2}N$. $L$ feeds into the number of control bits $t$ needed. Without going too far into details (which are discussed in subsection 5.3.1 of 'Quantum Computation and Quantum Information'), the general idea is that $2L + 1$ bits are needed for the continued fractions algorithm to approximate $\frac{s}{r}$, and this slots into the precision for phase estimation for an error rate of $\epsilon$ to result in $$t=2L+ 1 + \text{log}_{2} \left( 2 + \frac{1}{\epsilon} \right) =2\log_2 (N)+ 1 + \text{log}_{2} \left( 2 + \frac{1}{\epsilon} \right) $$

The following simplified quantum circuit represents the quantum order-finding process we discussed in this subsection:

To summarise, the following steps are taken during quantum order-finding

1. Initialise $L$ ($\geq log_{2} N$) qubits in the first register of control qubits to $\ket{0}$ and $t$ ($=2L+ 1 + \text{log}_{2} \left( 2 + \frac{1}{\epsilon} \right)$) qubits in the second register of target qubits $\ket{0}$, and pick an integer $x$ co-prime to $N$
2. Create a superposition of states in the first register through Hadamard gates, $\frac{1}{\sqrt{2^{n}}} \left( \sum^{2^{n}-1}_{k=0} \ket{k} \right) \ket{1} $
3. Apply control-$U^{k}$ gates to the second register using the first register as control qubits, to prepare $\frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \ket{k} \ket{x^{k} \ (\mathrm{mod}\ N) }$, which can also be presented as $\frac{1}{\sqrt{2^{n}}} \sum^{2^{n}-1}_{k=0} \frac{1}{\sqrt{r}} \sum^{r-1}_{s=0} e^{2\pi i k \frac{s}{r}} \ket{k} \ket{u_s}$
4. Apply the inverse Fourier Transform, and measure the first register to record an estimate of the phase $\frac{s}{r}$
5. Apply the continued fractions algorithm to extract order $r$ from the denominator of the phase estimate

Factoring

Factoring a large number $N$ can make use of order finding when $r$ is an even number. This can be done because $x^{r+s}=r^{s}\ (\mathrm{mod}\ N)$ (for $x$ co-prime to $N$) means that $x^{r} = 1 + kN$ (for an integer k), and through using the relation $y^{2} = (y-1)(y+1)$ for $y=x^{r/2}$, which together combines to form $$kN = (x^{r/2} - 1) (x^{r/2} + 1)$$ This all means that for non-trivial values of $y=x^{r/2}$, meaning $y\neq N-1 \ (\mathrm{mod}\ N) = -1 \ (\mathrm{mod}\ N))$, then both $(x^{r/2} - 1)$ and $(x^{r/2} + 1)$ will be factors of $N$ which are not $N$ or $1$. For a composite number $N$ composed of $m$ distinct factors, the probability that $r$ is even and that $x^{r/2} \neq -1 \ (\mathrm{mod}\ N))$ is at least $1 + \frac{1}{2^{m}}$ (with the derivation in section A4.3 in 'Quantum Computation and Quantum Information' for further reading), so it is likely that this is indeed the case. Assuming this condition is met, then the next step from there is to find the greatest common divider between these factors $(x^{r/2} - 1)$ and $(x^{r/2} + 1)$ and $N$, such as through Euclid's algorithm as explained below.

Euclid's algorithm works iteratively to find the greatest common divisor (or factor) between two integers $a$ and $b$, otherwise expressed as gcd$(a,b)$. The starting two steps are \begin{align} a &= b k_{1} + r_{1} \\ b &= r_{1} k_{2} + r_{2} \end{align} Where $k_{n}$ is the quotient, or integer multiplier of $b$, and $r_{n}$ is the remainder for a given step $n$. For further steps $n \geq 3$, the iteration formula instead becomes $$r_{n-2} = r_{n-1}k_{n} + r_{n}$$ Which continues until the condition $r_{n}=0$ is reached, whereupon $r_{n-1}$ is determined to be the greatest common divisor of $a$ and $b$. Note here, that if gcd$(a,b)=1$ then $a$ and $b$ are co-prime, so this is a useful test for such.

To grasp Euclid's algorithm, it's worth going through an example such as to find the greatest common divisor for $a=1068$ and $b=464$.

For $n=1$, we use $a = b k_{1} + r_{1}$ \begin{align} a &= b k_{1}+ r_{1}\\ 1068 &= 464 k_{1} + r_{1}\\ 1068 &= 464 \times 2 + 140 \\ k_{1} &= 2, \quad r_{1} = 140 \end{align} For $n=2$, we use $b = r_{1} k_{2} + r_{2} $ \begin{align} b &= r_{1} k_{2} + r_{2} \\ 464 &= 140 k_{2} + r_{2} \\ 464 &= 140 \times 3 + 44 \\ k_2 &= 3, \quad r_2 = 44 \end{align} For $n=3$, we introduce the iterative rule $r_{n-2} = r_{n-1}k_{n} + r_{n}$ \begin{align} r_{1} &= r_{2} k_{3} + r_{3}\\ 140 &= 44 k_{3} + r_{3} \\ 140 &= 44 \times 3 + 8 \\ k_{3} &= 3, \quad r_{3} = 8 \end{align} For $n=4$, \begin{align} r_{2} &= r_{3} k_{4} + r_{4}\\ 44 &= 8 k_{4} + r_{4} \\ 44 &= 8 \times 5 + 4 \\ k_{4} &= 5, \quad r_{4} = 4 \end{align} For $n=5$, \begin{align} r_{3} &= r_{4} k_{5} + r_{5}\\ 8 &= 4 k_{5} + r_{5} \\ 8 &= 4 \times 2 + 0 \\ k_{5} &= 2, \quad r_{5} = 0 \end{align} Here $r_{5}=0$ so the algorithm ends and we take the penultimate remainder $r_{4}=4$ as the greatest common divisor of $1068$ and $464$.

For factorisation of a number $N$, before implementing order-finding, realistically it is helpful to classically screen for if the factor $2$ when $N$ is even, and whether $N$ is a power of a positive integer $a$ (as $a^{b}$ for $b\geq 2$). The implementation of order finding would also likely have to be repeated multiple times, both to account for failing to meet the conditions ($r$ is even and that $x^{r/2} \neq -1 \ (\mathrm{mod}\ N))$ ) and to capture all the factors of $N$.

Bring everything together, the steps are to repeat the following until all the factors of $N$ are collected

1. Check if $N$ is even, if so return $2$ as a factor
2. Check if $N$ can be expressed as $a^{b}$ for positive integers $a$ and $b$, and $b \geq 2$. If this matches, return $a$ and focus on factoring $a$
3. Randomly choose an integer $x$ which is $1 \lt x \lt N$ and find gcd$(x,N)$. If gcd$(x,N) \geq 1$ then return it as a factor and restart this step, otherwise proceed to the next step with co-prime $x$ and $N$
4. Use quantum order-finding to find the order $r$ in $x^{r}=1 \ (\mathrm{mod}\ N)$
5. If $r$ is even and $x^{r/2} \neq -1 \ (\mathrm{mod}\ N))$, then compute and return gcd$(x^{r/2} - 1,N)$ and gcd$(x^{r/2} + 1,N)$. Otherwise, this algorithm fails and return to step 3

Shor's Algorithm

This subsection is currently a work in progress.

Risk to RSA

This subsection is currently a work in progress.

General Applications of QFT

This subsection is currently a work in progress.

Grover's Quantum Search Algorithm

This section is currently a work in progress.

Further Comprehension

To internalise the concepts presented here, it helps to write out derivations with pen-and-paper, and to consult multiple sources to see how each addresses the same concepts from slightly different points of view.

The YouTube playlist of the Qiskit 2021 Summer School course is freely available, along with lecture notes, and definitely worth checking out if you prefer learning in a lecture format. Specifically, lectures 2.1 and 2.2 cover the Deutsch-Jozsa and Grover's search algorithms.

Another recommended resource which I've consulted when writing this, is the seminal textbook 'Quantum Computation and Quantum Information' (10th Anniversary Edition) by Michael Nielsen and Isaac Chuang. While this is unfortunately not freely available, it is worth going through this or another textbook or similar text resource (like Ronald de Wolf's Quantum Computing Lecture Notes) and solve questions with pen-and-paper to build a more rigorous understanding after going through this tutorial. In the case of 'Quantum Computation and Quantum Information', I would recommend going through chapters 1, 5 and 6.