**Quantum Fourier Transform, Quantum Phase Estimation and Shor’s Algorithm**

**Using the Power of Quantum Computing to Solve Hard Problems — Part 4**

One of the most promising aspect of quantum computing is that it will, one day, have the ability to factor integers in polynomial time. Since the current classical cryptographic protocols, like RSA, rely on factoring large prime numbers being almost impossible for classical computers, Shor’s algorithm has the potential to disrupt our widely used crytographic protocols.

Before we get to Shor’s algorithm, we need to understand Quantum Fourier Transformation (QFT) and Quantum Phase Estimation (QPE).

**Quantum Fourier Transform**

Similar to classical fourier transform, the QFT is the quantum implementation of the discrete Fourier transform over the amplitudes of a wavefunction. The QFT will be used in QPE and also in Shor’s algorithm. You can simply think of QFT transforming a qubit from its computational basis of |0⟩ and |1⟩ to the states in Fourier basis |+⟩ and |-⟩. We will get into the implemention using Qiskit shortly.

**Quantum Phase Estimation**

QPE aims to estimate the phase θ associated with an eigenvalue e^(2πiθ) of a unitary operator U. Since U is unitary, all of its eigenvalues have a norm of 1. QPE is one of the most important and famous algorithms. It is a key subroutine of Shor’s factoring algorithm and for quantum simulation algorithms.

The quantum phase estimation algorithm uses phase kickback to write the phase of U, in the Fourier basis, to the t qubits in the counting register. We then use the inverse QFT to translate this from the Fourier basis into computational basis.

Let’s dive into implementing it using Qiskit. Qiskit has built in QPE but we will define our own using the labs from QGSS 2023. For the purpose of clarity, when defining the initial circuit, I will list out each gate separately. We will then use our knowledge of QPE and QFT to implement Shor’s algorithm.

Importing required libraries:

`from qiskit import QuantumCircuit, Aer, execute`

import numpy as np

from qiskit.visualization import plot_histogram

import matplotlib.pyplot as plt

from math import gcd

Defining QFT and QFT-dagger functions:

`#QFT Circuit`

def qft(n):

"""Creates an n-qubit QFT circuit"""

circuit = QuantumCircuit(n)

def swap_registers(circuit, n):

