Implementing Shor’s Factoring Algorithm and Comparing to Classical Factoring Algorithms

  1. The following gate can be used to factor the number of 15 with all possible values of a. This gate is presented withouth proof and taken from the Qiskit Textbook.

def c_amod15(a, power):
    """Controlled multiplication by a mod 15"""
    if a not in [2,4,7,8,11,13]:
        raise ValueError("'a' must be 2,4,7,8,11 or 13")
    U = QuantumCircuit(4)
    for _iteration in range(power):
        if a in [2,13]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a in [7,8]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [4, 11]:
            U.swap(1,3)
            U.swap(0,2)
        if a in [7,11,13]:
            for q in range(4):
                U.x(q)
    U = U.to_gate()
    U.name = f"{a}^{power} mod 15"
    c_U = U.control()
    return c_U
  1. (5 pts.) Explain why a = 2, 4, 7, 8, 11, and 13 are valid choices for factoring N = 15 but other numbers are not.

The requirements for a are that a is larger than 1 but smaller than the number chosen to be factored. In this case the number to factor is 15. a also must not share any common factors with the number to be factored. So a = 2, 4, 7, 8, 11, and 13 are valid options for a when the number to be factored is 15. 3, 5, 6, 9, 12 are not valid has they have factors in common with 15. a = 14 would be a valid option but is not implemented in the curent gate.

  1. (25 pts.) Implement this gate into a quantum phase estimation (QPE) algorithm and set up the appropriate post-processing to determing the factors of N = 15.

  2. (10 pts.) Does there appear to be a value of a which gives the best results (consider the most common output only)?

All values of a do work (check them in the below code) though a = 4 returns only one possible result during the simulation, and that is 3 and 5. One could argue that this is the best option.

  1. (10 pts.) What is the minimim number of counting qubits needed to accurately predict the factors of N = 15 for the a value you choose in part c? For a = 4 the period is 2 so this can be approximated with only a single counting qubit. Results will vary for other values of a.
##############################
##         IMPORTS          ##
##############################
# greatest common divisor
from math import gcd 
# Needed for pi and a few other things
import numpy as np
# For the continued fractions
from fractions import Fraction
# Needed to build the quantum circuit
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
# To simulate running the quantum circuit, transpile any built gates, and 
# visualuze the results of the simulation as a histogram
from qiskit_aer import AerSimulator
from qiskit import transpile
from qiskit.visualization import plot_histogram
# Needed to get the inverse QFT 
from qiskit.circuit.library import QFT
##############################
##  CREATE QUANTUM CIRCUIT  ##
##############################
# Specify variables
counting_qubits = 1  # number of counting qubits
N = 15 # The number of factor
a = 4 # A number such that 1<a<N and gcd(a,N) = 1


# Create QuantumCircuit with specified counting qubits
# plus 4 qubits for U to act on (the U gate needs 4 qubits)
q = QuantumRegister(counting_qubits + 4)
c = ClassicalRegister(counting_qubits)
qc = QuantumCircuit(q,c)

# Initialize counting qubits
# in state |+>
for q in range(counting_qubits):
    qc.h(q)

# And auxiliary register in state |1>
qc.x(counting_qubits)

# Do controlled-U operations
for q in range(counting_qubits):
                                # Control qubit + Target qubits (the last for qubits)
    qc.append(c_amod15(a, 2**q), [q] + [i+counting_qubits for i in range(4)])


# Do inverse-QFT
qc = qc.compose(QFT(counting_qubits, inverse=True), range(counting_qubits))

# Measure circuit
qc.measure(range(counting_qubits), range(counting_qubits))

##############################
## SIMULATE QUANTUM CIRCUIT ##
##############################
# Define the simulator, transpile the circuit, simulate the circuit, and then print
# a histogram of the results
aer_sim = AerSimulator()
t_qc = transpile(qc, aer_sim)
results = aer_sim.run(t_qc).result().get_counts()

##############################
##    CONVERT TO DECIMAL    ##
##############################
def convert_to_decimal(b):
    """
        Converts a binary string to an decimal using the quantum phase estimation
        algorithm
        Inputs:
            b (a binary string)
        Returns:
            sum (a float): the binary string converted to a decimal
    """
    sum = 0
    for i in range(len(b)):
        sum += int(b[i])/2**(i+1)
    return sum

# Note that unlike previous iterations of this function, we do not want to convert
# the result of the QPE algorithm to an angle, so we do not multiply by 2*pi

##########################################
## CONVERT SIMULATION RESULTS TO ANSWER ##
##########################################
# Sort the returned dictionary in order of counts of each result
sorted_results = sorted(results, key=results.get, reverse=True)

# For each possible result, convert it to a decimal
s_over_r_decimal = [convert_to_decimal(i) for i in sorted_results]

# For each possible decimal approximation of s/r, use continued fractions to
# get a more usable representation. Note that the limit_denominator function
# will limit the denomonator of the returned fraction to equal to or below the
# given number (r must be less than 15).
s_over_r = [Fraction(i).limit_denominator(N) for i in s_over_r_decimal]

# Extract the denominator (or r) from all of the fraction approximations to
# s.r
r_guess = [f.denominator for f in s_over_r]

# r must be even, so keep only the even values
r_guess = [r for r in r_guess if r % 2 == 0]

# Define a^(r/2), int because the gcd function must have an int input
a_r_2 = [int(a**(r/2)) for r in r_guess]

# The possible factors of N can be found using a, r, and the gcd function
factor1 = [gcd(num+1,15) for num in a_r_2]
factor2 = [gcd(num-1,15) for num in a_r_2]

# Print the guessed r values as well as the factors found with each r value
# and with each equation
print(r_guess)
print(factor1)
print(factor2)
[2]
[5]
[3]