#  Shor’s Factoring Algorithm and
Modular Exponentiation Operators

Robert L Singleton Jr† Quantum Division
SavantX Research Center
Santa Fe, New Mexico, USA

(8/7/2023)

###### Abstract

This manuscript provides a pedagogical presentation of Shor’s factoring
algorithm, which is a quantum algorithm for factoring very large numbers (of
order of hundreds to thousands of bits) in polynomial time. In contrast, all
known classical algorithms for the factoring problem take an exponential time
to factor such large numbers. Shor’s algorithm therefore has profound
implication for public-key encryption such as RSA and Diffie-Hellman key
exchange. In these notes, we assume no prior knowledge of Shor’s algorithm
beyond a basic familiarity with the circuit model of quantum computing. Shor’s
algorithm contains a number of moving parts, and can be rather daunting at
first. The literature is replete with derivations and expositions of Shor’s
algorithm, but most of them seem to be lacking in essential details, and none
of them provide a pedagogical presentation. They require a thicket of
appendices and assume a knowledge of quantum algorithms and classical
mathematics with which the reader might not be familiar. We therefore start
with first principle derivations of the quantum Fourier transform (QFT) and
quantum phase estimation (QPE), which are the essential building blocks of
Shor’s algorithm. We then go on to develop the theory of modular
exponentiation (ME) operators, one of the fundamental components of Shor’s
algorithm, and the place where most of the quantum resources are deployed. We
also delve into the number theory that establishes the link between
factorization and the period of the modular exponential function. We then
apply the QPE algorithm to obtain Shor’s factoring algorithm. We also discuss
the post-quantum processing and the method of continued fractions, which is
used to extract the exact period of the modular exponential function from the
approximately measured phase angles of the ME operator. The manuscript then
moves on to a series of examples. We first verify the formalism by factoring
$N=15$, the smallest number accessible to Shor’s algorithm. We then proceed to
factor larger integers, developing a systematic procedure that will find the
ME operators for any semi-prime $N=p\times q$ (where $q$ and $p$ are prime).
Finally, we factor the composite numbers $N=21,33,35,143,247$ using the Qiskit
simulator. It is observed that the ME operators are somewhat forgiving, and
truncated approximate forms are able to extract factors just as well as the
exact operators. This is because the method of continued fractions only
requires an approximate phase value for its input, which suggests that
implementing Shor’s algorithm might not be as difficult as first suspected.

$\dagger$ corresponding email<EMAIL_ADDRESS>or
<EMAIL_ADDRESS>
###### Contents

  1. I Introduction
  2. II The Quantum Fourier Transform
    1. II.1 General Definitions
    2. II.2 Qubit Ordering and the QFT Circuit
      1. II.2.1 Convention 2: Standard Physics/Mathematics
      2. II.2.2 Convention 1: OpenQASM/Qiskit
  3. III Quantum Phase Estimation
    1. III.1 Convention 2: Standard Physics/Mathematics
    2. III.2 Convention 1: OpenQASM/Qiskit
    3. III.3 General Phase Angles
    4. III.4 Generalized Input States
  4. IV Continued Fractions
  5. V Factoring with Shor’s Algorithm
    1. V.1 Basic Observation
    2. V.2 Period Finding and Factorization
    3. V.3 The Factorization Circuit
    4. V.4 Extracting the Exact Period from the Measured Phase
  6. VI Verifying the Theory: Factoring N = 15
    1. VI.1 Modular Exponentiation Operators
    2. VI.2 The Factorization Circuit
    3. VI.3 The Spectrum for a = 4
    4. VI.4 The Spectrum for a = 8
    5. VI.5 Analysis of the Phase Histogram
  7. VII Further Examples
    1. VII.1 Factoring Larger Numbers: N = 21 = 3 $\times$ 7, a = 2, r = 6
    2. VII.2 More Factoring
      1. VII.2.1 $N=35=5\times 7$, $a=4$, $r=6$
      2. VII.2.2 $N=33=3\times 11$, $a=7$, $r=10$
      3. VII.2.3 $N=143=13\times 11$, $a=5$, $r=20$
      4. VII.2.4 $N=247=13\times 19$, $a=2$, $r=36$
  8. VIII Conclusions and Outlook
  9. A Modular Exponentiation Operators
    1. A.1 Composite ME operators for N = 143, a = 5, r = 20, m = 8, u_ver=1
    2. A.2 Composite ME operators for N = 247, a = 2, r = 36, m = 9, u_ver=1

##  I Introduction

In this manuscript, we present a pedagogical construction of Shor’s factoring
algorithm [1], which can factor exponentially large integers in polynomial
time. All known classical algorithms for factorization require an exponential
time because they work by brute force, essentially testing all (or most)
numbers less than the number being factored. In contrast, Shor’s algorithm
exploits the massive parallelism inherent in quantum mechanics, so that all
numbers can be tested simultaneously rather than sequentially, thereby
providing for a polynomial factorization process. Since Shor’s algorithm can
factor massively large numbers very quickly, it has major implications for the
security of encryption standards such as RSA [2] and Diffie-Hellman [3, 4] key
exchange, rendering these encryption methods severely compromised. As most
internet communication is based on such public key encryption schemes, Shor’s
algorithm has profound implications for digital security.

Shor’s algorithm rests upon two fundamental quantum algorithms, the quantum
Fourier transform (QFT) and quantum phase estimation (QPE). The QFT, as the
name suggests, implements the discrete Fourier transform on a gated quantum
computer. Like the classical Fourier transform, it extracts frequency signals
from an input source, except that the QFT works by manipulating quantum bits
or qubits (two state quantum systems) on a gated quantum computer. The QPE
algorithm, in contrast, finds the complex phases or the Eigenvalues of an
arbitrary unitary linear operator. Shor’s algorithm elegantly combines the QFT
and QPE to construct a powerful quantum algorithm for factoring very large
integers. More precisely, by employing a specific and well chosen unitary
operator called the modular exponentiation (ME) operator, quantum parallelism
allows the QPE to extract the factors of exponentially large numbers in
polynomial time. One might say that the QPE algorithm is the workhorse of
Shor’s algorithm [5], and we shall spend most of our time on the associated ME
operators. The mathematics behind Shor’s algorithm is based on a simple but
profound result from number theory, which maps the factoring problem onto
another mathematical problem that finds the period of the modular exponential
function. The period of this function is directly related to the factors of
the number in question, and the QPE extracts this period using the method of
continued fractions, thereby providing the sought after factors.

Historically, Shor’s algorithm was motivated by Simon’s algorithm [6], but we
shall not discuss these (interesting) details. We assume no prior knowledge of
Shor’s algorithm beyond a basic familiarity with the circuit model of quantum
computing. Shor’s algorithm has a number of moving parts, and it is rather
complex. It contains a pre-processing phase that happens on a classical
computer, then the QPE is conducted on a quantum computer (using the
associated ME operators), and finally there is a post-processing phase that
takes place on a classical computer (employing the method of continued
fractions). The literature is thick with expositions of Shor’s algorithm, but
most of them seem to be lacking in some essential respect, and none of them
provide a satisfying pedagogical presentation. Consequently, this manuscript
is an attempt to derive Shor’s algorithm from first principles in a self-
contained manner assuming minimal familiarity with the requisite quantum
computing machinery. In Sections II and III, we therefore provide complete
derivations the QFT and the QPE algorithms, respectively. As we have
emphasized, these algorithms form the essential building blocks of Shor’s
algorithm. For the classical post-processing stage, one must utilize the
theory of continued fractions, and in Section IV we provide a brief
introduction to the subject, proving a number of fundamental theorems. In
Section V, we are finally ready to address Shor’s algorithm, which involves a
rigorous construction of the appropriate modular exponentiation operator $U$.
Shor’s algorithm then follows by applying the QPE algorithm to this operator.
After this, we move on to discuss the post-quantum processing in more detail,
where we apply the theory of continued fractions from Section IV to extract
the exact period of the modular exponential function $f(x)$ from the
approximately measured phase of the ME operator $U$.

In Section VI we apply the formalism to factor $N=15$, the smallest number
accessible to Shor’s algorithm, and we use this section to develop an all
purpose factoring script. The difficulty in factoring a number with Shor’s
algorithm does not lie in the magnitude of the number itself, but in the size
of the period $r$ of the modular exponential function $f(x)$ [7]. In Section
VII, we apply the formalism to the composites $N=21,33,35,143,247$, which have
periods ranging from $r=2$ to $r=36$. One might think that we have
accomplished nothing, since knowing the exact ME operator is equivalent to
knowing the period $r$ of the function $f(x)$, and Shor’s algorithm would
therefore be unnecessary. However, it turns out that we do not require the
exact ME operators! It is observed that the ME operators are somewhat
forgiving, and truncated approximate forms are able to extract factors just as
well as the exact operators. This is because the method of continued fractions
only requires an approximate phase value for its input, which suggests that
implementing Shor’s algorithm might not be as difficult as first suspected.
Finally, Section VIII provides some conclusions and closing remarks.

##  II The Quantum Fourier Transform

###  II.1 General Definitions

In this section we formulate of the quantum Fourier transform (QFT), where our
primary references are Refs. [8] and [9]. Given an $M$-vector of complex
numbers $\psi=(\psi_{0},\psi_{1},\cdots,\psi_{{\scriptscriptstyle M}-1})$, the
discrete Fourier transform
$\tilde{\psi}=(\tilde{\psi}_{0},\tilde{\psi}_{1},\cdots,\tilde{\psi}_{{\scriptscriptstyle
M}-1})$ is defined by

$\displaystyle\tilde{\psi}_{\ell}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi i\,\ell
k/M}\,\psi_{k}\ ,$ (1)

and the discrete inverse Fourier transform is therefore given by

$\displaystyle\psi_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{\ell=0}^{M-1}e^{-2\pi i\,\ell
k/M}\,\tilde{\psi}_{\ell}\ ,$ (2)

where the indices $\ell,k\in\\{0,1,\cdots,M-1\\}$. We wish to implement the
Fourier transform using an $m$-qubit quantum system, where $M=2^{m}$ is the
number of possible quantum states. The corresponding quantum Fourier transform
will be a linear unitary operator on the $m$-qubit Hilbert space, denoted by
$QFT$, whose action on the computational basis elements reproduces the
classical transform (1),

$\displaystyle QFT\,|\ell\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi i\,\ell
k/M}\,|k\rangle\ .$ (3)

Recall that a linear operator defined only on the basis states is sufficient
to give the operator on any state in the Hilbert space. We can now express the
$QFT$ operator in a very useful basis-dependent form,

$\displaystyle
QFT=QFT\cdot\mathbb{1}=QFT\cdot\underbrace{\sum_{\ell=0}^{M-1}|\ell\rangle\langle\ell|}_{\mathbb{1}}=\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}\sum_{\ell=0}^{M-1}e^{2\pi
i\,k\ell/M}\,|k\rangle\langle\ell|\ ,$ (4)

where we have used the decomposition of unity
$\mathbb{1}=\sum_{\ell=0}^{{\scriptscriptstyle M}-1}|\ell\rangle\langle\ell|$.
Since the quantum Fourier transform is unitary (an easy proof), that is to say
$QFT\cdot QFT^{\dagger}=\mathbb{1}$, then the inverse quantum Fourier
transform is given by

$\displaystyle
QFT^{-1}=QFT^{\dagger}=\frac{1}{\sqrt{M}}\sum_{\ell=0}^{M-1}\sum_{k=0}^{M-1}e^{-2\pi
i\,k\ell/M}\,|\ell\rangle\langle k|\ .$ (5)

Our aim in this Section is to construct a quantum circuit to implement the
$QFT$ operator and its inverse.

###  II.2 Qubit Ordering and the QFT Circuit

Figure 1: An $m$-bit binary integer $k$ can be encoded on a gated quantum
computer in either of the conventions described in the text. Convention 1 is
called the OpenQASM/Qiskit convention, while Convention 2 is called the
Physics/Mathematics convention. Convention 1 labels the top qubit by 0, and
works its way down to the last qubit labeled by $m-1$. Binary integers are
expressed with the standard bit encoding $k=k_{m-1}\cdots k_{1}k_{0}$, which
places the lowest order bit $k_{0}$ at the top of the circuit. Convention 2
numbers the qubits from 1 to $m$ running from top to bottom, and the bit-
ordering of integers is flipped to $k=k_{1}k_{2}\cdots k_{m}$, placing the
lowest order bit $k_{m}$ at the bottom of the circuit.

A quantum circuit has four distinct qubit ordering conventions of which we
must be cognizant, as they are all present in the literature. We can order the
qubits on the quantum circuit in two ways, and we can order the bits of a
binary integer in two ways, thereby giving four possible conventions. We will
be concerned with two such conventions. In computer science, one represents an
$m$-bit integer $k$ by $m$ binary digits $k_{r}\in\\{0,1\\}$ using the
standard notation $k=k_{m-1}\cdots k_{1}k_{0}$. The bit-ordering convention is
that $k_{0}$ is the lowest-order bit, so that the value of the integer is
given by $k=k_{0}\,2^{0}+k_{1}\,2^{1}+\cdots+k_{m-1}\,2^{m-1}$. We must also
label the qubits of the quantum circuit, which gives two more possible
conventions. In OpenQASM/Qiskit [10], the upper qubit of an $m$-qubit circuit
is labeled by 0, working its way down the circuit and ending with qubit $m-1$
at the bottom. For a general $m$-bit binary integer $k=k_{m-1}\cdots
k_{1},k_{0}$, the least significant bit $k_{0}$ is therefore placed at the top
of the circuit, encoded by the computational basis state $|k_{0}\rangle$ of
the upper qubit. The integers $k$ therefore correspond to the computational
basis states $|k\rangle\equiv|k_{m-1}\cdots
k_{0}\rangle\equiv|k_{0}\rangle\otimes|k_{1}\rangle\otimes\cdots\otimes|k_{m-2}\rangle\otimes|k_{m-1}\rangle$.
Note that the qubit ordering is opposite to the bit-string ordering. This will
be one of our primary conventions. Alternatively, we can number the qubits on
the circuit from $1$ to $m$, and express binary numbers by $m$-bit strings of
the form $k=k_{1}k_{2}\cdots k_{m-1}k_{m}$. In this convention, the lowest
order bit $k_{m}$ is placed at the bottom of the circuit, and the value of the
index integer is $k=k_{m}\,2^{0}+k_{m-1}\,2^{1}+\cdots+k_{1}\,2^{m-1}$. The
qubit ordering is the same as the bit-string ordering, and the integers $k$
correspond to the computational basis states $|k\rangle\equiv|k_{1}\cdots
k_{m}\rangle\equiv|k_{1}\rangle\otimes|k_{2}\rangle\otimes\cdots\otimes|k_{m-1}\rangle\otimes|k_{m}\rangle$.
This is the standard physics and mathematics convention. We shall use both
conventions, which are illustrated in Fig. 1. Quantum circuits will be
inverted horizontally between these two conventions, so it is important to
keep track of which convention is in use.

####  II.2.1 Convention 2: Standard Physics/Mathematics

We first work through the details of the QFT for Convention 2, the physics and
mathematics convention. We consider an $m$-qubit system in which the qubits
are ordered from top to bottom, starting with qubit-$1$ in the upper position
of the circuit and qubit-$m$ at the bottom of the circuit. The quantum system
has $M=2^{m}$ distinct states that can be indexed by an integer
$k\in\\{0,1,\cdots,M-1\\}$. We can express this integer in the binary form,

$\displaystyle k$ $\displaystyle=$ $\displaystyle k_{1}k_{2}\cdots
k_{m-1}k_{m}$ (6) $\displaystyle=$ $\displaystyle
2^{m-1}k_{1}+2^{m-2}k_{2}+\cdots+2^{1}\,k_{m-1}+2^{0}\,k_{m}\ ,$

where $k_{m}$ is the least significant bit. The computational basis elements
are then defined by

$\displaystyle|k\rangle$ $\displaystyle=$ $\displaystyle|k_{1}k_{2}\cdots
k_{m-1}k_{m}\rangle$ (7) $\displaystyle=$
$\displaystyle|k_{1}\rangle\otimes|k_{2}\rangle\otimes\cdots\otimes|k_{m-1}\rangle\otimes|k_{m}\rangle~{}~{}\text{where}~{}~{}k_{r}\in\\{0,1\\}\
.$

For example, the state labeled by $k=1$ is represented by

$\displaystyle|1\rangle=|0\cdots
01\rangle=|0\rangle\otimes\cdots\otimes|0\rangle\otimes|1\rangle\ ,$ (8)

with $k_{m}=1$ and all other bits $k_{r}=0$. This state will play a critical
role in Shor’s algorithm. Note that we are using a slightly ambiguous notation
in which $|1\rangle$ is used in different senses on the left- and right-hand
sides of equation (8). However, the meaning of the state $|1\rangle$ will be
clear from context, so this should cause no problems. Furthermore, this bit
convention implies the useful relation

$\displaystyle\frac{k}{M}=\frac{k_{1}}{2^{1}}+\frac{k_{2}}{2^{2}}+\cdots+\frac{k_{m-1}}{2^{m-1}}+\frac{k_{m}}{2^{m}}\
.$ (9)

Note that a sum over the index $k\in\\{0,1,\cdots,M-1\\}$ can be converted
into $m$ sums over the binary components $k_{r}\in\\{0,1\\}$ for
$r\in\\{1,2,\cdots,m\\}$,

$\displaystyle\sum_{k=0}^{M-1}~{}=\sum_{k_{m}\in\\{0,1\\}}\cdots\sum_{k_{2}\in\\{0,1\\}}\sum_{k_{1}\in\\{0,1\\}}\
.$ (10)

This allows us to express the quantum Fourier transform (3) in the form

$\displaystyle QFT\,|\ell\rangle$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi i\,\ell
k/M}\,|k\rangle$ (12) $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{m}}\cdots\sum_{k_{2}}\sum_{k_{1}}e^{2\pi
i\,\ell\big{(}2^{m-1}k_{1}+2^{m-2}k_{2}+\cdots+2^{1}\,k_{m-1}+2^{0}\,k_{m}\big{)}\big{/}2^{m}}\,|k_{1}k_{2}\cdots
k_{m-1}k_{m}\rangle~{}~{}~{}~{}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{1}=0,1}e^{2\pi i\,\ell
k_{1}/2^{1}}|k_{1}\rangle\otimes\sum_{k_{2}=0,1}e^{2\pi i\,\ell
k_{2}/2^{2}}|k_{2}\rangle\otimes\cdots\otimes\sum_{k_{m}=0,1}e^{2\pi i\,\ell
k_{m}/2^{m}}|k_{m}\rangle~{}~{}~{}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{1}}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{2}}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$ $\displaystyle\hskip
56.9055pt\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m-1}}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m}}\,|1\rangle\Big{)}_{m}\ .~{}~{}~{}$

In our current convention, the quantum state indexed by the integer $\ell$ is
given by

$\displaystyle|\ell\rangle$ $\displaystyle=$
$\displaystyle|\ell_{1}\ell_{2}\,\cdots\,\ell_{m-1}\ell_{m}\rangle~{}~{}~{}\text{with}~{}~{}\ell_{r}\in\\{0,1\\}$
(13) $\displaystyle=$
$\displaystyle|\ell_{1}\rangle\otimes|\ell_{2}\rangle\otimes\cdots\otimes|\ell_{m-1}\rangle\otimes|\ell_{m}\rangle\
,$

