Skip to main content
IBM Quantum Platform

Shor's algorithm

For this Qiskit in Classrooms module, students must have a working Python environment with the following packages installed:

  • qiskit v2.1.0 or newer
  • qiskit-ibm-runtime v0.40.1 or newer
  • qiskit-aer v0.17.0 or newer
  • qiskit.visualization
  • numpy
  • pylatexenc

To set up and install the packages above, see the Install Qiskit guide. In order to run jobs on real quantum computers, students will need to set up an account with IBM Quantum® by following the steps in the Set up your IBM Cloud account guide.

This module was tested and used three seconds of QPU time. This is an estimate only. Your actual usage may vary.

# Uncomment and modify this line as needed to install dependencies
#!pip install 'qiskit>=2.1.0' 'qiskit-ibm-runtime>=0.40.1' 'qiskit-aer>=0.17.0' 'numpy' 'pylatexenc'

Intro

In the early 1990s, there was growing excitement around the potential of quantum computers to solve problems that were hard for classical computers. A few talented computer scientists had devised algorithms that demonstrated the power of quantum computing for some niche, contrived problems, but nobody had found a single, "killer app" of quantum computing that was sure to revolutionize the field. That was, until 1994, when Peter Shor came up with what is now called Shor's algorithm for factoring large numbers.

It was well known at the time that finding the prime factors of a large number was extremely difficult for a classical computer. In fact, internet security protocols relied on this difficulty. Shor found a way to find these factors exponentially more efficiently by offloading some of the more challenging steps onto a theoretical, future quantum computer.

In this module, we will explore Shor's algorithm. First, we'll give a bit more context to the algorithm, formalizing the problem it solves and explaining the relevance to cyber security. Next, we'll give a primer on modular mathematics and how to apply it to the factoring problem, showing how factoring reduces to another problem called "order-finding." We'll show how the Quantum Fourier Transform and Quantum Phase Estimation that we learned about in a previous module come into play, and how to use them to solve the order-finding problem.

Finally, we'll run Shor's algorithm on a real quantum computer! Keep in mind, though, that this algorithm will really only be useful when we have a large, fault-tolerant quantum computer, which is still some years away. So, we'll just factor a small number to demonstrate how the algorithm works.


The factoring problem

The goal of the factoring problem is to find the prime factors of a number NN. For some numbers NN, this is pretty easy. For example, if NN is even, one of its prime factors will be 2. If NN is a prime power, meaning N=pkN=p^k for some prime number pp, it is also fairly easy to find pp: we just approximate the kextthk^ ext{th} root of NN and look for nearby primes that could be pp.

Where classical computers struggle, though, is when NN is odd and not a prime power. This is the case Shor's algorithm deals with. The algorithm finds two factors pp and qq such that N=pqN=pq. It can be applied recursively until all factors are prime. In the next sections we'll see how this problem is tackled.

Relevance to cyber security

Many cryptographic schemes have been built based on the fact that factoring large numbers is hard, including one that is commonly used today, called RSA. In RSA, a public key is created by multiplying two large prime numbers together to get N=pqN = p\cdot q. Then, anyone can use this public key to encrypt data. But only somebody with the private key, pp and qq, can decrypt that data.

If NN were easy to factor, then anyone would be able to determine what pp and qq are and break the encryption. But it's not. This is a famously difficult problem. In fact, the prime factors of a number called RSA1024, which is 1024 binary digits and 309 decimal digits long, has still not been found, despite a $100,000 prize being offered for its factoring way back in 1991.


Shor's solution

In 1994, Peter Shor realized that a quantum computer could factor a large number exponentially more efficiently than a classical computer could. His insight relied on the relationship between this factoring problem and modular arithmetic. We'll go through a brief primer on modular arithmetic, then we'll see how we can use this to factor NN.

Modular arithmetic

Modular arithmetic is a system of counting that is cyclic, meaning that while counting starts in the usual way, with integers 0, 1, 2, etc., at some point, after a period NN, the counting starts over again. Let's see how this works with an example. Say our period is 5. Then, as we're counting, where we would normally hit 5, we instead start over at 0:

0,1,2,3,4,0,1,2,3,4,0,1,2,...0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, ...

This is because in the "modulo-5" world, 5 is equivalent to 0. We say that 5mod5 =05\bmod 5 \ = 0. In fact, all multiples of 5 will be equivalent to 0mod50\bmod 5.

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

Use modular arithmetic to solve the following problem:

You depart on a long, transcontinental train ride at 8 AM. The train ride is 60 hours long. What time is it when you arrive?

Answer:

The period is 24, since there are 24 hours in a day. So, this problem can be written in modular arithmetic as:

(8+60)mod(24)=20(8+60)\text{mod}(24) = 20

So, you would arrive at your destination at 20:00, or 8 PM.

ZN\mathbb{Z}_N and ZN\mathbb{Z}_N^*

It's often useful to introduce two sets, ZN\mathbb{Z}_N and ZN\mathbb{Z}_N^*. ZN\mathbb{Z}_N is simply the set of numbers that exist in a "modulo-NN" world. For example, when we were counting modulo-5, the set would be Z5={0,1,2,3,4}\mathbb{Z}_5=\{0,1,2,3,4\}. Another example: Z15={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14}\mathbb{Z}_{15} = \{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14\}. We can perform addition and multiplication (modulo NN) on the elements in ZN\mathbb{Z}_N, and the result of each of these operations also is an element in ZN\mathbb{Z}_N, making ZN\mathbb{Z}_N a mathematical object called a ring.

There's a special subset of ZN\mathbb{Z}_N that is of particular interest to us for Shor's algorithm. That is the subset of numbers in ZN\mathbb{Z}_N such that the greatest common divisor between each element and NN is 1, so each element is "co-prime" to NN. If we take the set of these numbers along with the modular multiplication operation, then this forms another mathematical object, called a group. We call this group ZN\mathbb{Z}_N^*. Turns out that with ZN\mathbb{Z}_N^* (and finite groups in general), if we pick any element aZNa \in \mathbb{Z}_N^*, and repeatedly multiply aa to itself, we'll always eventually get the number 11. The minimum number of times one must multiply aa to itself in order to get 11 is called the order of aa. This fact will be very important to our discussion of how to factor numbers below.

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

What is Z15\mathbb{Z}_{15}^*?

Answer:

Z15={1,2,4,7,8,11,13,14}\mathbb{Z}_{15}^* = \{1,2,4,7,8,11,13,14\}

We excluded the following numbers:

3:GCD(3,15)=35:GCD(5,15)=56:GCD(6,15)=39:GCD(9,15)=310:GCD(10,15)=512:GCD(12,15)=3\begin{aligned} 3: GCD(3,15)=3 \\ 5: GCD(5,15)=5 \\ 6: GCD(6,15)=3 \\ 9: GCD(9,15)=3 \\ 10: GCD(10,15)=5 \\ 12: GCD(12,15)=3 \\ \end{aligned}

What is the order of each of the elements in Z15\mathbb{Z}_{15}^*?

Answer:

The order rr is the lowest number such that armod(15)=1a^r\text{mod}(15)=1 for each element aa.

11mod(15)=1,r=124mod(15)=1,r=442mod(15)=1,r=274mod(15)=1,r=484mod(15)=1,r=4112mod(15)=1,r=2134mod(15)=1,r=4142mod(15)=1,r=2\begin{aligned} 1^1\text{mod}(15) = 1, r=1 \\ 2^4\text{mod}(15) = 1, r=4 \\ 4^2\text{mod}(15) = 1, r=2 \\ 7^4\text{mod}(15) = 1, r=4 \\ 8^4\text{mod}(15) = 1, r=4 \\ 11^2\text{mod}(15) = 1, r=2 \\ 13^4\text{mod}(15) = 1, r=4 \\ 14^2\text{mod}(15) = 1, r=2 \\ \end{aligned}

Note that, while we were able to find the order of the numbers in Z15\mathbb{Z}_{15}^*, this is NOT an easy task in general, for larger NN. This is the crux of the factoring problem and why we need a quantum computer. We'll see why as we work through the rest of the notebook.

Apply modular arithmetic to the factoring problem

The key to finding factors pp and qq such that N=pqN=pq comes down to finding some other integer xx such that

x21modNx^2 \equiv 1 \bmod N and x≢±1modN.x \not\equiv \pm 1 \bmod N.

How does finding xx help us find the factors pp and qq? Let's go through the argument now. Since x21modNx^2 \equiv 1 \bmod N, that means that x210modNx^2 - 1 \equiv 0 \bmod N . In other words, x21x^2 - 1 is a multiple of NN. So, for some integer ll,

x21=lNx^2 - 1 = l N

We can factor x21x^2 - 1 to get:

(x+1)(x1)=lN(x+1)(x-1) = l N

From our initial assumptions we know that x≢±1modNx \not\equiv \pm 1 \bmod N, so NN doesn't divide evenly into either x+1x+1 or x1x-1. So, the two factors of NN, pp and qq must each divide into x1x-1 and x+1x+1. Either pp is a factor of x1x-1 and qq is a factor of x+1x+1, or vice versa. Therefore, if we compute the greatest common divisors (GCDs) between NN and both x1x-1 and x+1x+1, that will give us the factors pp and qq. Calculating the GCD between two numbers is a classically easy task that can be accomplished, for example, using Euclid's algorithm.

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

