# This code has been adapted and modified from IBM Qiskit 2021 and 
# uses the constant optimized modular exponentiation circuit for mod 15 as contained
# in https://arxiv.org/abs/1202.6614.


import numpy as np
import math
from decimal import *
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.visualization import plot_histogram
from qiskit import Aer, transpile, assemble
import pandas as pd
from fractions import Fraction
# 
# import math
# from math import gcd
# from numpy.random import randint
# 
# 
# from decimal import *
print("Imports Successful")

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"


def my_mod(a,n):
    getcontext().prec = 27
    return round((Decimal(a)/Decimal(n) - Decimal(a)//Decimal(n) ) * n)

def constant_optimized_modular_exponentation_modulus15(a, power):
    if a not in [2,7,8,11,13]:
        raise ValueError("'a' must be 2,7,8,11 or 13")
    U = QuantumCircuit(4)        
    for iteration in range(power):
        if a in [2,13]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
        if a in [7,8]:
            U.swap(2,3)
            U.swap(1,2)
            U.swap(0,1)
        if a == 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 = "%i^%i mod 15" % (a, power)
    control_U = U.control()
    return control_U

def inverse_qft(n):
    circuit = QuantumCircuit(n)
    
    for i in range(n//2):
        circuit.swap(i, n-1-i)
        
    for j in range(n):
        for m in range(j):
            circuit.cp(-np.pi/float(2**(j-m)), m, j)
        circuit.h(j)
    circuit.name = "QFT†"
    return circuit

N = 15
a = 7
n_count = 8
counting_register = QuantumRegister(size = n_count, name = "counting_register")
acting_register = QuantumRegister(size = 4, name="acting_register")
classic_register = ClassicalRegister(size = n_count, name="classic_register")
qc = QuantumCircuit(counting_register, acting_register ,classic_register)
initial_state = [1,0]

for q in range(8):
    qc.initialize(initial_state, q)
qc.draw(output = 'mpl', filename = "Step0")

for q in range(n_count):
    qc.h(q)

qc.draw(output = 'mpl', filename = "Step1")

qc.x(3+n_count)

qc.draw(output = 'mpl', filename = "Step1b")

for q in range(n_count):
    qc.append(constant_optimized_modular_exponentation_modulus15(a, 2**q), 
             [q] + [i+n_count for i in range(4)])
    

qc.measure(range(n_count,n_count + 4), range(4))
qc.barrier()
qc.draw(output = 'mpl', filename = "Step2")

qc.append(inverse_qft(n_count), range(n_count))
qc.draw(output = 'mpl', filename = "Step3")

# Measure circuit
qc.measure(range(n_count), range(n_count))
qc.draw(output = 'mpl', filename = "Step4")

qasm_sim = Aer.get_backend('qasm_simulator')
t_qc = transpile(qc, qasm_sim)
qobj = assemble(t_qc)
results = qasm_sim.run(qobj).result()
counts = results.get_counts()

plot_histogram(counts)

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2)  
    phase = decimal/(2**n_count)  
    measured_phases.append(phase)
    
    rows.append([f"{output}(bin) = {decimal:>3}(dec)", 
                 f"{decimal}/{2**n_count} = {phase:.2f}"])
headers=["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
df

rows = []
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(15)
    rows.append([phase, f"{frac.numerator}/{frac.denominator}", frac.denominator])

headers=["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
my_period_r = max(df["Guess for r"])
print("My period (r) is %i" % my_period_r)

# Confirm that the period is 4
xvals = np.arange(N)
xvals = [x.item() for x in xvals]
yvals = [my_mod(a**x, N) for x in xvals]

fig, ax = plt.subplots();
ax.plot(xvals, yvals, linewidth=1, linestyle='dotted', marker='x');
ax.set(xlabel='$x$', ylabel='$%i^x$ mod $%i$' % (a, N),
       title="Example of Periodic Function in Shor's Algorithm");
try: 
    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 ValueError:
    print('Could not find a period')

first_shared_factor = math.gcd((7**(int(my_period_r/2)) + 1), 15)
first_shared_factor
second_shared_factor = math.gcd((7**(int(my_period_r/2)) - 1), 15)
second_shared_factor

%qiskit_copyright