where $\ell$ takes the binary form

$\displaystyle\ell$ $\displaystyle=$
$\displaystyle\ell_{1}\ell_{2}\cdots\ell_{m-1}\ell_{m}=2^{m-1}\ell_{1}+2^{m-2}\ell_{2}+\cdots+2^{1}\,\ell_{m-1}+2^{0}\,\ell_{m}\
.$ (14)

We shall also introduce the notion of binary fractions corresponding to non-
negative $m$-bit phase angles,

$\displaystyle\Omega\equiv
0.\ell_{1}\ell_{2}\cdots\ell_{m-1}\ell_{m}\equiv\frac{\ell_{1}}{2^{1}}+\frac{\ell_{2}}{2^{2}}+\cdots+\frac{\ell_{m-1}}{2^{m-1}}+\frac{\ell_{m}}{2^{m}}\
.$ (15)

Note that the phase $\Omega=0.\ell_{1}\ell_{2}\cdots\ell_{m}$ and the
corresponding integer index $\ell=\ell_{1}\ell_{2}\cdots\ell_{m}$ are related
by

$\displaystyle\ell=M\Omega=2^{m}\,\Omega\ ,$ (16)

an expression we shall employ throughout the sequel. We now rewrite the
exponential terms in (12) as follows, working slowly through the algebra,
starting with qubit-1:

$\displaystyle 2\pi i\,\frac{\ell}{2^{1}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{1}}\Big{[}\Big{(}2^{m-1}\ell_{1}+2^{m-2}\ell_{2}+\cdots+2^{1}\,\ell_{m-1}\Big{)}+2^{0}\,\ell_{m}\Big{]}$
(17) $\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\Big{(}2^{m-2}\ell_{1}+2^{m-3}\ell_{2}+\cdots+2^{0}\,\ell_{m-1}\Big{)}+\Omega_{1}\Big{]}$
$\displaystyle 2\pi i\,\frac{\ell}{2^{2}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{2}}\Big{[}\Big{(}2^{m-1}\ell_{1}+2^{m-2}\ell_{2}+\cdots
2^{2}\ell_{m-2}\Big{)}+2^{1}\,\ell_{m-1}+2^{0}\,\ell_{m}\Big{]}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\Big{(}2^{m-3}\ell_{1}+2^{m-4}\ell_{2}+\cdots
2^{0}\ell_{m-2}\Big{)}+\Omega_{2}\Big{]}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\cdots$ $\displaystyle 2\pi i\,\frac{\ell}{2^{r}}$
$\displaystyle=$ $\displaystyle\frac{2\pi
i}{2^{r}}\Big{[}\Big{(}2^{m-1}\ell_{1}+2^{m-2}\ell_{2}+\cdots+2^{r}\,\ell_{m-r}\Big{)}+2^{r-1}\,\ell_{m-r+1}$
$\displaystyle\hskip
170.71652pt+\cdots+2^{1}\,\ell_{m-1}+2^{0}\,\ell_{m}\Big{]}$ $\displaystyle=$
$\displaystyle 2\pi
i\Big{[}\Big{(}2^{m-r-1}\ell_{1}+2^{m-r-2}\ell_{2}+\cdots+2^{0}\,\ell_{m-r}\Big{)}+\Omega_{r}\Big{]}$
$\displaystyle\cdots$ $\displaystyle 2\pi i\,\frac{\ell}{2^{m-1}}$
$\displaystyle=$ $\displaystyle\frac{2\pi
i}{2^{m-1}}\Big{[}\big{(}2^{m-1}\ell_{1}\big{)}+2^{m-2}\ell_{2}+\cdots+2^{1}\,\ell_{m-1}+2^{0}\,\ell_{m}\Big{]}$
(20) $\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\big{(}\ell_{1}\big{)}+\Omega_{m-1}\Big{]}$ $\displaystyle 2\pi
i\,\frac{\ell}{2^{m}}$ $\displaystyle=$ $\displaystyle\frac{2\pi
i}{2^{m}}\Big{[}2^{m-1}\ell_{1}+2^{m-2}\ell_{2}+\cdots+2^{1}\,\ell_{m-1}+2^{0}\,\ell_{m}\Big{]}$
(21) $\displaystyle=$ $\displaystyle 2\pi i\Big{[}\Omega_{m}\Big{]}\ ,$

where the partial-phase angles are defined by

$\displaystyle\Omega_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m}}{2^{1}}=0.\ell_{m}$ (22)
$\displaystyle\Omega_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m-1}}{2^{1}}+\frac{\ell_{m}}{2^{2}}=0.\ell_{m-1}\ell_{m}$
$\displaystyle\cdots$ $\displaystyle\Omega_{r}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m-r+1}}{2^{1}}+\cdots+\frac{\ell_{m-1}}{2^{r-1}}+\frac{\ell_{m}}{2^{r}}=0.\ell_{m-r+1}\cdots\ell_{m-1}\ell_{m}$
$\displaystyle\cdots$ $\displaystyle\Omega_{m-1}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{2}}{2^{1}}+\frac{\ell_{3}}{2^{2}}+\cdots+\frac{\ell_{m-2}}{2^{m-1}}+\frac{\ell_{m}}{2^{m-1}}=0.\ell_{2}\cdots\ell_{m-1}\ell_{m}$
(25) $\displaystyle\Omega_{m}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{1}}{2^{1}}+\frac{\ell_{2}}{2^{2}}+\cdots+\frac{\ell_{m-1}}{2^{m-1}}+\frac{\ell_{m}}{2^{m}}=0.\ell_{1}\cdots\ell_{m-1}\ell_{m}\
.$ (26)

Note that integer multiples of $2\pi i$ in the parentheses of equations
(17)–(21) do not contribute, as $e^{2\pi i\,n}=1$ for any integer $n$, thereby
permitting us to express the $QFT$ operation only in terms of the partial
phases

$\displaystyle\Omega_{r}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{r}\,\frac{\ell_{m-r+k}}{2^{k}}~{}~{}~{}\text{for}~{}~{}r\in\\{1,2,\cdots,m\\}\
.$ (27)

The quantum Fourier transform (12) of the $\ell$-state can therefore be
written in any one of three useful forms:

$\displaystyle QFT\,|\ell\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{1}}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{2}}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$ $\displaystyle\hskip
56.9055pt\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m-1}}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m}}\,|1\rangle\Big{)}_{m}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{m}}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{m-1}\ell_{m}}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$ (30)
$\displaystyle\hskip 42.67912pt\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{2}\cdots\ell_{m-1}\ell_{m}}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{1}\ell_{2}\cdots\ell_{m-1}\ell_{m}}\,|1\rangle\Big{)}_{m}$
$\displaystyle=$ $\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{2}}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$ $\displaystyle\hskip
71.13188pt\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-1}}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m}}\,|1\rangle\Big{)}_{m}\ .$

We next reverse the order of the qubits with a string of SWAP gates to form
the state

$\displaystyle|\psi_{\rm rev}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m}}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-1}}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$
$\displaystyle\hskip 136.5733pt\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{2}}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}}\,|1\rangle\Big{)}_{m}\ .$

The state $|\psi_{\rm rev}\rangle$ can be represented quite readily by a
quantum circuit. To see this, let us start with qubit-$m$ of (II.2.1). Since
$e^{2\pi i\,\Omega_{1}}=e^{2\pi i\,\ell_{m}/2}=(-1)^{\ell_{m}}$, and since
$\ell_{m}$ takes the binary values $0$ and $1$, we have

$\displaystyle\left.\begin{matrix}\displaystyle{\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}(\ell_{m}=0)}\,|1\rangle\Big{)}_{m}=\frac{|0\rangle_{m}+|1\rangle_{m}}{\sqrt{2}}=H|0\rangle_{m}}\\\\[15.0pt]
\displaystyle{\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}(\ell_{m}=1)}\,|1\rangle\Big{)}_{m}=\frac{|0\rangle_{m}-|1\rangle_{m}}{\sqrt{2}}=H|1\rangle_{m}}\end{matrix}~{}\right\\}=H|\ell_{m}\rangle~{}~{}~{}~{}~{}~{}~{}~{}$
(32)

for $\ell_{m}\in\\{0,1\\}$, where the single-qubit Hadamard gate is defined by

$\displaystyle H=\frac{1}{\sqrt{2}}\left[\begin{matrix}1&~{}1\\\
1&-1\end{matrix}\right]~{}~{}~{}\text{with}~{}~{}|0\rangle=\left[\begin{matrix}~{}1&\\\
~{}0&\end{matrix}\right]~{}~{}\text{and}~{}~{}|1\rangle=\left[\begin{matrix}~{}0&\\\
~{}1&\end{matrix}\right]\ .$ (33)

Figure 2: Convention 2 for the $QFT$: The standard physics/mathematics
convention.

We will often place subscripts on the single-qubit basis states such as
$|\ell_{m}\rangle_{m}$ to explicitly indicate their qubit position in the
quantum circuit. We will also denote the qubit upon which the Hadamard gate
acts by a superscript, e.g. $H^{m}$ explicitly states that $H$ acts on the
$m$-th qubit. Therefore, we can express (32) in the form

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}}\,|1\rangle\Big{)}_{m}$ $\displaystyle=$ $\displaystyle
H^{m}|\ell_{m}\rangle_{m}~{}~{}\text{with}~{}~{}\ell_{m}\in\\{0,1\\}\ .$ (34)

Moving on to the next qubit, $m-1$, we find

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{2}}\,|1\rangle\Big{)}_{m-1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,\Big{(}|0\rangle+e^{2\pi
i\,\ell_{m-1}/2^{1}}\cdot e^{2\pi i\,\ell_{m}/2^{2}}\,|1\rangle\Big{)}_{m-1}\
.$ (35)

The first exponential $e^{2\pi i\,\ell_{m-1}/2}=e^{\pi
i\,\ell_{m-1}}=(-1)^{\ell_{m-1}}$ gives a Hadamard gate acting on
$|\ell_{m-1}\rangle_{m-1}$, and the second exponential $e^{2\pi
i\,\ell_{m}/2^{2}}$ produces a controlled phase gate with angle
$\theta=2\pi/2^{2}$ (call it $CR_{2}$) acting on the target qubit
$|\ell_{m-1}\rangle_{m-1}$ with the control qubit $|\ell_{m}\rangle_{m}$, so
that

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{2}}\,|1\rangle\Big{)}_{m-1}$ $\displaystyle=$ $\displaystyle
H^{m-1}\cdot
C^{m}R_{2}^{m-1}\,|\ell_{m-1}\rangle_{m-1}~{}~{}\text{where}~{}~{}\ell_{m-1}\in\\{0,1\\}\
.~{}~{}~{}$ (36)

We have used superscripts on the controlled-$R$ gate, writing $C^{c}R_{2}^{t}$
to explicitly indicate the control qubit $c$ and the target qubit $t$.
Finally, let us examine the $1$-st qubit, where we find

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m}}\,|1\rangle\Big{)}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi i\,\ell_{1}/2}\cdot
e^{2\pi i\,\ell_{2}/2^{2}}\cdots e^{2\pi i\,\ell_{m-1}/2^{m-1}}\cdot e^{2\pi
i\,\ell_{m}/2^{m}}\,|1\rangle\Big{)}_{1}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(37) $\displaystyle=$ $\displaystyle H^{1}\cdot
C^{2}R_{2}^{1}\,\cdots\,C^{m-1}R_{m-1}^{1}\cdot
C^{m}R_{m}^{1}\,|\ell_{1}\rangle_{1}~{}~{}~{}\text{with}~{}~{}\ell_{1}\in\\{0,1\\}\
,$ (38)

where the single-qubit phase operator is defined by

$\displaystyle R_{n}$ $\displaystyle=$
$\displaystyle\left[\begin{matrix}1&0\\\\[-5.0pt] 0&e^{2\pi
i/2^{n}}\end{matrix}\right]\ .$ (39)

In terms of the standard phase gate $P(\theta)$, we can express the phase
rotation by

$\displaystyle R_{n}=P(\theta_{n})=\left[\begin{matrix}1&0\\\\[-5.0pt]
0&e^{i\theta_{n}}\end{matrix}\right]~{}~{}~{}\text{where}~{}~{}~{}\theta_{n}=\frac{2\pi}{2^{n}}=\frac{\pi}{2^{n-1}}\
.$ (40)

Figure 2 reproduces the corresponding $QFT$ circuit. Note that SWAP gates are
required at the end of the circuit to place the states back into their
original order. The inverse $QFT^{\dagger}$ is given by reading the circuit
backwards from right to left, starting with the SWAP gates, and inverting all
phase angles (replacing $R_{n}$ by $R_{n}^{\dagger}$). The circuit uses
$1+2+\cdots+m=\frac{1}{2}m(m+1)={\cal O}(m^{2})$ distinct gates, plus ${\cal
O}(m/2)$ SWAP gates.

What would have happened if we had not used the SWAP gates to reverse the
qubit order in state (II.2.1), but instead appealed directly to (30)? We would
have found terms like $|\ell_{m}\rangle_{1}$, and the index of $\ell$ would
not have paired properly with the associated qubit. By performing the SWAP
operations, we only find states of the form $|\ell_{k}\rangle_{k}$. It should
therefore cause no confusion if we henceforth drop the subscript on the basis
states and simply write $|\ell_{k}\rangle$.

####  II.2.2 Convention 1: OpenQASM/Qiskit

We now look at the convention used by OpenQASM/Qiskit. Consider an $m$ qubit
system in which the circuit for the computational basis states start with
qubit $0$ in the upper position and qubit $m-1$ in the lower position. There
are $M=2^{m}$ quantum states in the system, and the index integer
$k\in\\{0,1,\cdots,M-1\\}$ can be expressed in a binary form where $k_{0}$ is
the least significant bit,

$\displaystyle k$ $\displaystyle=$ $\displaystyle k_{m-1}k_{m-2}\cdots
k_{1}k_{0}$ (41) $\displaystyle=$ $\displaystyle
2^{m-1}k_{m-1}+2^{m-2}k_{m-2}+\cdots+2^{1}\,k_{1}+2^{0}\,k_{0}\ .$

The integer $k$ can then be used to label the computational basis states,

$\displaystyle|k\rangle$ $\displaystyle=$ $\displaystyle|k_{m-1}k_{m-2}\cdots
k_{1}k_{0}\rangle$ (42) $\displaystyle=$
$\displaystyle|k_{0}\rangle\otimes|k_{1}\rangle\otimes\cdots\otimes|k_{m-2}\rangle\otimes|k_{m-1}\rangle\
.$

In this convention, the state labeled by $k=1$ is given by

$\displaystyle|1\rangle=|0\cdots
01\rangle=|1\rangle\otimes|0\rangle\otimes\cdots\otimes|0\rangle\ ,$ (43)

with $k_{0}=1$ and all other bits $k_{r}=0$. As previously mentioned, this
state will play a critical role in Shor’s algorithm. We also record here the
convenient relation

$\displaystyle\frac{k}{M}=\frac{k_{m-1}}{2^{1}}+\frac{k_{m-2}}{2^{2}}+\cdots+\frac{k_{1}}{2^{m-1}}+\frac{k_{0}}{2^{m}}\
.$ (44)

As before, we can replace a sum over the index integer $k$ by $m$ sums over
the binary components $k_{r}\in\\{0,1\\}$ for $r\in\\{0,1,\cdots,m-1\\}$,

$\displaystyle\sum_{k=0}^{M-1}~{}=\sum_{k_{m-1}\in\\{0,1\\}}\cdots\sum_{k_{1}\in\\{0,1\\}}\sum_{k_{0}\in\\{0,1\\}}\
,$ (45)

so that expression (3) for the quantum Fourier transform becomes

$\displaystyle QFT\,|\ell\rangle$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi i\,\ell
k/M}\,|k\rangle$ (47) $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{m-1}}\cdots\sum_{k_{1}}\sum_{k_{0}}e^{2\pi
i\,\ell\big{(}2^{0}\,k_{0}+2^{1}\,k_{1}+\cdots+2^{m-2}k_{m-2}+2^{m-1}k_{m-1}\big{)}\big{/}2^{m}}\,|k_{m-1}k_{m-2}\cdots
k_{1}k_{0}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{0}=0,1}e^{2\pi i\,\ell
k_{0}/2^{m}}|k_{0}\rangle\otimes\sum_{k_{1}=0,1}e^{2\pi i\,\ell
k_{1}/2^{m-1}}|k_{1}\rangle\otimes\cdots\otimes\\!\\!\sum_{k_{m-1}=0,1}e^{2\pi
i\,\ell k_{m-1}/2^{1}}|k_{m-1}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m}}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m-1}}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$
$\displaystyle\hskip 56.9055pt\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{2}}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{1}}\,|1\rangle\Big{)}_{m-1}\ .~{}~{}~{}$

Next we use the form of the $\ell$-state given by our OpenQASM convention,

$\displaystyle|\ell\rangle$ $\displaystyle=$
$\displaystyle|\ell_{m-1}\ell_{m-2}\cdots\ell_{1}\ell_{0}\,\rangle~{}~{}~{}\text{with}~{}~{}\ell_{r}\in\\{0,1\\}$
(48) $\displaystyle=$
$\displaystyle|\ell_{0}\rangle\otimes|\ell_{1}\rangle\otimes\cdots\otimes|\ell_{m-2}\rangle\otimes|\ell_{m-1}\rangle\
,$

where the $m$-bit integer $\ell$ takes the binary form

$\displaystyle\ell$ $\displaystyle=$
$\displaystyle\ell_{m-1}\ell_{m-2}\cdots\ell_{1}\ell_{0}$ (49)
$\displaystyle=$ $\displaystyle
2^{m-1}\ell_{m-1}+2^{m-2}\ell_{m-2}+\cdots+2^{1}\,\ell_{1}+2^{0}\,\ell_{0}\ ,$

while $m$-bit binary fractions can be expressed by

$\displaystyle\Omega=0.\ell_{m-1}\ell_{m-2}\cdots\ell_{1}\ell_{0}\equiv\frac{\ell_{m-1}}{2^{1}}+\frac{\ell_{m-2}}{2^{2}}+\cdots+\frac{\ell_{1}}{2^{m-1}}+\frac{\ell_{0}}{2^{m}}\
.$ (50)

Similarly to the previous case, the exponential terms in (47) can be written