It might be tricky to understand each step of logic above, so try working through it with an example. Use N=15N=15 and x=11x=11. First, verify that x21mod(N)x^2 \equiv 1 \text{mod}(N) and x≢±1modNx \not\equiv \pm 1 \bmod N. Then continue to verify each step. Finally, calculate GCD(11±1,15)\text{GCD}(11\pm1,15) and verify that they are the factors of 1515.

Answer:

112=12111^2 = 121, which is 158+115*8 + 1, so 112mod15=111^2\bmod 15 = 1. \checkmark

111=10 11 - 1 = 10, which is not equivalent to 0mod150\bmod 15. \checkmark

11+1=12 11 + 1 = 12, which is not equivalent to 0mod150\bmod 15. \checkmark

Now, we know that (x+1)(x1)=lN(x+1)(x-1) = l N for some integer ll. This is verified when we plug in xx and NN: (12)(10)=l15(12)(10) = l 15 when l=8l = 8. \checkmark

Now, we need to calculate GCD(12,15)\text{GCD}(12,15) and GCD(10,15)\text{GCD}(10,15).

GCD(12,15)=3GCD(10,15)=5\begin{aligned} \text{GCD}(12,15) = 3 \\ \text{GCD}(10,15) = 5 \end{aligned}

So, we found our factors of 1515!

The algorithm

Now that we've seen how finding some integer xx such that x21modNx^2 \equiv 1\bmod N helps us to factor NN, we can go through Shor's algorithm. It essentially boils down to finding xx:

  1. Pick a random integer Choose a random integer aa such that 1<a<N1 < a < N.
  • Compute GCD(a,N)\text{GCD}(a, N) classically.
    • If GCD(a,N)>1\text{GCD}(a, N) > 1, you've already found a factor. Stop.
    • Otherwise, continue.
  1. Find order rr of aa modulo NN Find the smallest positive integer rr that satisfies ar1(modN)a^r \equiv 1 \pmod N.

  2. Check if the order is even

  • If rr is odd, go back to step 1 and pick a new aa.
  • If rr is even, continue to step 4.
  1. Compute x=ar/2modNx = a^{r/2} \bmod N
  • Check that x≢1(modN)x \not\equiv 1 \pmod N and x≢1(modN)x \not\equiv -1 \pmod N.
    • If x±1(modN)x \equiv \pm 1 \pmod N, go back to step 1 and pick a new aa.
  • Otherwise, compute the gcds to extract the factors:
p=GCD(x1,N),q=GCD(x+1,N)p = \text{GCD}(x-1, N), \quad q = \text{GCD}(x+1, N)

These will be nontrivial factors of NN.

  1. Recursively factor if needed
  • If pp and/or qq are not prime, apply the algorithm recursively to factor them completely.
  • Once all factors are prime, the factoring is complete.

Based on this procedure, it might not be obvious why a quantum computer is needed to complete this task. It's necessary because step 2, finding the order of aa modulo NN, is classically a very difficult problem. The complexity scales exponentially in the number NN. But with a quantum computer, we just have to make use of Quantum Phase Estimation to solve it. Step 4, finding the GCD of two integers, is actually a pretty easy thing to do classically. So, the only step that actually needs the power of a quantum computer is the order-finding step. We say that the factoring problem "reduces" to the order-finding problem.

The hard part: order finding

Now, we'll go through how we can use a quantum computer in order finding. First, let's clarify what we mean by "order." Of course, I already told you what the order means mathematically: it's the first non-zero integer rr such that ar=1(modN).a^r = 1 \pmod N. But let's see if we can gain a little more intuition for this concept.

For small enough NN, we can just determine the order by calculating each power of aa, taking the modulus NN of that number, then stopping when we find the power rr that satisfies ar=1mod(N)a^r = 1 \text{mod}(N). That's what we did with our example, N=15N=15, above. Let's take a look at some graphs of these modular powers for some sample values of aa and NN:

Value of a to the power of k modulo N versus the power k, where a=2 and N=15. We see that as k increases, a repeating pattern emerges, showing that a^k modulo N is periodic in k. Value of a to the power of k modulo N versus the power k, where a=5 and N=21. We see that as k increases, a repeating pattern emerges, showing that a^k modulo N is periodic in k.

Notice anything? These are periodic functions! And the order rr is the same as the period! So, order-finding is equivalent to period-finding.

Quantum computers are very well suited to finding the period of functions. For this, we can use an algorithmic subroutine called Quantum Phase Estimation. We discussed QPE and its relationship to the Quantum Fourier Transform in the previous module. For a detailed refresher, go to the QFT module or John Watrous' lesson on Quantum Phase Estimation in his Quantum Algorithms course. We'll go through the gist of the procedure now:

In Quantum Phase Estimation (QPE), we start with a unitary UU and an eigenstate of that unitary ψ|\psi\rangle. Then, we use QPE to approximate the corresponding eigenvalue, which, since the operator is unitary, will be of the form e2πiθe^{2\pi i \theta}. So, finding the eigenvalue is equivalent to finding the value of θ\theta in the periodic function. The circuit looks like this:

Circuit diagram of the Quantum phase estimation procedure. The top m control qubits are prepared in superpositions with Hadamard gates, then controlled-unitary gates are applied to the bottom qubits, which are in an eigenstate of the unitary. Finally, an inverse quantum Fourier transform is applied to the top qubits and they are measured.

where the number of control qubits (the top mm qubits in the figure above) determines the precision of the approximation.

In Shor’s algorithm, we use QPE on the unitary operator MaM_a:

MayaymodN. M_a|y\rangle \equiv |ay \mod N \rangle .

Here, y|y\rangle denotes a computational basis state of the multi-qubit register, where the binary value of the qubits corresponds to the integer yy. For instance, if N=15N=15 and y=2y = 2, then y|y\rangle is represented by the four-qubit basis state 0010|0010\rangle, since four qubits are required to encode numbers up to 15. (If this concept is unfamiliar, see the introductory Qiskit in classrooms module for a refresher on binary encoding of quantum states.)

Now, we need to figure out an eigenstate of this unitary. If we started in the state 1|1\rangle, we can see that each successive application of UU will multiply the state of our register by a(modN)a \pmod N, and after rr applications we will arrive at the state 1|1\rangle again. For example with a=3a = 3 and N=35N = 35:

M31=3M321=9M331=27M3(r1)1=12M3r1=1\begin{aligned} M_3|1\rangle &= |3\rangle & \\ M_3^2|1\rangle &= |9\rangle \\ M_3^3|1\rangle &= |27\rangle \\ & \vdots \\ M_3^{(r-1)}|1\rangle &= |12\rangle \\ M_3^r|1\rangle &= |1\rangle \end{aligned}

So superpositions of the states in this cycle (ψj|\psi_j\rangle) of the form:

ψj=1rk=0r1e2πijkrak|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}

are all eigenstates of MaM_a. (There are more eigenstates than just these. But we only care about the ones of the form above.)

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

Find an eigenstate of the unitary corresponding to a=2a=2 and N=15N = 15.

Answer:

M21=2M221=4M231=8M241=1\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2^2|1\rangle &= |4\rangle \\ M_2^3|1\rangle &= |8\rangle \\ M_2^4|1\rangle &= |1\rangle \\ \end{aligned}

So, the order r=4r=4. The eigenstates we're interested in will be an equal superposition of all of the states that were cycled through above, with various phases:

ψ0=12(1+2+4+8)ψ1=12(e2πi041+e2πi142+e2πi244+e2πi348)=12(1+i24i8)ψ2=12(e2πi041+e2πi242+e2πi444+e2πi648)=12(12+48)ψ3=12(e2πi041+e2πi342+e2πi644+e2πi948)=12(1i24+i8)\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{1}{4}}|2\rangle+e^{2 \pi i \frac{2}{4}}|4\rangle+e^{2 \pi i \frac{3}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{2}{4}}|2\rangle+e^{2 \pi i \frac{4}{4}}|4\rangle+e^{2 \pi i \frac{6}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(e^{2 \pi i \frac{0}{4}}|1\rangle+e^{2 \pi i \frac{3}{4}}|2\rangle+e^{2 \pi i \frac{6}{4}}|4\rangle+e^{2 \pi i \frac{9}{4}}|8\rangle) \\ &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}

Let's say we were able to initialize our qubit state into one of these eigenstates (spoiler — we're not. Or, at least not easily. We'll explain why and what we can do instead shortly). Then we could use QPE to estimate the corresponding eigenvalue, ωj=e2πiθj\omega_j = e^{2 \pi i \theta_j} where θj=jr\theta_j = \frac{j}{r}. Then, we'll be able to determine the order rr by the simple equation:

r=jθj.r = \frac{j}{\theta_j}.

But remember, I said that QPE estimates θj\theta_j — it doesn't give us an exact value. We need the estimate to be good enough to differentiate between rr and r+1r+1. The more control qubits mm, we have, the better the estimate will be. In the problems at the end of the lesson, you'll be asked to determine the minimum mm needed to factor a number NN.

Now, we have to reconcile a problem. All of the above explanation of how to find rr begins with preparing the eigenstate ψj=1rk=0r1e2πijkrak|\psi_j\rangle = \tfrac{1}{\sqrt{r}}\sum_{k=0}^{r-1}{e^{\frac{2 \pi i j k}{r}} |a^k \rangle}. But we don't know how to do that without already knowing what rr is. The logic is circular. We need a way to estimate the eigenvalue without initializing the eigenstate.

Instead of starting with an eigenstate of MaM_a, we can prepare the initial state into the nn-qubit state corresponding to 1|1\rangle in binary (as in, 000...01|000...01\rangle). Although this state itself is obviously not an eigenstate of MaM_a, it is a superposition over all of the eigenstates ψk|\psi_k\rangle:

1=1rk=0r1ψk|1\rangle = \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle}

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

Verify that 1|1\rangle is equivalent to the superposition over the eigenstates you found for N=15N=15 and a=2a=2 in the previous check-in question.

Answer:

The four eigenstates were:

ψ0=12(1+2+4+8)ψ1=12(1+i24i8)ψ2=12(12+48)ψ3=12(1i24+i8)\begin{aligned} |\psi_0\rangle &= \frac{1}{2}(|1\rangle+|2\rangle+|4\rangle+|8\rangle) \\ |\psi_1\rangle &= \frac{1}{2}(|1\rangle+i|2\rangle-|4\rangle-i|8\rangle) \\ |\psi_2\rangle &= \frac{1}{2}(|1\rangle-|2\rangle+|4\rangle-|8\rangle) \\ |\psi_3\rangle &= \frac{1}{2}(|1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ \end{aligned}

So,

1rk=0r1ψk=12(ψ0+ψ1+ψ2+ψ3)=14(1+2+4+8+1+i24i8+12+48+1i24+i8)=14(41)=1\begin{aligned} \frac{1}{\sqrt{r}} \sum\limits_{k=0}^{r-1}{|\psi_k\rangle} &= \frac{1}{2}(|\psi_0\rangle + |\psi_1\rangle + |\psi_2\rangle + |\psi_3\rangle ) \\ &= \frac{1}{4}(|1\rangle+|2\rangle+|4\rangle+|8\rangle+|1\rangle+i|2\rangle-|4\rangle-i|8\rangle+|1\rangle-|2\rangle+|4\rangle-|8\rangle + |1\rangle-i|2\rangle-|4\rangle+i|8\rangle) \\ &= \frac{1}{4}(4|1\rangle) = |1\rangle \end{aligned}

How does this let us find the order rr? Since the starting state is a superposition over all of the eigenstates of the form listed above, then the QPE algorithm simultaneously estimates each of the θk\theta_k corresponding to these eigenstates. So, the measurement of the mm control qubits at the end will yield an approximation to the value k/rk/r where k{0,1,2,...,r1}k \in \{0,1,2,...,r-1\} is one of the eigenvalues randomly chosen. If we repeat this circuit a few times and get a few samples with different values of kk, we will quickly be able to deduce rr.


Implement in Qiskit

As we mentioned earlier, our hardware isn't at the point where we can factor huge numbers like RSA1024. We'll just factor a small number to demonstrate how the algorithm works. For this demo, we will use a simplified version of the code presented in the Shor's algorithm tutorial. If you want more details, please visit the tutorial.

We'll run the algorithm using our standard framework for solving quantum problems, called the Qiskit patterns framework. This consists of four steps:

  1. Mapping your problem to a quantum circuit
  2. Optimize the circuit to be run on quantum hardware
  3. Execute your circuit on the quantum computer
  4. Post-process the measurements

1. Map

Let's factor N=15N=15, selecting a=2a=2 to be our co-prime integer.

First, we need to construct the circuit that will implement the modular multiplication unitary, MaM_a. This is actually the trickiest part of the whole implementation, and can be very computationally expensive, depending on how it's done. For this, we'll cheat a little bit: We know that we're starting in the state 1|1\rangle, and from an earlier check-in question,

M21=2M22=4M24=8M28=1\begin{aligned} M_2|1\rangle &= |2\rangle & \\ M_2|2\rangle &= |4\rangle \\ M_2|4\rangle &= |8\rangle \\ M_2|8\rangle &= |1\rangle \\ \end{aligned}

So, we will construct a unitary that does the correct operations on these four states, but leaves all other states alone. This is cheating because we're using our knowledge of the order of 2mod152\bmod 15 to simplify the unitary. If we were actually trying to factor a number whose factors were unknown to us, we wouldn't be able to do this.

Check your understanding

Read the question(s) below, think about your answer, then click the triangle to reveal the solution.

With your knowledge of how the M2M_2 operator transforms the states above, construct the operator from a series of SWAP gates, which exchange the states of two qubits. (Hint: writing each state i|i\rangle in binary will help.)

Answer:

Let's rewrite the action of M2M_2 on the states in binary:

M20001=0010M20010=0100M20100=1000M21000=0001\begin{aligned} M_2|0001\rangle &= |0010\rangle \\ M_2|0010\rangle &= |0100\rangle \\ M_2|0100\rangle &= |1000\rangle \\ M_2|1000\rangle &= |0001\rangle \\ \end{aligned}

Each of these actions can be accomplished with a simple SWAP. M20001M_2|0001\rangle is achieved by swapping the states of qubit 00 and 11. M20010M_2|0010\rangle is achieved by swapping the states of qubit 11 and 22. And so on. So, we can deconstruct the M2M_2 matrix into the following series of SWAP gates:

M2=SWAP(0,1)SWAP(1,2)SWAP(2,3)M_2 = SWAP(0,1)SWAP(1,2)SWAP(2,3)

Remembering that operators act right to left, let's check that this has the effect that we want on each of the states:

M20001=SWAP(0,1)SWAP(1,2)SWAP(2,3)0001=SWAP(0,1)SWAP(1,2)0001=SWAP(0,1)0001=0010M20010=SWAP(0,1)SWAP(1,2)SWAP(2,3)0010=SWAP(0,1)SWAP(1,2)0010=SWAP(0,1)0100=0100M20100=SWAP(0,1)SWAP(1,2)SWAP(2,3)0100=SWAP(0,1)SWAP(1,2)1000=SWAP(0,1)1000=1000M21000=SWAP(0,1)SWAP(1,2)SWAP(2,3)1000=SWAP(0,1)SWAP(1,2)0100=SWAP(0,1)0010=0001\begin{aligned} M_2|0001\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0001\rangle \\ &= SWAP(0,1)SWAP(1,2)|0001\rangle \\ &= SWAP(0,1)|0001\rangle \\ &=|0010\rangle \checkmark \\ M_2|0010\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0010\rangle \\ &= SWAP(0,1)SWAP(1,2)|0010\rangle \\ &= SWAP(0,1)|0100\rangle \\ &=|0100\rangle \checkmark \\ M_2|0100\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|0100\rangle \\ &= SWAP(0,1)SWAP(1,2)|1000\rangle \\ &= SWAP(0,1)|1000\rangle \\ &=|1000\rangle \checkmark \\ M_2|1000\rangle &= SWAP(0,1)SWAP(1,2)SWAP(2,3)|1000\rangle \\ &= SWAP(0,1)SWAP(1,2)|0100\rangle \\ &= SWAP(0,1)|0010\rangle \\ &=|0001\rangle \checkmark \\ \end{aligned}

We can now code the circuit that is equivalent to this operator in Qiskit.

First, we import the necessary packages:

# Import necessary packages
 
import numpy as np
from fractions import Fraction
from math import floor, gcd, log
 
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit.library import QFTGate
from qiskit.transpiler import generate_preset_pass_manager
from qiskit.visualization import plot_histogram
 
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler

Then, we make the M2M_2 operator:

def M2mod15():
    """
    M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)
 
    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)
 
    U = U.to_gate()
    U.name = f"M_{b}"
 
    return U
# Get the M2 operator
M2 = M2mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M2, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

Output of the previous code cell

The QPE algorithm uses a controlled-UU gate. So, now that we have an M2M_2 circuit, we need to make it a controlled-M2M_2 circuit:

def controlled_M2mod15():
    """
    Controlled M2 (mod 15)
    """
    b = 2
    U = QuantumCircuit(4)
 
    U.swap(2, 3)
    U.swap(1, 2)
    U.swap(0, 1)
 
    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()
 
    return c_U
# Get the controlled-M2 operator
controlled_M2 = controlled_M2mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M2, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Output of the previous code cell

Now we have our controlled-UU gate. But in order to run the Quantum Phase Estimation algorithm, we'll need controlled-U2U^2, controlled-U4U^4, up to controlled-U2m1U^{2^{m-1}}, where mm is the number of qubits used to estimate the phase. The more qubits, the more precise the phase estimation will be. We'll use m=8m=8 control qubits for our phase estimation procedure. So, we need:

Ma2kya2kymodNM_{a^{2^k}}|y\rangle \equiv |a^{2^k} y \bmod N \rangle

where the index kk, with 0km1=70 \le k \le m-1 = 7, corresponds to the control qubit. Now let's calculate a2kmodNa^{2^k}\bmod N for each value of kk:

def a2kmodN(a, k, N):
    """Compute a^{2^k} (mod N) by repeated squaring"""
    for _ in range(k):
        a = int(np.mod(a**2, N))
    return a
k_list = range(8)
b_list = [a2kmodN(2, k, 15) for k in k_list]
 
print(b_list)

Output:

[2, 4, 1, 1, 1, 1, 1, 1]

Since a2kmodN=1a^{2^k} \bmod N = 1 for k2k \ge 2, all of the corresponding operators (M8M_8 and up) are equivalent to the identity. So, we only need to construct one more matrix, M4.M_4.

Note: This simplification only works here because the order of 2mod152 \bmod 15 is 44. Once k=2k=2 (so 2k=42^k = 4), every subsequent power of the operator is the identity. In general, for larger numbers NN or different choices of aa, you cannot skip constructing the higher powers. This is one reason why this is considered a toy example: the small numbers allow shortcuts that wouldn’t work for larger cases.

def M4mod15():
    """
    M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)
 
    U.swap(1, 3)
    U.swap(0, 2)
 
    U = U.to_gate()
    U.name = f"M_{b}"
 
    return U
# Get the M4 operator
M4 = M4mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(4)
circ.compose(M4, inplace=True)
circ.decompose(reps=2).draw(output="mpl", fold=-1)

Output:

Output of the previous code cell

And as before, we make it a controlled-M4M_4 operator:

def controlled_M4mod15():
    """
    Controlled M4 (mod 15)
    """
    b = 4
    U = QuantumCircuit(4)
 
    U.swap(1, 3)
    U.swap(0, 2)
 
    U = U.to_gate()
    U.name = f"M_{b}"
    c_U = U.control()
 
    return c_U
# Get the controlled-M4 operator
controlled_M4 = controlled_M4mod15()
 
# Add it to a circuit and plot
circ = QuantumCircuit(5)
circ.compose(controlled_M4, inplace=True)
circ.decompose(reps=1).draw(output="mpl", fold=-1)

Output:

Output of the previous code cell

Now, we can put it all together to find the order of 2mod152\bmod 15 with a quantum circuit, using phase estimation:

# Order finding problem for N = 15 with a = 2
N = 15
a = 2
 
# Number of qubits
num_target = floor(log(N - 1, 2)) + 1  # for modular exponentiation operators
num_control = 2 * num_target  # for enough precision of estimation
 
# List of M_b operators in order
k_list = range(num_control)
b_list = [a2kmodN(2, k, 15) for k in k_list]
 
# Initialize the circuit
control = QuantumRegister(num_control, name="C")
target = QuantumRegister(num_target, name="T")
output = ClassicalRegister(num_control, name="out")
circuit = QuantumCircuit(control, target, output)
 
# Initialize the target register to the state |1>
circuit.x(num_control)
 
# Add the Hadamard gates and controlled versions of the
# multiplication gates
for k, qubit in enumerate(control):
    circuit.h(k)
    b = b_list[k]
    if b == 2:
        circuit.compose(
            M2mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    elif b == 4:
        circuit.compose(
            M4mod15().control(), qubits=[qubit] + list(target), inplace=True
        )
    else:
        continue  # M1 is the identity operator
 
# Apply the inverse QFT to the control register
circuit.compose(QFTGate(num_control).inverse(), qubits=control, inplace=True)
 
# Measure the control register
circuit.measure(control, output)
 
circuit.draw("mpl", fold=-1)

Output:

Output of the previous code cell

2. Optimize

Now that we've mapped our circuit, the next step is to optimize it to be run on a particular quantum computer. First we need to load the backend.

service = QiskitRuntimeService()
 
backend = service.backend("ibm_marrakesh")

If you don't have time available on your account or want to use a simulator for any reason, you can run the cell below to set up a simulator that will mimic the quantum device we selected above:

pm = generate_preset_pass_manager(optimization_level=2, backend=backend)
 
transpiled_circuit = pm.run(circuit)
 
print(f"2q-depth: {transpiled_circuit.depth(lambda x: x.operation.num_qubits==2)}")
print(f"2q-size: {transpiled_circuit.size(lambda x: x.operation.num_qubits==2)}")
print(f"Operator counts: {transpiled_circuit.count_ops()}")
transpiled_circuit.draw(output="mpl", fold=-1, style="clifford", idle_wires=False)

Output:

2q-depth: 188
2q-size: 281
Operator counts: OrderedDict({'sx': 548, 'rz': 380, 'cz': 281, 'measure': 8, 'x': 6})
Output of the previous code cell

3. Execute

# Sampler primitive to obtain the probability distribution
sampler = Sampler(backend)
 
# Turn on dynamical decoupling with sequence XpXm
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XpXm"
# Enable gate twirling
sampler.options.twirling.enable_gates = True
 
pub = transpiled_circuit
job = sampler.run([pub], shots=1024)
result = job.result()[0]
counts = result.data["out"].get_counts()
plot_histogram(counts, figsize=(35, 5))

Output:

Output of the previous code cell

We see four clear peaks at 00000000, 01000000, 10000000, and 11000000, with some counts in other bitstrings due to noise in the quantum computer. We'll ignore these and keep only the dominant four by imposing a threshold: only counts above this threshold are considered a true signal above the noise.

# Dictionary of bitstrings and their counts to keep
counts_keep = {}
# Threshold to filter
threshold = np.max(list(counts.values())) / 2
 
for key, value in counts.items():
    if value > threshold:
        counts_keep[key] = value
 
print(counts_keep)

4. Post-process

For Shor's algorithm, quite a lot of the algorithm is performed classically. So, we'll put the rest of it in the "post-processing" step, after we've gotten our measurements from the quantum computer. Each of the measurements above can be converted into integers, which, after we divide by 2m2^m, are our approximations for kr\frac{k}{r}, where kk is random each time.

a = 2
N = 15
 
FACTOR_FOUND = False
num_attempt = 0
 
while not FACTOR_FOUND:
    print(f"\nATTEMPT {num_attempt}:")
    # Here, we get the bitstring by iterating over outcomes
    # of a previous hardware run with multiple shots.
    # Instead, we can also perform a single-shot measurement
    # here in the loop.
    bitstring = list(counts_keep.keys())[num_attempt]
    num_attempt += 1
    # Find the phase from measurement
    decimal = int(bitstring, 2)
    phase = decimal / (2**num_control)  # phase = k / r
    print(f"Phase: theta = {phase}")
 
    # Guess the order from phase
    frac = Fraction(phase).limit_denominator(N)
    r = frac.denominator  # order = r
    print(f"Order of {a} modulo {N} estimated as: r = {r}")
 
    if phase != 0:
        # Guesses for factors are gcd(a^{r / 2} ± 1, 15)
        if r % 2 == 0:
            x = pow(a, r // 2, N) - 1
            d = gcd(x, N)
            if d > 1:
                FACTOR_FOUND = True
                print(f"*** Non-trivial factor found: {x} ***")

Output:


ATTEMPT 0:
Phase: theta = 0.0
Order of 2 modulo 15 estimated as: r = 1

ATTEMPT 1:
Phase: theta = 0.75
Order of 2 modulo 15 estimated as: r = 4
*** Non-trivial factor found: 3 ***

Conclusion

After going through the module, you may be struck with a new appreciation of the brilliance of Peter Shor to have come up with such a clever algorithm. But hopefully you've also reached a new level of understanding of its deceptive simplicity. Although the algorithm may seem impressively (or intimidatingly) complex, if you break it down into each step of logic and go through it slowly, you, too, will be able to run Shor's algorithm.

While we're far from using this algorithm to factor numbers like RSA1024, our quantum computers are getting better every day, and once a threshold called fault tolerance is reached, then algorithms such as these will be soon to follow. It's an exciting time to be learning about quantum computing!


Problems

Critical concepts:

  • Modern cryptographic systems rely on the classical difficulty of factoring large integers.
  • Modular arithmetic — including the structures ZN\mathbb{Z}_N and ZN\mathbb{Z}_N^* — provides the mathematical foundation for Shor’s algorithm.
  • The problem of factoring an integer NN can be reduced to the problem of finding the order of a number modulo NN.
  • Quantum order finding uses quantum phase estimation techniques to determine the period of the function axmodNa^x \mod N.
  • Shor’s algorithm consists of a classical–quantum hybrid workflow that selects a base, performs quantum order finding, and then classically computes the factors from the result.

True/False:

  1. T/F The efficiency of Shor's algorithm threatens the security of RSA encryption.
  2. T/F Shor's algorithm can be executed efficiently on any modern-day quantum computer.
  3. T/F Shor’s algorithm uses quantum phase estimation (QPE) as a key subroutine.
  4. T/F The classical part of Shor’s algorithm involves calculating the greatest common divisor (GCD).
  5. T/F Shor’s algorithm only works for factoring even numbers.
  6. T/F A successful run of Shor’s algorithm always guarantees the correct factors.

Short answer:

  1. Why is Shor’s algorithm considered a potential future threat to RSA encryption?
  2. Why is finding the period, or order, of a modular exponential function helpful for factoring a number in Shor’s algorithm?

Challenge problems:

  1. How many control qubits mm do we need for a given number NN that we're trying to factor to obtain the precision in the QPE needed to find the correct value of the order rr?

  2. Following the procedure we outlined here to factor 15, now try to factor 21.

Was this page helpful?
Report a bug, typo, or request content on GitHub.