import numpy as np

# Importing standard Qiskit libraries
from qiskit import QuantumCircuit, transpile, Aer, IBMQ
from qiskit.tools.jupyter import *
from qiskit.visualization import *
from ibm_quantum_widgets import *
from qiskit.providers.aer import QasmSimulator
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, Aer, transpile, assemble
from qiskit.visualization import plot_histogram
from math import gcd
from numpy.random import randint
import pandas as pd
from fractions import Fraction
print("Imports Successful")

# Loading your IBM Quantum account(s)
provider = IBMQ.load_account()

#defining the U gates for mod 63
def c_amod63(a, power):
    """Controlled multiplication by a mod 63"""
    if a not in [2, 4, 5, 8, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 26, 29, 31, 32, 34, 37, 38, 40, 41, 43, 44, 46, 47, 50, 52, 53, 55, 58, 59, 61, 62]:
        raise ValueError("'a' must be 2, 4, 5, 8, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 26, 29, 31, 32, 34, 37, 38, 40, 41, 43, 44, 46, 47, 50, 52, 53, 55, 58, 59, 61, or 62")
    #the number of qubits used is 6 
    U = QuantumCircuit(6) 
    
    #implementing the swap gates for rotation 
    #we implement every number and see common behaivor between the numbers 
    for iteration in range(power):
        
        if a in [2,61]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)

        if a in [4, 59]:
            U.swap(0,2)
            U.swap(1,3)
            U.swap(2,4)
            U.swap(3,5)
            
        if a in [5, 58]:
            U.swap(0,2)
            U.swap(1,3)
            U.swap(0,4)
            U.swap(1,5)
            
        if a in [8, 55]:
            U.swap(0,3)
            U.swap(1,4)
            U.swap(2,5)
            
        if a in [10, 53]:
            U.swap(0,2)
            U.swap(1,3)
            U.swap(0,4)
            U.swap(1,5)
            U.swap(0,2)
            U.swap(1,3)
            U.swap(0,4)
            U.swap(1,5)
            
        if a in [11, 52]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            
        if a in [13, 50]:
            U.swap(0,5)
            U.swap(1,5)
            U.swap(2,5)
            U.swap(3,5)
            U.swap(4,5)
            
        if a in [16, 47]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            
       # if a in [17, 46]:
       # if a in [19, 44]:
        if a in [20, 43]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            
        if a in [22, 41]:
            U.swap(0,5)
            U.swap(1,5)
            U.swap(2,5)
            U.swap(3,5)
            U.swap(4,5)
        
        if a in [23, 40]: 
            U.swap(0,2)
            U.swap(1,3)
            U.swap(2,4)
            U.swap(3,5)
            
        if a in [25, 38]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
       # if a in [26, 37]:
    
        if a in [29, 34]:
            U.swap(0,1)
            U.swap(1,2)
            U.swap(2,3)
            U.swap(3,4)
            U.swap(4,5)
            
        if a in [31, 32]:
            U.swap(0,5)
            U.swap(1,5)
            U.swap(2,5)
            U.swap(3,5)
            U.swap(4,5)
    
        if a in [62]:
            U.swap(0,3)
            U.swap(1,4)
            U.swap(2,5)
            
        if a in [5, 11, 13, 17, 19, 23, 25, 29, 31, 37, 41, 43, 47, 53, 55, 59, 61]:
            for q in range(6):
                U.x(q)
    U = U.to_gate()
    U.name = "%i^%i mod 63" % (a, power)
    c_U = U.control()
    return c_U

# Specify variables
n_count = 10  # number of counting qubits
a = 13

def qft_dagger(n):
    """n-qubit QFTdagger the first n qubits in circ"""
    qc = QuantumCircuit(n)
    # Don't forget the Swaps!
    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)
    qc.name = "QFT†"
    return qc

# Create QuantumCircuit with n_count counting qubits
# plus 4 qubits for U to act on
qc = QuantumCircuit(n_count + 10, n_count)

# Initialize counting qubits
# in state |+>
for q in range(n_count):
    qc.h(q)
    
# And auxiliary register in state |1>
qc.x(3+n_count)

# Do controlled-U operations
for q in range(n_count):
    qc.append(c_amod63(a, 2**q), 
             [q] + [i+n_count for i in range(6)])

# Do inverse-QFT
qc.append(qft_dagger(n_count), range(n_count))

# Measure circuit
qc.measure(range(n_count), range(n_count))
qc.draw(fold=-1)  # -1 means 'do not fold' 

aer_sim = Aer.get_backend('aer_simulator')
t_qc = transpile(qc, aer_sim)
qobj = assemble(t_qc)
results = aer_sim.run(qobj).result()
counts = results.get_counts()
plot_histogram(counts)

rows, measured_phases = [], []
for output in counts:
    decimal = int(output, 2)  # Convert (base 2) string to decimal
    phase = decimal/(2**n_count)  # Find corresponding eigenvalue
    measured_phases.append(phase)
    # Add these values to the rows in our table:
    rows.append([f"{output}(bin) = {decimal:>3}(dec)", 
                 f"{decimal}/{2**n_count} = {phase:.2f}"])
# Print the rows in a table
headers=["Register Output", "Phase"]
df = pd.DataFrame(rows, columns=headers)
print(df)

Fraction(0.666)

# Get fraction that most closely resembles 0.666
# with denominator < 63
Fraction(0.666).limit_denominator(63)


rows = []
for phase in measured_phases:
    frac = Fraction(phase).limit_denominator(63)
    rows.append([phase, f"{frac.numerator}/{frac.denominator}", frac.denominator])
# Print as a table
headers=["Phase", "Fraction", "Guess for r"]
df = pd.DataFrame(rows, columns=headers)
print(df)

