1 COMPUTATION AND COMPLEXITY CLASSES

Shor’s Algorithm and the Quantum Fourier Transform

Fang Xi Lin
McGill University
fangxi.lin@mail.mcgill.ca

Abstract

Large numbers have traditionally been believed to be difficult to factor efficiently on a classical
computer. Shor’s quantum algorithm gives a way to factor integers in polynomial time using a
quantum computer. In addition, the algorithm also allows the computation of discrete logarithms
in polynomial time. The algorithm relies in a crucial way on the quantum Fourier transform.
We will briefly introduce quantum mechanics and quantum computation, then describe both the
quantum Fourier transform and Shor’s algorithm in detail.

Introduction
The problem of how to factor a large integer efficiently has been studied extensively in number
theory. It is generally believed that factorization of a number n is hard to do in a efficient way. That
is, it cannot be done in a number of steps which is polynomial in the length of the integer we’re
trying to factor1 . The RSA cryptosystem, among others, relies on the presumed difficulty of this
task. Classically, the fastest known algorithm is the General Number Field Sieve (GNFS) algorithm,
which works in super-polynomial, but sub-exponential time.
In 1994, Peter Shor discovered an algorithm that can factor numbers in polynomial time using a
quantum computer[10] , a drastic improvement over the GNFS. Shor’s algorithm consists of a classical
and a quantum part. The classical part involves modular exponentiation via repeated squaring,
which can be performed quickly. The quantum part involves a “quantum Fourier transform”. We
will prove that in certain cases, the quantum Fourier transform can be constructed in polynomial
time, wherein lies the efficiency of the algorithm.
Although Shor’s factoring algorithm is much more publicized, Shor’s ideas will allow us to
compute discrete logarithms, which is also believed to be a hard task for a classical computer in the
same sense that factoring numbers is.
There have been several successful experimental demonstrations of the factoring algorithm. In
2001, A group from IBM was able to factor 15 using a quantum computer of 7 qubits implemented
using nuclear magnetic resonance[13] . Since then, 15 and 21 were factored using photonic qubits.
However, there were concerns that some of these experimentations were not true compilations of
Shor’s algorithm[12] .

1 Computation and Complexity Classes

The ability to compute is limited by two resources: space (memory) and time. The difficulty of
computing allows problems to be categorized into different complexity classes. This is the subject of
1
Note that a number of size d has input length log d.

1
 2 QUANTUM MECHANICS AND QUANTUM COMPUTATION

study for a computation complexity theorist. We will, however, try to give some intuitive insight
into the theory of complexity classes.
Consider an algorithm which takes in an input of length n (for example, the number of digits
in a number). We call this a polynomial time algorithm if it doesn’t take more than C nk steps,
for some fixed C, k > 0, to compute the answer. We denote this by O(nk ). These are considered
efficient algorithms (although n1000 is not really efficient in practice). The class of problems solved
by these algorithms is called P.
Another important complexity class is called NP. This is the class of problems whose solutions
can be verified in polynomial time. For example, once we find the factorization of some number
N = pq, we can efficiently verify that pq = N . Indeed, we have P ⊆ NP (the reverse inclusion is a
open problem, one of the seven Millennium problems).
Both complexity classes presented above are bounded by time. There are also a number of
complexity classes bounded by space. PSPACE is such a class, which contains problems that can be
solved with a polynomial amount of bits in input size.
There are two more complexity classes which are important to this paper. The first is the BPP,
which are problems that can be solved with a bounded probability of error in polynomial time. The
second one is the BQP, which is essentially the same thing on a quantum machine. Factoring numbers
using Shor’s algorithm is BQP.
The known relationship between these complexity classes is:

P ⊆ BPP, NP, BQP ⊆ PSPACE (1)

In addition, we also have BPP ⊆ BQP. The relationship between BPP, BQP and NP is unknown.

2 Quantum Mechanics and Quantum Computation

Shor’s algorithm is a quantum algorithm. That is, it exploits the fact that we can have superpositions
of quantum states. We will follow Nielsen and Chuang[8] to build up some concepts in quantum
mechanics, quantum computation, and quantum circuits in the following section.

2.1 Quantum Mechanics

There are multiple formulations of quantum mechanics. The one which will be useful to us is the
matrix mechanics (due to Heisenberg) formulation. Without digging into subtle details, objects
in quantum mechanics have essentially a one-to-one correspondence to objects in linear algebra.
The quantum mechanical vector spaces are called Hilbert spaces, which are normed complex vector
spaces. We will assume that the Hilbert spaces are finite dimensional in this paper.
The standard quantum mechanical notation2 for a (column) vector in a vector space is a |ψ" (a
ket). It is used to represent a quantum state. The notation for a dual (row) vector is #φ| (a bra).
The inner product between two vectors is then denoted #φ|ψ" (a bra-ket), which yields in complex
number in general. In general, we require that the quantum states be normalized, i.e. | #ψ|ψ" |2 = 1.
Physically, | #φ|ψ" |2 is the probability of observing |ψ" in the state |φ". We will call the zero vector
0 instead of |0" since the latter will be used for something else.
The matrices are called operators in quantum mechanics. We will require that all our operators
to be unitary, i.e., given an operator A, we need A† := (AT )∗ = A−1 . We will use #φ| A |ψ" to denote
the inner product between #φ| and A |ψ".3
2
Due to Paul Dirac
3
Equivalently, between !φ| A† and |ψ"

