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 gcdDefining 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 guessesAfter 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:
