10
From Circuits to Algorithms

In the last chapter we became comfortable with putting together gates to create circuits for simple algorithms. We’re now ready to look at more advanced quantum algorithms and considerations on how and when to use them.

Our target in this chapter is Peter Shor’s 1995 algorithm for factoring large integers almost exponentially faster than classical methods. To get there we need more tools, such as phase estimation, the Quantum Fourier Transform, and function period finding. These are important techniques in their own rights, but are necessary in combination for quantum factoring.

We also return to the idea of complexity that we first saw for classical algorithms in section 2.8. This allows us to understand what ‘‘almost exponentially faster’’ means.

This chapter contains more mathematics and equations for quantum computing than we have encountered previously. I recommend that you take the time to understand the linear algebra and complex number computations throughout. While they may appear daunting at first, the techniques are used frequently in quantum computing algorithms.

Topics covered in this chapter

10.1 Quantum Fourier Transform

The Quantum Fourier Transform (QFT) is widely used in quantum computing, notably in Shor’s factorization algorithm in section 10.6 . If that weren’t enough, the Hadamard H is the 1-qubit QFT and we’ve seen many examples of its use.

Most treatments of the QFT start by comparing it to the classical Discrete Fourier Transform and then the Fast Fourier Transform. If you don’t know either of these, don’t worry. I’m presenting the QFT in detail for its own sake in quantum computing. Should you know about or read up about the classical analogs, the similarities should be clear.

10.1.1 Roots of unity

We are all familiar with square roots. For example, √4 is equal to either 2 or −2. We can also write √2 = 2^1/2 and say that there are two ‘‘2nd-roots of 2.’’ Similarly, 5 is a cube root, or ‘‘3rd-root,’’ of 125. In general, we talk about an ‘‘Nth-root’’ for some natural number N. When we consider the complex numbers, we have a rich collection of such Nth-roots for 1.

For N = 1, there is only one first root of unity, and that is 1 itself.

When N = 2, there are two second roots of unity: 1 and −1. −1 is primitive.

With N = 3, we can start to see a pattern. For this, remember Euler’s formula

and |e^{ϕ i}| = 1. If ϕ = 2π, we go all the way around the unit circle and back to 1. If ϕ = 2π/3 then we go only one-third of the way. Rotating another 2π/3 radians we get to two-thirds around. The third roots of unity are

and ω₁ and ω₂ are primitive.

For N = 4 we can do something similar, though the situation is much simpler.

Notice something else here: for the first time, we see overlaps in the collection of roots of unity for different N. This happens when the greatest common divisor of the different N is greater than 1.

For N = 5 and above we continue in this way.

When there is a fraction in the exponent of such expressions, it is common to pull the i into the numerator, as I have done here.

For some N as we have considered, let ω be an Nth root of unity. What is ω⁻¹? This one is easy because

Complex conjugation to the rescue again!

The polynomial x^N − 1 can be factored as

Every Nth root of unity satisfies x^N − 1 = 0. The first factor (x − 1) is 0 only for ω = 1. Therefore, every other Nth root of unity ω ≠ 1 satisfies We can re-express this as Taking it a step further, let . Then with the second equality holding because we are using ω instead of ω as the Nth root of unity. Thought of another way, we are adding together the same set of roots of unity, just in a different order. In these equivalent versions, this is called the summation formula. To see this, let This is a primitive Nth root of unity and for any 0 < j < N, is another Nth root of unity ≠ 1. From the summation formula, When j = 0, we are simply adding 1 to itself N times.

10.1.2 The formula

As we have seen, a general quantum state |ϕ⟩ on n qubits can be written

for N = 2ⁿ. There are N amplitudes a_j corresponding to the N standard basis kets |j⟩. For a fixed |ϕ⟩, we get a complex-valued function where a(j) = a_j. We also know .

Definition

The quantum Fourier Transform of |ϕ⟩ is

where We can simplify this by letting ω = e^{2π i/N}, which is a primitive Nth root of unity. When we have 1 qubit, N = 2 and ω = −1. So For |ψ⟩ = |0⟩, a₀ = 1 and a₁ = 0. So and . For |ψ⟩ = |1⟩, a₀ = 0 and a₁ = 1. So and . All together,

This is none other than H!

So QFT₁ = H. Is QFT_n = H^⊗n?

No.

Matrix

From its definition, you can see that QFT_n has matrix

Since ω is an Nth root of unity, we can simplify several of the exponents. For example,

• ω^{(N−1)(N−1)} = ω^{N2 − 2N + 1} = ω because ω^N2 = (ω^N)^N = 1^N = 1.
• ω^kN = ω^N = 1 for k in Z.
• ω^2(N−1) = ω^N−2.
• ω⁻¹ = ω^N−1.

Applying rules such as these, we get

For n = 1, this is as previously noted. But when n = 2, N = 4 and ω = i, So QFT_n ≠ H^⊗n, in general.

Recursive matrix

The higher Hadamard matrices H^⊗n are defined recursively in terms of the lower ones by

Is there a similar decomposition for QFT_n? Yes, but there are other matrix factors. Namely, where