2
2 QUANTUM MECHANICS AND QUANTUM COMPUTATION 2.2 Quantum Computation

There is a useful way of representing linear operators called the outer product representation.
The outer product is a ket-bra, which is a linear operator defined by (|ψ" #φ|) |φ# " := #φ# |φ" |ψ".
Now let A be an operator
! and let |j" be any orthonormal basis. Then any quantum state |ψ"
can be written as |ψ" = j αj |j" and #j|ψ" = αj . Therefore we have
 
$ $ $
 |j" #j| |ψ" = #j|ψ" |j" = αj |j" = |ψ" (2)
j j j
!
for any arbitrary |ψ". Hence, we have I = j |j" #j|. This allows us to represent any operator A as
a linear combination of outer products since
$ $
A = IAI = |j" #j| A |k" #k| = #j| A |k" |j" #k| (3)
j,k j,k

so that #j| A |k" is the j, k th entry of the matrix representation of A.

We will make use of the tensor product, ⊗, extensively. In quantum mechanics, tensor products
are used to describe multi-particle systems (in our case, multi-qubit systems).
Definition 1. Let HA be an N dimensional Hilbert space with orthonormal basis |j" and let HB
be an M dimensional Hilbert space with orthonormal basis |k". Then the tensor product HA ⊗ HB
is the N M dimensional system spanned by the pair |j" ⊗ |k" where
' # (
#j | ⊗# k # | (|j" ⊗ |k") = #j # |j" #k # |k" = δjj ! δkk! . (4)
! !
To every vector
! |ψA " = j αj |j" ∈ HA and |ψB " = k βk |k" ∈ HB we can associate a vector
|ψA " ⊗ |ψB " = jk αj βk |j" ⊗ |k" ∈ HA ⊗ HB . We call the resulting state a product state. We
sometimes suppress the ⊗ symbol and write |ψA " ⊗ |ψB " as |ψA " |ψB " or even |ψA ψB ". Not all states
in HA ⊗ HB can be written as |ψA " ⊗ |ψB " for some |ψA " ∈ HA and |ψB " ∈ HB (otherwise we
would have an N + M dimensional system!). The states which are not product states are said to be
entangled.
Similarly, for every operator A acting on |ψA " and operator B acting on |ψB ", we can define an
operator A ⊗ B acting on |ψA " ⊗ |ψB " as (A ⊗ B)(|ψA " ⊗ |ψB ") = A |ψA " ⊗ B |ψB ".

2.2 Quantum Computation

Just as classical computers use bits to compute, quantum computers use quantum bits or qubits.
A qubit is a vector |ψ" = α0 |0" + α1 |1". Operators which act on qubits are two-by-two unitary
matrices. We call these operators quantum gates, since they are analogous to the logic gates in
classical computer science. In general, a circuit of quantum gates, or a quantum circuit, looks like

Figure 1: A quantum circuit. Time goes from left to right, wires represent qubits, and rectangles represent
unitary operators/quantum gates.

3
 3 CLASSICAL PART: REDUCTION TO ORDER FINDING

The most important gate to us will be the Hadamard gate.

Definition 2. Let {|0" , |1"} be a basis. Then the Hadamard gate is
) *
1 1 1
H→√ (5)
2 1 −1
In particular, it sends the basis vectors to a uniform superposition of the basis vectors. That is, we
get 12 chance of measuring either |0" or |1". In a quantum circuit, it is represented by

H
Figure 2: A Hadamard gate

As we will see in section 4, the quantum Fourier transform can be represented with Hadamard
gates. We will also need the notion of a controlled gate.
Definition 3. If U is a single qubit unitary operation, a controlled-U is a two qubit operation on a
control and a target qubit such that if the control qubit is set, then the gate will act on the target
qubit. If not, the target qubit is left alone. That is, we have


|0" |0" )→ |0" |0"


|0" |1" )→ |0" |1"
Uc : (6)

 |1" |0" →
) |1" U |0"


|1" |1" )→ |1" U |1"

where the left qubit is the control qubit and the right qubit is the target qubit. In a circuit, it
looks like

U
Figure 3: A controlled-U gate. On the top is the control qubit, on the bottom is the target qubit.

3 Classical Part: Reduction to Order Finding