$\displaystyle 2\pi i\,\frac{\ell}{2^{m}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{m}}\Big{[}2^{m-1}\ell_{m-1}+2^{m-2}\ell_{m-2}+\cdots+2^{1}\,\ell_{1}+2^{0}\,\ell_{0}\Big{]}$
(51) $\displaystyle=$ $\displaystyle 2\pi i\Big{[}\Omega_{0}\Big{]}$
$\displaystyle 2\pi i\,\frac{\ell}{2^{m-1}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{m-1}}\Big{[}\big{(}2^{m-1}\ell_{m-1}\big{)}+2^{m-2}\ell_{m-2}+\cdots+2^{1}\,\ell_{1}+2^{0}\,\ell_{0}\Big{]}$
$\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\big{(}\ell_{m-1}\big{)}+\Omega_{1}\Big{]}$ $\displaystyle\cdots$
$\displaystyle 2\pi i\,\frac{\ell}{2^{m-r}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{m-r}}\Big{[}\Big{(}2^{m-1}\ell_{m-1}+2^{m-2}\ell_{m-2}+\cdots+2^{m-r}\,\ell_{m-r}\Big{)}+2^{m-r-1}\,\ell_{m-r-1}$
$\displaystyle\hskip 207.7052pt+\cdots+2^{1}\,\ell_{1}+2^{0}\,\ell_{0}\Big{]}$
$\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\Big{(}2^{r-1}\ell_{m-1}+2^{r-2}\ell_{m-2}+\cdots+2^{0}\,\ell_{m-r}\Big{)}+\Omega_{r}\Big{]}$
$\displaystyle\cdots$ $\displaystyle 2\pi i\,\frac{\ell}{2^{2}}$
$\displaystyle=$ $\displaystyle\frac{2\pi
i}{2^{2}}\Big{[}\Big{(}2^{m-1}\ell_{m-1}+2^{m-2}\ell_{m-2}+\cdots
2^{2}\ell_{2}\Big{)}+2^{1}\,\ell_{1}+2^{0}\,\ell_{0}\Big{]}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(54) $\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\Big{(}2^{m-3}\ell_{m-1}+2^{m-4}\ell_{m-2}+\cdots
2^{0}\ell_{2}\Big{)}+\Omega_{m-2}\Big{]}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle 2\pi i\,\frac{\ell}{2^{1}}$ $\displaystyle=$
$\displaystyle\frac{2\pi
i}{2^{1}}\Big{[}\Big{(}2^{m-1}\ell_{m-1}+2^{m-2}\ell_{m-2}+\cdots+2^{1}\,\ell_{1}\Big{)}+2^{0}\,\ell_{0}\Big{]}$
(55) $\displaystyle=$ $\displaystyle 2\pi
i\Big{[}\Big{(}2^{m-2}\ell_{m-1}+2^{m-3}\ell_{m-2}+\cdots+2^{0}\,\ell_{1}\Big{)}+\Omega_{m-1}\Big{]}\
,$

where the partial phases are now defined by

$\displaystyle\Omega_{0}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m-1}}{2^{1}}+\frac{\ell_{m-2}}{2^{2}}+\cdots+\frac{\ell_{1}}{2^{m-1}}+\frac{\ell_{0}}{2^{m}}=0.\ell_{m-1}\cdots\ell_{1}\ell_{0}$
(56) $\displaystyle\Omega_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m-2}}{2^{1}}+\frac{\ell_{m-3}}{2^{2}}+\cdots+\frac{\ell_{1}}{2^{m-2}}+\frac{\ell_{0}}{2^{m-1}}=0.\ell_{m-2}\cdots\ell_{1}\ell_{0}$
$\displaystyle\cdots$ $\displaystyle\Omega_{r}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{m-r-1}}{2^{1}}+\frac{\ell_{m-r-2}}{2^{2}}+\cdots+\frac{\ell_{1}}{2^{m-r-1}}+\frac{\ell_{0}}{2^{m-r}}=0.\ell_{m-r-1}\cdots\ell_{1}\ell_{0}$
$\displaystyle\cdots$ $\displaystyle\Omega_{m-2}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{1}}{2^{1}}+\frac{\ell_{0}}{2^{2}}=0.\ell_{1}\ell_{0}$
(59) $\displaystyle\Omega_{m-1}$ $\displaystyle\equiv$
$\displaystyle\frac{\ell_{0}}{2^{1}}=0.\ell_{0}\ .$ (60)

The general phase takes the form

$\displaystyle\Omega_{r}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{m-r}\,\frac{\ell_{m-r-k}}{2^{k}}~{}~{}~{}\text{for}~{}~{}r\in\\{0,1,\cdots,m-1\\}\
,$ (61)

and the quantum Fourier transform (47) can now be expressed in one of three
equivalent ways:

$\displaystyle QFT\,|\ell\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m}}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{m-1}}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$
$\displaystyle\hskip 56.9055pt\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{2}}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\ell/2^{1}}\,|1\rangle\Big{)}_{m-1}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{0}}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{1}\ell_{0}}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$ (64)
$\displaystyle\hskip 51.21504pt\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{m-2}\cdots\ell_{1}\ell_{0}}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,0.\ell_{m-1}\cdots\ell_{1}\ell_{0}}\,|1\rangle\Big{)}_{m-1}$
$\displaystyle=$ $\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{0}}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$ $\displaystyle\hskip
71.13188pt\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-2}}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-1}}\,|1\rangle\Big{)}_{m-1}\ .$

As before, we invert the qubits with SWAP gates to form the state

$\displaystyle|\psi_{\rm rev}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-1}}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-2}}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$
$\displaystyle\hskip 142.26378pt\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{1}}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{0}}\,|1\rangle\Big{)}_{m-1}\ ,$

which can be expressed in terms of basic gates to give the following circuit.
Starting with qubit-$0$ of (II.2.2), and using $e^{2\pi
i\,\Omega_{m-1}}=e^{2\pi i\,\ell_{0}/2}=(-1)^{\ell_{0}}$, we find a state
similar to the previous case,

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{m-1}}\,|1\rangle\Big{)}_{0}$ $\displaystyle=$ $\displaystyle
H^{0}|\ell_{0}\rangle~{}~{}~{}\text{where}~{}~{}\ell_{0}\in\\{0,1\\}\ .$ (66)

The other qubits give corresponding results, so we move on to the bottom of
the circuit, expressing qubit $m-1$ by

$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
i\,\Omega_{0}}\,|1\rangle\Big{)}_{m-1}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi i\,\ell_{m-1}/2}\cdot
e^{2\pi i\,\ell_{m-2}/2^{2}}\cdots e^{2\pi i\,\ell_{1}/2^{m-1}}\cdot e^{2\pi
i\,\ell_{0}/2^{m}}\,|1\rangle\Big{)}_{m-1}$ (67) $\displaystyle=$
$\displaystyle H^{m-1}\cdot C^{m-2}R_{2}^{m-1}\cdots C^{1}R_{m-1}^{m-1}\cdot
C^{0}R_{m}^{m-1}\,|\ell_{m-1}\rangle$

with $\ell_{m-1}\in\\{0,1\\}$, where the phase operators are defined by

$\displaystyle R_{n}$ $\displaystyle=$
$\displaystyle\left[\begin{matrix}1&0\\\\[-5.0pt] 0&e^{2\pi
i/2^{n}}\end{matrix}\right]=P\big{(}\pi/2^{n-1}\big{)}\ .$ (68)

Figure 3 reproduces the commensurate QFT circuit in the OpenQASM/Qiskit
convention. The inverse $QFT^{\dagger}$ is given by reading the circuit in
reverse order from right to left, starting with the SWAP gates, and inverting
the sign of the phase gates. Note that the QFT circuits are reversed between
the two conventions.

Figure 3: Convention 1 for the $QFT$: OpenQASM/Qiskit.

##  III Quantum Phase Estimation

We now turn to quantum phase estimation (QPE), which is the workhorse of
Shor’s algorithm. Our primary references for this section are Refs. [8] and
[9]. We shall work through the calculations in both Conventions 1 and 2 as
outlined in Section II.2, corresponding to the OpenQASM and the standard
physics conventions, respectively.

###  III.1 Convention 2: Standard Physics/Mathematics

We first examine Convention 2, the traditional physics and mathematics
convention. Consider a linear unitary operator $U$ with Eigenvalue $e^{2\pi
i\theta}$ and Eigenstate $|u\rangle$, where $\theta$ is a real phase such that
$0\leq\theta<1$:

$\displaystyle U|u\rangle=e^{2\pi i\theta}\,|u\rangle\ .$ (69)

Since $U$ is unitary, its Eigenvalues have a norm of unity. We wish to build a
QPE circuit that will output an approximate ($m$-bit binary) value for the
phase angle $\theta$. The circuit will consist of a front-end and a back-end.
There are two registers in the QPE front-end circuit: (i) a control register
consisting of $m$ qubits and (ii) a work register containing $n$ qubits. We
store the Eigenstate $|u\rangle$ in the work register. To construct the front-
end circuit, we first apply Hadamard gates to every control qubit, forming the
state

$\displaystyle|\psi_{1}\rangle=H|0\rangle\otimes\cdots\otimes
H|0\rangle\otimes|u\rangle=\frac{1}{\sqrt{M}}\sum_{k-0}^{M-1}|k\rangle\otimes|u\rangle\
,$ (70)

where $M=2^{m}$ is the total number of computational basis states. Given the
unitary operator $U$, we assume that we are able to build a family of $m$
controlled-$U^{p}$ gates for $p\in\\{2^{0},2^{1},\cdots,2^{m-1}\\}$ that
operate on the work register containing the state $|u\rangle$. Note that the
action of a $CU^{p}$ operator for a single control qubit takes the form

$\displaystyle CU^{p}H|0\rangle\otimes|u\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\,CU^{p}\,\Big{(}|0\rangle+|1\rangle\Big{)}\otimes|u\rangle=\frac{1}{\sqrt{2}}\,|0\rangle\otimes|u\rangle+\frac{1}{\sqrt{2}}\,|1\rangle\otimes
U^{p}|u\rangle$ (71) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}|0\rangle\otimes|u\rangle+|1\rangle\otimes
e^{2\pi ip\,\theta}|u\rangle\Big{)}=\frac{1}{\sqrt{2}}\Big{(}|0\rangle+e^{2\pi
ip\,\theta}\,|1\rangle\Big{)}\otimes|u\rangle\ .~{}~{}~{}~{}$ (72)

This is an example of phase kickback, in which the phase operation in the
target register makes its way back into the control register. We now string
these gates together to form the front-end of the circuit composed of the
gates $C^{n}U^{2^{n}}$ for $n\in\\{0,1,\cdots,m-1\\}$, as illustrated in Fig.
4, with the least significant power $p=2^{0}$ attached to the least
significant $m$-th target qubit (as Convention 2 dictates). We see that the
output state of the front-end becomes

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,2^{m-1}\theta}\,|1\rangle\Big{)}_{1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,2^{m-2}\theta}\,|1\rangle\Big{)}_{2}\otimes\cdots\otimes$
$\displaystyle\hskip 142.26378pt\Big{(}|0\rangle+e^{2\pi
i\,2^{1}\theta}\,|1\rangle\Big{)}_{m-1}\otimes\Big{(}|0\rangle+e^{2\pi
i\,2^{0}\theta}\,|1\rangle\Big{)}_{m}\otimes|u\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{1}=0,1}e^{2\pi
i\,\theta\,2^{m-1}k_{1}}|k_{1}\rangle\otimes\sum_{k_{2}=0,1}e^{2\pi
i\,\theta\,2^{m-2}k_{2}}|k_{2}\rangle\otimes\cdots\otimes\sum_{k_{m}=0,1}e^{2\pi
i\,\theta\,2^{0}k_{m}}|k_{m}\rangle\otimes|u\rangle~{}~{}~{}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\sum_{k_{1}}\sum_{k_{2}}\cdots\sum_{k_{m}}e^{2\pi
i\,\theta\big{(}2^{m-1}k_{1}+2^{m-2}k_{2}+\cdots+2^{0}k_{m}\big{)}}|k_{1}\rangle\otimes|k_{2}\rangle\otimes\dots\otimes|k_{m}\rangle\otimes|u\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
ik\,\theta}\,|k\rangle\otimes|u\rangle\ ,$ (74)

where $M=2^{m}$ and $k=2^{m-1}k_{1}+2^{m-2}k_{2}+\cdots
2^{1}k_{m-1}+2^{0}k_{m}$, with $k_{r}\in\\{0,1\\}$.

Figure 4: QPE front-end: Convention 2 (physics and mathematics).

This result is valid for a general phase angle $\theta$.

For simplicity, let us first suppose that the binary form of the phase angle
terminates after exactly $m$ bits, so that

$\displaystyle\theta$ $\displaystyle=$ $\displaystyle
0.\theta_{1}\,\theta_{2}\,\cdots\theta_{m-1}\theta_{m}~{}~{}~{}\text{where}~{}~{}\theta_{r}\in\\{0,1\\}$
(75) $\displaystyle=$
$\displaystyle\frac{\theta_{1}}{2}+\frac{\theta_{2}}{2^{2}}+\cdots+\frac{\theta_{m-1}}{2^{m-1}}+\frac{\theta_{m}}{2^{m}}\
.$ (76)

We will shortly extend the argument to general phase angles that do not not
terminate. We see that (76) implies that $M\theta$ can be written as a binary
integer,

$\displaystyle\ell_{\theta}\equiv M\theta=2^{m}\theta$ $\displaystyle=$
$\displaystyle
2^{m-1}\,\theta_{1}+2^{m-2}\,\theta_{2}+\cdots+2^{0}\,\theta_{m}$ (77)
$\displaystyle=$
$\displaystyle\theta_{1}\,\theta_{2}\,\cdots\,\theta_{m}\in\\{0,1,\cdots,M-1\\}\
.$ (78)

Figure 5: QPE: Convention 2. The QPE circuit with the inverse quantum Fourier
transform $QFT^{\dagger}$ represented abstractly by a box.

Upon using the relation $\theta=\ell_{\theta}/M$, we can now express the
output state of the front-end as

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
i\,k\ell_{\theta}/M}|k\rangle\otimes|u\rangle$ (79) $\displaystyle=$
$\displaystyle QFT\,|\ell_{\theta}\rangle\otimes|u\rangle\ ,$ (80)

where we have used the definition of the $QFT$ operator (3). Therefore, the
back-end of the QPE circuit will consist of an inverse $QFT$ acting on the
control register, as illustrated in Fig. 5. In Fig. 6 we expand the
$QFT^{\dagger}$ circuit explicitly. In either case, the final output state is
given by

$\displaystyle|\psi_{3}\rangle=QFT^{\dagger}\,|\psi_{2}\rangle=|\ell_{\theta}\rangle\otimes|u\rangle\
,$ (81)

where $\ell_{\theta}\equiv 2^{m}\theta\in\\{0,1,\cdots,M-1\\}$. Some authors
denote the state $|\ell_{\theta}\rangle$ by $|2^{m}\theta\rangle$. Upon
measuring the control register, we will find
$\ell_{\theta}=\theta_{1}\cdots\theta_{m}$ for $\theta_{r}\in\\{0,1\\}$, and
the corresponding phase is then given exactly by
$\theta=\ell_{\theta}/2^{m}=0.\theta_{1}\cdots\theta_{m}$, in agreement with
(75). We see that the QPE circuit in Fig. 5 does indeed extract the correct
phase.

Figure 6: The $QFT^{\dagger}$ circuit has been expanded to illustrate the full
complexity of the QPE circuit.

###  III.2 Convention 1: OpenQASM/Qiskit

We now examine the OpenQASM/Qiskit convention. As before, we string the
$CU^{p}$ operators together to form the front-end of the circuit, with the
smallest power $p$ being attached to the $0$-th qubit (the upper and least
significant qubit of the circuit). The corresponding front-end is illustrated
in Fig. 7.

Figure 7: QPE front-end: Convention 1 (OpenQASM/Qiskit).

Similarly to (74), the output of the front-end can be expressed by

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2^{m/2}}\Big{(}|0\rangle+e^{2\pi
i\,2^{0}\theta}\,|1\rangle\Big{)}_{0}\otimes\Big{(}|0\rangle+e^{2\pi
i\,2^{1}\theta}\,|1\rangle\Big{)}_{1}\otimes\cdots\otimes$
$\displaystyle\hskip 25.6073pt\Big{(}|0\rangle+e^{2\pi
i\,2^{m-2}\theta}\,|1\rangle\Big{)}_{m-2}\otimes\Big{(}|0\rangle+e^{2\pi
i\,2^{m-1}\theta}\,|1\rangle\Big{)}_{m-1}\otimes|u\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
ik\,\theta}\,|k\rangle\otimes|u\rangle\ ,$ (83)

where $M=2^{m}$ and $k=2^{m-1}k_{m-1}+2^{m-2}k_{m-2}+\cdots
2^{1}k_{1}+2^{0}k_{0}$, with $k_{r}\in\\{0,1\\}$. We first consider the simple
case of an $m$-bit phase angle,

$\displaystyle\theta$ $\displaystyle=$ $\displaystyle
0.\theta_{m-1}\,\theta_{m-2}\,\cdots\,\theta_{1}\,\theta_{0}~{}~{}~{}\text{where}~{}~{}\theta_{r}\in\\{0,1\\}$
(84) $\displaystyle=$
$\displaystyle\frac{\theta_{m-1}}{2^{1}}+\frac{\theta_{m-2}}{2^{2}}+\cdots+\frac{\theta_{1}}{2^{m-1}}+\frac{\theta_{0}}{2^{m}}\
.$

As before, this expression implies that $M\theta$ is a binary integer between
$0$ and $M-1$:

$\displaystyle\ell_{\theta}\equiv M\theta=2^{m}\theta$ $\displaystyle=$
$\displaystyle
2^{m-1}\,\theta_{m-1}+2^{m-2}\,\theta_{m-2}+\cdots+2^{1}\,\theta_{1}+2^{0}\,\theta_{0}$
(85) $\displaystyle=$
$\displaystyle\theta_{m-1}\,\theta_{m-2}\,\cdots\,\theta_{1}\,\theta_{0}\in\\{0,1,\cdots,M-1\\}\
.$ (86)

Again, the relation $\theta=\ell_{\theta}/M$ and definition (3) of the QFT
imply that the output of the front-end becomes

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
ik\,\theta}|k\rangle\otimes|u\rangle=\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
i\,k\ell_{\theta}/M}|k\rangle\otimes|u\rangle$ (87) $\displaystyle=$
$\displaystyle QFT\,|\ell_{\theta}\rangle\otimes|u\rangle\ .$

Upon taking the inverse Fourier transform, the final state of the QPE circuit
is

$\displaystyle|\psi_{3}\rangle=QFT^{\dagger}\,|\psi_{2}\rangle=|\ell_{\theta}\rangle\otimes|u\rangle\
.$ (88)

The back-end of the QPE circuit therefore consists of an inverse $QFT$
operator, as illustrated in Fig. 8. We also give the full circuit for the
inverse $QFT$ in Fig. 9. We see that a measurement of the control register
gives the integer $\ell_{\theta}=\theta_{m-1}\cdots\theta_{0}$, which in turn
provides the correct measured phase
$\theta=\ell_{\theta}/2^{m}=0.\theta_{m-1}\cdots\theta_{0}$, in agreement with
(84).

Figure 8: QPE front- and back-end: Convention 1 (OpenQASM/Qiskit)  Figure 9:
QPE: Convention 1. Expanded $QFT^{\dagger}$.

###  III.3 General Phase Angles

We now turn to the case in which the phase angle $\theta$ in (69) is a general
real number between 0 and 1. The analysis here applies for both Conventions.
Recall that the output of the front-end in Fig. 4 or Fig. 7 is

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}e^{2\pi
ik\,\theta}\,|k\rangle\otimes|u\rangle\ ,$ (89)

where this result holds for general phase angles $\theta$. The output of the
QPE circuit is then given by
$|\psi_{3}\rangle=QFT^{\dagger}\,|\psi_{2}\rangle$, where the inverse of the
$QFT$ operator takes the form

$\displaystyle
QFT^{\dagger}=\frac{1}{\sqrt{M}}\sum_{\ell=0}^{M-1}\sum_{k=0}^{M-1}e^{-2\pi
i\,k\ell/M}\,|\ell\rangle\langle k|\ .$ (90)

