import numpy as np
from numpy import pi
from qiskit import QuantumCircuit, Aer, transpile
from qiskit.visualization import plot_histogram, plot_bloch_multivector
from qiskit.circuit import QuantumRegister, ClassicalRegister

from qiskit.circuit.library import QFT

sim = Aer.get_backend('statevector_simulator')
init_state = '101'
qc1 = QuantumCircuit(3)
qc1.initialize(init_state)
qc1.append(QFT(num_qubits=3), [0, 1, 2])
display(qc1.draw())
sv1 = sim.run(transpile(qc1, sim)).result().get_statevector()
plot_bloch_multivector(sv1)


sim = Aer.get_backend('statevector_simulator')
init_state = '101'

qc2 = QuantumCircuit(3)
qc2.initialize(init_state)
qc2.barrier()

qc2.h(2)
qc2.cp(pi/2, 1, 2)
qc2.cp(pi/4, 0, 2)

qc2.h(1)
qc2.cp(pi/2, 0, 1)

qc2.h(0)

qc2.swap(0, 2)

display(qc2.draw())

sv2 = sim.run(qc2).result().get_statevector()
plot_bloch_multivector(sv2)

# implement general QFT

def create_qft(n, draw=False):
    qc = QuantumCircuit(n)
    for i in reversed(range(n)):
        qc.h(i)
        for e, j in enumerate(reversed(range(i))):
            qc.cp(pi/2**(e+1), j, i)
    for i in range(n//2):
        qc.swap(i, n-i-1)
    if draw:
        display(qc.draw())
    return qc.to_gate(label='myqft{}'.format(n))

sim = Aer.get_backend('statevector_simulator')
init_state = '101'

qc3 = QuantumCircuit(3)
qc3.initialize(init_state)
qc3.append(create_qft(3, True), [0, 1, 2])
display(qc3.draw())
sv3 = sim.run(transpile(qc3, sim)).result().get_statevector()
plot_bloch_multivector(sv3)

def create_carrier_orig(name=None, draw=False):
    c0 = QuantumRegister(1, 'c0')
    a0 = QuantumRegister(1, 'a0')
    b0 = QuantumRegister(1, 'b0')
    c1 = QuantumRegister(1, 'c1')
    qc = QuantumCircuit(c0, a0, b0, c1)
    qc.toffoli(1, 2, 3)
    qc.cx(1, 2)
    qc.toffoli(0, 2, 3)
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'carrier_orig')

def create_carrier_orig_inv(name=None):
    qc = QuantumCircuit(4)
    qc.append(create_carrier_orig().inverse(), range(4))
    return qc.to_gate(label=name or 'carrier_orig_inv')

create_carrier_orig('carrier', True)

def create_carrier(name=None, draw=False):
    c0 = QuantumRegister(1, 'c0')
    a0 = QuantumRegister(1, 'a0')
    b0 = QuantumRegister(1, 'b0')
    c1 = QuantumRegister(1, 'c1')
    qc = QuantumCircuit(c0, a0, b0, c1)
    qc.toffoli(1, 2, 3)
    qc.cx(1, 2)
    qc.toffoli(0, 2, 3)
    qc.cx(1, 2)
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'carrier')

def create_carrier_inv(name=None):
    qc = QuantumCircuit(4)
    qc.append(create_carrier().inverse(), range(4))
    return qc.to_gate(label=name or 'carrier_inv')

create_carrier('carrier', True)

def create_sum(name=None, draw=False):
    c0 = QuantumRegister(1, 'c0')
    a0 = QuantumRegister(1, 'a0')
    b0 = QuantumRegister(1, 'b0')
    qc = QuantumCircuit(c0, a0, b0)
    qc.cx(a0, b0)
    qc.cx(c0, b0)
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'sum')

def create_sum_inv(name=None):
    qc = QuantumCircuit(3)
    qc.append(create_sum().inverse(), range(3))
    return qc.to_gate(label=name or 'sum_inv')

create_sum('sum', True)

a = QuantumRegister(3, 'a')
b = QuantumRegister(4, 'b')
c = QuantumRegister(3, 'c')
creg = ClassicalRegister(4, 'creg')
qc = QuantumCircuit(a, b, c, creg)
qc.x(a[0])  # a = 101(2)
qc.x(a[2])
qc.x(b[1])  # b = 110(2)
qc.x(b[2])
qc.append(create_carrier_orig(), [c[0], a[0], b[0], c[1]])
qc.append(create_carrier_orig(), [c[1], a[1], b[1], c[2]])
qc.append(create_carrier_orig(), [c[2], a[2], b[2], b[3]])
qc.cx(a[2], b[2])
qc.append(create_sum(), [c[2], a[2], b[2]])
qc.append(create_carrier_orig_inv(), [c[1], a[1], b[1], c[2]])
qc.append(create_sum(), [c[1], a[1], b[1]])
qc.append(create_carrier_orig_inv(), [c[0], a[0], b[0], c[1]])
qc.append(create_sum(), [c[0], a[0], b[0]])
qc.barrier()
qc.measure(b, creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)

a = QuantumRegister(3, 'a')
b = QuantumRegister(4, 'b')
c = QuantumRegister(3, 'c')
creg = ClassicalRegister(4, 'creg')
qc = QuantumCircuit(a, b, c, creg)
qc.x(a[0])  # a = 101(2)
qc.x(a[2])
qc.x(b[1])  # b = 110(2)
qc.x(b[2])
qc.append(create_carrier(), [c[0], a[0], b[0], c[1]])
qc.append(create_carrier(), [c[1], a[1], b[1], c[2]])
qc.append(create_carrier(), [c[2], a[2], b[2], b[3]])
# qc.cx(a[2], b[2])     # this was put into create_carrier()
qc.append(create_sum(), [c[2], a[2], b[2]])
qc.append(create_carrier_inv(), [c[1], a[1], b[1], c[2]])
qc.append(create_sum(), [c[1], a[1], b[1]])
qc.append(create_carrier_inv(), [c[0], a[0], b[0], c[1]])
qc.append(create_sum(), [c[0], a[0], b[0]])
qc.barrier()
qc.measure(b, creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)

def create_adder(n, name=None, draw=False):
    """create adder with n bits a, n+1 bits b and n bits carrier"""
    a = QuantumRegister(n, 'a')
    b = QuantumRegister(n+1, 'b')
    c = QuantumRegister(n, 'c')
    qc = QuantumCircuit(a, b, c)
    for i in range(n-1):
        qc.append(create_carrier(), [c[i], a[i], b[i], c[i+1]])
    qc.append(create_carrier(), [c[n-1], a[n-1], b[n-1], b[n]])
    
    qc.append(create_sum(), [c[n-1], a[n-1], b[n-1]])
    for i in reversed(range(n-1)):
        qc.append(create_carrier_inv(), [c[i], a[i], b[i], c[i+1]])
        qc.append(create_sum(), [c[i], a[i], b[i]])
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'adder')

def create_adder_inv(n, name=None):
    qc = QuantumCircuit(3*n+1)
    qc.append(create_adder(n).inverse(), range(3*n+1))
    return qc.to_gate(label=name or 'adder_inv')


create_adder(3, 'adder', True)

a = QuantumRegister(3, 'a')
b = QuantumRegister(4, 'b')
c = QuantumRegister(3, 'c')
creg = ClassicalRegister(4, 'creg')
qc = QuantumCircuit(a, b, c, creg)
qc.x(a[0])  # a = 5
qc.x(a[2])
qc.x(b[1])  # b = 6
qc.x(b[2])
qc.append(create_adder(3), range(10))
qc.measure(b, creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)

def create_adder_mod(n, bigN, name=None, draw=False):
    """
    :param n: number of bits for a, c and bigN; b has n+1 bits
    :param bigN: mod number
    :return: Gate of 4n+2 qubits
        0 to n-1 (n qubits) for a
        n to 2n (n+1 qubits) for b
        2n+1 to 3n (n qubits) for c (init to 0)
        3n+1 to 4n (n qubits) for bigN (init to bigN when use)
        4n+1 (1 qubit) for temp (init to 0)
    """
    a = QuantumRegister(n, 'a')
    b = QuantumRegister(n+1, 'b')
    c = QuantumRegister(n, 'c')
    bN = QuantumRegister(n, 'bigN')
    t = QuantumRegister(1, 't')
    qc = QuantumCircuit(a, b, c, bN, t)
    # block 1
    qc.append(create_adder(n), list(a)+list(b)+list(c))
    for i in range(n):
        qc.swap(a[i], bN[i])
    qc.append(create_adder_inv(n), list(a)+list(b)+list(c))
    qc.x(b[n])
    qc.cx(b[n], t[0])
    qc.x(b[n])
    tempN = bigN
    i = 0
    while tempN != 0:
        if tempN % 2 != 0:
            qc.cx(t[0], a[i])
        i = i + 1
        tempN = tempN // 2
    qc.append(create_adder(n), list(a)+list(b)+list(c))
    tempN = bigN
    i = 0
    while tempN != 0:
        if tempN % 2 != 0:
            qc.cx(t[0], a[i])
        i = i + 1
        tempN = tempN // 2
    for i in range(n):
        qc.swap(a[i], bN[i])
    # block 2
    qc.append(create_adder_inv(n), list(a)+list(b)+list(c))
    qc.cx(b[n], t[0])
    qc.append(create_adder(n), list(a)+list(b)+list(c))
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'adder_mod')


create_adder_mod(3, 5, 'adder_mod', True)

n = 3
bigN = 5
#qc = QuantumCircuit(4*n+2)
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bigN')
t = QuantumRegister(1, 't')
creg = ClassicalRegister(4, 'creg')
qc = QuantumCircuit(a, b, c, bN, t, creg)
qc.x(a[2]) # a = 4
qc.x(b[2]) # b = 4
qc.x(bN[0]) # bN = 5
qc.x(bN[2]) 
qc.append(create_adder_mod(n, bigN), list(a) + list(b) + list(c) + list(bN) + list(t))
qc.measure(b, creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)

def create_ctrl_mult_mod(n, bigN, m, name=None, draw=False):
    """
    Calculate zm mod N and output to b

    :param n: number of bits for a, c and bigN; b has n+1 bits
    :param bigN: mod number
    :param m: multiplier
    :return: Gate of 4n+2 qubits
        0 (1 qubit) for x
        1 to n (n qubits) for z
        n+1 to 2n (n qubits) for a
        2n+1 to 3n+1 (n+1 qubits) for b
        3n+2 to 4n+1 (n qubits) for c (init to 0)
        4n+2 to 5n+1 (n qubits) for bigN (init to bigN when use)
        5n+2 (1 qubit) for temp (init to 0)
    """
    x = QuantumRegister(1, 'x')
    z = QuantumRegister(n, 'z')
    a = QuantumRegister(n, 'a')
    b = QuantumRegister(n+1, 'b')
    c = QuantumRegister(n, 'c')
    bN = QuantumRegister(n, 'bN')
    t = QuantumRegister(1, 't')
    qc = QuantumCircuit(x, z, a, b, c, bN, t)
    next_mod = m # m 2^0 mod N
    for j in range(n):
        temp_mod = next_mod
        i = 0
        while temp_mod != 0:
            if temp_mod % 2 != 0:
                qc.ccx(x[0], z[j], a[i])
            i = i + 1
            temp_mod = temp_mod // 2

        qc.append(create_adder_mod(n, bigN), list(a) + list(b) + list(c) + list(bN) + list(t))

        temp_mod = next_mod
        i = 0
        while temp_mod != 0:
            if temp_mod % 2 != 0:
                qc.ccx(x[0], z[j], a[i])
            i = i + 1
            temp_mod = temp_mod // 2
    
        next_mod = (next_mod * 2) % bigN # update for m 2^i+1 mod N 
    
    # seems no need to copy z to b if x=0.  may try to remove this block and test
    # qc.x(x[0])
    # for j in range(n):
    #     qc.ccx(x[0], z[j], b[j])
    # qc.x(x[0])

    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'ctrl_mult_mod')

def create_ctrl_mult_mod_inv(n, bigN, m):
    qc = QuantumCircuit(5*n+3)
    qc.append(create_ctrl_mult_mod(n, bigN, m).inverse(), range(5*n+3))
    return qc.to_gate(label='ctrl_mult_mod_inv')


create_ctrl_mult_mod(3, 5, 3, 'ctrl_mult_mod', True)

n = 3
bigN = 5
m = 3
x = QuantumRegister(1, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
creg = ClassicalRegister(n, 'creg')
qc = QuantumCircuit(x, z, a, b, c, bN, t, creg)
qc.x(x[0]) # x = 1
qc.x(z[0]) # z = 3
qc.x(z[1])
qc.x(bN[0]) # bN = 5
qc.x(bN[2])
qc.append(create_ctrl_mult_mod(n, bigN, m), list(x) + z[:] + list(a) + list(b) + list(c) + list(bN) + list(t))
# for i in range(n):
#     qc.measure(b[i], creg[i])
qc.measure(b[0:n], creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)

from sympy import mod_inverse

def create_mod_exp(n, bigN, y, nx, name=None, draw=False):
    """
    Calculate y^x mod N and output to z

    :param n: number of bits for z, a, c and bigN; b has n+1 bits
    :param bigN: mod number
    :param y: multiplier
    :param nx: number of bits for x
    :return: Gate of 4n+2 qubits
        0 to nx-1 (n qubit) for x
        nx to nx+n-1 (n qubits) for z (init to 1 when use)
        nx+n to nx+2n-1 (n qubits) for a (init to 0)
        nx+2n to nx+3n (n+1 qubits) for b (init to 0)
        nx+3n+1 to nx+4n (n qubits) for c (init to 0)
        nx+4n+1 to nx+5n (n qubits) for bigN (init to bigN when use)
        nx+5n+1 (1 qubit) for temp (init to 0)
    """
    x = QuantumRegister(nx, 'x')
    z = QuantumRegister(n, 'z')
    a = QuantumRegister(n, 'a')
    b = QuantumRegister(n+1, 'b')
    c = QuantumRegister(n, 'c')
    bN = QuantumRegister(n, 'bN')
    t = QuantumRegister(1, 't')
    qc = QuantumCircuit(x, z, a, b, c, bN, t)
    m = y
    for i in range(nx):
        qc.append(create_ctrl_mult_mod(n, bigN, m), 
            [x[i]] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
        for j in range(n):
            qc.cswap(x[i], z[j], b[j])
        m_inv = mod_inverse(m, bigN)
        qc.append(create_ctrl_mult_mod_inv(n, bigN, m_inv), 
            [x[i]] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
        m = (m * m) % bigN # update for next cycle
    if draw:
        display(qc.draw())
    return qc.to_gate(label=name or 'mod_exp')


create_mod_exp(3, 5, 3, 3, 'mod_exp', True)

n = 3
bigN = 5
y = 3
nx = 3
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
creg = ClassicalRegister(n, 'creg')
qc = QuantumCircuit(x, z, a, b, c, bN, t, creg)
qc.x(x[0]) # x = 3
qc.x(x[1])
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 5
qc.x(bN[2])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.measure(z, creg)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)   # 3^3 mod 5 = 2 or 010(2)

n = 5
bigN = 15
y = 7
nx = 8
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
cregx = ClassicalRegister(nx, 'cregx')
qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
qc.h(x)
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 15 = 1111(2)
qc.x(bN[1])
qc.x(bN[2])
qc.x(bN[3])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.append(create_qft(nx), x)
qc.measure(x, cregx)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)


from fractions import Fraction
Fraction(192/2**8).limit_denominator(221)

import math

print(math.gcd(7**2+1, 15))
print(math.gcd(7**2-1, 15))

n = 6
bigN = 55
y = 7
nx = 12
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
cregx = ClassicalRegister(nx, 'cregx')
qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
qc.h(x)
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 55 = 110111(2)
qc.x(bN[1])
qc.x(bN[2])
qc.x(bN[4])
qc.x(bN[5])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.append(create_qft(nx), x)
qc.measure(x, cregx)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)


import matplotlib.pyplot as plt
import math
from fractions import Fraction

def is_solution(k, nx, a, N, p, q):
    r = Fraction(k/2**nx).limit_denominator(N).denominator
    if r % 2 == 1:
        return False
    if math.gcd(a**(r//2)+1, N) == p and math.gcd(a**(r//2)-1, N) == q:
        return True
    if math.gcd(a**(r//2)+1, N) == q and math.gcd(a**(r//2)-1, N) == p:
        return True
    return False

def analyze_result(counts, nx, a, N, p, q):
    xycorrect = {}
    xywrong = {}
    for k,v in counts.items():
        if is_solution(int(k, 2), nx, a, N, p, q):
            xycorrect[int(k, 2)] = v
        else:
            xywrong[int(k, 2)] = v
            
    plt.title('Factor {}'.format(N))
    plt.xlabel('measured values on {} qubits'.format(nx))
    plt.ylabel('counts')
    plt.scatter(xycorrect.keys(), xycorrect.values(), s=1, c='g')
    plt.scatter(xywrong.keys(), xywrong.values(), s=1, c='r')
    plt.legend(['correct', 'wrong'])

    print('correct count {}, wrong count {}, total {}, correct rate {}'.format(
        sum(xycorrect.values()),
        sum(xywrong.values()),
        sum(xycorrect.values()) + sum(xywrong.values()),
        float(sum(xycorrect.values())) / float(sum(xycorrect.values()) + sum(xywrong.values()),)))


analyze_result(counts, 12, 7, 55, 5, 11)

from fractions import Fraction
Fraction(3891/2**12).limit_denominator(55)


import math

print(math.gcd(7**10+1, 55))
print(math.gcd(7**10-1, 55))

n = 8
bigN = 221
y = 7
nx = 16
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
cregx = ClassicalRegister(nx, 'cregx')
qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
qc.h(x)
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 221 = 11011101(2)
qc.x(bN[2])
qc.x(bN[3])
qc.x(bN[4])
qc.x(bN[6])
qc.x(bN[7])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.append(create_qft(nx), x)
qc.measure(x, cregx)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)


analyze_result(counts, 16, 7, 221, 13, 17)

from fractions import Fraction
Fraction(1365/2**16).limit_denominator(221)

import math

print(math.gcd(3**24+1, 221))
print(math.gcd(3**24-1, 221))

17*13

n = 9
bigN = 437
y = 7
nx = 18
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
cregx = ClassicalRegister(nx, 'cregx')
qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
qc.h(x)
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 437 = 110110101(2)
qc.x(bN[2])
qc.x(bN[4])
qc.x(bN[5])
qc.x(bN[7])
qc.x(bN[8])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.append(create_qft(nx), x)
qc.measure(x, cregx)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)


analyze_result(counts, 18, 7, 437, 19, 23)

from fractions import Fraction
Fraction(123128/2**18).limit_denominator(437)

import math

print(math.gcd(7**33+1, 437))
print(math.gcd(7**33-1, 437))

23*19

qc_t = transpile(qc, backend)
print('circuit depth {}, operator count {}'.format(qc_t.depth(), qc_t.count_ops()))

n = 10
bigN = 851
y = 7
nx = 20
x = QuantumRegister(nx, 'x')
z = QuantumRegister(n, 'z')
a = QuantumRegister(n, 'a')
b = QuantumRegister(n+1, 'b')
c = QuantumRegister(n, 'c')
bN = QuantumRegister(n, 'bN')
t = QuantumRegister(1, 't')
cregx = ClassicalRegister(nx, 'cregx')
qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
qc.h(x)
qc.x(z[0]) # z = 1 init
qc.x(bN[0]) # bN = 851 = 1101010011(2)
qc.x(bN[1])
qc.x(bN[4])
qc.x(bN[6])
qc.x(bN[8])
qc.x(bN[9])

qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
qc.append(create_qft(nx), x)
qc.measure(x, cregx)
display(qc.draw())
backend = Aer.get_backend('qasm_simulator')
job = backend.run(transpile(qc, backend), shots=1024)
counts = job.result().get_counts()
print(counts)


analyze_result(counts, 20, 7, 851, 23, 37)

from fractions import Fraction
Fraction(0/2**20).limit_denominator(851)

def gen_qc(n, bigN, y, nx):
    x = QuantumRegister(nx, 'x')
    z = QuantumRegister(n, 'z')
    a = QuantumRegister(n, 'a')
    b = QuantumRegister(n+1, 'b')
    c = QuantumRegister(n, 'c')
    bN = QuantumRegister(n, 'bN')
    t = QuantumRegister(1, 't')
    cregx = ClassicalRegister(nx, 'cregx')
    qc = QuantumCircuit(x, z, a, b, c, bN, t, cregx)
    qc.h(x)
    qc.x(z[0]) # z = 1 init
    tmp_bigN = bigN
    i = 0
    while tmp_bigN != 0:
        if tmp_bigN % 2 == 1:
            qc.x(bN[i])
        tmp_bigN = tmp_bigN // 2
        i = i + 1

    qc.append(create_mod_exp(n, bigN, y, nx), x[:] + z[:] + a[:] + b[:] + c[:] + bN[:] + list(t))
    qc.append(create_qft(nx), x)
    qc.measure(x, cregx)
    backend = Aer.get_backend('qasm_simulator')
    qc_t = transpile(qc, backend)
    return qc_t

qc_15  = gen_qc(4, 15, 7, 8)
qc_55  = gen_qc(6, 55, 7, 12)
qc_221 = gen_qc(8, 221, 7, 16)
qc_437 = gen_qc(9, 437, 7, 18)
qc_851 = gen_qc(10, 851, 7, 20)



print(qc_15.depth(), qc_15.count_ops())
print(qc_851.depth(), qc_851.count_ops())

n_array = [4, 6, 8, 9, 10]
qc_array = [qc_15, qc_55, qc_221, qc_437, qc_851]
depth_array = []
ops_array = []
for qc_i in qc_array:
    depth_array.append(qc_i.depth())  
    ops_array.append(sum(qc_i.count_ops().values()))
print(depth_array)
print(ops_array)

import matplotlib.pyplot as plt
import numpy

opsfit = numpy.polyfit(n_array, ops_array, 3)
print(opsfit)
opsmodel = numpy.poly1d(opsfit)
opsline = numpy.linspace(1, 10, 10)

plt.xlabel('n qubits for N')
plt.ylabel('depth and count_ops')
plt.scatter(n_array, depth_array, marker='o', c='blue')
plt.scatter(n_array, ops_array, marker='x', c='orange')
plt.plot(opsline, opsmodel(opsline), c='orange')
plt.legend(['depth', 'count_ops', 'cout_ops poly deg 3 fit'])
plt.show()

opsline2 = numpy.linspace(1, 100, 100)
plt.xlabel('n qubits for N')
plt.ylabel('count_ops')
plt.scatter(n_array, ops_array, marker='x', c='orange')
plt.plot(opsline2, opsmodel(opsline2), c='orange')
plt.legend(['count_ops', 'cout_ops poly deg 3 fit'])
plt.show()

qc_7663 = gen_qc(13, 7663, 7, 26)


qc_35263 = gen_qc(16, 35263, 7, 32)


opsline3 = numpy.linspace(1, 20, 20)
plt.xlabel('n qubits for N')
plt.ylabel('count_ops')
plt.scatter(n_array, ops_array, marker='x', c='orange')
plt.scatter([13, 16], [sum(qc_7663.count_ops().values()), sum(qc_35263.count_ops().values())])
plt.plot(opsline3, opsmodel(opsline3), c='orange')
plt.legend(['count_ops', 'count_ops for 7663 and 35263', 'cout_ops poly deg 3 fit'])
plt.xticks(range(0, 22, 2))
plt.show()

import qiskit.tools.jupyter
%qiskit_version_table