Shor’s algorithm does not allow us to factor a number directly. Instead, it allows us to find the order
of an element a modulo n in polynomial time. We will show that the problem of finding a non-trivial
factor to n can be reduced (efficiently) to finding the order of a non-trivial element in Z/nZ[6] .
Lemma 4. Given a composite number n, and x non-trivial square root of 1 modulo n (i.e. x2 = 1
(mod n) but x is neither 1 nor −1 mod n, then either gcd(x − 1, n) or gcd(x + 1, n) is a non-trivial
factor of n.
Proof. Since x2 ≡ 1 (mod n), we have x2 − 1 ≡ 0 (mod n). Factoring, we get (x − 1)(x + 1) ≡ 0
(mod n). This implies that n is a factor of (x + 1)(x − 1). Since (x ± 1) +≡ 0 (mod n), n has a
non-trivial factor with x + 1 or x − 1. To find this common factor efficiently, we apply Euclid’s
algorithm to get gcd(x − 1, n) or gcd(x + 1, n).

4
4 FOURIER TRANSFORMS

Example 5. Let n = 55 = 5 · 11. We find that 34 is a square root of 1 mod n since 342 = 1156 =
1 + 21 · 55. Computing, we get gcd(33, 55) = 11 and gcd(35, 55) = 5.

Lemma 6. Let n be odd, then at least half the elements in (Z/nZ)× have even order.

Proof. Suppose ord(x) = r is odd. Then (−x)r = (−1)r xr = (−1)r = −1 (mod n). Hence −x must
have order 2r, which is even. Therefore, at least half the elements in (Z/nZ)× have even order.

Equipped with these tools, we will proceed to prove the main result that allows us to reduce
factorization of n to finding the order of an element in Z/nZ.

Theorem 7. Let n be an odd integer and let n = pe11 pe22 . . . pekk be the prime factorization of n. Then
the probability that a uniform randomly chosen x ∈ Z/nZ has even order r and xr/2 +≡ −1 (mod n)
k−1
is at least 1 − 12 .

Proof. By the Chinese Remainder Theorem, choosing x ∈ (Z/nZ)× (uniform) randomly is equivalent
to choosing xi ∈ (Z/pei i Z)× for each pi randomly. Let r be the order of x and let ri be the order of
xi . In particular, xr/2 is never 1 modulo n. We want to show that the probability of either r being
k−1
odd or xr/2 ≡ −1 (mod n) is at most 12 .
Note that r = lcm(r1 , r2 , . . . , rk ) (where lcm denotes the least common multiple). To see this,
xr ≡ 1 (mod n) ⇒ xr ≡ 1 (mod pei i ), hence r is a multiple of each ri . It is the least such number
and hence the least common multiple of the ri ’s.
Suppose that r is odd. This happens only if all of the ri ’s are odd. ri is odd with probability at
k
most one-half by lemma 6. Hence, r is odd with probability at most 12 .
Now suppose that r is even. We still have to worry about the possibility that xr/2 ≡ ±1
(mod n). By the Chinese Remainder Theorem, this happens only if xr/2 ≡ ±1 (mod pei i ) for every
pi . We need to avoid these cases since ≡ +1 means r wasn’t the order, and ≡ −1 doesn’t yield a
useful factorization. The probability of choosing an x such that one of these two cases happens is
2 · 2−k = 2−k+1 .
Combining the probabilities, we get a success probability of at least (1 − 2−k )(1 − 2−k+1 ) ≥
1 − 3 · 2−k .

By lemma 4 and theorem 7, given a composite number n and the order r of some x ∈ Z/nZ, we
can compute gcd(xrx /2 ± 1, n) efficiently using Euclid’s algorithm. This gives a non trivial factor of
n unless r is odd or xrx /2 ≡ −1 (mod n). In particular, if n is a semi-prime, i.e. it is a product of
two primes p, q, then theorem 7 implies that n will be factored with probability 12 .

4 Fourier Transforms
Since finding a factor of n given the order of some element in Z/nZ can be done efficiently even
on a classical computer, it still remains to be shown that we can find the order of the element
efficiently. It is unknown how to quickly find the order of a given element on a classical computer,
but Shor’s order finding algorithm will allow us to do so by employing a quantum computer. The
order finding algorithm relies crucially on a unitary operator Fn , the “quantum Fourier transform”
(QFT) operator, which acts like a discrete Fourier transform. We assume the knowledge of the usual
Fourier transform for the following section.
Before discussing the quantum Fourier transform, we will talk a bit about the discrete Fourier
transform (DFT) as well as the Fast Fourier Transform (FFT) algorithm. The QFT will be
constructed to be essentially the equivalent of the FFT with a quantum circuit[11] .

5
4.1 Discrete Fourier Transform 4 FOURIER TRANSFORMS

4.1 Discrete Fourier Transform

Definition 8. 4 Let f = (f0 , f1 , . . . , fN −1 ) be a vector in CN . The discrete Fourier transform is a
map

F :CN → CN
f )→ f˜

defined by
N −1
˜ 1 $ −jk
fj := √ ζ fk (7)
N k=0

where ζ = exp( 2πi th root of unity.

N ) is the N

We will use f˜ to denote the DFT of f . As in the case of the usual Fourier transform, there’s an
inverse Fourier transform given by the expected formula.

Lemma 9. The inverse discrete Fourier transform, F −1 , is given by

n−1
1 $ jk ˜
fj = √ ζ fk (8)
N k=0

We can check how the Fourier transform acts on the standard basis. Let {e1 , e2 , . . . , eN } be the
standard basis of CN , where el denotes the vector has has 1 at the lth component and 0 elsewhere
(i.e. elj = δjl . Then the DFT of el is given by

N −1
1 $ −jk 1
ẽlj = √ ζ δjl = √ ζ −jl (9)
N k=0 n

The matrix representation of F in the standard basis is

 
1 1 1 ... 1
1 ζ −1 ζ −2 · · · ζ −(N −1) 
1 
 ζ −2 ζ −4

· · · ζ −2(N −1) 
F → √ 1  (10)
N . .. .. .. .. 
 .. . . . . 
2
1 ζ −(N −1) ζ −2(N −1) · · · ζ −(N −1)

Example 10. Consider the N = 2 DFT. We have ζ −1 = exp(− 2πi

2 ) = −1 and so
) *
1 1 1
F→√ (11)
2 1 −1

which is a rotation by 45◦ , and also looks like a Hadamard gate! We will see more of this when we
discuss the QFT.
4
We use the physicists and mathematicians’ convention to define the DFT and everything that follows. Computer
scientists usually have the sign on the exponent of ζ reversed.

6
4 FOURIER TRANSFORMS 4.2 Fast Fourier Transform

Since performing the DFT on a vector f is like a matrix multiplication, it takes N (complex)
multiplications and N − 1 additions for each component and there are N components so we have
N 2 multiplications and N (N − 1) additions. Since additions can be efficiently computed, the speed
is limited by the N 2 multiplications. Hence, we shall only consider the number of multiplications.
As we shall see shortly, the FFT will allow us to perform DFT in fewer operations by exploiting
some symmetries of the DFT.
Many of the properties of the DFT are analogous to the properties of the FT. For example, we
claim the DFT convolution theorem holds[14] .
Definition 11. A circular or cyclic convolution of f and g, denoted by f ∗g, is a map CN ×CN → CN
given component-wise by
N
$ −1
(f ∗ g)j = fk gj−k (12)
k=0

where g−a := gN −a .
Theorem 12. (Circular Convolution Theorem) Let h = f ∗ g, then h̃j , the j th component of h̃, is
given by f˜j · g̃j where · denotes the usual multiplication.

4.2 Fast Fourier Transform

For the following, we assume that N = 2m .
Consider the roots of unity, ζ. Observe that
) *
j+N/2 2πi · (j + N/2)
ζ = exp
N
) *
2πij (13)
= exp exp (πi)
N
= −ζ j
and similarly
ζ j+N = ζ j . (14)
This suggests that we may be able to split f into smaller parts.
Since f has N = 2m components, we can divide f into an even and an odd part. More precisely,
we define feven = (f0 , f2 , . . . , fN −2 ) and fodd = (f1 , f3 , . . . , fN −1 ). Now, we can rewrite equation (7)
as

N/2−1 N/2−1
1 $ −j·2k 1 $ −j·(2k+1)
f˜j = √ ζ fevenk + √ ζ foddk
N k=0 N k=0
N/2−1 N/2−1
1 $ −j·2k ζ −j $ −j·2k
=√ ζ fevenk + √ ζ foddk
N k=0 N k=0
  (15)
N/2−1 N/2−1
1 1 $ ζ −j $
= √ 1 ζ −j·2k fevenk + 1 ζ −j·2k foddk 
2 N/2 k=0 N/2 k=0
1 2 3
= √ f˜evenj + ζ −j f˜oddj
2

7
4.2 Fast Fourier Transform 4 FOURIER TRANSFORMS

where in the second step we pull out a factor of ζ −j , called the twiddle factor, out of the fodd term,
and in the fourth step we apply the definition of the DFT to feven and fodd .
Since the DFT is periodic with period n (i.e. f˜j+N = f˜j ), feven and fodd are periodic with period
N/2. Hence, combining (13), (14), (15), we get

√
2f˜ = f˜evenj + ζ −j f˜oddj if 0 ≤ j ≤ N/2 − 1 (16)
√
2f˜ = f˜evenj − ζ −j f˜oddj if N/2 ≤ j ≤ N − 1 (17)

This gives us a way to compute the DFT of a vector of size N in terms of smaller vectors of
size N/2 = 2m−1 . To compute feven and fodd , we require (N/2)2 + (N/2) √ = N /2 multiplications.
2 2
−j ˜
To compute ζ fodd , we require another multiplications. For the 2, we can “collect” them
√ N/2 √
and then multiply the total factor of 2m = N in at the end of the recursion as N additional
multiplications (instead of multiplying in every step). Hence, we require N 2 /2 + N/2 multiplications
in total, which is about a factor of two faster than the N 2 multiplications in the original DFT.
Since the smaller vectors still have length divisible by 2, we can apply this procedure recursively
until we get a DFT of size 2, which are all “classical Hadamard transforms” as in example 10. In
general, the FFT allows us to compute the DFT in O(N log N ) = O(2m log 2m ) operations, which is
still exponential in m.

Example 13. (Fast multiplication of two polynomials[5] ) Let p(x) = α0 + α1 x + . . . + αN −1 xN −1

and q(x) = β0 + β1 x + . . . + βN −1 xN −1 be two polynomials with complex coefficients. Then
4N −1 5 N −1 
$ $
p(x)q(x) = α i xi  βj xj 
i=0 j=0
42N −2 5 (18)
$
= λk x k
k=0

where
N
$ −1
λk = αl βk−l (19)
l=0

Therefore, computing p(x)q(x) and hence computing the λk ’s directly takes N 2 multiplications.
However, equation (19) looks very much like a discrete convolution. Let us view p(x), q(x) as vectors
with their coefficients as components. That is, p → (α0 , α1 , . . . , αN −1 ), q → (β0 , β1 , . . . , βN −1 ). Now,
we can append 0 as necessary to p and q to make them 2N dimensional vectors (because we want p
and q to have the same dimension as p ∗ q). We can express the λk ’s as

2N
$ −1
λk = αl mod 2N · βk−l mod 2N (20)
l=0

which gives us the circular convolution of p and q! Hence, we can compute λk ’s by first performing
the DFT on the p and q vectors which takes N log N multiplications. Then, we multiply the resulting
vector component-wise which takes 2N multiplications. Finally, we take the inverse DFT and obtain
the coefficients λk .

8
4 FOURIER TRANSFORMS 4.3 Quantum Fourier Transform

4.3 Quantum Fourier Transform

Finally, we move onto the main topic of this section, the quantum Fourier transform. We will follow
our steps for constructing the DFT and FFT.

! −114. Let {|0" , |1" , . . . , |N − 1"} be an orthonormal basis for a quantum system and let
Definition
|φ" = N j=0 |j" be a quantum state. Then the quantum Fourier transform FN is a map defined by

N
$ −1 N
$ −1 N −1
1 $ −jk
|φ" = |j" )→ √ ζ |k" (21)
j=0 j=0
N k=0

In particular, the basis states transform as

N −1
1 $ −jk
|j" )→ √ ζ |k" (22)
N k=0

And hence we get an representation for FN

N −1
1 $ −jk
FN =√ ζ |k" #j| (23)
N j,k=0

Note that since ζ̄ = ζ −1 , we have

N −1
1 $ jk
Fn† = √ ζ |k" #j| (24)
N j,k=0

We can easily check that the quantum Fourier transform is unitary.

N −1 N −1
1 $ −jk $
FN FN† = ζ |k" #j| ζ rs |r" #s|
N
j,k=0 r,s=0
N
$ −1
1
= ζ rs−jk |k" #j|r" #s|
N
j,k,r,s=0
N
$ −1
1
= ζ rs−jk δjr |k" #s|
N
j,k,r,s=0
4N −1 5 (25)
N −1
1 $ $
= ζ r(s−k) |k" #s|
N
k,s=0 r=0
N
$ −1
= δks |k" #s|
k,s=0
N
$ −1
= |k" #k| = I
k=0

!N −1 2 3
2πir(s−k)
where we used the fact that r=0 exp N = N δsk .

9
4.3 Quantum Fourier Transform 4 FOURIER TRANSFORMS

As in the case of DFT, constructing FN naïvely is not very efficient. We will try to implement
the QFT as a quantum circuit efficiently. Recall that we assumed N = 2m . Consider the Fourier
transformed state (22), if we express j in binary as j1 j2 . . . jm ∈ {0, 1}m and j = j1 2m−1 + j2 2m−2 +
. . . + jm 20 , we will see that it is in fact a product state.
Theorem 15. |j" is a product state which can be written as a product of m qubits[3]

1 2 3 2 3 2 3
|j" → √ |0" + e−2πi(0.jm ) |1" ⊗ |0" + e−2πi(0.jm−1 jm ) |1" ⊗ . . . ⊗ |0" + e−2πi(0.j1 j2 ...jm ) |1"
N
(26)

where (0.j1 j2 . . . jm ) = j1 2−1 + j2 2−2 + . . . + jm 2−m denotes the binary fraction.

Proof. To see this, write out the binary expansion of |k".
n−1
1 $ −jk 1 $ !m
2m−r kr
|j" → √ ζ |k" = √ ζ −j r=1 |k1 " |k2 " . . . |km "
N k=0 N k1 ,k2 ,...,km ∈{0,1}
$ m
6
1 m−r k
=√ ζ −j2 r
|kr "
N k1 ,k2 ,...,km ∈{0,1} r=1
 
6m $
1  m−r k
(27)
=√ ζ −j2 r
|kr "
N r=1 k ∈{0,1}
r
m 2 3 1 62 3
6 m
1 m−r m−r /2m
=√ |0" + ζ −j2 |1" = √ |0" + e−2πij2 |1"
N r=1
n r=1
1 62 3
m
−r
=√ |0" + e−2πij2 |1"
N r=1

where in the second step we expanded the exponential as product and regrouped terms. Notice the
similarities between this procedure and the FFT (we essentially did the whole FFT recursion in one
fell swoop).
Expanding the j in the “twiddle factor” in binary, we get

4 m
5 4 m
5
$ $
exp −2πi 2m−l jl /2r = exp −2πi 2m−r−l jl
l=1 l=1
(28)
= exp (−2πi(0.jm−r+1 jm−r+2 . . . jm ))
so that (27) gives

1 62 3
m
−2πi(0.jm−r+1 jm−r+2 ...jm )
|j" → √ |0" + e |1" (29)
n r=1

as was required.

Note that the value of e−2πi(0.jm−r+1 jm−r+2 ...jm ) is either 1 or −1, like a Hadamard transformed
qubit. Moreover, note that the last qubit depends on all the input qubits but the dependence
decreases as we go further. We can use this to construct a quantum circuit. We first need a new
quantum gate.

10
4 FOURIER TRANSFORMS 4.3 Quantum Fourier Transform

Definition 16. A rotation gate is a unitary operator defined as

) *
1 '0 (
Rs := (30)
0 exp −2πi
2s

Now consider the following circuit which almost gives us the desired transformation

| j1 › H R2 ... Rm ...
| j2 › ... H ... Rm–1 ...
...

| jm–1› ... H R2
| jm› ... H
Figure 4: A quantum circuit for efficient quantum Fourier transform

Applying H to |j1 ", the first qubit of |j" = |j1 " |j2 " . . . |jm ", we get
1
√ (|0" + e−2πi(0.j1 ) |1") ⊗ |j2 " . . . |jm "
2
Applying the controlled R2 to this, we get
1
√ (|0" + e−2πi(0.j1 j2 ) |1") ⊗ |j2 " . . . |jm "
2
Keep going through the circuit until we get
1
√ (|0" + e−2πi(0.j1 j2 ...jm ) |1") ⊗ |j2 " . . . |jm "
2
for the first qubit. For the second qubit, we do the same thing and get
1
√ (|0" + e−2πi(0.j1 j2 ...jm ) |1") ⊗ (|0" + e−2πi(0.j2 ...jm ) |1") ⊗ |j3 " . . . |jm "
22
and so on until the mth qubit, after which we have
1 2 3 2 3 2 3
√ |0" + e−2πi(0.j1 j2 ...jm ) |1" ⊗ |0" + e−2πi(0.j2 ...jm ) |1" ⊗ . . . ⊗ |0" + e−2πi(0.jm ) |1"
2m
which is almost what we wanted, except in the reverse order! To remedy this, we add 1m/22 swap
gates at the end of the circuit.
!mWe can count the number of gates in the circuit. From bottom up, we have 1 + 2 + . . . + m =
j=1 = m(m + 1)/2 Hadamard gates and controlled rotations gates. In addition, we have the 1m/22
swap gates we put in at the end. Hence, the circuit is polynomial in m. This is an exponential speed
up over the classical FFT! However, since the QFT acts on quantum states, we can’t just apply the
QFT to data sets as with the DFT. Moreover, we could only construct this when we had N = 2m .
In general, we can construct QFT in polynomial time only if N is smooth. There are ways to get
around this[7] , but we won’t cover them.

11
 5 QUANTUM PART: ORDER FINDING ALGORITHM

5 Quantum Part: Order Finding Algorithm

The results of section 3 reduced the problem of factoring n to finding the order r of an element
a in (Z/nZ)× . Shor provides such an algorithm by employing the QFT. Although there’s an
equivalent way of doing this using a quantum circuit via phase estimation5 , we will be following
Shor’s approach[11] . Readers interested in the alternative approach can consult our reference for the
phase estimation algorithm[2] .
For the following section, we will assume that n is an composite odd integer which is not a power
of prime (the algorithm fails otherwise). If n is even, we can just factor out all the powers of 2 until
we get an odd integer, then run the algorithm on the resulting integer. We can test whether n is a
prime efficiently using classical primality tests such as the Miller-Rabin test[6],[9] and the AKS test[1] .
We can also test if n is a power of prime efficiently by taking the k th root of n until n1/k < 2.
Given n, we choose N = 2m such that n2 ≤ N < 2n2 (i.e. choose the unique power of 2 in that
range). We will be working with two registers (two arrays of qubits). The first one will hold a
number x (mod N ), the second one will hold a number mod n. Each of them holds m qubits. At
first, the registers are
|0" ⊗ |0"
We put the!first register in the uniform superposition of numbers x (mod N ) by using the QFT,
−1 −0x
|0" )→ √1N Nx=0 ζ |x". We get

N −1
1 $
√ |x" ⊗ |0" (31)
N x=0

Now suppose f (x) = ax (mod N ). Note that the period of f is the same as ord(a) = r. Given
some base a, Can we compute f (x) efficiently? The answer is yes, we can just exponentiate by
repeated squaring!
Example 17. Let a = 2, N = 15. Suppose we want to compute f (10). Naïvely, this requires 10
multiplications. However, we can apply repeated squaring

210 = (22 )5
(22 )5 = (2(22 )2 )2

and thus

22 = 2 · 2
24 = 22 · 22
25 = 24 · 2
210 = 25 · 25

