Implementing and Testing the Deutsch-Jozsa Algorithm

# Needed to set up the quantum circuit
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
# Needed to simulate running a quantum computer
from qiskit_aer import AerSimulator
# Neded to visualize the results of running a quantum computer
from qiskit.visualization import plot_histogram
# Needed for vectors and stuff
import numpy as np
# This import let's us visualize the Qiskit qubits as state vectors (bra-ket notation)
from qiskit.quantum_info import Statevector
# This import let's us visualize 
from qiskit.visualization import plot_bloch_multivector
## Can also perform QFT using the built-in function
from qiskit.circuit.library import QFT
# importing Qiskit transpile
from qiskit import transpile
  1. (15 pts.) Using the quantum phase estimation (QPE) algorithm we created in class, with the controlled U gate being the controlled phase gate with \(\phi=\frac{\pi}{3}\), determine how many counting qubits are needed to match the eigenvalue to one decimal place, two decimal places, and three decimal places. From these results, comment on the expected number of qubits you would need to achieve the computational precision. Here we will define this as the precision held by the standard float (about 10 decimal places).
# Write a little function to convert the output of the quantum circuit
# to an angle

##############################
##     CONVERT TO ANGLE     ##
##############################
def convert_to_angle(b):
    """
        Converts a binary string to an angle using the quantum phase estimation
        algorithm
        Inputs:
            b (a binary string)
        Returns:
            angle (a float): the corresponding angle
    """
    sum = 0
    for i in range(len(b)):
        sum += int(b[i])/2**(i+1)
    angle = sum*2*np.pi
    return angle

## WARNING: This starts to really chug above 10 qubits 
## (remember this is a simulated quantum computer,
# not a real one)
def qpe_cp (phi, counting_qubits):
    
    # Need one more qubit than bit
    q = QuantumRegister(counting_qubits + 1)
    c = ClassicalRegister(counting_qubits)
    qc = QuantumCircuit(q,c)
    
    ## First we need to make sure the last qubit is in an eigenstate of the U gate. 
    ## The up and down states are both possible eigenvectors but up is trivial 
    ## (has an eigenvalue of 1) so let's go with down
    qc.x(counting_qubits)
    qc.barrier()
    
    # Apply H-Gates to counting qubits:
    for qubit in range(counting_qubits):
        qc.h(qubit)
    qc.barrier()
        
    # Do the controlled-U operations:
    repetitions = 1
    for control in range(counting_qubits):
        for i in range(repetitions):
            qc.cp(phi, control, counting_qubits)
        qc.barrier()
        repetitions *= 2
    
    # Do the inverse QFT:
    qc = qc.compose(QFT(counting_qubits, inverse=True), range(counting_qubits))
    qc.barrier()
    
    # Measure of course!
    for n in range(counting_qubits):
       qc.measure(n,n)
    
    return qc

def simulate_and_sort (qc, shots=2048):
    # Run the simulation
    simulator = AerSimulator()
    t_qpe = transpile(qc, simulator)
    results = simulator.run(t_qpe, shots=shots).result()
    answer = results.get_counts()
    # Sort the returned dictionary in order of counts of each result
    sorted_answer = sorted(answer, key=answer.get, reverse=True)
    return sorted_answer

def compare_decimal_places (true, approx):
    # Convert floats to strings
    true = str(true)
    approx = str(approx)
    # Determine the minimum length to avoid index out of range errors
    min_length = min(len(true), len(approx))

    max_match = None
    
    # Iterate over the characters of both strings
    for i in range(min_length):
        if true[i] != approx[i]:
            max_match = i
            break
  
    true_dec = true.find(".")
    approx_dec = approx.find(".")

    if max_match == None:
        #print("The first digits do not match.")
        return 0
    elif true_dec != approx_dec:
        #print("The numbers before the decimal are of different lengths")
        return 0
    elif max_match < true_dec:
        #print("The numbers are not matched before the decimal")
        return 0
    else:
        #print("The numbers match to", max_match-1-true_dec, "places")
        return max_match-1-true_dec
        
phi = np.pi/3
one_match = False
two_match = False
three_match = False
for num_q in range(1,13):
    qc = qpe_cp(phi, num_q)
    sorted_answer = simulate_and_sort(qc)
    best_answer = convert_to_angle(sorted_answer[0])
    match = compare_decimal_places (phi, best_answer)
    if not one_match and match==1:
        print("The minimum number of counting qubits for a one decimal match is:", num_q)
        one_match = True
    if not two_match and match==2:
        print("The minimum number of counting qubits for a two decimal match is:", num_q)
        two_match = True
    if not three_match and match==3:
        print("The minimum number of counting qubits for a three decimal match is:", num_q)
        three_match = True        
    
    
The minimum number of counting qubits for a one decimal match is: 6
The minimum number of counting qubits for a two decimal match is: 9
The minimum number of counting qubits for a three decimal match is: 12

2. (15 pts.) Rewrite the QPE algorithm from class to use any other gate for U (that is not the phase gate). Set up your algorithm and prove that you can get (at least approximately) both eigenvalues of your new U gate.

## WARNING: This starts to really chug above 10 qubits 
## (remember this is a simulated quantum computer,
# not a real one)
def qpe_cz (counting_qubits, eigenvector):
    
    # Need one more qubit than bit
    q = QuantumRegister(counting_qubits + 1)
    c = ClassicalRegister(counting_qubits)
    qc = QuantumCircuit(q,c)
    
    ## First we need to make sure the last qubit is in an eigenstate of the U gate. 
    ## The up and down states are both possible eigenvectors but up is trivial 
    ## (has an eigenvalue of 1) so let's go with down
    if eigenvector=="down":
        qc.x(counting_qubits)
        qc.barrier()
    
    # Apply H-Gates to counting qubits:
    for qubit in range(counting_qubits):
        qc.h(qubit)
    qc.barrier()
        
    # Do the controlled-U operations:
    repetitions = 1
    for control in range(counting_qubits):
        for i in range(repetitions):
            qc.cz(control, counting_qubits)
        qc.barrier()
        repetitions *= 2
    
    # Do the inverse QFT:
    qc = qc.compose(QFT(counting_qubits, inverse=True), range(counting_qubits))
    qc.barrier()
    
    # Measure of course!
    for n in range(counting_qubits):
       qc.measure(n,n)
    
    return qc

def simulate_and_sort (qc, shots=2048):
    # Run the simulation
    simulator = AerSimulator()
    t_qpe = transpile(qc, simulator)
    results = simulator.run(t_qpe, shots=shots).result()
    answer = results.get_counts()
    # Sort the returned dictionary in order of counts of each result
    sorted_answer = sorted(answer, key=answer.get, reverse=True)
    return sorted_answer

def convert_to_eigenvalue (angle):
    return np.exp(1.0j*angle)
num_q = 10
qc = qpe_cz(num_q, "down")
sorted_answer = simulate_and_sort(qc)
best_answer = convert_to_eigenvalue(convert_to_angle(sorted_answer[0]))
print("The eigenvalue should be -1")
print("The eigenvalue is", np.round(best_answer, 3))
The eigenvalue should be -1
The eigenvalue is (-1+0j)
num_q = 10
qc = qpe_cz(num_q, "up")
sorted_answer = simulate_and_sort(qc)
best_answer = convert_to_eigenvalue(convert_to_angle(sorted_answer[0]))
print("The eigenvalue should be 1")
print("The eigenvalue is", np.round(best_answer, 3))
The eigenvalue should be 1
The eigenvalue is (1+0j)