• N = 2ⁿ,
• I_N is the N by N identity matrix,
• ω is the primitive Nth root of unity e^{2π i/N},
• Ω_N is the diagonal matrix
  
  and
• P_2ⁿ⁺¹ is a shuffle transform defined by setting its j, k entry via the formula
  
  Remember that the first row and column index in the matrix is 1.

We do not derive this here but its effect is to move the vector entries with odd numbered indices to the front, followed by the even ones. For example,

where Given this recursive decomposition, we can break a QFT down into smaller and smaller gates. This is seen in the relatively simple structure of the circuit.

10.1.3 The circuit

The QFT₃ circuit can be constructed entirely from H, S, T, and SWAP gates. [20, Section 5.1]

The last SWAP is necessary to reverse the order of the qubit quantum states. We can rewrite this with explicit gates.

An = S gate is a rotation given by the matrix

An = T gate is a rotation given by the matrix Define = . So = and = . Rewriting the above circuit in terms of gives

The pattern becomes clear when we extend this to 4 qubits and QFT₄.

In this section we have discussed the Quantum Fourier Transform QFT_n on n qubits. In some circuits we need the inverse Quantum Fourier Transform QFT_n⁻¹. We get this by running QFT_n backwards. This is possible because all gates ultimately correspond to unitary operations, which are invertible and hence reversible.

10.2 Factoring

It is hard to do a web search and find ‘‘practical’’ applications of factoring integers. Many of the results say things like ‘‘factoring integers is a tool in for factoring polynomials,’’ and ‘‘factoring integers is useful for solving some differential equations.’’ These examples are fine, but it seems like math to do more math.

One area where factoring comes into play is cryptography. There are cryptographic protocols that involve not being able to easily factor large numbers, in order to keep your information safe while it travels across the Internet or is stored in databases.

In this section, we examine some of the mathematics of factoring and focus in on one famous algorithm. This is Shor’s algorithm, named after its discoverer, mathematician Peter Shor. It can factor large integers almost exponentially faster than classical methods given a sufficiently large and powerful quantum computer.

If one does not understand what ‘‘sufficiently large and powerful’’ means, it is easy to make poor predictions of whether, or when, quantum computing can cause problems for encryption.

10.2.1 The factoring problem

Most people don’t worry about factoring integers once they graduate from high school or are no longer teenagers. In this book, we first encountered factoring when we discussed integers. In

the −1 is a unit and 2, 3, and 5 are examples of prime numbers. They are irreducible in the sense that their only factors are 1 and themselves.

As of this writing, the largest confirmed prime number has 24,862,048 digits. [21] Determining whether a number is prime or factoring that number are different but related problems. The part of mathematics that studies these topics, as well as many others, is called number theory.

While it can be hard to factor large integers, it is not simply the size that makes it difficult. The number

is bigger than the above prime but is trivially factorable.

When we factor an N in Z greater than 1 we express it as a product

for some integer n ≥ 1 with all the p_j prime numbers. We sort the primes from smallest to largest, which means p_j ≤ p_k if j < k. Some primes may be repeated more than once. If, for example, 11 appears in the factorization five times then we say 11 has multiplicity 5. The notation we use to show that a p_j is a factor of N is p_j | N. Read this as ‘‘p_j divides N.’’

If n = 1 then N itself is prime. While the goal is to find all the p_j factors, we really want to start by finding one p_j and then continue by factoring the smaller integer N / p_j. As we explore factoring, we focus on that first prime factor and then re-use our methods, basic or advanced, to do what’s left.

10.2.2 Big integers

On a 64 bit classical processor, the largest natural number you can represent is

For positive and negative integers, one bit is reserved for the sign and so the largest number we can represent is These may look big, but from a factoring perspective they are small. We know all the primes that go well above these values.

What would you do if you wanted to factor 18,446,744,073,709,551,617, which is (2⁶⁴ − 1) + 2? The processor cannot hold the number directly in its hardware.

We fix this by creating ‘‘bignum’’ or ‘‘big integer’’ objects in software. These can be as large as computer memory permits. Essentially, we divide up the big integer across multiple hardware integers and then write software routines for negation, addition, subtraction, multiplication, division, and (related) quotients and remainders. Of these, the division operations are the trickiest to code correctly, though it is well known how to do it.

Scientific applications often support big integers. These include Maple, Mathematica, MATLAB^®, and cryptographic software. Python, Erlang, Haskell, Julia, Ruby, and many variants of Lisp and Scheme are programming languages which provide built-in big integers. Other languages provide the functionality via library extensions.

The GNU Multiple Precision Arithmetic Library implements highly optimized low-level routines for arbitrary precision integers and other number types. It is distributed under the dual GNU LGPL v3 and GNU GPL v2 licenses. [9]

Let’s do a quick review of classical ways of factoring integers and then move on to Shor’s approach.

10.2.3 Classical factoring: basic methods

Classical factoring techniques can be broken down into basic attempts to break a number N down into prime factors that require arithmetic, and more advanced versions that require sophisticated mathematics. We go into detail on the former and survey the latter.

If N = 1 or N = −1, we are done. If N is negative, we note the unit is −1 and assume that N is positive.

There is one and only one even prime number and that is 2. If N is even, meaning that the last digit is one of 0, 2, 4, 6, or 8, then p₁ = 2. Continue with N = N / 2. If this N is even, then p₂ = 2, and reset N = N / 2. Keep going until N is not prime. If N ever becomes 1 in the factoring processes, we are done. At this point N is odd and we may have some initial p_j factors equal to 2. To see if 3 is a factor we use a trick you may have learned while young.

We can tell if 3 is a factor of N if 3 divides the much smaller integer that is the sum of the digits. To prove this, let N be as above. Then

Convince yourself of this. If N ≡ 0 mod 3 then N is divisible by 3. By modular arithmetic, which we discussed in section 3.7,

because every 10 ^j ≡ 1 mod 3!

There are no special tricks for 4 because we already factored out any 2s. 5 | N if the last digit is 0 or 5. Because 6 is 2 × 3, there is nothing to do.

For the prime 7 we can do trial division: divide N by 7 and if there is no remainder, 7 is a factor. There is nothing to do for 8, 9, or 10. For the prime 11 it is trial division again.

What I have described here is one mode of attack: start with a list of sorted primes and see if N is divisible by each in turn. You don’t need to look at all primes < N or even ≤ N/2, you need to look at those ≤ √N.

This is not a practical method of finding factors once N gets larger because division is computationally expensive and you have to do many of them. It’s not difficult to find what is called the integer square root of N, the largest integer s such that s² ≤ N. This is ⌊√N⌋. If s² = N then N is a square.

How would we go about creating a list of primes to use in trial division? A classical way (I really mean ‘‘classical’’ – it’s from ancient Greece) is the Sieve of Eratosthenes. Let’s see how this works to find all primes less than or equal to 30. We ‘‘sift’’ through the numbers and eliminate the composite values.

We begin by creating a list of all numbers in N less than or equal to 30. We underline all the numbers and selectively remove the underlines for the non-primes. We ‘‘mark’’ a number by removing its underline. When we are done, only the primes are underlined. Since 1 is not a prime, we mark it.

The first number that is not marked is 2. It must be prime. Now mark each of 2 + 2, 2 + 2 + 2, 2 + 2 + 2 + 2 and so on. The numbers we are marking are all multiples of 2 and so are not prime. We get:

After 2, the next number that is not marked is 3. It must be prime. We mark 3 + 3, 3 + 3 + 3, etc.

It is efficient to simply add 3 to the last number we marked. We can use multiplication, but that’s more computation than we need.

After continuing in this manner, we end up with

The integers that remain are prime. This is not a terrible way to create a list of a few hundred or thousand primes, but you would not want to do this to find primes with millions of digits.

A simple Python implementation of the sieve is in Listing 10.1 . The output for the example above and a larger set of primes is  

 
print(list_of_primes(30)) 
 
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29] 
 
print(list_of_primes(300)) 
 
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 
149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 
227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293]
  

For a given number N, there are number theoretic algorithms such as the Miller-Rabin test if N is prime, or at least prime with high probability.

A further technique we can try is testing whether N is a power of another integer. Using an efficient algorithm called Newton’s method, we can determine whether N is an mth power. The most efficient method to do this test does so in O(log(N)) time. [1]

Which m do we need to try?

If 2^m > N, then any other positive integer raised to the mth power is also too big. Therefore we need 2^m ≤ N. By taking base 2 logarithms and noting that both sides are greater than one, we need to look at all m in Z such that

For N = 10⁵⁰, for example, a maximum of m = 167 suffices.

Another basic method is due to mathematician Pierre de Fermat. We can always represent N as the difference of two squares N = u² − v². The importance of this is that we then have

If we can find choices of u and v so that neither u + v nor u − v are 1, then we have a factorization.

---

---

Another way of saying the above is

Since N is odd, N − 1 and N + 1 are both even and so one possible pair of choices is u = (N + 1)/2 and v = (N − 1)/2. Suppose N = 87. Then setting yields u + v = 87 and u − v = 1. This choice of u and v does not help us. These particular choices for finding a pair never works because In any case, since 87 is prime we cannot find a good pair which works.

So is there a better alternative or approach?

Start by computing the integer square root s of N. If s² = N we have found a factor. Otherwise, set u = s + 1 and consider u² − N. If this is a square, then it is v². If it is not, increment u by 1 and try again. Since N = (u + v)(u − v), we can stop when u + v ≥ n. (In fact, we can stop a little sooner, but this shows there is point after which we need not continue.)

Let’s try N = 143. s = 11 and since s² ≠ 144, set u = s + 1 = 12. u² − N = 144 − 143 = 1 and so set v = 1. v is a perfect square and so

This works well when a u + v and u − v are both close to √N. When they are not, we need to do many iterations and computations.

10.2.4 Classical factoring: advanced methods

In Fermat’s method, we tried to find u and v such that

This is called a congruence of square modulo N. If we had such a congruence, we could create the factorization This does not say (u + v)(u − v) = N but rather that there is a non-zero integer c such that (u + v)(u − v) = c N.

Given this, let g = gcd(u + v, N). Then g and N/g is a possible factorization of N. If either of these is 1 then we have failed and we have to find new candidates for u and v.

The key question is: how efficiently can we find u and v that are likely to work?

In Fermat’s method we work outwards from the integer square root of N, hoping to find a factor near it. Let’s briefly look at more advanced number theoretic approaches.

The following is known as Dixon’s method after its discoverer, mathematician John D. Dixon. For a given N to factor, choose B and hence P. Let n be the number of primes in P. We want to find n + 1 distinct positive integers f_j such that j² mod N is B-smooth. Said another way,

for non-negative integer exponents e_j,k. How do we find these f_j? We choose them randomly between 1 + the integer square root of N and N. Suppose N = 17. The the integer square root of N is 4. In Python, the code

 
#import random 
 
random.randint(5, 17)
  

returns a random integer greater than or equal to 5 and less than 17.

Now that we have our collection of f_j, look at the product

By linear algebraic methods on the matrix we can find replacement values for the exponents so that the congruence still holds and e_{j, 1} +···+ e_{j, n} is even. This means the right-hand side of is a square, as is the left. This equation has the form and we can proceed to test if we have a factor. If we don’t we can

• look for different f_j, or
• increase B and keep trying.

Full treatments of the algorithm put bounds on the value of B. Can we find better congruences of squares, faster? The f_j we use are smooth numbers. Can we sift through all the possibilities better so that we need not consider as many bad candidates?

10.3 How hard can that be, again

In section 2.8 we first saw and used the O( ) notation for sorting and searching. Bubble sort runs in O(n²) time and merge sort is O(n log(n)). A brute force search is O(n) but adding sorting and random access allows a binary search to be O(log(n)). All algorithms are of polynomial time because we can bound them, in this case, by O(n²). More precisely, there is a hierarchy of time complexities. For the examples above,

Polynomial time is higher than all of them but we can say each runs in at least polynomial time.

We are concerned with this because there is a special distinction between polynomial time and exponential time. The latter describes exponential growth and can cause problems to become quickly intractable. All these descriptions apply to large n. Even though an exponential running time sounds bad, the algorithm may still be feasible when n is small.

Something running in exponential time has its execution time proportional to 2 ^f(n). Here f(n) is a polynomial in n such as n² − n + 1. Double exponential time is even worse: the time is proportional to 2 ^{2 f(n)}.

Subexponential is an improvement over exponential. It means for every b > 1 in R, the running time is less than bⁿ. Set b = 2 to compare with exponential time. [13]

For example, an algorithm that runs in is subexponential.

There are many subtleties concerning classical and quantum complexity, and a full discussion involves traditional and quantum Turing machines. As our need is only to understand where quantum computing can be better than classical computing, I don’t delve into a full and formal description of complexity theory here.

Fermat’s method of factoring integers is of exponential time. The general number field sieve is subexponential. There are no known polynomial time classical algorithms for factoring.

When you learned to add and multiply, you learned deterministic algorithms for each. Given the same input numbers, you executed a series of steps that always led to the same answer. No randomness was involved. If you encounter random numbers as you continue through this book, chances are you are looking at a non-deterministic algorithm.

A non-deterministic algorithm may involve probability and choices that could ultimately produce different answers. Even if the answers are the same, the choices made inside the algorithm may cause the algorithm to run quite quickly or excruciatingly slowly. More formally, a bad choice can cause an exponential time run while a good one can enable polynomial time execution.

The complexity class P refers to the collection of problems that can be solved in polynomial time using a deterministic algorithm.

The class NP is the collection of problems whose solutions can be checked in polynomial time. We do not know if P = NP but you could become wealthy via industry prizes if you could prove it or its contrary. Checking a solution by asking a question like ‘‘Is p × q equal to N?’’ is an example of a decision problem.

Consider integer factorization: we do not know if there is a classical way of factoring a composite integer in polynomial time but, once done, we can easily multiply together the factors to see if the solution is correct. Otherwise said, we do not know if classical factorization is in P but we know it is in NP.

If P ≠ NP we can represent the set of problems as shown below. Not all NP problems are equally difficult. They are defined by their ability to check solutions in polynomial time.

Two new terms are shown: NP-Hard and NP-Complete. The first is easier to explain: the problems in NP-Hard are at least as hard as the most difficult ones in NP.

There is a valuable notion of reducing one problem to another. If A and B are two problems, we can ‘‘reduce A to B in polynomial time’’ if by a series of manipulations we can express A as an instance of B.

As an analogy, suppose one problem (A) is ‘‘travel to the Empire State building in New York City.’’ For problem B, we express it as ‘‘drive to any location in New York City.’’ B in general is hard because of the need to get to any location and the difficulty in driving to it. A can be reduced to B because we need to get to one location somehow.

A problem is in NP-Hard if any NP problem can be reduced to it. As you can see from the diagram, the class of NP-Hard problems intersects NP but is larger. The intersection is the NP-Complete class, the hardest problems in NP to which all others can be reduced. An NP-Hard problem that is not in NP cannot be solved or verified in polynomial time.

There are many, many (many, many, …) complexity classes that you may find interesting enough to learn more about beyond what we have discussed. BPP is the class of decision problems that can execute in polynomial time but with a caveat. BPP stands for ‘‘bounded-error probabilistic polynomial time.’’ Since probability is involved, an algorithm that solves this problem may not always return the correct answer. The ‘‘bounded’’ part of the name is made concrete by insisting the correct answer is gotten at least two-thirds of the time.

The class P is in BPP but it is not known if all problems in NP are in BPP.

The above classifications refer to problems solved on a classical computer. For quantum computers, we consider the BQP class, those solvable in bounded-error polynomial time on a quantum computer. BQP is short for ‘‘bounded-error quantum polynomial time.’’ Loosely speaking, BQP corresponds to quantum computers as BPP does for classical computers.

However, a stronger relationship exists. Any problem in BPP is also in BQP. This follows from our being able to perform classical operations on a quantum computer, though it might be extremely slow. We examined this in section 9.3 and subsection 9.3.1 when we saw that nand is universal for classical computing and we can create a nand on a quantum computer via a quantum Toffoli gate.

10.4 Phase estimation

Let U be an N by N square matrix with complex entries. From section 5.9 , the solutions λ of the equation

are the eigenvalues {λ₁, λ₂, …, λ_N} of U. Some of the λ_j may be equal. If a particular eigenvalue λ shows up k times among the N, we say λ has multiplicity k.

Each eigenvalue λ_j corresponds to an eigenvector v_j so that

can take each v_j to be a unit vector.

When U is unitary, we can say even more: each λ_j has absolute value 1 and so can be represented as

This is slightly different notation than we’ve used before, but it is common. It is the product 2π ϕ that is the full radian measure of the rotation.

Now we are ready to pose the question whose solution we outline in this section:

Some authors call the eigenvector |ψ⟩ an eigenket. I don’t in this book.

We are looking for a fractional approximation to ϕ. By increasing the number of qubits we use in the circuit, we can get a closer fit to ϕ.

Suppose U, |ψ⟩, and e^{2π ϕ i} are as above. By applying U multiple times

If k = 2 ^j then .

We use this in the circuit for phase estimation in the following way.

This is a controlled-U^k gate. The indicates |ψ⟩ is represented on an n qubit register and U^k takes n quantum state inputs and outputs. This U^k is only run if the input on the wire above it is |1⟩.

Let U be a single qubit gate that operates like

Then up to differing only by a global phase, these two circuits are the same:

Note how the value of the control qubit is modified.

This extends to larger matrices acting as a unit multiple of the identity and is exactly the case for how U operates on its eigenvector |ψ⟩.

Let m be the number of bits we need to get us as close as we need to the original phase. [20, Section 5.2] We use m qubits on the upper portion of the circuit, the ‘‘upper register.’’

When m = 2, our phase estimation circuit is

In general, we keep adding gates on the bottom up to U^2m−1 as is shown in Figure 10.1. By the nature of application of the phase shifts described above, these modify the quantum states in the upper register. We then take the inverse Quantum Fourier Transform.

Continuing with our 2-qubit case, the state of the upper register before applying QFT₂⁻¹ is

I’ve highlighted two of the terms in the exponents to show the pattern.

In the 2-qubit case,

because ω = e^(2πi)/22 = i is a primitive fourth root of unity. Since this matrix is unitary, its inverse is its conjugate transpose For N = 2² and its inverse Quantum Fourier Transform is for ω = e^2πi/N, a primitive Nth root of unity. The Greek letter ‘‘χ’’ is pronounced ‘‘chı’ In our circuit, a_k = (1/2) e^{2πϕk i}. Putting all this together, the final upper register before measurement is

Obvious, right? Let me state why we have some of the elements in these messy expressions that we do:

• ϕ is the phase we want to estimate.
• The |j⟩₂ are the standard basis kets.
• The sums involving the |j⟩₂ are from superpositions.
• The exponential forms involving e are roots of unity or other complex numbers with absolute value 1.
• The QFT_n⁻¹ introduces the minus sign into the exponent while using QFT_n would not have.

While the expressions look complicated, because they are, each part is there because of some gate action.

The general form of the last expression for m qubits instead of 2 is

Including the top register and the bottom register (which is |ψ⟩), the full pre-measurement state at the end of the circuit is The probability amplitudes are for the expression to the left of ⊗ |ψ⟩. Now let’s get an estimate for ϕ. We do so by finding a good rational approximation between 0 and 1. Starting with ϕ, multiply it by 2^m and choose the closest c in Z that is within 1/2 of 2^m ϕ. To be concrete, choose Let Then 0 ≤ 2^m d ≤ 1/2. We are looking for c/2^m to be a good approximation to ϕ. In fact, we want the measured state of the upper register to be |c⟩. By increasing m, the number of qubits in each of the registers, we can get a more accurate approximation. This is the crux of what the algorithm does, find |c⟩ as the measured standard basis ket in the upper register as a good approximation to 2^m ϕ.

Given c and d, we rewrite the upper register quantum state as

The probability of getting |c⟩ is the square of the absolute value of the probability amplitude when j = c. That is, If d = 0 the sum on the right is exactly 1 and the measured result of the upper register is |c⟩ with probability 1. This happens when ϕ = c/2^m is a rational number. All is not lost if d ≠ 0. In that case, the circuit returns the correct answer with probability 4/π² ≈ 0.405. [18, Section 3.7]

This means we must run the circuit multiple times to get the correct answer. If you run the circuit 28 or more times, the probability of never having gotten the correct answer is less than 10⁻⁶, as we saw in section 6.3 .

10.5 Order and period finding

Consider the function a^k on whole numbers k for a fixed a in N greater than 1. For example, if a = 3 then the first 12 values are

As we look at larger exponents k, the values of 3^k will just get larger and larger. If we instead use modular arithmetic as we saw in section 3.7 , 3^k cannot get arbitrarily large. For example, modulo M = 13, the values we get are Working modulo M = 16 yields Finally, for modulo M = 22 we get In each case, the sequence starts repeating. That is, the sequences, and hence the functions, are periodic. If we define f_a(x) = a^x mod M for a coprime to M, then the smallest positive integer r such that f_a(x) = f_a(x + r) for all x is called the period of f_a.

For M = 13 in the first example, the period r = 3. For M = 16, r = 4. In the final example with M = 22, the period r = 5.

Given such an a as above, we can look at all a¹, a², a³, … modulo M and ask: what is the smallest r in N, if it exists, such that a^r ≡ 1 mod M? If there is such an r, it is called the order of a mod M. Then a^x+r ≡ a^x mod M by multiplication by a^x. Thus r is also the period of f_a. For this reason, the ‘‘period finding problem’’ is equivalent to the ‘‘order finding problem.’’

In the rest of this section we develop a hybrid quantum-classical algorithm to find the order of such an a, with a few conditions.

This is an important algorithm because one of its applications is integer factorization. With such an r in hand, and assuming it is even,

Given a good a, a good even r, and Euclid’s algorithm, we might be able to find a factorization of M.

We use phase estimation with some post-processing using continued fractions to compute the order r of a modulo M. For 0 ≤ j ≤ r−1, we find approximations to the phases ϕ_j = j/r accurate to 2ℓ_bits + 1 bits with probability ≥ (1−ε)/r. The smaller our ε, the larger our ℓ_ε is.

I’ve chosen to use ℓ to remind us that ℓ_bits and ℓ_ε are lengths of quantum registers.

10.5.1 Modular exponentiation

Since we are finding periods of functions or orders of numbers modulo M, we need to be able to compute a^x mod M in a quantum way.

For a binary string y of length ℓ_bits, that is, y is in {0, 1}^ℓ_bits, define

Remember that a is coprime to M, which means they share no non-trivial factors. This is the same as saying gcd(a,M) = 1 and so there exists b and c in Z so that ab + Mc = 1. Looking at this modulo M, Thus a is invertible modulo M with a⁻¹ = b.

When I write an expression like |j mod M⟩ it is shorthand for |j mod M⟩_ℓbits.

The |y⟩ are the computational basis vectors in a 2^ℓ_bits-dimensional vector space over C. U is a 2^ℓ_bits by 2^ℓ_bits square matrix. The lower right 2^ℓ_bits − M by 2^ℓ_bits − M submatrix is the identity matrix. The upper left M by M submatrix is a permutation matrix because a is invertible modulo M. All other matrix entries are zero. Hence all of U is a permutation matrix consisting of 0s and 1s, and and so is unitary.

If we apply U multiple times, we get multiplication and hence exponentiation by repeated squaring:

For example,

As we saw at the end of section 9.4, for any non-negative integer z with binary representation z_k−1 z_k−2 ··· z₁ z₀ with z₀ the least significant bit,

We can define In practice we will allow k to be as large as ℓ_ε.

With these observations and calculations, we know how to do quantum modular exponentiation. Like many such quantum subroutines, we may need additional qubits for the computations outside the main circuit description. As in section 9.4, we can do modular exponentiation in O(ℓ_bits³) gates using algorithms similar to the simple classical versions.

We have one last thing to observe about U before we move on to the circuit for order finding and its analysis. With r the order of a modulo M, the kets

are eigenvectors of U for 0 ≤ j < r. To test this, we apply U to each of these expressions: Even though we use r explicitly in the above equations, remember that we do not know what it is yet! We have shown only that the eigenvalues corresponding to the eigenvectors |w_j⟩ are as j goes from 0 to r − 1. We use phase estimation to get at those eigenvalues and hence r.

10.5.2 The circuit

To set up our circuit, we need two quantum registers. Since we use phase estimation, the first register needs enough qubits to get us the accuracy we need. This is ℓ_ε. The number of qubits in the second register is ℓ_bits.

The circuit for the quantum portion of the order finding algorithm is in Figure 10.2 . The part of the circuit that I have labeled U_PE is the core of the setup for the phase estimation. Referring back to the general phase estimation circuit in Figure 10.1, U_PE needs to prepare an eigenvector in the second register and handle the controlled U^2j gates.

Unfortunately, there is nothing that allows us to create

directly. The phases ϕ_j we wish to estimate are j/r. When in a quandary in a quantum algorithm, create a superposition! Consider The 1/√r out front is the common probability amplitude of all the kets. Squaring this and multiplying by r gives us the required value of 1. Expanding the expression, where the equality involving δ_{k, 0} follows from the extended summation formula for roots of unity in subsection 10.1.1.

The value |1⟩_ℓbits is not |111···111⟩_ℓbits but is the ket |000···001⟩_ℓbits, assuming the right-most bit is the least significant.

So while we cannot prepare an individual eigenvector |w_j⟩, it is simple to prepare |1⟩_ℓbits as the circuit below.

This is the first part of U_PE for the lower quantum register. The rest of U_PE is the sequence of controlled-U^2j gates from Figure 10.1 .

Unlike the phase estimation algorithm, we are using a superposition of eigenvectors instead of a single one as the input to the lower quantum register. We’re not done, because what comes out of our circuit after applying the inverse QFT requires more work to get to r. In the next section, we use continued fractions to get us closer to the factorization.

10.5.3 The continued fraction part

We chose ℓ_ε so that for 0 ≤ j ≤ r−1 we can find approximations ϕ_j to the phases j/r of U’s eigenvalues accurate to 2ℓ_bits + 1 bits with probability ≥ (1−ε)/r. j satisfies 0 ≤ j < r where each j is as probable as another. This is another way of saying that the j we get are uniformly distributed.

The output of the phase estimation algorithm after measurement is

We’ve done a lot of quantum work and we have in hand a rational number c_j/2^ε approximation to j/r. We want to get to j/r and then r.

If gcd(c_j, 2^ε) ≠ 1 then we try again because the fraction is not in reduced form. If c_j and 2^ε are coprime, we test if

If this succeeds, we take r = 2^ε and we are done. c_j and 2^ε will be coprime if and only if c_j is odd, and there are many such candidates less than 2^ε.

We might just get a bad result with c_j/2^ε being, essentially, garbage. In this case, we can try again or decrease ε, thereby increasing ℓ_ε.

We now employ this result about approximations and continued fractions:

This statement is quite powerful: given two rational numbers, if the second is sufficiently close to the first, then the second is a convergent of the first. Reworking this for our situation gives the following result.

If c_j/2^ε is a good approximation to j/r then it is accurate to 2ℓ_bits + 1 bits with probability ≥ (1−ε)/r. What does it mean for some x to be accurate to some y for some number of bits b? Simply that

Therefore, by our definition of ℓ_bits and noting that r ≤ M. Now we start computing the convergents for the known reduced fraction c_j/2^ε. Among these will be j/r in reduced form. We test the denominators of the convergents as candidates for r. When we find one that works, we are done.

The complexity of the overall algorithm is dominated by the quantum modular exponentiation and so is O(ℓ_bits³).

10.6 Shor’s algorithm

We now have the tools we need to sketch Shor’s algorithm for factoring integers in polynomial time on a sufficiently large quantum computer.

The complete algorithm has both classical and quantum components. Work is done on both kinds of machines to get to the answer. It is the quantum portion that drops us down to polynomial complexity in the number of gates by use of phase estimation, order finding, modular exponentiation, and the Quantum Fourier Transform.

Let odd M in Z be greater than 3 for which you have already tried the basic tricks from subsection 10.2.3 to check it is not a multiple of 3, 5, 7, and so on. So that you don’t waste your time, you should also try trial division using a small list of primes, though this is not necessary. It is necessary, however, to make sure that M is not a power of a prime number, and you can use Newton’s method to test this.

So M is an odd positive number in Z that is not a power of a prime. It has a reasonable chance of being composite.

The following is the general approach to Shor’s algorithm given M as above:

1. Choose a random number a such that 1 < a < M. Keep track of these values since we might need to repeat this step again.
2. Check if gcd(a, M) = 1. If not, we have found a factor of M and we are done. This is pretty unlikely but now we know that a and M are coprime: they have no integer factors in common.
3. Now find the non-zero order r of a mod M. This means that a^r ≡ 1 mod M. If r is odd, go back to step 1 and try again with a different a.
4. If r is even, we have

   ar ≡ 1 mod M	⇒ ar − 1 ≡ 0 mod M
   	⇒ (ar/2 − 1) (ar/2 + 1) ≡ 0 mod M

5. Now look at gcd(a^r/2 − 1, M) and, if necessary, gcd(a^r/2 + 1, M). If either of these does not equal 1, we have found a factor of M and we have succeeded.
6. If both of these greatest common divisors are 1, we repeat all the above from step 1 with another random a. We continue to do this until we find a factor.