We therefore obtain the final state

$\displaystyle|\psi_{3}\rangle$ $\displaystyle=$ $\displaystyle
QFT^{\dagger}|\psi_{2}\rangle=\frac{1}{M}\sum_{\ell=0}^{M-1}\sum_{k=0}^{M-1}e^{-2\pi
i\,k\ell/M}\,e^{2\pi ik\,\theta}\,|\ell\rangle\otimes|u\rangle$ (91)
$\displaystyle=$
$\displaystyle\frac{1}{M}\sum_{\ell=0}^{M-1}\sum_{k=0}^{M-1}\big{[}e^{2\pi
i\,(\theta-\ell/M)}\big{]}^{k}\,|\ell\rangle\otimes|u\rangle=\frac{1}{M}\sum_{\ell=0}^{M-1}\frac{1-e^{2\pi
i\,(\theta-\ell/M)M}}{1-e^{2\pi
i\,(\theta-\ell/M)}}\,|\ell\rangle\otimes|u\rangle\ ,~{}~{}~{}~{}~{}~{}$ (92)

where we have performed an exact finite geometric sum. Since this is such an
important result, we summarize it below:

$\displaystyle|\psi_{3}\rangle=\sum_{\ell=0}^{M-1}A_{\ell}(\theta)\,|\ell\rangle\otimes|u\rangle$
(93)

with amplitudes

$\displaystyle A_{\ell}(\theta)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{M}\,\frac{1-e^{2\pi i\,(\theta-\ell/M)M}}{1-e^{2\pi
i\,(\theta-\ell/M)}}$ (94) $\displaystyle=$
$\displaystyle\frac{1}{M}\,\frac{1-e^{2\pi
i\,(\theta-\theta_{\ell})M}}{1-e^{2\pi
i\,(\theta-\theta_{\ell})}}=\frac{1}{M}\,\frac{1-e^{2\pi
i\,(\ell_{\theta}-\ell)}}{1-e^{2\pi i\,(\ell_{\theta}-\ell)/M}}\ .$ (95)

In (94) we have expressed the amplitude in terms of the fundamental quantities
$\theta$ and $\ell$, while (95) expresses the amplitude in terms of the
$m$-bit measured phase $\theta_{\ell}\equiv\ell/M$ and the “mode number”
$\ell_{\theta}\equiv\theta M$ (which might or might not be an integer). We now
find that the probability of measuring the $\ell$-th state is

$\displaystyle P_{\ell}(\theta)=\big{|}A_{\ell}(\theta)\big{|}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{M^{2}}\,\frac{\sin^{2}\left[\pi\left(\theta-\displaystyle\frac{\ell}{M}\right)M\right]}{\sin^{2}\left[\pi\left(\theta-\displaystyle\frac{\ell}{M}\right)\right]}$
(96) $\displaystyle=$
$\displaystyle\frac{1}{M^{2}}\,\frac{\sin^{2}\left[\pi\big{(}\theta-\theta_{\ell}\big{)}M\right]}{\sin^{2}\left[\pi\big{(}\theta-\theta_{\ell}\big{)}\right]}=\frac{1}{M^{2}}\,\frac{\sin^{2}\left[\pi\big{(}\ell_{\theta}-\ell\big{)}\right]}{\sin^{2}\left[\pi\big{(}\ell_{\theta}-\ell\big{)}/M\right]}\
.~{}~{}~{}~{}~{}$ (97)

The probability $P_{\ell}(\theta)$ is maximum for the state
$\ell\in\\{0,1,\cdots,M-1\\}$ for which
$\delta=\theta-\ell/M=\theta-\theta_{\ell}$ is minimum. For large values of
$M=2^{m}$, the probability $P_{\ell}(\theta)$ is sharply peaked about
$\theta$, as illustrated in Fig. 10. We note that the value of
$\ell_{\theta}\equiv M\theta$ need not be an $m$-bit integer (nor indeed, an
integer at all). The point to be emphasized is that when $\ell_{\theta}$ is
not an integer, then the QPE circuit becomes probabilistic in nature. This
will turn out to be a key feature of Shor’s factoring circuit.

Figure 10: The probability $P_{\ell}(\theta)=|A_{\ell}(\theta)|^{2}$ as a
function of the measured phase angle $\theta_{\ell}\equiv\ell/M$ for
$M=2^{4}=16$, where $\theta=1/3$. Note that $A_{\ell}(\theta)\approx 1$ for
$\theta_{\ell}\approx\theta$.

In fact, this analysis implicitly assumes that $\ell_{\theta}=M\theta$ in
amplitude (94) is not an integer, as the expression for $A_{\ell}(\theta)$
then requires a more delicate treatment. To see this, let us write the
amplitude in terms of the quantity $\ell_{\theta}$, so that

$\displaystyle A_{\ell,\ell_{\theta}}\equiv
A_{\ell}(\ell_{\theta}/M)=\frac{1}{M}\,\frac{1-e^{2\pi
i\,(\ell_{\theta}-\ell)}}{1-e^{2\pi i\,(\ell_{\theta}-\ell)/M}}\ .$ (98)

We have already considered the case in which $\ell_{\theta}$ is a non-negative
integer in the previous section, and consistency demands that amplitude (98),
or equivalently (94), must collapses to the Kronecker-delta form

$\displaystyle A_{\ell}(\theta)=\delta_{\ell,\ell_{\theta}}$ (99)

when $\ell_{\theta}$ becomes an integer. We will now show how the expression
for $A_{\ell}(\theta)$ in (98) reduces to this simpler form. Note that the
numerator in (98) vanishes for every integer $\ell_{\theta}$ (and therefore
for every value of $\ell_{\theta}-\ell$, as $\ell$ is itself an integer).
Consequently, the only way one can obtain a non-zero probability is when the
denominator also vanishes, so that we have the indeterminate form $0/0$.
However, the denominator vanishes only for $\ell_{\theta}-\ell$ such that

$\displaystyle(\ell_{\theta}-\ell)/M=n~{}~{}\text{for any}~{}n\in\mathbb{Z}\
,$ (100)

or equivalently,

$\displaystyle\ell=\ell_{\theta}-nM\ .$ (101)

The index $\ell\in\\{0,1,\cdots,M-1\\}$ is consequently out of range for every
value of $n\in\mathbb{Z}$ except $n=0$; therefore, $A_{\ell}(\theta)=0$ except
when

$\displaystyle\ell=\ell_{\theta}\equiv M\theta\ .$ (102)

The value of the amplitude for $\ell=\ell_{\theta}$ must of course be unity
(up to an arbitrary phase), and indeed it is. As previously mentioned, when
$\ell=\ell_{\theta}$ in (98), we obtain the indeterminate form $0/0$. We
therefore replace $\ell_{\theta}-\ell$ by a small displacement $\varepsilon$
and the take the limit $\varepsilon\to 0$, in which case the associated
amplitude becomes

$\displaystyle A=\lim_{\varepsilon\to 0}\frac{1}{M}\,\frac{1-e^{2\pi
i\,\varepsilon}}{1-e^{2\pi i\,\varepsilon/M}}=1\ ,$ (103)

as expected. This establishes expression (99). In this case, the final state
(93) becomes the Eigenstate

$\displaystyle|\psi_{3}\rangle=|\ell_{\theta}\rangle\otimes|u\rangle\
,~{}~{}\text{where}~{}\ell_{\theta}=M\theta\in\\{0,1,\cdots,M-1\\}\ ,$ (104)

as we obtained in the previous section.

###  III.4 Generalized Input States

Suppose now that we choose the work state $|u\rangle$ to be a linear sum over
all Eigenstates of $U$. Then the output $|\psi_{2}\rangle$ in equation (89)
will become a corresponding sum over these Eigenstates. That is to say, for
the Eigenstates

$\displaystyle U|u_{s}\rangle=e^{2\pi i\,\theta_{s}}\,|u_{s}\rangle\ ,$ (105)

let us populate the work register with the linear combination of states

$\displaystyle|u\rangle={\sum}_{s}a_{s}\,|u_{s}\rangle\ .$ (106)

We now see that the output of the front-end (89) becomes

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\,\sum_{k=0}^{M-1}{\sum}_{s}\,a_{s}\,e^{2\pi
ik\,\theta_{s}}\,|k\rangle\otimes|u_{s}\rangle\ ,$ (107)

and that the control and work registers are now entangled. The output of QPE
circuit is given by $|\psi_{3}\rangle=QFT^{\dagger}|\psi_{2}\rangle$, and
using the form of $QFT^{\dagger}$ in (90), we find

$\displaystyle|\psi_{3}\rangle$ $\displaystyle=$ $\displaystyle
QFT^{\dagger}|\psi_{2}\rangle=\sum_{\ell=0}^{M-1}{\sum}_{s}A_{\ell}(\theta_{s})\,|\ell\rangle\otimes|u_{s}\rangle\
,$ (108)

where the amplitudes $A_{\ell}(\theta_{s})$ are defined by

$\displaystyle A_{\ell}(\theta_{s})$ $\displaystyle\equiv$
$\displaystyle\frac{a_{s}}{M}\,\frac{1-e^{2\pi
i\,(\theta_{s}-\ell/M)M}}{1-e^{2\pi i\,(\theta_{s}-\ell/M)}}$ (109)
$\displaystyle=$ $\displaystyle\frac{a_{s}}{M}\,\frac{1-e^{2\pi
i\,(\theta_{s}-\theta_{\ell})M}}{1-e^{2\pi
i\,(\theta_{s}-\theta_{\ell})}}=\frac{a_{s}}{M}\,\frac{1-e^{2\pi
i\,(\ell_{s}-\ell)}}{1-e^{2\pi i\,(\ell_{s}-\ell)/M}}$ (110)

when $\ell_{s}\equiv M\theta_{s}$ is not an integer. In contrast, for integer
values $\ell_{s}\in\\{0,1,\cdots,M-1\\}$, then expression (109) reduces to

$\displaystyle A_{\ell}(\theta_{s})\equiv a_{s}\,\delta_{\ell,\ell_{s}}\ .$
(111)

For general values of $\theta_{s}$, we therefore find that the probability of
measuring a specific $\ell$-$s$ state is given by

$\displaystyle P_{\ell}(\theta_{s})=\big{|}A_{\ell}(\theta_{s})\big{|}^{2}$
$\displaystyle=$
$\displaystyle\frac{|a_{s}|^{2}}{M^{2}}~{}\frac{\sin^{2}\left[\pi\left(\theta_{s}-\displaystyle\frac{\ell}{M}\right)M\right]}{\sin^{2}\left[\pi\left(\theta_{s}-\displaystyle\frac{\ell}{M}\right)\right]}$
(112) $\displaystyle=$
$\displaystyle\frac{|a_{s}|^{2}}{M^{2}}~{}\frac{\sin^{2}\left[\pi\big{(}\theta_{s}-\theta_{\ell}\big{)}M\right]}{\sin^{2}\left[\pi\big{(}\theta_{s}-\theta_{\ell}\big{)}\right]}=\frac{|a_{s}|^{2}}{M^{2}}~{}\frac{\sin^{2}\left[\pi\big{(}\ell_{s}-\ell\big{)}\right]}{\sin^{2}\left[\pi\big{(}\ell_{s}-\ell\big{)}/M\right]}\
.~{}~{}~{}~{}~{}$ (113)

As before, for large values of $M=2^{m}$, the probability is sharply peaked
about the control states for which $\ell\approx\ell_{s}\equiv M\theta_{s}$,
while the amplitudes $a_{s}$ determine the most likely values of $\theta_{s}$.

When all of the angles $\theta_{s}$ become $m$-bit rational numbers, this
analysis simplifies considerably. In this case, the parameters

$\displaystyle\ell_{s}\equiv M\theta_{s}=2^{m}\theta_{s}$ (114)

are all integers, and the output of the front-end (107) can be written

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\,\sum_{k=0}^{M-1}{\sum}_{s}\,a_{s}\,e^{2\pi
i\,k\theta_{s}}\,|k\rangle\otimes|u_{s}\rangle$ (115) $\displaystyle=$
$\displaystyle{\sum}_{s}\,a_{s}\,\frac{1}{\sqrt{M}}\,\sum_{k=0}^{M-1}\,e^{2\pi
i\,k\ell_{s}/M}\,|k\rangle\otimes|u_{s}\rangle$ (116) $\displaystyle=$
$\displaystyle{\sum}_{s}\,a_{s}\,QFT\,\,|\ell_{s}\rangle\otimes|u_{s}\rangle\
,$ (117)

and therefore

$\displaystyle|\psi_{3}\rangle=QFT^{\dagger}\,|\psi_{2}\rangle={\sum}_{s}a_{s}\,|\ell_{s}\rangle\otimes|u_{s}\rangle\
.$ (118)

A measurement of the system now yields one of the states
$|\ell_{s}\rangle\otimes|u_{s}\rangle$ with probability $P_{s}=|a_{s}|^{2}$.
All other states have vanishing probability! This result can also be obtained
directly from the general final state (108). In this case, when all phase
angles $\theta_{s}$ becomes $m$-bit fractions, then the amplitudes (109)
reduce to $A_{\ell}(\theta_{s})=a_{s}\,\delta_{\ell,\ell_{s}}$, so that only
the terms for which $\ell=\ell_{s}\equiv M\theta_{s}$ for some value of $s$
will contribute, and the general expression (108) collapses to (118).

##  IV Continued Fractions

In an effort to render this manuscript self-contained, in this section we take
a mathematical digression to briefly introduce the theory of continued
fractions, a topic with which many readers might not be entirely familiar. In
Section V.4, we will employ continued fractions in the post-quantum processing
stage of Shor’s algorithm to extract the exact phase $\phi_{s}=s/r$ of the
modular exponentiation operator $U_{a,{\scriptscriptstyle N}}$ from the
measured (and approximate) $m$-bit phase value $\tilde{\phi}$. As we have
emphasized, this will provide the requisite period $r$ from which the factors
of $N$ can be inferred. Continued fractions, however, are interesting in their
own right, and they provide a number of fascinating connections between the
integers and the real numbers. In classical mathematics, for example,
continued fractions were employed to find rational approximations to many
irrational numbers of interest (something quite useful before the advent of
modern computers and calculators). We introduce the subject by considering the
following infinite continued fraction

$\displaystyle
x=1+\displaystyle\frac{1}{1+\displaystyle\frac{1}{1+\displaystyle\frac{1}{1+\displaystyle\frac{1}{1+\,\ddots}}}}\
,$ (119)

where the fraction $x$ telescopes downward without end. Note that the
denominator of the fractional piece after the initial 1 takes the same form as
the continued fraction itself, so that we can express (119) by the equation

$\displaystyle x=1+\frac{1}{x}\ .$ (120)

At the risk of introducing a spurious (negative) solution, we multiply (120)
by $x$ to obtain the quadratic equation $x^{2}=x+1$. The positive solution to
this equation is $x=(1+\sqrt{5})/2$, which we recognize as the golden mean!
Continued fractions are interesting indeed. In fact, any real number can be
expressed as a continued fraction consisting of a sequence of well-chosen
integers using a simple and easily executed algorithm.

We shall concentrate on continued fractions that terminate after a finite
number of iterations, thereby producing a rational number. We define a finite
continued fraction as a number of the form

$\displaystyle
x=a_{0}+\displaystyle\frac{1}{a_{1}+\displaystyle\frac{1}{a_{2}+\displaystyle\frac{1}{a_{3}+\ddots+\frac{1}{a_{n}}}}}\
,$ (121)

where $a_{0}$ is an integer (positive or negative) and
$a_{1},a_{2},\cdots,a_{n}$ are all positive integers. We will denote a
continued fraction by enumerating its integer coefficients in square brackets,
so that $x=[a_{0};a_{1},a_{2},\cdots,a_{n}]$. It is traditional to offset the
first integer $a_{0}$ with a semicolon. For our purposes, the most important
attributes of continued fractions are their so-called convergents, whose
definition is formalized below.

###### Definition 1

Suppose $x=[a_{0};a_{1},a_{2},\cdots,a_{n}]$ is a continued fraction. Any
continued fraction of the form $[a_{0};a_{1},a_{2},\cdots,a_{m}]$ for $m\leq
n$ is called a convergent of the original continued fraction for $x$.

For example,

$\displaystyle[a_{0};a_{1},a_{2},a_{3}]=a_{0}+\displaystyle\frac{1}{a_{1}+\displaystyle\frac{1}{a_{2}+\displaystyle\frac{1}{a_{3}}}}$
(122)

is a convergent of the continued fraction
$x=[a_{0};a_{1},a_{2},a_{3},\cdots,a_{n}]$. In fact, we can regard the
convergents of a continued fraction of $x$ as systematically improving
rational approximations to $x$. As we have emphasized, any finite continued
fraction gives a rational number. It turns out that the converse is also true,
namely, that any rational number can be represented as a finite continued
fraction with integer coefficients. There exists a simple and efficient
(polynomial time) algorithm for determining the associated continued fraction
of a rational number. The algorithm is best explained through an example, so
let us consider the following rational approximation to $\pi$:

$\displaystyle 3.1415$ $\displaystyle=$
$\displaystyle\frac{31415}{10000}=\frac{6283}{2000}=3+\frac{283}{2000}=3+\displaystyle\frac{1}{\displaystyle\frac{~{}2000~{}}{~{}283~{}}}$
(123) $\displaystyle=$ $\displaystyle
3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{19}{283}~{}}=3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{\displaystyle\frac{283}{19}}~{}}=3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{14+\displaystyle\frac{17}{19}}~{}}$
(124) $\displaystyle=$ $\displaystyle
3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{14+\displaystyle\frac{1}{\displaystyle\frac{19}{17}}}~{}}=3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{14+\displaystyle\frac{1}{1+\displaystyle\frac{2}{17}}}~{}}$
(125) $\displaystyle=$ $\displaystyle
3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{14+\displaystyle\frac{1}{1+\displaystyle\frac{1}{\displaystyle\frac{17}{2}}}}~{}}=3+\displaystyle\frac{1}{~{}7+\displaystyle\frac{1}{14+\displaystyle\frac{1}{1+\displaystyle\frac{1}{8+\displaystyle\frac{1}{2}}}}~{}}~{}~{}~{},$
(126)

and therefore $3.1415=[3;7,14,1,8,2]$. This algorithm consists of successive
inversions of rational numbers, followed by splitting the inverted form into
an integer plus a rational piece, and continuing this process again. For
example, in line (123) the rational contribution $283/2000$ is inverted to
form the equivalent number $1/(2000/283)$, and in the first line of (124), the
denominator $2000/283>1$ is then split into its equivalent form $7+19/283$. We
then invert the rational piece to give $1/(283/19)$, and we split $283/19>1$
into an equivalent form $14+17/19$. We continue in this fashion until the
procedure terminates.

In a certain sense, continued fractions are a more natural representation of
real numbers than their decimal counterparts. This is because the continued
fraction representation of rational numbers always terminates after a finite
number of iterations, i.e. $x\in\mathbb{Q}$ iff $x=[a_{0};a_{1},\cdots,a_{n}]$
for some finite sequence of integer coefficients $a_{\ell}$. In contrast,
decimal representations of rational numbers need not be finite, e.g. the
rational number $2/3=0.666\cdots$ has an infinite number of digits, whereas
the continued fraction expansion $2/3=[0;1,2]$ has only two non-zero
coefficients.

Note that we have required the coefficients of continued fractions to be
integers. This is because the infinite continued fraction expansion of an
irrational number is then unique. Furthermore, the continued fraction
expansion of a rational number is almost unique. It turns out that there are
only two possible continued fraction expansions for any given rational number,
provided the coefficients $a_{\ell}$ are integers (we call this semi-
uniqueness). To see this, suppose that $x=[a_{0};a_{1},\cdots,a_{n}]$ with
$a_{n}>1$ is a continued fraction expansion with integer coefficients. We can
rewrite the last coefficient as $a_{n}=(a_{n}-1)+1/1$, and consequently we can
also express the rational number by the continued fraction
$x=[a_{0};a_{1},\cdots,a_{n}-1,1]$. This means that, without loss of
generality, we may take a continued fraction representation of a rational
number to have either an even or an odd number of terms, providing the
coefficients are integers (and we shall use this fact in proving Theorem 4
below). We will, however, sometimes find it convenient to generalize the
notion of a continued fraction to allow for rational (or even real)
coefficients $a_{\ell}$. But we pay a price for doing so, as we can no longer
be assured of the semi-uniqueness of the associated continued fraction.

We now prove three essential theorems concerning continued fractions. The
first two establish that the convergents of a continued fraction take the form
$p_{n}/q_{n}$, where $p_{n}$ and $q_{n}$ are special sequences of relatively
prime integers. The third theorem can be used to relate the period $r$ to the
sequence of denominators $q_{n}$, and as we shall see in the next section, it
will be essential in extracting the period from the measured phase. Our
primary references for this section are Refs. [8] and [9].

###### Theorem 2

Let $[a_{0};a_{1},\cdots,a_{m}]$ be a continued fraction, where the
coefficients $a_{\ell}$ can be either rational numbers or integers. The
convergents $x_{n}\equiv[a_{0};a_{1},\cdots,a_{n}]$ for $n\leq m$ are equal to
the ratio $x_{n}=p_{n}/q_{n}$, where $p_{n}$ and $q_{n}$ are defined through
the sequence

$\displaystyle p_{n}$ $\displaystyle=$ $\displaystyle a_{n}p_{n-1}+p_{n-2}$
$\displaystyle q_{n}$ $\displaystyle=$ $\displaystyle a_{n}q_{n-1}+q_{n-2}$
(127)

for $2\leq n\leq m$, with the seed values

$\displaystyle p_{0}$ $\displaystyle=$ $\displaystyle a_{0}\hskip
42.67912ptq_{0}=1$ (128) $\displaystyle p_{1}$ $\displaystyle=$ $\displaystyle
a_{1}a_{0}+1~{}~{}~{}q_{1}=a_{1}\ .$ (129)

Furthermore, if the coefficients $a_{\ell}$ are positive integers, then
$p_{n}$ and $q_{n}$ are also positive integers (and they strictly increase in
magnitude).

PROOF. The proof will be through induction on $n$. For $n=0$, the convergent
is just $x_{0}\equiv[a_{0}]$, which corresponds to $x_{0}=p_{0}/q_{0}$ for
$p_{0}=a_{0}$ and $q_{0}=1$, thereby validating (128). For $n=1$, we have the
convergent

$\displaystyle
x_{1}\equiv[a_{0};a_{1}]=a_{0}+\frac{1}{a_{1}}=\frac{a_{1}a_{0}+1}{a_{1}}\ ,$
(130)

so that $x_{1}=p_{1}/q_{1}$ with $p_{1}=a_{1}a_{0}+1$ and $q_{1}=a_{1}$,
thereby validating (129). This takes care of the initial seeding. To provide a
bit of intuition, let us explicitly verify the $n=2$ case for
$x_{2}\equiv[a_{0};a_{1},a_{2}]$. From (127) we have

$\displaystyle p_{2}$ $\displaystyle=$ $\displaystyle
a_{2}p_{1}+p_{0}=a_{2}a_{1}a_{0}+a_{2}+a_{0}$ (131) $\displaystyle q_{2}$
$\displaystyle=$ $\displaystyle a_{2}q_{1}+q_{0}=a_{2}a_{1}+1\ ,$ (132)

and we can therefore express the convergent $x_{2}$ in the form:

$\displaystyle x_{2}$ $\displaystyle\equiv$
$\displaystyle[a_{0};a_{1},a_{2}]=a_{0}+\displaystyle\frac{1}{a_{1}+\displaystyle\frac{1}{a_{2}}}$
(133) $\displaystyle=$ $\displaystyle
a_{0}+\frac{a_{2}}{a_{2}a_{1}+1}=\frac{a_{2}a_{1}a_{0}+a_{0}+a_{2}}{a_{2}a_{1}+1}=\frac{p_{2}}{q_{2}}\
,$ (134)

thereby validating the theorem for $n=2$. We now assume (127) holds for some
$n\geq 3$ with $x_{n}=[a_{0};a_{1},a_{2},\cdots,a_{n}]$, and we wish to prove
that it continues to hold for $n+1$. Note that any convergent
$x_{n+1}=[a_{0};a_{1},\cdots,a_{n-1},a_{n},a_{n+1}]$ (which has $n+1$
coefficients) may be expressed in the alternative form

$\displaystyle
x_{n+1}=\underbrace{[a_{0};a_{1},\cdots,a_{n-1},a_{n}+1/a_{n+1}]}_{n~{}\text{coefficients}}\
,$ (135)

which contains only $n$ coefficients, albeit rational coefficients. We can
therefore apply the induction hypothesis to (135). To this end, let
$\tilde{p}_{\ell}/\tilde{q}_{\ell}$ be the sequence of convergents associated
with the second form of the continued fraction for $x_{n+1}$. The induction
hypothesis now gives

$\displaystyle
x_{n+1}=[a_{0};a_{1},\cdots,a_{n-1}\,,\,\underbrace{a_{n}+1/a_{n+1}}_{\tilde{a}_{n}}]=\frac{\tilde{p}_{n}}{\tilde{q}_{n}}\
,$ (136)

where

$\displaystyle\frac{\tilde{p}_{n}}{\tilde{q}_{n}}=\frac{\tilde{a}_{n}\,\tilde{p}_{n-1}+\tilde{p}_{n-2}}{\tilde{a}_{n}\,\tilde{q}_{n-1}+\tilde{q}_{n-2}}\
.$ (137)

It is clear that $\tilde{p}_{n-2}=p_{n-2}$, $\tilde{p}_{n-1}=p_{n-1}$ and
$\tilde{q}_{n-2}=q_{n-2}$, $\tilde{q}_{n-1}=q_{n-1}$, and we therefore find

$\displaystyle x_{n+1}$ $\displaystyle=$
$\displaystyle\frac{\displaystyle{\Big{(}a_{n}+\frac{1}{a_{n+1}}\Big{)}}p_{n-1}+p_{n-2}}{\displaystyle{\Big{(}a_{n}+\frac{1}{a_{n+1}}\Big{)}}q_{n-1}+q_{n-2}}=\frac{(a_{n}a_{n+1}+1)p_{n-1}+p_{n-2}a_{n+1}}{(a_{n}a_{n+1}+1)p_{n-1}+p_{n-2}a_{n+1}}$
(138) $\displaystyle=$
$\displaystyle\frac{a_{n+1}(a_{n}p_{n-1}+p_{n-2})+p_{n-1}}{a_{n+1}(a_{n}q_{n-1}+q_{n-2})+q_{n-1}}~{}~{}~{}\Leftarrow~{}~{}\text{use}~{}(\ref{eq_pm_qm})$
(139) $\displaystyle=$
$\displaystyle\frac{a_{n+1}p_{n}+p_{n-1}}{a_{n+1}q_{n}+q_{n-1}}\equiv\frac{p_{n+1}}{q_{n+1}}\
.$ (140)

Thus, the theorem is true for $n+1$. It is obvious that if the coefficients
$a_{\ell}$ are positive integers, then $p_{n}$ and $q_{n}$ are as well. This
completes the proof. $\blacksquare$

###### Theorem 3

If the coefficients of the continued fraction
$[a_{0};a_{1},a_{2},\cdots,a_{m}]$ are integers, then the integers $p_{n}$ and
$q_{n}$ of Theorem 2 are relatively prime, and satisfy the relation

$\displaystyle q_{n}p_{n-1}-p_{n}q_{n-1}=(-1)^{n}$ (141)

for $n\geq 1$.

PROOF. We prove (141) by induction on $n$. From (128) and (129) we find

$\displaystyle n=1:~{}~{}q_{1}p_{0}-p_{1}q_{0}$ $\displaystyle=$
$\displaystyle a_{1}a_{0}-(a_{1}a_{0}+1)\cdot 1=-1=(-1)^{1}\ ,$ (142)

so that (141) holds for $n=1$. Similarly, expressions (131) and (132) imply
that

$\displaystyle n=2:~{}~{}q_{2}p_{1}-p_{2}q_{1}$ $\displaystyle=$
$\displaystyle(a_{2}a_{1}+1)(a_{1}a_{0}+1)-(a_{2}a_{1}a_{0}+a_{2}+a_{0})a_{1}=1=(-1)^{2}\
,~{}~{}~{}~{}~{}$ (143)

so that (141) also holds for $n=2$. Let us now assume that (141) holds for
some $n\geq 3$, and let us prove that it continues to hold for $n+1$. Taking
$n\to n+1$ in (127), the induction hypothesis gives

$\displaystyle n+1:~{}~{}q_{n+1}p_{n}-p_{n+1}q_{n}$ $\displaystyle=$
$\displaystyle(a_{n+1}q_{n}+q_{n-1})p_{n}-(a_{n+1}p_{n}+p_{n-1})q_{n}$ (144)
$\displaystyle=$
$\displaystyle-\big{(}q_{n}p_{n-1}-p_{n}q_{n-1}\big{)}=-(-1)^{n}=(-1)^{n+1}\
.$

Therefore, (127) holds for $n+1$, and this completes the proof of the first
part of the theorem.

We must now show that $p_{n}$ and $q_{n}$ have no common factors other than
unity. Let us therefore assume that $k_{n}\geq 1$ is a common factor of
$p_{n}$ and $q_{n}$, so that

$\displaystyle
p_{n}=k_{n}\tilde{p}_{n}~{}~{}\text{and}~{}~{}q_{n}=k_{n}\tilde{q}_{n}~{}~{}\text{with}~{}~{}\frac{p_{n}}{q_{n}}=\frac{\tilde{p}_{n}}{\tilde{q}_{n}}\
.$ (145)

Let us also assume that $k_{n-1}\geq 1$ is a common factor of $p_{n-1}$ and
$q_{n-1}$, so that

$\displaystyle
p_{n-1}=k_{n-1}\tilde{p}_{n-1}~{}~{}\text{and}~{}~{}q_{n-1}=k_{n-1}\tilde{q}_{n-1}~{}~{}\text{with}~{}~{}\frac{p_{n-1}}{q_{n-1}}=\frac{\tilde{p}_{n-1}}{\tilde{q}_{n-1}}\
.$ (146)

Then (141) now takes the form

$\displaystyle
k_{n}k_{n-1}\big{(}\tilde{q}_{n}\tilde{p}_{n-1}-\tilde{p}_{n}\tilde{q}_{n-1}\big{)}=(-1)^{n}\
.$ (147)

It is obvious that $\tilde{p}_{n}/\tilde{q}_{n}$ has the same continued
fraction representation as $p_{n}/q_{n}$ for all $n\geq 1$, so that
$\tilde{q}_{n}\tilde{p}_{n-1}-\tilde{p}_{n}\tilde{q}_{n-1}=(-1)^{n}$, from
which it follows that

$\displaystyle k_{n}k_{n-1}=1\ .$ (148)

For $n=1$, we have $k_{1}k_{0}=1$. Since $1=q_{0}=k_{0}\,\tilde{q}_{0}$, and
since $k_{0}$ and $q_{0}$ are both integers, we must have $k_{0}=1$ (and
$\tilde{q}_{0}=1$). Therefore, $k_{1}=1$. For $n=2$, expression (148) becomes
$k_{2}k_{1}=1$, and therefore $k_{2}=1$. Continuing in this fashion, we find
$k_{n}=1$ for all $n\geq 1$, and hence $p_{n}$ and $q_{n}$ are relatively
prime. $\blacksquare$

We move on to our primary result, the theorem that will allow us to extract
the desired period from an approximately measured phase angle.

###### Theorem 4

Let $x$ be a rational number. If two relatively prime integers $p$ and $q$
satisfy

$\displaystyle\left|\frac{p}{q}-x\right|\leq\frac{1}{2q^{2}}\ ,$ (149)

then $p/q$ is necessarily a convergent of $x$.

PROOF. Let $p/q=[a_{0};a_{1},\cdots,a_{n}]$ be the continued fraction
representation for $p/q$, and define the convergents $p_{\ell}/q_{\ell}$ for
$\ell=0,1,\cdots,n$ as in Theorem 2, so that $p_{n}/q_{n}=p/q$. The object of
the proof is to construct the continued fraction representation for $x$, and
this will explicitly show that $p/q$ is one of its convergents.

Let us first define the error $\delta$ for $p/q=p_{n}/q_{n}$ by

$\displaystyle x-\frac{p_{n}}{q_{n}}\equiv\frac{\delta}{2q_{n}^{2}}\ .$ (150)

Note that inequality (149) gives $0\leq\delta\leq 1$. We now define the
parameter

$\displaystyle\lambda\equiv
2\,\frac{q_{n}p_{n-1}-p_{n}q_{n-1}}{\delta}-\frac{q_{n-1}}{q_{n}}\ ,$ (151)

and with some algebra we can show that equations (150) and (151) imply

$\displaystyle x=\frac{\lambda p_{n}+p_{n-1}}{\lambda q_{n}+q_{n-1}}\ .$ (152)

To see this, note that (151) allows us to write

$\displaystyle\frac{2}{\delta}\,\big{(}q_{n}p_{n-1}-p_{n}q_{n-1}\big{)}=\lambda+\frac{q_{n-1}}{q_{n}}=\frac{\lambda
q_{n}+q_{n-1}}{q_{n}}~{}~{}~{}\Rightarrow~{}~{}$ (153)
$\displaystyle\frac{\delta}{2}=\frac{q_{n}}{\lambda
q_{n}+q_{n-1}}\big{(}q_{n}p_{n-1}-p_{n}q_{n-1}\big{)}\ .$ (154)

Using this result in equation (150) gives

$\displaystyle x$ $\displaystyle=$
$\displaystyle\frac{p_{n}}{q_{n}}+\frac{\delta}{2}\,\frac{1}{q_{n}^{2}}=\frac{p_{n}}{q_{n}}+\frac{q_{n}p_{n-1}-p_{n}q_{n-1}}{q_{n}\big{(}\lambda
q_{n}+q_{n-1}\big{)}}$ (155) $\displaystyle=$
$\displaystyle\frac{p_{n}\big{(}\lambda
q_{n}+q_{n-1}\big{)}+\big{(}q_{n}p_{n-1}-p_{n}q_{n-1}\big{)}}{q_{n}\big{(}\lambda
q_{n}+q_{n-1})}$ (156) $\displaystyle=$ $\displaystyle\frac{\big{(}\lambda
p_{n}+p_{n}q_{n-1}/q_{n}\big{)}+\big{(}p_{n-1}-p_{n}q_{n-1}/q_{n}\big{)}}{\lambda
q_{n}+q_{n-1}}=\frac{\lambda p_{n}+p_{n-1}}{\lambda q_{n}+q_{n-1}}\ .$ (157)

Now that we have established (152), Theorem 2 implies $x$ is a continued
fraction with coefficients $a_{0},a_{1},\cdots,a_{n},a_{n+1}$, where
$a_{n+1}=\lambda$, i.e.

$\displaystyle
x=[a_{0};a_{1},\cdots,a_{n},a_{n+1}]=[a_{0};a_{1},\cdots,a_{n},\lambda]\ .$
(158)

Without loss of generality, we can assume that $n$ is even, so that Theorem 3
gives

$\displaystyle\lambda=\frac{2}{\delta}-\frac{q_{n-1}}{q_{n}}>2-\frac{q_{n-1}}{q_{n}}>2-1>1\
.$ (159)

Thus, $\lambda$ is a rational number greater one, and it therefore has a
continued fraction expansion of the form $\lambda=[b_{0};b_{1},\cdots,b_{m}]$
(since $\lambda>0$ we must have $b_{0}>0$, and we henceforth drop the
semicolon after the initial coefficient $b_{0}$). Therefore, we find

$\displaystyle x=[a_{0};a_{1},\cdots,a_{n},b_{0},\cdots,b_{m}]\ ,$ (160)

which shows that $x$ is a finite continued fraction with $p/q=p_{n}/q_{n}$ as
one of its convergents. $\blacksquare$

##  V Factoring with Shor’s Algorithm

###  V.1 Basic Observation

We now address an essential observation on our way to building Shor’s
algorithm. Let $N$ be the positive integer we wish to factor. We assume that
$N$ is not even, and not a power of a prime number (otherwise we can find a
factor quickly). We say that two integers $a,b\in\mathbb{Z}$ are congruent
modulo $N$ provided the difference $a-b$ is divisible by $N$, and we express
this by writing

$\displaystyle a=b~{}({\rm mod}~{}N)\ .$ (161)

That is to say, $a$ and $b$ are congruent modulo $N$ provided there exists
another integer $m\in\mathbb{Z}$ such that

$\displaystyle a-b=mN\ .$ (162)

Shor’s factoring algorithm relies on the following observation. Suppose we can
find a non-trivial or proper square root of unity modulo $N$. In other words,
suppose that we have found an integer $b$ such that

$\displaystyle b^{2}=1~{}({\rm mod}~{}N)\ .$ (163)

Then $b$ is a modular square root of unity. By a proper square root, we mean
that

$\displaystyle b\neq\pm 1~{}({\rm mod}~{}N)\ .$ (164)

Equation (163) implies that there exists an integer $m\in\mathbb{Z}$ such that

$\displaystyle b^{2}-1=mN\ .$ (165)

The latter equation is the key to the factoring algorithm, as it can be
expressed as

$\displaystyle(b+1)(b-1)=mN\ ,$ (166)

and therefore, we see that the greatest common divisors $d_{\pm}\equiv{\rm
gcd}(b\pm 1,N)$ are factors of $N$. We note that finding the greatest common
divisor of two integers can be performed very quickly (in polynomial time) on
a classical computer. It can be shown that (163) and (164) indeed lead to non-
trivial or proper factors, in that $d_{\pm}\neq 1,N$. This is formalized in
the following theorem.

###### Theorem 5

Let $b\in\mathbb{Z}$ be a proper square root of unity modulo $N$. That is to
say, let  $b^{2}=1~{}({\rm mod}~{}N)$ and $b\neq\pm 1~{}({\rm mod}~{}N)$. Then
${\rm gcd}(b+1,N)$ and ${\rm gcd}(b-1,N)$ are proper factors of $N$.

PROOF. First consider $d={\rm gcd}(b-1,N)$, which is indeed a factor of $N$.
We will show that $d\neq 1$ and $d\neq N$. The proof will be by contradiction.

\- First assume $d={\rm gcd}(b-1,N)=N$. Then $N$ divides $b-1$, so that
$b-1=mN$ for some $m\in\mathbb{Z}$, or equivalently, $b=1~{}({\rm mod}~{}N)$.
This contradicts the fact that $b$ is a proper root of unity.

\- Now assume $d={\rm gcd}(b-1,N)=1$. Since $b-1$ and $N$ are relatively
prime, there exists integers $u,v\in\mathbb{Z}$ such that

$\displaystyle(b-1)u+Nv=1\ .$ (167)

Multiplying both sides by $b+1$ gives the expression

$\displaystyle b+1=(b^{2}-1)u+(b+1)Nv\ .$ (168)

Let us divide both sides of this equation by $N$, and employ equation (165)
for the first term $(b^{2}-1)u$ on the right-hand side of (168):

$\displaystyle\frac{b+1}{N}=\underbrace{~{}\frac{b^{2}-1}{N}~{}}_{\text{integer}~{}\equiv~{}m}\\!\cdot\,u+(b+1)\cdot
v=m\cdot u+(b+1)\cdot v\in\mathbb{Z}\ ,$ (169)

where $u,v,m,b$ are all integers. Thus, $N$ divides $b+1$, so that
$b=-1~{}({\rm mod}~{}N)$. Again, this contradicts the fact that $b$ is a
proper root of unity. Similar reasoning holds for $d={\rm gcd}(b+1,N)$.
$\blacksquare$

###  V.2 Period Finding and Factorization

Consider two integers $a$ and $N$ such that $a<N$. The integer $N$ is the
number we wish to factor, while $a$ is an initial “guess” for one of the
factors. We will usually refer to $a$ as the base. In fact, we can randomly
choose the base from $\\{2,3,\cdots,N-1\\}$, provided that $a$ and $N$ are
relatively prime, so that ${\rm gcd}(a,N)=1$ (otherwise we have found a non-
trivial factor of $N$). Let us now define the order of $a$ modulo $N$ as the
least positive integer $r$ such that

$\displaystyle a^{r}=1~{}({\rm mod}~{}N)\ .$ (170)

Therefore, if $r$ is even, then $b=a^{r/2}$ is an integer, and it is a square
root of unity. If it is also non-trivial, we can perform the factoring
algorithm outlined above based on this value of $b$. If the order $r$ is odd
and $a$ is not a perfect square, then we must try a new base $a$. On the other
hand, if $r$ is odd and $a$ is a perfect square, then $b=a^{r/2}$ is still an
integer square root of unity, and it can be used in the algorithm [12]. For
these cases of $r$ and $a$, we see that $b=a^{r/2}$ is a square root of unity,
and provided that it is non-trivial, then the factors of $N$ are given by
${\rm gcd}(a^{r/2}\pm 1,N)$.

There is an equivalent way of looking at this that employs the periodic
modular exponential function defined by

$\displaystyle f_{a\,{\scriptscriptstyle N}}(x)=a^{x}~{}({\rm mod}~{}N)\ ,$
(171)

Figure 11: The function $f(x)=a^{x}~{}({\rm mod}~{}N)$ for $N=15$ and the
bases $a\in\\{4,8,11,2\\}$.

where $a,x$, and $N$ are non-negative integers. We will usually drop the
subscripts and simply write $f(x)$. The modular order $r$ is nothing more than
the period of $f(x)$. To see this, note that $f(0)=1$, and since $r$ is the
smallest integer such that $f(r)=a^{r}~{}({\rm mod}~{}N)=1$, we see that
$f(r)=1$. In fact, for any argument $x$, we have $f(x+r)=f(x)$, and thus the
order $r$ is the period of $f(x)$. As an example, let us take $N=15$. The base
$a$ and the number $N$ cannot have any non-trivial common factors, and the
base must satisfy $1<a<N$, which limits the allowed values to
$a\in\\{2,4,7,8,11,13,14\\}$. Figures 11 and 12 illustrate the functions
$f(x)$ for $N=15$ for these values of $a$ (note that $a=14$ gives a trivial
square root with no factors, so we do not bother to provide a plot). We
summarize below the factorization algorithm from the last section based on the
periods $r$ for the cases specified in the Figures. For the bases $a$ in Fig.
11 we find:

$\bullet~{}a=4\Rightarrow r=2$: $\bullet~{}a=2\Rightarrow r=4$:
\- $a^{r/2}-1=4^{1}-1=3~{}\Rightarrow~{}{\rm gcd}(3,15)=3$ \-
$a^{r/2}-1=2^{2}-1=4-1=3~{}\Rightarrow~{}{\rm gcd}(3,15)=3$
\- $a^{r/2}+1=4^{1}+1=5~{}\Rightarrow~{}{\rm gcd}(5,15)=5$ \-
$a^{r/2}+1=2^{2}+1=4+1=5~{}\Rightarrow~{}{\rm gcd}(5,15)=5$

$\bullet~{}a=11\Rightarrow r=2$:
\- $a^{r/2}-1=11^{1}-1=10~{}\Rightarrow~{}{\rm gcd}(10,15)=5$
\- $a^{r/2}+1=11^{1}+1=12~{}\Rightarrow~{}{\rm gcd}(12,15)=3$

$\bullet~{}a=8\Rightarrow r=4$: $\bullet~{}a=14\Rightarrow r=2$:
\- $a^{r/2}-1=8^{2}-1=64-1=63~{}\Rightarrow~{}{\rm gcd}(63,15)=3$ \-
$a^{r/2}-1=13$
\- $a^{r/2}+1=8^{2}+1=64+1=65~{}\Rightarrow~{}{\rm gcd}(65,15)=5$ \-
$a^{r/2}+1=15\Leftarrow$ trivial

Figure 12: The function $f(x)=a^{x}~{}({\rm mod}~{}N)$ for $N=15$ and the
bases $a\in\\{7,13\\}$.

And for the bases $a$ in Fig. 12 we find:
$\bullet~{}a=7\Rightarrow r=4$:
\- $a^{r/2}-1=7^{2}-1=49-1=48~{}\Rightarrow~{}{\rm gcd}(48,15)=3$
\- $a^{r/2}+1=7^{2}+1=49+1=50~{}\Rightarrow~{}{\rm gcd}(50,15)=5$

$\bullet~{}a=13\Rightarrow r=4$:
\- $a^{r/2}-1=13^{2}-1=169-1=168~{}\Rightarrow~{}{\rm gcd}(168,15)=3$
\- $a^{r/2}+1=13^{2}+1=169+1=170~{}\Rightarrow~{}{\rm gcd}(170,15)=5$

###  V.3 The Factorization Circuit

Our primary references for this section are Refs. [8, 9, 11]. We will show
that one can find the period $r$ of the modular exponential function
$f_{a\,{\scriptscriptstyle N}}(x)=a^{x}~{}({\rm mod}~{}N)$ by exploiting the
quantum phase estimation (QPE) algorithm developed in Section III. The crucial
step in this procedure is to define a unitary operator
$U_{a\,{\scriptscriptstyle N}}$ whose Eigen-phases contain information about
the period $r$. We then employ the QPE algorithm to find the phases of
$U_{a\,{\scriptscriptstyle N}}$, thereby permitting us to determine the exact
value of $r$. Recall that the QPE algorithm consists of two quantum registers,
a control register and a work register. The control register contains $m$
qubits that dictate the resolution of the measured output phase of
$U_{a\,{\scriptscriptstyle N}}$, while the work register encodes information
about the number $N$, and we therefore take the number of work qubits to be
$n=\lceil\log_{2}N\rceil$ (the binary length of $N$). We shall define a linear
unitary operator $U_{a\,{\scriptscriptstyle N}}$ by its action on the
computational basis states of the work register,

$\displaystyle U_{a\,{\scriptscriptstyle N}}\,|w\rangle=|a\cdot w~{}({\rm
mod}~{}N)\rangle\ .$ (172)

We will usually drop the $N$ and $a$ subscripts and write
$U=U_{a\,{\scriptscriptstyle N}}$ for simplicity. We shall refer to $U$ as
either the phase operator or the modular exponentiation (ME) operator. To
continue, let us now solve the Eigenvalue problem for the ME operator $U$.
This is accomplished by the simple observation that

$\displaystyle U\,|a^{x}~{}({\rm mod}~{}N)\rangle$ $\displaystyle=$
$\displaystyle|a^{x+1}~{}({\rm mod}~{}N)\rangle$ (173)

for any non-negative integer $x$. Let us now define the $r$ states

$\displaystyle|u_{s}\rangle=\frac{1}{\sqrt{r}}\sum_{k=0}^{r-1}\,e^{-2\pi
ik\,s/r}\,|a^{k}~{}({\rm
mod}~{}N)\rangle~{}~{}~{}\text{for}~{}~{}s\in\\{0,1,\cdots,r-1\\}\ ,$ (174)

from which equation (173) gives

$\displaystyle U|u_{s}\rangle=e^{2\pi
i\,\phi_{s}}\,|u_{s}\rangle~{}~{}~{}\text{with}~{}~{}\phi_{s}=\frac{s}{r}\ .$
(175)

The possible phases of the ME operator $U$ are therefore $\phi_{s}=s/r$ for
$s\in\\{0,1,\cdots,r-1\\}$, where $r$ is the period of the function
$f(x)=a^{x}~{}({\rm mod}~{}N)$. To prove this result, it is instructive to
expand the states in (174) term-by-term,

$\displaystyle|u_{s}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{r}}\,\Big{[}\,|a^{0}~{}({\rm
mod}~{}N)\rangle+e^{-2\pi i\,s/r}\,|a^{1}~{}({\rm mod}~{}N)\rangle+e^{-2\pi
i\,2s/r}\,|a^{2}~{}({\rm mod}~{}N)\rangle+\cdots+~{}~{}~{}$ (176)
$\displaystyle\hskip 28.45274pte^{-2\pi i\,(r-2)s/r}\,|a^{r-2}~{}({\rm
mod}~{}N)\rangle+e^{-2\pi i\,(r-1)s/r}\,|a^{r-1}~{}({\rm
mod}~{}N)\rangle\Big{]}\ .$

The series terminates after $r$ terms because $a^{r}~{}({\rm mod}~{}N)=1$,
which leads us back to the first term $|1\rangle=|a^{0}~{}({\rm
mod}~{}N)\rangle$. From relation (173), we can now easily prove expression
(175):

$\displaystyle U|u_{s}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{r}}\Big{[}\,|a^{1}~{}({\rm
mod}~{}N)\rangle+e^{-2\pi i\,s/r}|\,a^{2}~{}({\rm mod}~{}N)\rangle+e^{-2\pi
i\,2s/r}|\,a^{3}~{}({\rm mod}~{}N)\rangle+\cdots+$ (178) $\displaystyle\hskip
28.45274pte^{-2\pi i\,(r-2)s/r}|\,a^{r-1}~{}({\rm mod}~{}N)\rangle+e^{-2\pi
i\,(r-1)s/r}|\,a^{r}~{}({\rm mod}~{}N)\rangle\Big{]}$ $\displaystyle=$
$\displaystyle e^{2\pi i\,s/r}\,\frac{1}{\sqrt{r}}\Big{[}\,e^{-2\pi
i\,s/r}\,|a^{1}~{}({\rm mod}~{}N)\rangle+e^{-2\pi i\,2s/r}\,|a^{2}~{}({\rm
mod}~{}N)\rangle+e^{-2\pi i\,3s/r}\,|a^{3}~{}({\rm mod}~{}N)\rangle$
$\displaystyle\hskip 28.45274pt+\cdots+e^{-2\pi i\,(r-1)s/r}|\,a^{r-1}~{}({\rm
mod}~{}N)\rangle\,+\,\underbrace{\,~{}e^{-2\pi
i\,s}~{}\,}_{1}\,|\,a^{0}~{}({\rm mod}~{}N)\rangle\Big{]}$ $\displaystyle=$
$\displaystyle e^{2\pi i\,s/r}\,|u_{s}\rangle\ .$ (179)

Also note that the phases $e^{2\pi ik\,s/r}$ sum to zero over $s$ for any
value of the non-zero integers $r$ and $k$. This is easy to prove, as the sum
is geometric and can be performed exactly:

$\displaystyle\sum_{s=0}^{r-1}e^{2\pi i\,ks/r}$ $\displaystyle=$
$\displaystyle\sum_{s=0}^{r-1}\big{[}e^{2\pi
i\,k/r}\big{]}^{s}=\frac{1-\big{[}e^{2\pi i\,k/r}\big{]}^{r}}{1-e^{2\pi
i\,k/r}}=\frac{1-e^{2\pi i\,k}}{1-e^{2\pi i\,k/r}}=0$ (180)

for any integer $k\neq 0$. The sum over $s$ in (176) therefore gives

$\displaystyle\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}|u_{s}\rangle=|1\rangle\ ,$
(181)

which will prove to be a key ingredient for Shor’s algorithm. Note that only
the first term $|a^{0}~{}({\rm mod}~{}N)\rangle=|1\rangle$ contributes to the
sum, as all other terms have phases that add to zero. In fact, if we multiply
(174) or (176) by $e^{2\pi i\,ks/r}$, we remove the phase of the term
$|a^{k}~{}({\rm mod}~{}N)\rangle$, and upon summing over $s$ we find the
generalized result

$\displaystyle\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}e^{2\pi
ik\,\phi_{s}}\,|u_{s}\rangle=|a^{k}~{}({\rm
mod}~{}N)\rangle=|f(k)\rangle~{}~{}~{}\text{where}~{}~{}\phi_{s}=s/r\ .$ (182)

Figure 13: First attempt at Shor’s factoring circuit for the physics and
mathematics Convention 2. The control register has $m$ qubits, and the work
register has $n=\lceil{\rm log}_{2}N\rceil$ qubits. We populate the work
register with one of the Eigenstates $|u_{s}\rangle$, where
$U|u_{s}\rangle=e^{2\pi i\phi_{s}}|u_{s}\rangle$ with phase $\phi_{s}=s/r$ for
some $s\in\\{0,1,\cdots,r-1\\}$. Measuring the control register projects the
wavefunction into an Eigenstate $|\tilde{\ell}\,\rangle\otimes|u_{s}\rangle$,
where the state is indexed by the $m$-bit integer
$\tilde{\ell}=\tilde{\phi}_{1}\cdots\tilde{\phi}_{m}$, from which we obtain
the measured phase
$\tilde{\phi}_{\ell}=\tilde{\ell}/2^{m}=0.\tilde{\phi}_{1}\dots\tilde{\phi}_{m}$
to $m$ bits of accuracy. We expect $\tilde{\phi}_{\ell}\approx\phi_{s}$,
thereby allowing us to determine $r$. The only problem with this reasoning is
that we do not know the state $|u_{s}\rangle$ in advance since we do not a
priori know the period $r$.

The central observation here is that every phase $\phi_{s}=s/r$, except
$\phi_{0}=0$, contains the period $r$. This is the basis of Shor’s algorithm:
by measuring a phase $\phi_{s}=s/r$ of the ME operator $U$ for which $s$ and
$r$ are relatively prime, we can infer the period $r$ of the function $f(x)$.
And from the corresponding factoring procedure outlined above, we can then
find non-trivial factors of $N$.

We can now address the quantum circuit for Shor’s algorithm, the first attempt
of which is illustrated in Fig. 13. We employ a QPE algorithm with the phase
operator $U$ as defined in (172). We emphasize that the QPE circuit consists
of a control register of $m$ qubits and a work register of $n$ qubits. As
stated above, the work register stores information about the number $N$.
Suppose that the work register is populated by a specific Eigenstate
$|u_{s}\rangle$ with phase $\phi_{s}=s/r$, so that $n=\lceil{\rm
log}_{2}N\rceil$. Note that there are only $r$ phase states $|u_{s}\rangle$
out of a possible $2^{n}$ work states, so that the states $|u_{s}\rangle$ are
very sparse indeed. Upon measuring the final state of the control register,
the wavefunction collapses into an Eigenstate
$|\tilde{\ell}\,\rangle\otimes|u_{s}\rangle$ for
$\tilde{\ell}\in\\{0,1,\cdots,M-1\\}$. The control register index takes the
measured value $\tilde{\ell}=\tilde{\phi}_{1}\cdots\tilde{\phi}_{m}$, where
each $\tilde{\phi}_{k}\in\\{0,1\\}$ is the measured outcome of qubit
$k\in\\{1,2,\cdots,m\\}$ of the control register. From this we can readily
infer the measured phase to be
$\tilde{\phi}_{\ell}=\tilde{\ell}/2^{m}=0.\tilde{\phi}_{1}\dots\tilde{\phi}_{m}$.
From here on, we will place a tilde over measured quantities. If the $m$-bit
resolution is sufficiently large, then we expect
$\tilde{\phi}_{\ell}\approx\phi_{s}$, from which we can extract the value of
$r$ (provided $s$ and $r$ are relatively prime). In fact, as we will establish
in the next section, we should choose $m=2n+1$ to ensure sufficient phase
resolution.

By populating the work register with $|u_{s}\rangle$, Fig. 13 uses the QPE
algorithm to find the phase $\phi_{s}=s/r$. The problem with this method,
however, is that we do not know the state $|u_{s}\rangle$ in advance, as it
depends upon the unknown period $r$. In fact, if we knew the Eigenstates
$|u_{s}\rangle$, then we should also know the value of the period $r$, and we
would have no need for Shor’s algorithm. We can circumvent this difficulty by
employing (181), which we repeat here for emphasis:

$\displaystyle|1\rangle=\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}|u_{s}\rangle\ .$
(183)

Figure 14: Shor’s factoring circuit. Convention 2. We populate the work
register with the state $|1\rangle$, which is just a uniform linear
combination of the Eigenstates $|u_{s}\rangle$. This requires no prior
knowledge of these Eigenstates. Through phase kickback, the control register
becomes a linear combination of terms involving the states $|u_{s}\rangle$ and
their phases $\phi_{s}$. Measuring the control register thereby projects the
system into a state $|\tilde{\ell}\,\rangle\otimes|\tilde{u}_{s}\rangle$ for
some $\tilde{\ell}\in\\{0,1,\cdots,M-1\\}$ and $s\in\\{0,1,\cdots,r-1\\}$. The
control state is indexed by the $m$-bit integer
$\tilde{\ell}=\tilde{\phi}_{1}\cdots\tilde{\phi}_{m}$, where each
$\tilde{\phi}_{k}\in\\{0,1\\}$ is the measured outcome of qubit
$k\in\\{1,2,\cdots,m\\}$ of the control register. We then evaluate the $m$-bit
measured phase
$\tilde{\phi}_{\ell}=\tilde{\ell}/2^{m}=0.\tilde{\phi}_{1}\cdots\tilde{\phi}_{m}$,
from which we can extract the exact phase $\phi_{s}=s/r$ using the method of
continued fractions (provided $s$ and $r$ are relatively prime).

This suggests that we populate the work register with the easily prepared
state $|1\rangle=|0\cdots 01\rangle$, which is just a uniform linear
combination of the phase Eigenstates $|u_{s}\rangle$. This is achieved by
initializing the lowest order qubit in the work register to the 1-state, and
all other work qubits to 0-states. This new circuit is illustrated in Fig. 14.
Populating the work register with a linear superposition of the states
$|u_{s}\rangle$ has the effect of rendering the state of the control register
(through phase kickback) as a linear combination of the Eigenstates
$|u_{s}\rangle$ and their corresponding phases $\phi_{s}$. More specifically,
the state of the quantum system right after the front-end of the circuit
(position 2 in the Figure) will be given by

$\displaystyle|\psi_{2}\rangle=\frac{1}{\sqrt{rM}}\sum_{k=0}^{M-1}\sum_{s=0}^{r-1}e^{2\pi
ik\,\phi_{s}}\,|k\rangle\otimes|u_{s}\rangle\ .$ (184)

Finally, the $QFT^{\dagger}$ operation (position 3 in the Figure) transforms
this state into

$\displaystyle|\psi_{3}\rangle=\sum_{\ell=0}^{M-1}\sum_{s=0}^{r-1}A_{\ell}(\phi_{s})\,|\ell\rangle\otimes|u_{s}\rangle~{}~{}~{}\text{with}~{}~{}\phi_{s}=s/r\
,$ (185)

where the amplitudes are defined by

$\displaystyle A_{\ell}(\phi_{s})$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{r}M}\,\frac{1-e^{2\pi
i\,(\phi_{s}-\ell/M)M}}{1-e^{2\pi i\,(\phi_{s}-\ell/M)}}$ (186)
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{r}M}\,\frac{1-e^{2\pi
i\,(\phi_{s}-\phi_{\ell})M}}{1-e^{2\pi
i\,(\phi_{s}-\phi_{\ell})}}=\frac{1}{\sqrt{r}M}\,\frac{1-e^{2\pi
i\,(\ell_{s}-\ell)}}{1-e^{2\pi i\,(\ell_{s}-\ell)}/M}\ ,$ (187)

where $\phi_{\ell}\equiv\ell/M$ and $\ell_{s}\equiv M\phi_{s}$. The control
register has therefore obtained knowledge of the phases $\phi_{s}=s/r$ and
their corresponding Eigenstates $|u_{s}\rangle$.

Finally, we must measure the control register (position 4 in the Figure with
the bold red line), which collapses the quantum state $|\psi_{3}\rangle$ of
(185) into a definite Eigenstate state

$\displaystyle|\tilde{\ell}\,\rangle\otimes|u_{\tilde{s}}\rangle\ .$ (188)

In other words, the control register collapses to a state
$|\tilde{\ell}\,\rangle$ labeled by the $m$-bit integer index
$\tilde{\ell}=\tilde{\phi}_{1}\,\tilde{\phi}_{2}\cdots\tilde{\phi}_{m}\in\\{0,1,\cdots,M-1\\}$,
where the measured phase is then given by
$\tilde{\phi}_{\ell}\equiv\tilde{\ell}/M=0.\tilde{\phi}_{1}\,\tilde{\phi}_{2}\cdots\tilde{\phi}_{m}$.
Likewise, the work register is projected into a state $|u_{\tilde{s}}\rangle$
determined by a random choice of the phase integer
$\tilde{s}\in\\{0,1,\cdots,r-1\\}$. The ME phase angle will be written
$\tilde{\phi}_{s}=\tilde{s}/r$, and we define the corresponding mode
$\tilde{\ell}_{s}=M\tilde{\phi}_{s}$. The probability of measuring the
Eigenstate state (188) is therefore

$\displaystyle P_{\tilde{\ell},\tilde{s}}$ $\displaystyle=$
$\displaystyle\big{|}A_{\tilde{\ell}}(\tilde{\phi}_{s})\big{|}^{2}=\frac{1}{rM^{2}}~{}\frac{\sin^{2}\left[\pi\left(\tilde{\phi}_{s}-\displaystyle\frac{\tilde{\ell}}{M}\right)M\right]}{\sin^{2}\left[\pi\left(\tilde{\phi}_{s}-\displaystyle\frac{\tilde{\ell}}{M}\right)\right]}$
(189) $\displaystyle=$
$\displaystyle\frac{1}{rM^{2}}~{}\frac{\sin^{2}\left[\pi\big{(}\tilde{\phi}_{s}-\tilde{\phi}_{\ell}\big{)}M\right]}{\sin^{2}\left[\pi\big{(}\tilde{\phi}_{s}-\tilde{\phi}_{\ell}\big{)}\right]}=\frac{1}{rM^{2}}~{}\frac{\sin^{2}\left[\pi\big{(}\tilde{\ell}_{s}-\tilde{\ell}\,\big{)}\right]}{\sin^{2}\left[\pi\big{(}\tilde{\ell}_{s}-\tilde{\ell}\,\big{)}/M\right]}\
.~{}~{}~{}~{}~{}~{}$ (190)

In applying Shor’s algorithm it will be critical that $\tilde{\phi}_{\ell}$ is
a rational number, and that $\tilde{s}$ and $r$ are relatively prime. We will
usually denote the state $|u_{\tilde{s}}\rangle$ by the simpler form
$|\tilde{u}_{s}\rangle$, in which case we will drop the tilde from the ME
phase and write $\phi_{s}=s/r$. Note that the control and work registers are
entangled, and that the most likely control-output states $\tilde{\ell}$ are
those for which $\tilde{\phi}_{\ell}\approx\phi_{s}$ for some
$s\in\\{0,1,\cdots,r-1\\}$. The rub is then to extract the exact phase
$\phi_{s}$ from the approximately measured $m$-bit phase
$\tilde{\phi}_{\ell}$, and in particular we must pull out the integer period
$r$ from the rational number $\tilde{\phi}_{\ell}$. This requires the
mathematical technique of continued fractions, which we shall apply in the
next section. If we obtain the zero-phase state $\phi_{0}=0$, or if $s$ and
$r$ have non-trivial common factors, then the measurement will fail, and we
must try again. Once we have a potential value for the period $r$, we must
explicitly check that $r$ is even, and that $a^{r}=1~{}({\rm mod}~{}N)$. We
must also ensure that $a^{r/2}\neq\pm 1~{}({\rm mod}~{}N)$. If these criteria
are not met, then we must try again. Thus, Shor’s algorithm is probabilistic
in nature, but the chance of a successful run is quite high, and usually
requires at most two or three attempts.

Let us now work through the circuit in Fig. 14 in more detail. First, we must
prepare the zero-state

$\displaystyle|0\rangle^{\otimes m}\otimes|0\rangle^{\otimes n}\ ,$ (191)

and then apply an $X$ operator to the lowest order qubit of the work register,
thereby producing the initial state

$\displaystyle|\psi_{0}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle^{\otimes m}\otimes|1\rangle\ .$ (192)

We are using a short-hand notation for the states in the $n$-qubit work
register, where we denote the computational basis states by their
corresponding integer index, rather than breaking them out into their
respective tensor products of single-particle qubit states. In the physics
Convention 2, the explicit form of the work register would be

$\displaystyle|1\rangle=|0\rangle_{1}\otimes\cdots\otimes|0\rangle_{n-1}\otimes|1\rangle_{n}\
,$ (193)

and in the Qiskit Convention 1 we would have

$\displaystyle|1\rangle=|1\rangle_{0}\otimes|0\rangle_{1}\otimes\cdots\otimes|0\rangle_{n-1}\
.$ (194)

We should also note that the phase operator $U$ is distributed over the work
register, and it takes different forms between Conventions 1 and 2.

Hadamard gates $H^{\otimes m}$ are then applied to the control register, which
splits the state $|0\rangle^{\otimes m}$ into a linear superposition of all
possible control states $|k\rangle$ for $k\in\\{0,1,\cdots,M-1\\}$ with
$M=2^{m}$. Upon using (181), we can therefore express the resulting state in
either of two forms,

$\displaystyle|\psi_{1}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes|1\rangle$
(195) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{rM}}\sum_{k=0}^{M-1}\sum_{s=0}^{r-1}\,|k\rangle\otimes|u_{s}\rangle\
,$ (196)

each of which reveals something essential about the quantum system. For
example, (196) tells us that the state $|\psi_{1}\rangle$ is in fact just a
uniform linear superposition involving the Eigenstates $|u_{s}\rangle$, which
was crucial to the above measurement analysis. Next, the circuit operates on
the state $|\psi_{1}\rangle$ with a sequence of $m$ controlled-phase operators
$CU^{p}$ for $p\in\\{2^{0},2^{1},\cdots,2^{m-1}\\}$, giving the state

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{r\,M}}\sum_{k=0}^{M-1}\sum_{s=0}^{r-1}e^{2\pi
ik\,\phi_{s}}\,|k\rangle\otimes\,|u_{s}\rangle~{}~{}~{}\text{for}~{}~{}\phi_{s}=\frac{s}{r}\
.$ (197)

By employing (182), we can express this state in the alternative form

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}\,|k\rangle\otimes\,\left(\frac{1}{\sqrt{r}}\sum_{s=0}^{r-1}e^{2\pi
i\,k\phi_{s}}|u_{s}\rangle\right)$ (198) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes|a^{k}~{}({\rm
mod}~{}N)\rangle=\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes|f(k)\rangle\
.$ (199)

Recall that the period of the function $f(x)=a^{x}~{}({\rm mod}~{}N)$ is $r$,
and we therefore find the following sequence of operations on the state
$|1\rangle$:

$\displaystyle U^{1}|1\rangle$ $\displaystyle=$ $\displaystyle|a~{}({\rm
mod}~{}N)\rangle$ $\displaystyle U^{2}|1\rangle$ $\displaystyle=$
$\displaystyle U|a~{}({\rm mod}~{}N)\rangle=|a^{2}~{}({\rm mod}~{}N)\rangle$
$\displaystyle\cdots$ $\displaystyle U^{k}|1\rangle$ $\displaystyle=$
$\displaystyle U|a^{k-1}~{}({\rm mod}~{}N)\rangle=|a^{k}~{}({\rm
mod}~{}N)\rangle$ $\displaystyle\cdots$ $\displaystyle U^{r-1}|1\rangle$
$\displaystyle=$ $\displaystyle U|a^{r-2}~{}({\rm
mod}~{}N)\rangle=|a^{r-1}~{}({\rm mod}~{}N)\rangle$ $\displaystyle
U^{r}|1\rangle$ $\displaystyle=$ $\displaystyle U|a^{r-1}~{}({\rm
mod}~{}N)\rangle=|a^{r}~{}({\rm mod}~{}N)\rangle=|1\rangle$ $\displaystyle
U^{r+1}|1\rangle$ $\displaystyle=$ $\displaystyle U|a^{r}~{}({\rm
mod}~{}N)\rangle=|a^{r+1}~{}({\rm mod}~{}N)\rangle=|a~{}({\rm
mod}~{}N)\rangle$ $\displaystyle U^{r+2}|1\rangle$ $\displaystyle=$
$\displaystyle U|a^{r+1}~{}({\rm mod}~{}N)\rangle=|a^{r+2}~{}({\rm
mod}~{}N)\rangle=|a^{2}~{}({\rm mod}~{}N)\rangle$ $\displaystyle\cdots$
$\displaystyle~{}~{}\ ,$

which can be summarized by

$\displaystyle
U^{x}|1\rangle=|f(x)\rangle~{}~{}\text{and}~{}~{}|f(x+1)\rangle=U\,|f(x)\rangle$
(201)

for any non-negative integer $x$. We shall use these relations later in the
manuscript when constructing the ME operators $U$. Note that we can express
the front-end output state from (199) in the form

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes
U^{k}\,|1\rangle\ ,$ (202)

which shows that the ME operators between point-1 and point-2 of Fig. 14 have
the effect of augmenting the work register by $U^{k}|1\rangle$ for every mode
$k\in\\{0,1,\cdots,M-1\\}$. We have almost finished analyzing the circuit. Now
that we have state $|\psi_{2}\rangle$ of (197) in hand, we apply the inverse
quantum Fourier transform $QFT^{\dagger}$ to the control register, producing
the final state

$\displaystyle|\psi_{3}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\ell=0}^{M-1}\sum_{s=0}^{r-1}\,A_{\ell}(\phi_{s})\,|\ell\rangle\otimes|u_{s}\rangle\
,$ (203)

where the amplitudes $A_{\ell}(\phi_{s})$ are given by (186). We now have a
linear superposition of the states $|\ell\rangle\otimes|u_{s}\rangle$, and as
described above, a measurement on the control register will produce any one of
them,

$\displaystyle\sum_{\ell=0}^{M-1}\sum_{s=0}^{r-1}\,A_{\ell}(\phi_{s})\,|\ell\rangle\otimes|u_{s}\rangle\to|\tilde{\ell}\rangle\otimes|\tilde{u}_{s}\rangle$
(204)

with probability $P_{\tilde{\ell},s}=|A_{\tilde{\ell}}(\phi_{s})|^{2}$. The
measured phase $\tilde{\phi}_{\ell}$ is given in terms of the measured output
index $\tilde{\ell}$ of the control register by
$\tilde{\phi}_{\ell}=\tilde{\ell}/2^{m}$, and Shor’s algorithm is designed to
give dominant peaks close to the exact phases $\phi_{s}=s/r$. We will have
more to say about the details of this procedure in later sections. For now, we
summarize the action of the Shor circuit in Table 1. Note that when the exact
phases $\phi_{s}=s/r$ can be represented by $m$-bit fractions, the phase
histogram simplifies considerably: There are exactly $r$ equally likely peaks,
each corresponding to one of the phases $\phi_{s}$. This situation is
relatively rare, and does not occur in most cases, thereby allowing for more
complex phase histograms.

Table 1: Quantum Period Finding

  1. 1.

Initialize the state:
Prepare the state $|0\rangle^{\otimes m}\otimes|0\rangle^{\otimes n}$, and
apply $X$ to the lowest order qubit of the work register to produce the
initial state

$\displaystyle|\psi_{0}\rangle$ $\displaystyle=$
$\displaystyle|0\rangle^{\otimes m}\otimes|1\rangle\ ,$ (205)

where the number of work qubits is given by $n=\lceil\log_{2}N\rceil$.

  2. 2.

Randomize the control register:
Apply $H^{\otimes m}$ to the control register to give

$\displaystyle|\psi_{1}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes|1\rangle$
(206) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{rM}}\sum_{k=0}^{M-1}\sum_{s=0}^{r-1}\,|k\rangle\otimes|u_{s}\rangle\
,$ (207)

where the total number of quantum states is $M=2^{m}$.

  3. 3\.

Modular exponentiation (ME):
Conditionally apply the ME operators $U^{p}$ for
$p\in\\{2^{0},2^{1},\cdots,2^{m-1}\\}$ successively to the control qubits to
produce the state

$\displaystyle|\psi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{r\,M}}\sum_{k=0}^{M-1}\sum_{s=0}^{r-1}e^{2\pi
ik\,\phi_{s}}\,|k\rangle\otimes\,|u_{s}\rangle~{}~{}~{}\text{where}~{}~{}\phi_{s}=\frac{s}{r}$
(208) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{M}}\sum_{k=0}^{M-1}|k\rangle\otimes|a^{k}~{}({\rm
mod}~{}N)\rangle\ .$ (209)

  4. 4\.

Perform the inverse Fourier transform:

$\displaystyle|\psi_{3}\rangle$ $\displaystyle=$ $\displaystyle
QFT^{\dagger}|\psi_{2}\rangle=\sum_{\ell=0}^{M-1}\sum_{s=0}^{r-1}\,A_{\ell}(\phi_{s})\,|\ell\rangle\otimes|u_{s}\rangle\
,$ (210)

where the amplitudes $A_{\ell}(\phi_{s})$ are given by (186).

  5. 5\.

Perform a measurement of the control register:
The state then collapses into an Eigenstate,

$\displaystyle\sum_{\ell=0}^{M-1}\sum_{s=0}^{r-1}\,A_{\ell}(\phi_{s})\,|\ell\rangle\otimes|u_{s}\rangle\to|\tilde{\ell}\,\rangle\otimes|\tilde{u}_{s}\rangle\
,$ (211)

with probability $P_{\tilde{\ell},s}=|A_{\tilde{\ell}}(\phi_{s})|^{2}$. The
probability peaks at values of $\tilde{\ell}$ close to the exact phases
$\phi_{s}=s/r$, and the outcomes of the measurement are equally distributed
between the values of $s\in\\{0,1,\cdots,r-1\\}$.

  6. 6\.

Apply the method of continued fractions to extract the exact phase
$\phi_{s}=s/r$ from the approximately measured $m$-bit phase
$\tilde{\phi}_{\ell}$. This provides the exact period $r$, and therefore the
factors of $N$.

###  V.4 Extracting the Exact Period from the Measured Phase

We turn now to the (non-trivial) task of extracting the exact phase
$\phi_{s}=s/r$ from the approximately measured phase $\tilde{\phi}_{\ell}$ of
the control resister using the method of continued fractions outlined in
Section IV. For ease of notation, we shall henceforth drop the superfluous
$\ell$-subscript from $\tilde{\phi}_{\ell}$ and simply write $\tilde{\phi}$.
Since the measured phase $\tilde{\phi}$ is a positive rational number less
than one, it can be expressed as a finite continued fraction of the form

$\displaystyle\tilde{\phi}=\displaystyle\frac{1}{a_{1}+\displaystyle\frac{1}{a_{2}+\displaystyle\frac{1}{a_{3}+\ddots+\frac{1}{a_{\scriptscriptstyle
R}}}}}$ (212)

for some integer $R$, where $a_{1},a_{2},\cdots,a_{\scriptscriptstyle R}$ are
all positive integer coefficients. We shall drop the zero digit $a_{0}=0$, and
simply write $\tilde{\phi}=[a_{1},a_{2},\cdots,a_{\scriptscriptstyle R}]$ to
denote the form of the continued fraction. Suppose now that $s$ and $r$ are
two relatively prime positive integers that satisfy the inequality

$\displaystyle\left|\frac{s}{r}-\tilde{\phi}\right|\leq\frac{1}{2r^{2}}\ .$
(213)

Then by Theorem 4 of the Section IV, the ratio $s/r$ is necessarily a
convergent of the continued fraction (212) for $\tilde{\phi}$, so that
$s/r=[a_{1},a_{2},\cdots,a_{q}]$ for some $q\leq R$. We shall exploit this
fact to extract the exact period $r$ of the modular exponential function
$f_{a,{\scriptscriptstyle N}}(x)=a^{x}~{}({\rm mod}~{}N)$.

Continued fractions and their convergents can be calculated efficiently using
a variation of the Euclidean algorithm for finding the greatest common divisor
of two integers. Let us call the convergents of $\tilde{\phi}$ by
$s_{0}/r_{0},s_{1}/r_{1},\cdots,s_{\ell}/r_{\ell},\cdots,s_{q}/r_{q}$. We then
cycle through these convergents, from the smallest to the largest values of
the $r$’s, checking to ensure that

$\displaystyle r\text{ is even}$ (214) $\displaystyle a^{r/2}\neq\pm 1~{}({\rm
mod}~{}N)$ (215) $\displaystyle a^{r}=1~{}({\rm mod}~{}N)$ (216)

for each $r=r_{\ell}$. Condition (215) simply means that $b=a^{r/2}$ is not a
trivial root of unity. In passing, we note that $r$ can in fact be odd,
provided that $a$ is a perfect square, so that $b=a^{r/2}$ is still an integer
[12]. Apart from this caveat, if any of the conditions (214)–(216) are not met
for $r_{\ell}$, then the trial fails, and we move on to the next convergent.
The special cases $s_{\ell}=0$ and $s_{\ell}=1$ correspond to the phases
$\phi_{0}=0/r_{\ell}=0$ and $\phi_{1}=1/r_{\ell}$. Technically, we cannot use
the method of continued fractions for these cases since $0/r_{\ell}$ and
$1/r_{\ell}$ are not the ratio of two primes (as 0 and 1 are not prime). In
the former case, $\phi_{0}=0$ yields no information, and we must move on to
the next iteration. However, $\phi_{1}=1/r_{\ell}$ yields an integer
$r_{\ell}$, which could in principle be the correct period that we are
seeking. Therefore we shall check the case $s_{\ell}=1/r_{\ell}$ just to make
sure, even though $1/r_{\ell}$ does not satisfy the conditions of Theorem 4.
In any event, a detailed analysis shows that the probability of achieving a
solution after just a few trials is quite high. The solution for the smallest
value of $r_{\ell}$ gives the period $r=r_{\ell}$ that we are seeking, and the
factors of $N$ are then given by ${\rm gcd}(a^{r/2}\pm 1,N)$.

In closing this section, we should call attention to an important feature of
the method. We must somehow ensure that inequality (213) always holds, as it
would be quite cumbersome if we had to check this by hand every time. However,
this requirement can be hard-wired into the algorithm itself by choosing an
appropriate number of qubits $m$ for the control register, one that is based
on the inequality (213) itself. Recall that the work register has
$n=\lceil\log_{2}N\rceil$ qubits, and thus for any phase $\phi_{s}=s/r$, we
have $r<N\leq 2^{n}$. This leads to the inequality

$\displaystyle\frac{1}{2^{2n+1}}\leq\frac{1}{2r^{2}}\ .$ (217)

Therefore, if we take the control register to have $m=2n+1$ qubits, then a
measurement of the phase $\tilde{\phi}$ will necessarily have sufficient
precision to ensure that

$\displaystyle\left|\frac{s}{r}-\tilde{\phi}\right|=\left|\phi_{s}-\tilde{\phi}\right|\leq\frac{1}{2^{2n+1}}\leq\frac{1}{2r^{2}}\
,$ (218)

and then inequality (213) is automatically satisfied. This is the source of
our previous requirement that $m=2n+1$. Of course we must use even more
control qubits to account for machine error. For a probability of success at
least as large as $1-\varepsilon$, where $\varepsilon>0$ is a small
probability of failure, we must use $m=2n+1+n_{\varepsilon}$ qubits in the
control register, where $n_{\varepsilon}\equiv\lceil{\rm
log}_{2}\big{(}2+1/(2\varepsilon)\big{)}\rceil$. With this choice of $m$, upon
taking a measurement of the control register we find: (i) inequality (213)
will always be satisfied, (ii) the exact phase will be of the form
$\phi_{s}=s/r$, where $s\in\\{0,1,\cdots,r-1\\}$ is randomly selected, and
(iii) the ratio $s/r$ will be a convergent of the continued fraction for
$\tilde{\phi}$ (provided that $s$ and $r$ are relatively prime). However, we
should emphasize that conditions (214)–(216) must also be satisfied. If they
are not, then the method will fail, and we must move on to another iteration
of Shor’s algorithm. However, a complete error analysis shows that the
probability of success is quite high, and typically only a few iterations will
be required before obtaining a factor. We summarize Shor’s algorithm in Table
2.

Table 2: Shor’s Factorization Algorithm

  1. 1.

Chose a random number $a\in\\{2,3,\cdots,N-1\\}$ for the base. If ${\rm
gcd}(a,N)=1$, then proceed to the next step (otherwise we have found a non-
trivial factor of $N$ as required).

  2. 2.

Use quantum phase estimation (QPE) to measure the phase $\tilde{\phi}$ of the
modular exponentiation operator $U_{a{\scriptscriptstyle N}}$ defined by

$U_{a{\scriptscriptstyle N}}|w\rangle=|a\cdot w~{}({\rm
mod}~{}N)\rangle~{}~{}\text{with}~{}~{}U_{a{\scriptscriptstyle
N}}|u_{s}\rangle=e^{2\pi i\,s/r}\,|u_{s}\rangle\ ,$

where $s\in\\{0,1,\cdots,r-1\\}$. The QPE is the only quantum component of
Shor’s algorithm. This is also the bottleneck of the algorithm, as (i) most of
the quantum resources are deployed here, and (ii) a different operator
$U_{a{\scriptscriptstyle N}}$ is required for every choice of $a$ and $N$.

  3. 3.

We then use the method of continued fractions to extract the exact period $r$
from the approximately measured phase $\tilde{\phi}$. To do this, we examine
all convergents $\phi_{s}=s/r$ of $\tilde{\phi}$ such that

$\left|\frac{s}{r}-\tilde{\phi}\right|\leq\frac{1}{2r^{2}}\ ,$

which is achieved by requiring the number of control qubits to be $m=2n+1$. We
then check the convergents from the smallest to the largest values of $r$. If
$r$ is odd, then return to step 1. If $a^{r/2}=\pm 1~{}({\rm mod}~{}N)$, then
return to step 1. If $a^{r}=1~{}({\rm mod}~{}N)$, then we have found a
solution for $r$, and we proceed to the next step; otherwise return to step 1.

  4. 4.

Factors of $N$ are now given by ${\rm gcd}(a^{r/2}\pm 1,N)$. Finding the
greatest common divisor can be done very efficiently on a classical computer
using Euclid’s algorithm.

##  VI Verifying the Theory: Factoring N = 15

To highlight the principal aspects of Shor’s algorithm, we now build the
computational machinery to factor the number $N=15$. In the next section we
will employ the scripts developed here to factor larger and more complex
numbers. We shall employ IBM’s circuit simulator Qiskit. This means that we
must use the Qiskit qubit ordering convention in which the upper 0-th qubit
corresponds to the lowest order bit. We developed the Shor factorization
circuit in Section V.3, which is illustrated in Fig. 14. However, this
analysis used the physics and mathematics ordering convention rather than the
Qiskit convention. Consequently, we must convert to the Qiskit ordering
displayed in Fig. 15. In this convention, the measurement of the control
register (position 4 in the Figure, and indicated by the bold red line across
the register) gives the output state
$|\tilde{\ell}\,\rangle=|\tilde{\phi}_{m-1}\cdots\tilde{\phi}_{1}\,\tilde{\phi}_{0}\rangle$.

Figure 15: Shor factorization algorithm for Qiskit convention 1. In
comparison, Fig. 14 uses the physics Convention 2.

Just as in the previous section, we will employ the convention in which the
measured ($m$-bit) phase in the control register is written with a tilde, and
where the output state is indexed by the binary integer

$\displaystyle\tilde{\ell}$ $\displaystyle\equiv$
$\displaystyle\tilde{\phi}_{m-1}\,\tilde{\phi}_{m-2}\,\cdots\,\tilde{\phi}_{1}\,\tilde{\phi}_{0}~{}~{}~{}\text{where}~{}~{}\tilde{\phi}_{k}\in\\{0,1\\}$
(219) $\displaystyle=$ $\displaystyle
2^{m-1}\,\tilde{\phi}_{m-1}+2^{m-2}\,\tilde{\phi}_{m-2}+\cdots+2^{1}\,\tilde{\phi}_{1}+2^{0}\,\tilde{\phi}_{0}\
.$ (220)

One must keep in mind that there is always an implicit $m$-bit resolution
associated with any measured quantity in the control register, and in
particular, the corresponding angular phase is given by the $m$-bit fraction

$\displaystyle\tilde{\phi}_{\ell}\equiv\frac{\tilde{\ell}}{2^{m}}$
$\displaystyle=$
$\displaystyle\frac{\tilde{\phi}_{m-1}}{2^{1}}+\frac{\tilde{\phi}_{m-2}}{2^{2}}+\cdots+\frac{\tilde{\phi}_{1}}{2^{m-1}}+\frac{\tilde{\phi}_{0}}{2^{m}}$
(221) $\displaystyle=$ $\displaystyle
0.\tilde{\phi}_{m-1}\,\tilde{\phi}_{m-2}\,\cdots\,\tilde{\phi}_{1}\,\tilde{\phi}_{0}\
.$ (222)

For ease of notation, we shall drop the $\ell$-subscript and denote the
measured phase by $\tilde{\phi}$.

Since we will employ both binary and decimal numbers in this section, we will
often denote binary numbers using a bracket with a 2-subscript, writing
$\tilde{\phi}=[0.\tilde{\phi}_{m-1}\,\tilde{\phi}_{m-2}\,\cdots\,\tilde{\phi}_{1}\,\tilde{\phi}_{0}]_{2}$.
We will always place a tilde over measured quantities, so that $\tilde{\phi}$
denotes the phase as determined by a measurement of the $m$-bit control
register, as opposed to the (as yet undetermined) exact phase $\phi_{s}=s/r$.
If the control register contains a sufficient number of qubits, then the
measurement should be quite accurate, and the measured phase will be very
close to the exact phase, $\tilde{\phi}\approx\phi_{s}=s/r$. Since
$\tilde{\phi}$ and $\phi_{s}=s/r$ are approximately equal, we can extract the
exact value of the integers $s$ and $r$ using the method of continued
fractions (provided that $s$ and $r$ are relatively prime). The integer $r$ is
the period of the modular exponential function $f_{a\,{\scriptscriptstyle
N}}(x)=a^{x}~{}({\rm mod}~{}N)$ that we seek.

Let us now concentrate on the specific example of $N=15$. As we have
established, the work space requirement is $n=\lceil\log_{2}15\rceil=4$
qubits, and the control register must therefore contain $m=2n+1=9$ qubits (for
simplicity we consider only perfect measurements in which $n_{\epsilon}=0$).
For a given integer $N$, we choose the base $a$ to be a random integer such
that $1<a<N$ and ${\rm gcd}(a,N)=1$. Therefore, for $N=15$ we can only choose
$a\in\\{2,4,7,8,11,13,14\\}$. It turns out that $a=14$ gives a trivial root of
unity, so we can neglect this choice. In fact, in this example we shall take
either $a=4$ or $a=8$, where the former gives the period $r=2$ and the latter
gives $r=4$. Figure 16 illustrates the Qiskit circuit for Shor’s algorithm
with the base $a=8$. The first 9 qubits comprise the control register and are
labeled by the index $c$, and the last 4 qubits are the work register and are
labeled by $w$. We will therefore denote quantum states of the work register
by $|w_{3}w_{2}w_{1}w_{0}\rangle$. In the Qiskit qubit convention, the work
register is initially populated by the state

$\displaystyle|1\rangle=|0001\rangle=|1\rangle_{0}\otimes|0\rangle_{1}\otimes|0\rangle_{2}\otimes|0\rangle_{3}\
,$ (223)

with the least significant bit being $w_{0}=1$. The controlled modular
exponentiation operators $CU_{8,15}^{p}$ for
$p\in\\{2^{0},2^{1},\cdots,2^{8}\\}$ are represented by purple boxes attached
to their respective control qubits, and the operator $QFT^{\dagger}$ is
represented by the large purple rectangle on the far right of the control
register. At the end of the circuit, the control register undergoes a
measurement on all $m$ qubits. The work register might or might not undergo a
measurement, and we shall return to this point later in the section.

Figure 16: Qiskit version of Shor’s factoring circuit for $N=15$ and $a=8$,
with $m=9$ qubits in the control register and $n=4$ qubits in the work
register.

###  VI.1 Modular Exponentiation Operators

Let us now explore the modular exponentiation (ME) operators
$U_{a\,{\scriptscriptstyle N}}$ in more detail. For every choice of the number
$N$ and the base $a$, we must design a separate implementation of the operator
$U_{a\,{\scriptscriptstyle N}}$, and this is in fact the real bottleneck of
Shor’s algorithm. Indeed, this bottleneck occurs in two senses: (i) the ME
operators consume the greatest majority of the quantum resources of the
algorithm, and in the general case this will be of order $72\,n^{3}$ gates
[5], and (ii) even specialized cases of the ME operators are often highly non-
trivial to construct. The ME operators for $N=15$ with $a=4$ and $a=8$ are
illustrated in Fig. 17,

Figure 17: Modular exponentiation operators $U_{a,15}$ for $a=4$ in the left
panel and $a=8$ in the right panel. Their action on the work space is given by
$U_{4,15}|w_{3}w_{2}w_{1}w_{0}\rangle=|w_{1}w_{0}w_{3}w_{2}\rangle$ and
$U_{8,15}|w_{3}w_{2}w_{1}w_{0}\rangle=|w_{0}w_{3}w_{2}w_{1}\rangle$, which can
be reproduced by a sequence of SWAP gates.

and Table 3 gives the corresponding operations $U_{a,15}|w\rangle=|a\cdot
w~{}({\rm mod}~{}15)\rangle$ for every basis element $|w\rangle$ in the work
register. Each ME operator $U_{a,15}$ has two columns in the Table: One for
the decimal representation of the basis elements $|w\rangle$, and another for
the corresponding binary representation $|w_{3}w_{2}w_{1}w_{0}\rangle$, where
in binary form $w=[w_{3}w_{2}w_{1}w_{0}]_{2}$. The operators $U_{a,15}$ in
Fig. 17 were determined through simple inspection of the results of the binary
columns in Table 3. For example, the $a=8$ operator performs a permutation of
the 4-bit binary work states,
$U_{8,15}|w_{3}w_{2}w_{1}w_{0}\rangle=|w_{0}w_{3}w_{2}w_{1}\rangle$, which can
be implemented by the three SWAP gates in the right panel of Fig. 17.
Similarly, the ME operator for $a=4$ performs two SWAP operations on the work
register, so that
$U_{4,15}|w_{3}w_{2}w_{1}w_{0}\rangle=|w_{1}w_{0}w_{3}w_{2}\rangle$. The other
ME operators for $a\in\\{2,7,11,13\\}$ act similarly, and are illustrated in
Fig. 18.

Table 3:  Modular Exponentiation Operators $U=U_{a,15}$ for $a=4$ and $a=8$.
$U_{4,15}|w\rangle=|4\cdot w~{}({\rm mod}~{}15)\rangle$ |  $U_{8,15}|w\rangle=|8\cdot w~{}({\rm mod}~{}15)\rangle$
---|---
$U|1\rangle~{}=|4\rangle$ |  $U|0001\rangle=|0100\rangle$ | $U|1\rangle=|8\rangle$ |  $U|0001\rangle=|1000\rangle$
$U|2\rangle~{}=|8\rangle$ |  $U|0010\rangle=|1000\rangle$ |  $U|2\rangle=|1\rangle$ |  $U|0010\rangle=|0001\rangle$
$U|3\rangle~{}=|12\rangle$ |  $U|0011\rangle=|1100\rangle$ | $U|3\rangle=|9\rangle$ |  $U|0011\rangle=|1001\rangle$
$U|4\rangle~{}=|1\rangle$ |  $U|0100\rangle=|0001\rangle$ | $U|4\rangle=|2\rangle$ |  $U|0100\rangle=|0010\rangle$
$U|5\rangle~{}=|5\rangle$ |  $U|0101\rangle=|0101\rangle$ |  $U|5\rangle=|10\rangle$ |  $U|0101\rangle=|1010\rangle$
$U|6\rangle~{}=|9\rangle$ |  $U|0110\rangle=|1001\rangle$ | $U|6\rangle=|3\rangle$ |  $U|0110\rangle=|0011\rangle$
$U|7\rangle~{}=|13\rangle$ |  $U|0111\rangle=|1101\rangle$ |  $U|7\rangle=|11\rangle$ |  $U|0111\rangle=|1011\rangle$
$U|8\rangle~{}=|2\rangle$ |  $U|1000\rangle=|0010\rangle$ | $U|8\rangle=|4\rangle$ |  $U|1000\rangle=|0100\rangle$
$U|9\rangle~{}=|6\rangle$ |  $U|1001\rangle=|0110\rangle$ |  $U|9\rangle~{}=|12\rangle$ |  $U|1001\rangle=|1100\rangle$
$U|10\rangle=|10\rangle$ |  $U|1010\rangle=|1010\rangle$ | $U|10\rangle=|5\rangle$ |  $U|1010\rangle=|0101\rangle$
$U|11\rangle=|14\rangle$ |  $U|1011\rangle=|1110\rangle$ |  $U|11\rangle=|13\rangle$ |  $U|1011\rangle=|1101\rangle$
$U|12\rangle~{}=|3\rangle$ |  $U|1100\rangle=|0011\rangle$ | $U|12\rangle=|6\rangle$ |  $U|1100\rangle=|0110\rangle$
$U|13\rangle~{}=|7\rangle$ |  $U|1101\rangle=|0111\rangle$ |  $U|13\rangle=|14\rangle$ |  $U|1101\rangle=|1110\rangle$
$U|14\rangle~{}=|11\rangle$ |  $U|1110\rangle=|1011\rangle$ | $U|14\rangle=|7\rangle$ |  $U|1110\rangle=|0111\rangle$

Figure 18: Modular exponentiation operators $U_{a,15}$ for $a=2,7,11,13$ as
labeled.

###  VI.2 The Factorization Circuit

Our next goal is to construct the Python script used in creating the Shor
circuit in Fig. 16, after which we shall run this code on the Qiskit simulator
to factor $N=15$. We will construct the requisite script slowly in sequences
of Python code segments, each one illustrating an essential element of the
algorithm. First, we must import the necessary python packages:

    # import basics
    import numpy as np
    from random import randint
    from math import gcd

    # import Qiskit tools
    from qiskit import Aer, transpile, assemble
    from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

    # import plot tools
    from qiskit.visualization import plot_histogram
    import matplotlib.pyplot as plt

The first few imports are for basic mathematical functionality. For example,
we can create random integers with ${\tt randint}()$, and find the greatest
common divisor with ${\tt gcd}()$. We have also imported packages to run
Qiskit, and for post-processing the Qiskit data. We now select $N=15$ and the
corresponding base $a$:

    # number to factor
    N = 15

    # random number a in [2,N-1] wtih gcd(a,N)=1
    n = 0
    while n == 0:
        a = randint(2, N-1)
        if gcd(a, N) == 1: n = 1
        print("**:", a, N, gcd(a, N))

This code segment chooses a random integer between $2$ and $N-1$ inclusive,
and makes sure that the choice does not contain a non-trivial factor in common
with $N$ (otherwise we have found a sought after factor of $N$). This piece of
code can be omitted if we wish to work only with a specific value of $a$
(provided of course that we set the values of $N$ and $a$ here).

The next code segment defines the ME operators $CU_{a,{\scriptscriptstyle
N}}^{p}$ for $N=15$ for all permissible choices of $a$. This is really the
heart of Shor’s algorithm. For a general value of $N$, we would not be able to
implement all values of $a$ (we do this here only for purposes of
illustration), as there are an exponentially large number of them. Finally, we
define a subroutine for the inverse Fourier transform $QFT^{\dagger}$. The
code is given below:

    # modular exponentiation gates: p = 2^0, 2^1, .... , 2^(m-1)
    def c_Uamod15(a, p):
        U = QuantumCircuit(4)
        # concatenate U-factors to form U^p
        for iteration in range(p):
            if a in [2,13]:
                U.swap(0,1)
                U.swap(1,2)
                U.swap(2,3)
            if a in [7,8]:
                U.swap(2,3)
                U.swap(1,2)
                U.swap(0,1)
            if a in [4, 11]:
                U.swap(1,3)
                U.swap(0,2)
            if a in [7,11,13]:
                for q in range(4):