from google.colab import drive
drive.mount('/content/drive')

!pip install pylatexenc
!pip install qiskit
 # the output will be cleared after installation
from IPython.display import clear_output
clear_output()

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt 
plt.style.use('seaborn-notebook')
from scipy.fft import fft

def periodfind(x,I,N):  #  find period a of function f(x)=x^a mod N
    arr=[]
    for i in range(N):
        f= (x**i)%(I)
        arr.append(f)
    return arr
N=16
I=15
a=11
y=np.array(periodfind(a,I,N))
x=np.arange(N)
print('f(x) = ',y)  # Values of f(x)

# plotting
plt.style.use('seaborn-whitegrid')
def plot_fig(xvals,yvals,a,I): # to make the plot nice
    yvals=list(yvals)
    fig, ax = plt.subplots()
    ax.plot(xvals, yvals, linewidth=1, linestyle='--', marker='o')
    ax.set(xlabel='$x$', ylabel='f(x)=${a}^x$(mod {I})'.format(a=a,I=I),
           title='Period of function f(x)=${a}^x$(mod {I})'.format(a=a,I=I))
    try: # plot r on the graph
        r = yvals[1:].index(1) +1 
        plt.annotate(s='', xy=(0,1), xytext=(r,1), arrowprops=dict(arrowstyle='<->'))
        plt.annotate(s='$r=%i$' % r, xy=(r/3,1.5))
    except:
        print('Could not find period, check a < N and have no common factors.')
    display(plt.show())
plot_fig(x,y,a,I)
    
z=fft(y)  #Apply Fast fourier Transform

plt.style.use('seaborn-whitegrid')
plt.title('Fast Fourier transform of function f(x)=$3^x$(mod 16)')
plt.bar(x,np.abs(z),width=0.5)
display(plt.show())
print('Period is 1/(a/N) is ',2)

%matplotlib inline
import matplotlib.pyplot as plt
from qiskit import IBMQ
# Importing standard Qiskit libraries and configuring account
from qiskit.compiler import transpile, assemble
# Loading your IBM Q account(s)
# provider = IBMQ.load_account() # uncomment to run on quantum computerr

import numpy as np
import math
import qiskit
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, Aer, execute
from qiskit.tools.monitor import job_monitor
from qiskit.visualization import plot_histogram, plot_bloch_multivector, plot_state_hinton
from qiskit import Aer

def measure(circuit,n,c):
    for i in range(n):
        circuit.measure(i,i)
    return circuit

pi = np.pi
qc_q = QuantumRegister(4, 'input') #Initialize the Working Register
qc_c = QuantumRegister(2, 'period') #Initialize the Control Register
c = ClassicalRegister(4) #Register which will eventually hold measurements
Shor = QuantumCircuit(qc_q, qc_c,c)

Shor.x(4) # It is common convention to add a not gate before Modular Exponentiation
Shor.h(qc_q[0:4]) #Initialize 4 Hadamard Gates to allow for all possibilities
Shor.barrier()
Shor.draw('mpl')

#This is a compiled version of modular exponentiation that will not generalize outside of fermat primes.
Shor.cx(qc_q[0],qc_c[0]) 
Shor.cx(qc_q[0],qc_c[1])
Shor.barrier()
Shor.draw(output='mpl')

#This can be ommited by classically swapping the order of measurments if noise is an issue.
Shor.swap(qc_q[0],qc_q[3])
Shor.swap(qc_q[1],qc_q[2])

#The Inverse Quantum Fourier Transform
Shor.h(0)
Shor.crz(-pi/2,qc_q[0],qc_q[1])
Shor.h(1)

Shor.crz(-pi/2,qc_q[1],qc_q[2])
Shor.crz(-pi/4,qc_q[0],qc_q[2])
Shor.h(2)

Shor.crz(-pi/2,qc_q[2],qc_q[3])
Shor.crz(-pi/4,qc_q[1],qc_q[3])
Shor.crz(-pi/8,qc_q[0],qc_q[3])
Shor.h(3)

Shor.barrier()
Shor.draw(output='mpl')

#Measure Working Register
Shor.measure(qc_q[3],3)
Shor.measure(qc_q[2],2)
Shor.measure(qc_q[1],1)
Shor.measure(qc_q[0],0)

Shor.draw(output='mpl')

import pandas as pd
df=pd.read_csv("/content/drive/My Drive/QC-project/SHor code/shor_format_15_11.csv")
print(df.head(16))
plt.plot(df["input_reg_value"],df["period_reg_value"])
plt.show()

# # Load our saved IBMQ accounts and get the least busy backend device with less than or equal to nqubits
# provider = IBMQ.get_provider(hub='ibm-q')
# from qiskit.providers.ibmq import least_busy
# backend = least_busy(provider.backends(filters=lambda x: x.configuration().n_qubits >= N 
#                                        and not x.configuration().simulator 
#                                        and x.status().operational==True))
# print("least busy backend: ", backend)

simulator = Aer.get_backend('qasm_simulator')
# backend = provider.get_backend("ibmq_16_melbourne") # to run in the actual quantum computer

# Replace simulator in execute function with 'backend' to run on a quantum computer
job = execute(Shor, simulator, shots=8192)
result = job.result()
counts = result.get_counts(Shor)
plot_histogram(counts)

output=[]
for keys in list(counts.keys()):
  output.append(int(keys,2))
print('The values with maximum probabilities are: {}'.format(output))

# Classical analysis of result(classical)
print(output)
y_over_N=[out/N for out in output[1:] ]
print('y over N is',y_over_N)

def find_period(output):
    y_over_N=[out/N for out in output[1:] ]
    from fractions import Fraction
    temp=[]
    num=len(y_over_N)
    for i in range(num):
        temp.append(Fraction(y_over_N[i]).limit_denominator(N))
        temp[i]=Fraction(temp[i].denominator,temp[i].numerator)
    return temp[i].numerator

a= find_period(output)
print('The period of the function is ',a)

# find factors by taking GCD
from math import gcd
factors=[]

if (a/2)-(a//2)==0: # test if (a/2) is integer
    p=int(a/2)
    factors.append(gcd(11**p+1,15))
    factors.append(gcd(11**p-1,15))                  
print(factors)

# General Implementation of Shor's Algorithm
# Using Qiskit library
from qiskit.aqua import QuantumInstance
from qiskit.aqua.algorithms import Shor

n=63  # n is the number to be factored.
factors= Shor(n)  # Shor() is a function in qiskit.aqua.algorithms that performs shor's algorithm

simulator = Aer.get_backend('qasm_simulator')
result_dict = factors.run(QuantumInstance(simulator))
result = result_dict['factors']
result


qc_shor=factors.construct_circuit()
qc_shor.draw('mpl')

import qiskit
qiskit.__qiskit_version__