which requires 4 multiplications instead of 10. Notice that if we were calculating f (20), we would
only need 5 instead of 20 multiplications.
We need to apply f to the contents of the first register and store the result of f (x) in the second
register. To do so, we can construct f as a quantum function[11] . It turns out that this is the
bottleneck of the algorithm since implementing f on a quantum computer requires a lot of quantum
gates. Still, Shor’s algorithm is much faster than factoring on a classical computer. We have then
5
I like this way much better, however, I couldn’t include it in the main text due to time constraints.

12
5 QUANTUM PART: ORDER FINDING ALGORITHM

N −1
1 $
√ |x" ⊗ |f (x)" (32)
N x=0

Apply the QFT to the first register, we get

N −1 N −1
FN ⊗ I $ 1 $
√ |x" ⊗ |f (x)" = √ (FN |x") ⊗ |f (x)"
N x=0 N x=0
N −1
(33)
1 $ −xy
= ζ |y" ⊗ |f (x)"
N
x,y=0

We perform a measurement and compute the probability that we get a particular state |y" |f (k)",
where 0 ≤ k < r. Summing over all the possibilities,
7 72 7 72
7 N −1 7 7 7
71 $ 7 7 1 $ 7
7 #y| #f (k)| ζ −xy 7 7
|y" |f (x)"7 = 7 ζ −xy 7
(34)
7N 7
7 x,y=0 7 7N 7 x:f (x)≡f (k)

The sum is over all x, 0 ≤ x < N such that f (x) ≡ f (k) (mod n) i.e. ax ≡ ak (mod n). Since
ord a = r, this is equivalent to summing over all x such that x = k mod r. Writing x = br + k, the
probability is then
7 72 7 72
7 )(N −k−1)/r* 7 7 )(N −k−1)/r* 7
71 $ 7 7 $ 7
7 ζ −(br+k)y 7
= 71 ζ −bry 7
(35)
7N 7 7N 7
7 b=0 7 7 7 b=0

since |ζ −ky |2 = 1. Moreover, since ζ −bry+N = ζ −bry , we can reduce ry modulo n. Replace ry by
{ry}, where N/2 ≤ {ry} ≤ N/2. We can approximate the sum inside by an integral. So
)(N −k−1)/r*
$ 8 )(N −k−1)/r* ) *
1 −b{ry} 1 −2πib{ry}/N 1
ζ 3 e db + O (36)
N N b=0 N
b=0

Let u = rb/N, du = dbr/N , we get

" #
8 r (N −k−1) ) *
1 N r
−2πiu{ry}/r 1
e du + O (37)
r u=0 N
Now since k < r, approximating the upper bound of the integral by 1 will only give us an error of
O( N1 ) and so we get
8 ) *
1 1 −2πiu{ry}/r 1
e du + O (38)
r u=0 N

Now suppose that − 12 ≤ {ry}r ≤ 2 . We can integrate the function and see that the integral is
1

minimized when {ry}r = ± 2 . The integral will evaluate to πr if this were the case. However, this
1 2

happens if there exists a constant d such that

7 7
r r 7y d 7 1
− ≤ ry − dN ≤ ⇔ 77 − 77 ≤ (39)
2 2 N r 2N

13
 6 DISCRETE LOGARITHMS

This looks familiar! It is in fact the error bound for the best approximation of Ny . That is, we
want to find a best approximation of dr such such r < n. There is at most one such fraction since
N > n2 . We can compute dr by computing the continued fraction of Ny and truncate where necessary.
If dr is in its lowest terms and gcd(d, r) = 1, then we get r and can use it for the rest of the algorithm,
which is done classically. If not, the algorithm fails.
There are φ(r) numbers relatively prime to r. Moreover, there are r values for ak since ord(a) = r.
Hence, there are rφ(r) states which allows us to obtain r, and each state occurs with probability of
at least 3r12 . Therefore, we will get r with probability at least φ(r) φ(r)
3r . Since r > log log r for some
C

constant C [4] , we can repeat the algorithm O(log log r) times and almost guarantee that we find r.

6 Discrete Logarithms
Just as the RSA cryptosystem is based off the presumed difficulty of factoring a number classically,
the Diffie-Hellman key exchange protocol is based off the presumed difficulty of computing the
discrete logarithms efficiently. We will consider how to apply the ideas we developed in the previous
sections to compute discrete logarithms. We will only treat the special case when p − 1 is smooth
(i.e. its prime factors are all less than logC p for fixed C) and refer the reader to Shor’s paper for the
general case[11] . The general case is a bit more technical but contains the same ideas.
Let x ≡ g r (mod p). We want to compute r given x, g, p. Note that f (a, b) = g 0 = 1 only if
a ≡ −rb (mod p − 1). We start out with three registers all initialized to |0".

|0" ⊗ |0" ⊗ |0"

We can apply Fp−1 to the first two registers (i.e. apply Fp−1 ⊗ Fp−1 ⊗ I) to obtain

p−2
1 $
|a" |b" |0" (40)
p−1
a,b=0

Now suppose f (a, b) = xa g −b (mod p). We can compute f (a, b) efficiently by repeated squaring
as before. Put the result in register 3 and we get

p−2
1 $
|a" |b" |xa g −b " (41)
p−1
a,b=0

Apply the QFT to the first two registers again, we get