for qubit in range(n//2):

circuit.swap(qubit, n-qubit-1)

return circuit

def qft_rotations(circuit, n):

"""Performs qft on the first n qubits in circuit (without swaps)"""

if n == 0:

return circuit

n -= 1

circuit.h(n)

for qubit in range(n):

circuit.cp(np.pi/2**(n-qubit), qubit, n)

qft_rotations(circuit, n)

qft_rotations(circuit, n)

swap_registers(circuit, n)

return circuit

#Inverse Quantum Fourier Transform

def qft_dagger(qc, n):

for qubit in range(n//2):

qc.swap(qubit, n-qubit-1)

for j in range(n):

for m in range(j):

qc.cp(-np.pi/float(2**(j-m)), m, j)

qc.h(j)

We will now setup QPE circuit with four counting qubits. We will use phase gate with θ = 1/3 and we will set ƛ = 2π/3. Let’s setup the circuit below, we will define gates manually as opposed to using the QFT-dagger function defined earlier:

phase_register_size = 4

qpe4 = QuantumCircuit(phase_register_size+1, phase_register_size)

import math

qpe4.h(0)

qpe4.h(1)

qpe4.h(2)

qpe4.h(3)

qpe4.x(4)

qpe4.cp(2*math.pi/3, 0, 4)

for _ in range(2**(1)):

qpe4.cp(2*math.pi/3, 1, 4)

for _ in range(2**(2)):

qpe4.cp(2*math.pi/3, 2, 4)

for _ in range(2**(3)):

qpe4.cp(2*math.pi/3, 3, 4)

qpe4.barrier()

qpe4.swap(0,3)

qpe4.h(0)

qpe4.swap(1,2)

qpe4.cp(-math.pi/4, 0,1)

qpe4.h(1)

qpe4.cp(-math.pi/4,0,2)

qpe4.cp(-math.pi/2,1,2)

qpe4.cp(-math.pi/8,0,3)

qpe4.h(2)

qpe4.cp(-math.pi/4,1,3)

qpe4.measure(0,0)

qpe4.cp(-math.pi/2,2,3)

qpe4.measure(1,1)

qpe4.h(3)

qpe4.measure(2,2)

qpe4.measure(3,3)

qpe4.draw("mpl")

Now we will use the AerSimulator to simulate this circuit and plot a histogram:

`sim = Aer.get_backend('aer_simulator')`

shots = 2000

count_qpe4 = execute(qpe4, sim, shots=shots).result().get_counts()

plot_histogram(count_qpe4, figsize=(9,5))

Let us now write a function to convert the bit string into the estimate of θ in decimal form.

`#Grab the highest probability measurement`

max_binary_counts = 0

max_binary_val = ''

sum_item = 0

for key, item in count_qpe4.items():

sum_item += item

if item > max_binary_counts:

max_binary_counts = item

max_binary_val = key

##Function to convert a binary string to decimal

def bin_to_dec(bin_string):

i = 1

dec_val = 0

for ch in bin_string:

dec_val = dec_val + (int(ch) / 2**i)

i = i +1

return dec_val

estimated_phase = bin_to_dec(max_binary_val)

phase_accuracy_window = 2 ** -(len(max_binary_val)) # highest power of 2 (i.e. smallest decimal) this circuit can estimate

print(estimated_phase, phase_accuracy_window)

You will see results of 0.3125 0.0625.

We were expecting θ to be 0.333…., and we estimated 0.3125 with a 6.25% accuracy, which is pretty good. To increase the precision, we can add more counting qubits or use Iterative Phase Estimation (IPE) which is a variant of QPE and requires only one auxiliary qubit. We will not be implementing IPE for the purpose of this article.

**Shor’s Algorithm**

We will now construct a set of functions to implement Shor’s algorithm. The goal of the Shor’s algorithm is to find the prime factors (p, q) of some large number N. We will define the factoring problem into a period finding problem to take advantage of quantum computer to solve this problem in polynomial time. This will comprise the following steps:

**Step 1: **Choose a co-prime a, where the greatest common divisor of a and N is 1.

**Step 2: **Find the order (periodicity) of a modulo N, i.e. the smallest integer r such that a^r mod N = 1.

**Step 3: **Obtain the factor of N by computing the greatest common divisor of a^(r/2) ± 1 and N.

We will implement a circuit to factor a very simple number N = 15 and we will start with a = 7. The function will take an argument for k (the number of qubits) and compose, run, and process Shor’s algorithm to guess the factors.

`#This function will return a ControlledGate object which repeats the action`

# of U, 2^k times

def cU_multi(k):

sys_register_size = 4

circ = QuantumCircuit(sys_register_size)

for _ in range(2**k):

circ.append(U, range(sys_register_size))

U_multi = circ.to_gate()

U_multi.name = '7Mod15_[2^{}]'.format(k)

cU_multi = U_multi.control()

return cU_multi

def shor_qpe(k):

a = 7

#Step 1. Begin a while loop until a nontrivial guess is found

N = 15

count = k

while True:

#Step 2a. Construct a QPE circuit with m phase count qubits

# to guess the phase phi = s/r using the function cU_multi()

qc = QuantumCircuit(count+4, count)

qc.x(count)

for i in range(count):

qc.h(i)

for i in range(count):

qc.append(cU_multi(i), [i,8,9,10,11])

qft_dagger(qc, count)

for n in range(count):

qc.measure(n,n)

#Step 2b. Run the QPE circuit with a single shot, record the results

# and convert the estimated phase bitstring to decimal

sim = Aer.get_backend('aer_simulator')

shots = 1

solution = execute(qc, sim, shots=shots, memory=True) #sim.run(transpile(qc, sim), shots, memory=True)

readings = solution.result().get_memory()

phase = int(readings[0], 2) / 2**count

#Step 3. Use the Fraction object to find the guess for r

frac = Fraction(phase).limit_denominator(N)

s, r = frac.numerator, frac.denominator

#Step 4. Now that r has been found, use the builtin greatest common deonominator

# function to determine the guesses for a factor of N

guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]

#Step 5. For each guess in guesses, check if at least one is a non-trivial factor

# i.e. (guess != 1 or N) and (N % guess == 0)

for guess in guesses:

if (guess != 1 or N) and (N % guess == 0):

break_out_flag = 1

break

if break_out_flag == 1:

break

#Step 6. If a nontrivial factor is found return the list 'guesses', otherwise

# continue the while loop

return guesses

After a few runs, the algorithm should be able to guess the non-trivial factors 3 and 5. This is an easy factoring problem, we are still a long way away from having efficient enough quantum computers to solve the hard factoring problems.

This is Part 4 in a series of articles on quantum computing. The links to all the parts are below: