import math
import random
import numpy as np
import qrisp
import qiskit
from qrisp.interface import VirtualQiskitBackend
from qiskit_alice_bob_provider.local.provider import AliceBobLocalProvider

provider = AliceBobLocalProvider()
backend = provider.get_backend("EMU:40Q:LOGICAL_TARGET")
vrtl_qasm_sim = VirtualQiskitBackend(backend)

def is_prime(n):
    if n <= 1:
        return False
    elif n <= 3:
        return True
    elif n % 2 == 0 or n % 3 == 0:
        return False
    i = 5
    while i * i <= n:
        if n % i == 0 or n % (i + 2) == 0:
            return False
        i += 6
    return True


def find_period_classical(g, N):
    # returns the period of g in Z_N
    e = 1
    while pow(g, e, N) != 1:
        e = e + 1
    return e


def classical_Shor(N):
    # returns the factors of N by using Shor's algorithm
    # with a classical brute-force period finding routine
    if N % 2 == 0:
        return (2, N // 2)
    if is_prime(N):
        return (1, N)

    while True:
        a = random.randint(2, N - 1)
        if math.gcd(a, N) != 1:
            return sorted([math.gcd(a, N), N // math.gcd(a, N)])

        r = find_period_classical(a, N)

        if r % 2 == 0:
            x = pow(a, r // 2, N)
            p = math.gcd(x - 1, N)
            q = math.gcd(x + 1, N)
            if p != 1 and p != N:
                return sorted([p, N // p])
            if q != 1 and q != N:
                return sorted([q, N // q])

classical_Shor(15)

from sympy import continued_fraction_convergents, continued_fraction_iterator, Rational


def get_r_candidates(approx):
    rationals = continued_fraction_convergents(
        continued_fraction_iterator(Rational(approx))
    )
    return [rat.q for rat in rationals]

def get_phase_candidates(a, N):
    qg = qrisp.QuantumModulus(N)
    qg[:] = 1
    qpe_res = qrisp.QuantumFloat(2 * qg.size + 1, exponent=-(2 * qg.size + 1))
    qrisp.h(qpe_res)
    x = a
    for i in range(len(qpe_res)):
        with qrisp.control(qpe_res[i]):
            qg *= x
            x = (x * x) % N
    qrisp.QFT(qpe_res, inv=True)
    qiskit_qc = qpe_res.qs.compile().to_qiskit()
    # Compile the circuit to qiskit for technical limitations
    # of the qrisp framework.
    cr = qiskit.ClassicalRegister(len(qpe_res))
    qiskit_qc.add_register(cr)
    qiskit_qc.measure(range(len(qg), len(qpe_res) + len(qg)), range(len(cr)))
    results = qiskit.execute(qiskit_qc, backend=backend, shots=5).result().get_counts()
    return results

def binary_to_float(binaries):
    # Convert binary string to floating-point decimal
    results = []
    for binary in binaries:
        decimal = 0
        for i, bit in enumerate(binary):
            if bit == "1":
                decimal += 2 ** -(i + 1)
        results.append(decimal)
    return results

def quantum_shor(a, N):
    meas_res = get_phase_candidates(a, N)
    r_candidates = sum(
        [get_r_candidates(approx) for approx in binary_to_float(meas_res.keys())], []
    )
    for cand in sorted(set(r_candidates)):
        if (a**cand) % N == 1:
            r = cand
            break
    else:
        raise Exception(
            "Could not find the right candidate for the period. Please try resampling."
        )
    if r % 2:
        raise Exception(
            "Found an odd period for the given base. Please retry with another value of a."
        )
    p = np.gcd(a ** (r // 2) + 1, N)
    q = int(N / p)
    if p > q:
        return (q, p)
    return (p, q)

a = 2
N = 15
# WARNING: The code below takes approximatively 10 minutes to run
quantum_shor(a, N)

from qiskit_alice_bob_provider.processor.logical_cat import (
    _logical_bit_flip_error,
    _logical_phase_flip_error,
)


def logical_error_rate(distance, average_nb_photons, kappa_1, kappa_2):
    # Returns the logical error rate for the given parameters.
    return _logical_bit_flip_error(
        d=distance, nbar=average_nb_photons
    ) + _logical_phase_flip_error(
        d=distance, nbar=average_nb_photons, k1=kappa_1, k2=kappa_2
    )

def error_rate_upper_bound(V, success_rate):
    if success_rate <= 0 or success_rate > 1:
        raise Exception("The success rate should be between 0 and 1")
    return 1 - np.exp(np.log(success_rate) / V)

a = 2
N = 15
qg = qrisp.QuantumModulus(N)
qg[:] = 1
qpe_res = qrisp.QuantumFloat(2 * qg.size + 1, exponent=-(2 * qg.size + 1))
qrisp.h(qpe_res)
x = a
for i in range(len(qpe_res)):
    with qrisp.control(qpe_res[i]):
        qg *= x
        x = (x * x) % N
qrisp.QFT(qpe_res, inv=True)

# Parameters of the `EMU:40Q:LOGICAL_TARGET` backend
distance = 15
n_qubits = 40
kappa_1 = 100
kappa_2 = 1e7
average_nb_photons = 19

qc_qiskit_TARGET = qpe_res.qs.compile().to_qiskit()

print("Circuit width before transpilation:", qc_qiskit_TARGET.width())
print("Circuit depth before transpilation:", qc_qiskit_TARGET.depth())

qc_qiskit_transpiled_TARGET = qiskit.transpile(qc_qiskit_TARGET, backend=backend)

print("Circuit width after transpilation:", qc_qiskit_transpiled_TARGET.width())
print("Circuit depth after transpilation:", qc_qiskit_transpiled_TARGET.depth())

q_volume_TARGET = qc_qiskit_transpiled_TARGET.depth() * (qc_qiskit_transpiled_TARGET.width())  #width * depth of the circuit
print('Quantum volume of the transpiled circuit:', q_volume_TARGET)

erub_TARGET = error_rate_upper_bound(q_volume_TARGET, 0.5)
print('Error rate upper bound:', f"{erub_TARGET:.3e}")
ler_TARGET = logical_error_rate(distance,average_nb_photons,kappa_1,kappa_2)
print('Backend logical error rate:', f"{ler_TARGET:.3e}")

backend = provider.get_backend("EMU:15Q:LOGICAL_EARLY")

# Parameters of the `EMU:15Q:LOGICAL_EARLY` backend
distance = 13
n_qubits = 15
kappa_1 = 100
kappa_2 = 100000
average_nb_photons = 7

qc_qiskit_EARLY = qpe_res.qs.compile().to_qiskit()

print("Circuit width before transpilation:", qc_qiskit_EARLY.width())
print("Circuit depth before transpilation:", qc_qiskit_EARLY.depth())

qc_qiskit_transpiled_EARLY = qiskit.transpile(qc_qiskit_EARLY, backend=backend)

print("Circuit width after transpilation:", qc_qiskit_transpiled_EARLY.width())
print("Circuit depth after transpilation:", qc_qiskit_transpiled_EARLY.depth())

q_volume_EARLY = qc_qiskit_transpiled_EARLY.depth() * (qc_qiskit_transpiled_EARLY.width())  #width * depth of the circuit
print('Quantum volume of the transpiled circuit:', q_volume_EARLY)

erub_EARLY = error_rate_upper_bound(q_volume_EARLY, 0.5)
print('Error rate upper bound:', f"{erub_EARLY:.3e}")
ler_EARLY = logical_error_rate(distance,average_nb_photons,kappa_1,kappa_2)
print('Backend logical error rate:', f"{ler_EARLY:.3e}")

# WARNING: The code below takes approximatively 30 minutes to run"
for _ in range(10):
     print(quantum_shor(a,N))

import pandas as pd

data = [
    [
        qc_qiskit_transpiled_EARLY.width(),
        qc_qiskit_transpiled_EARLY.depth(),
        q_volume_EARLY,
        f"{ler_EARLY:.3e}",
        f"{erub_EARLY:.3e}",
    ],
    [
        qc_qiskit_transpiled_TARGET.width(),
        qc_qiskit_transpiled_TARGET.depth(),
        q_volume_TARGET,
        f"{ler_TARGET:.3e}",
        f"{erub_TARGET:.3e}",
    ],
]
pd.DataFrame(
    data,
    columns=["Width", "Depth", "Volume", "Logical Error Rate", "Upper bound"],
    index=["EMU:40Q:LOGICAL_TARGET", "EMU:15Q:LOGICAL_EARLY"],
)