p−2
$
1
|ψ" = ζ −ac ζ −bd |c" |d" |xa g −b " (42)
(p − 1)2
a,b,c,d=0

2 3
where ζ = exp 2πi
p−1
We perform the measurement and compute the probabilities that we get the a particular state

14
REFERENCES

|c" |d" |g k ".

7 72
7 72 77 1 $
7
7
7 7
7#c| #d| #g |ψ"7 = 77
k
2
ζ −(ac+bd) 7
7 (43)
7 (p − 1) a,b:a−rk=b 7
7 7
7 1
p−2
$ 72
7 −((b+rk)c+bd) 7
=7 ζ 7 (44)
7 (p − 1)2 7
b=0
7 72
7 1
p−2
$ 7
7 −(brc+bd) 7
=7 ζ 7 (45)
7 (p − 1)2 7
b=0

if d + rc +≡ 0 (mod p − 1), the sum is over all the (p − 1)st roots of unity and hence 0. If d + rc ≡ 0
(mod p − 1), we get
7 72 7 7
p−2 72
7 1
p−2
$ 7 7 1 $
7 7 7 7
7 2
ζ −(brc+bd) 7 = 7 2
17 (46)
7 (p − 1) 7 7 (p − 1) 7
b=0 b=0
1
= (47)
(p − 1)2

Hence, we only need to measure pairs (c, d) such that rc + d ≡ 0 (mod p − 1). Then, we can
then recover r by computing r ≡ −c−1 d (mod p − 1). The algorithm will fail unless a and p − 1 are
relatively prime. The probability of success is, as with factoring numbers, φ(p−1)
p−1 > log log p−1 .
C

7 Conclusion
We showed that Shor’s algorithm allows us to factor numbers much more quickly than classical
algorithms. It runs in O(log3 n) where n is the number we are trying to factor. With Shor’s algorithm,
factoring becomes a BQP problem since we have a bounded probability of failure on each run of the
algorithm. By applying the algorithm multiple times, we can be more and more sure to factor n.
The main bottleneck of the algorithm is implementing modular exponentiation using a quantum
circuit. The algorithm employs the QFT, which is basically a quantum version of the FFT, in a
crucial way. In addition to factoring numbers, Shor’s ideas also allowed us to compute discrete
logarithms in polynomial time. Shor’s algorithm is a real quantum algorithm which allows us to test
the abilities of quantum computers. So far, we’ve only been able to factor numbers up to 21 using
the algorithm since it requires coherent control of many qubits.

References
[1] Manindra Agrawal, Neeraj Kayal, and Nitin Saxena. “PRIMES is in P”. In: Ann. of Math 2
(2002), pp. 781–793.
[2] Donny Cheung. “Using Generalized Quantum Fourier Transforms in Quantum Phase Estimation
Algorithms”. MA thesis. Waterloo, Ontario, Canada: University of Waterloo, 2003.
[3] R. Cleve et al. “Quantum algorithms revisited”. In: Proceedings of The Royal Society (1998).
[4] G. H. Hardy and E. M. Wright. An introduction to the theory of numbers (Sixth ed.) Oxford
University Press, 2008.

15
REFERENCES REFERENCES

[5] A. Ignjatović. “Polynomial Multiplication and The Fast Fourier Transform (FFT)”. 2013. url:
http://www.cse.unsw.edu.au/~cs3121/Lectures/Topic3.pdf.
[6] Gary L. Miller. “Riemann’s Hypothesis and Tests for Primality”. In: Journal of Computer and
System Sciences 13.3 (1976), pp. 300–317.
[7] Michele Mosca and Christof Zalka. “Exact quantum Fourier transforms and discrete logarithm
algorithms”. In: Symposium on Foundations of Computer Science (1994). arXiv:quant-ph/
0301093 [quant-ph].
[8] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information.
Cambridge University Press, 2000.
[9] Michael O Rabin. “Probabilistic algorithm for testing primality”. In: Journal of Number Theory
12.1 (1980), pp. 128 –138. issn: 0022-314X. doi: 10.1016/0022- 314X(80)90084- 0. url:
http://www.sciencedirect.com/science/article/pii/0022314X80900840.
[10] Peter W. Shor. “Algorithm for Quantum Computation: Discrete Logarithms and Factoring”.
In: Symposium on Foundations of Computer Science (1994).
[11] Peter W. Shor. “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms
on a Quantum Computer”. In: SIAM J.Sci.Statist.Comput. (1997).
[12] John A Smolin, Graeme Smith, and Alex Vargo. “Pretending to factor large numbers on a
quantum computer”. In: Pre-Print (2013). arXiv:1301.7007 [quant-ph].
[13] Lieven M.K Vandersypen et al. “Experimental Realization of Shor’s Quantum Factoring
Algorithm Using Nuclear Magnetic Resonance”. In: Letters to Nature (2001).
[14] Ruye Wang. Convolution theorem for Discrete Periodic Signal. Apr. 2013. url: http : / /
fourier.eng.hmc.edu/e180/e101.1/e101/Fourier_Transform_D/node9.html.

16