In 2001, IBM Research scientists demonstrated factoring the number 15 via Shor’s algorithm on a 7-qubit NMR quantum computer. [28]

10.7 Summary

In this chapter we did some hard math to ultimately understand how to factor integers much faster than we can classically. Along the way we went deeply into several non-trivial quantum algorithms that are used in other quantum applications. These algorithms include the Quantum Fourier Transform, phase estimation, and order finding. These form a good basis for you to understand other quantum algorithms and their circuits.

We next turn our attention to the connections between the slightly abstract concepts we have seen so far and the physical quantum computers that we can build today.

References

[1]

Daniel J. Bernstein. ‘‘Detecting Perfect Powers in Essentially Linear Time’’. In: Mathematics of Computation 67.223 (July 1998), pp. 1253–1283.

[2]

David M. Bressoud. Factorization and primality testing. Undergraduate Texts in Mathematics. Springer-Verlag New York, 1989.

[3]

Patrick J. Coles et al. Quantum Algorithm Implementations for Beginners. 2018. url: https://arxiv.org/abs/1804.03719.

[4]

Thomas H. Cormen et al. Introduction to Algorithms. 3rd ed. The MIT Press, 2009.

[5]

R. Crandall and C. Pomerance. Prime Numbers: a Computational Approach. 2nd ed. Springer, 2005.

[6]

Sanjoy Dasgupta, Christos H. Papadimitriou, and Umesh Vazirani. Algorithms. McGraw-Hill, Inc., 2008.

[7]

Thomas G. Draper. Addition on a Quantum Computer. 2000. url: https://arxiv.org/abs/quant-ph/0008033.

[8]

D. S. Dummit and R. M. Foote. Abstract Algebra. 3rd ed. Wiley, 2004.

[9]

Free Software Foundation. The GNU Multiple Precision Arithmetic Library. url: https://gmplib.org/.

[10]

R.W. Hamming. Numerical Methods for Scientists and Engineers. Dover Books on Engineering. Dover, 1986.

[11]

Jeffrey Hoffstein, Jill Pipher, and Joseph H. Silverman. An Introduction to Mathematical Cryptography. 2nd ed. Undergraduate Texts in Mathematics 152. Springer Publishing Company, Incorporated, 2014.

[12]

Kenneth Ireland and Michael Rosen. A Classical Introduction to Modern Number Theory. 2nd ed. Graduate Texts in Mathematics 84. Springer-Verlag New York, 1990.

[13]

Burt Kaliski. ‘‘Subexponential Time’’. In: ed. by Henk C. A. van Tilborg and Sushil Jajodia. Springer US, 2011, pp. 1267–1267. url: https://doi.org/10.1007/978-1-4419-5906-5_436.

[14]

A. Ya. Khinchin. Continued Fractions. Revised. Dover Books on Mathematics. Dover Publications, 1997.

[15]

Neal Koblitz. A Course in Number Theory and Cryptography. 2nd ed. Graduate Texts in Mathematics 114. Springer-Verlag, 1994.

[16]

S. Lang. Algebra. 3rd ed. Graduate Texts in Mathematics 211. Springer-Verlag, 2002.

[17]

Richard J. Lipton and Kenneth W. Regan. Quantum Algorithms via Linear Algebra: A Primer. The MIT Press, 2014.

[18]

N. David Mermin. Quantum Computer Science: An Introduction. Cambridge University Press, 2007.

[19]

S.J. Miller et al. An Invitation to Modern Number Theory. Princeton University Press, 2006.

[20]

Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information. 10th ed. Cambridge University Press, 2011.

[21]

Joe Palca. The World Has A New Largest-Known Prime Number. url: https://www.npr.org/2018/12/21/679207604/the-world-has-a-new-largest-known-prime-number.

[22]

A.O. Pittenger. An Introduction to Quantum Computing Algorithms. Progress in Computer Science and Applied Logic. Birkhäuser Boston, 2012.

[23]

Eleanor Rieffel and Wolfgang Polak. Quantum Computing: A Gentle Introduction. 1st ed. The MIT Press, 2011.

[24]

Rich Rines and Isaac Chuang. High Performance Quantum Modular Multipliers. 2018. url: https://arxiv.org/abs/1801.01081.

[25]

Lidia Ruiz-Perez and Juan Carlos Garcia-Escartin. Quantum arithmetic with the Quantum Fourier Transform. May 2017. url: https://arxiv.org/abs/1411.5949.

[26]

Robert Sedgewick and Kevin Wayne. Algorithms. 4th ed. Addison-Wesley Professional, 2011.

[27]

Peter W. Shor. ‘‘Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer’’. In: SIAM J. Comput. 26.5 (Oct. 1997), pp. 1484–1509.

[28]

Lieven M. K. Vandersypen et al. ‘‘Experimental realization of Shor’s quantum factoring algorithm using nuclear magnetic resonance’’. In: Nature 414.6866 (2001), pp. 883–887.

[29]

John Watrous. CPSC 519/619: Introduction to Quantum Computing. url: https://cs.uwaterloo.ca/~watrous/LectureNotes/CPSC519.Winter2006/all.pdf.