import unittest

import hypothesis.strategies as st
from hypothesis import given, settings, note

import matplotlib.pyplot as plt
import numpy as np
import math
from qiskit import QuantumCircuit, Aer, transpile, assemble, execute
from qiskit.visualization import plot_histogram
from qiskit.circuit.library import CCXGate, CXGate, CSwapGate, HGate, SwapGate, CPhaseGate
from numpy.random import randint
import pandas as pd
from fractions import Fraction
from math import gcd # greatest common divisor

"""Function to compute the elements of Z_n."""
def multiplicative_group(n):
    """Returns the multiplicative group modulo n.

    Args:
        n: Modulus of the multiplicative group.
    """
    print("multiplicative group")
    assert n > 1
    group = [1]
    for x in range(2, n):
        if math.gcd(x, n) == 1:
            group.append(x)
    return group

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

qft_dagger(8)

def general_modular_exp(a, N, power):
    mult_group = multiplicative_group(N)
    if a not in mult_group:
        raise ValueError("'a' must be in " + str(mult_group))
    print(mult_group)

general_modular_exp(7, 15, 2)

def VBE_Adder(qubit_size):
    #outputs a + b
    #inverse outputs b - a (but in two's complement)
    #move to carry
    qc = QuantumCircuit(3*qubit_size + 1)
    for i in range(qubit_size-1):
        qc.ccx(i, i+qubit_size, 2+i+(2*qubit_size))
    
    #toffoli to second register
    qc.ccx(qubit_size-1, (2*qubit_size)-1, 2*qubit_size)
    
    #add a to b
    for i in range(qubit_size):
        qc.cx(i, i+qubit_size)
    
    #second register carry and last b
    for i in range(qubit_size-1):
        qc.ccx(i+qubit_size, 1+i+(2*qubit_size), 2+i+(2*qubit_size))
        
    qc.ccx((2*qubit_size)-1, 3*qubit_size, 2*qubit_size)
    
    qc.cx(3*qubit_size, (2*qubit_size)-1)
    
    #adder overflow
    for i in range(qubit_size-1):
        qc.ccx((2*qubit_size) - 2 - i, (3*qubit_size) - i - 1, (3*qubit_size) - i) 
        qc.cx((qubit_size) - 2 - i, (2*qubit_size) - 2 - i)
        qc.ccx((qubit_size) - 2 - i, (2*qubit_size) - 2 - i, (3*qubit_size) - i) 
        qc.cx((qubit_size) - 2 - i, (2*qubit_size) - 2 - i) 
        qc.cx((3*qubit_size) - i - 1, (2*qubit_size) - 2 - i) 
    
    #print(qc.draw(output='text'))
    qc.name = "VBE ADDER"
    #qc.to_gate()
    return qc    

#print(VBE_Adder(4))

def Modular_adder(qubit_size, N):
    binN = bin(N)[2:].zfill(qubit_size)[::-1]
        
    #move to carry
    qc = QuantumCircuit(4*qubit_size + 2)
    
    qc.append(VBE_Adder(qubit_size), [i for i in range(3*qubit_size+1)])
        
    for i in range(qubit_size):
        qc.swap(i, i + 1 + 3*qubit_size)
    
    qc.append(VBE_Adder(qubit_size).inverse(), [i for i in range(3*qubit_size+1)])
    
    qc.x(2*qubit_size)
    qc.cx(2*qubit_size, 4*qubit_size+1)
    qc.x(2*qubit_size)
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.cx(4*qubit_size + 1, i) 
    
    qc.append(VBE_Adder(qubit_size), [i for i in range(3*qubit_size+1)])
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.cx(4*qubit_size + 1, i) 
    
    for i in range(qubit_size):
        qc.swap(i, i + 1 + 3*qubit_size)
    
    qc.append(VBE_Adder(qubit_size).inverse(), [i for i in range(3*qubit_size+1)])
    
    qc.cx(2*qubit_size, 4*qubit_size+1)
    
    qc.append(VBE_Adder(qubit_size), [i for i in range(3*qubit_size+1)])
    
    qc.name = "MODULAR ADDER"
    
    return qc  

Modular_adder(4,11)

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:60% !important; }</style>"))

size = 4
N = 7

qc = QuantumCircuit(4*size + 2)

# a = 3,2,1,0

#qc.x(3)
qc.x(2)
qc.x(1)
qc.x(0)


# b = 8,7,6,5,4

#qc.x(7)
qc.x(6)
qc.x(5)
qc.x(4)

# carry = 12,11,10,9

# modulus N = 16,15,14,13
binN = bin(N)[2:].zfill(size)[::-1]
for i in range(size):
    if binN[i] == "1":
        qc.x(3*size + 1 + i) 

# temporary carry = 17 

print(qc)

qc.measure_all()

qc.append(Modular_adder(size,N), [i for i in range(4*size+2)])

qc.measure_all()

backend = Aer.get_backend('aer_simulator') 
job = execute(qc, backend, shots=1, memory=True)
readings = job.result().get_memory()
second_set = readings[0][0:4*size+2]
first_set = readings[0][(4*size)+3:(8*size)+6]
print(readings)
print(first_set)
print(second_set)

print('a     =  ' + str(first_set[3*size + 2: 4*size + 2]))
print('b     = ' + str(first_set[2*size + 1: 3*size + 2]))
print('b out = ' + str(second_set[2*size + 1: 3*size + 2]))
#print('carry =  ' + str(second_set[0:size]))


def Modular_multiplier(qubit_size, N, a):        
    qc = QuantumCircuit(5*qubit_size + 3)
    
    for i in range(qubit_size):
        mod = bin(a*(2**i) % N)[2:].zfill(qubit_size)[::-1]
        
        for j in range(qubit_size):
            if mod[j] == "1":
                qc.ccx(0, i+1, qubit_size + 1 + j) 

        #print([qubit_size + i + 1 for i in range(4*qubit_size+2)])
        qc.append(Modular_adder(qubit_size, N), [qubit_size + i + 1 for i in range(4*qubit_size+2)])
             
        #can just repeat the above method, but this looks symmetric 
        #the circuit diagram is more intuitive
        for j in range(qubit_size):
            if mod[::-1][j] == "1":                
                qc.ccx(0, i+1, 2*qubit_size - j) 
                
    #set x if control is false
    qc.x(0)
    for i in range(qubit_size):
        qc.ccx(0,i+1,2*qubit_size + 1 + i) 
    qc.x(0)
        
    qc.name = "MODULAR MULTIPLIER " + str(a) + "x mod " + str(N)  
        
    #print(qc.draw(output='text'))
    return qc  

Modular_multiplier(4, 15, 13)

size = 2
N = 2
a = 1

qc = QuantumCircuit(5*size + 3)

#control qubit
qc.x(0)

# x = 4,3,2,1
#qc.x(4)
#qc.x(3)
#qc.x(2)
#qc.x(1)

# N = 21,20,19,18
binN = bin(N)[2:].zfill(size)[::-1]
for i in range(size):
    if binN[i] == "1":
        qc.x(4*size + 2 + i) 

# temporary carry = 17 

#print(qc)

qc.measure_all()

qc.append(Modular_multiplier(size, N, a), [i for i in range(5*size+3)])

#print(qc)

qc.measure_all()

backend = Aer.get_backend('aer_simulator') 
job = execute(qc, backend, shots=1, memory=True)
readings = job.result().get_memory()
second_set = readings[0][0:5*size+3]
first_set = readings[0][(5*size)+4:(10*size)+8]
print(readings)
print(first_set)
print(second_set)

print('x     =  ' + str(first_set[4*size + 2: 5*size + 2]))
print('b     = ' + str(first_set[2*size + 1: 3*size + 2]))
print('b out = ' + str(second_set[2*size + 1: 3*size + 2]))
#print('carry =  ' + str(second_set[0:size]))


def Modular_exponentiation(qubit_size, N, a):            
    qc = QuantumCircuit(6*qubit_size + 3)
    
    #need to set x to 1
    qc.x(qubit_size+1)
        
    for i in range(qubit_size):
    #for i in range(1):
        qc.cx(i,qubit_size)
        print([i + qubit_size for i in range(5*qubit_size+3)])
        qc.append(Modular_multiplier(qubit_size, N, a**(2**i)%N), [i + qubit_size for i in range(5*qubit_size+3)])
        for j in range(qubit_size):
            qc.swap(qubit_size+1+j, 3*qubit_size+1+j)
        #we dont use A for this multiplier because we need the modular multiplicative inverse
        #print(a**(2**i)%N)
        #print(pow(a**(2**i)%N, -1, N))
        qc.append(Modular_multiplier(qubit_size, N, pow(a**(2**i)%N, -1, N)).inverse(), [i + qubit_size for i in range(5*qubit_size+3)])
        qc.cx(i,qubit_size)
                            
    qc.name = "MODULAR EXPONENTIATION"    
        
    #print(qc.draw(output='text'))
    return qc  

Modular_exponentiation(3, 7, 3)

qubit_size = 2
N = 3
a = 1

qc = QuantumCircuit(6*qubit_size + 3)

# x = 3,2,1,0
#qc.x(3)
qc.x(2)
#qc.x(1)
#qc.x(0)

# N = 21,20,19,18
binN = bin(N)[2:].zfill(qubit_size)[::-1]

print(binN)

for i in range(qubit_size):
    if binN[i] == "1":
        qc.x(5*qubit_size + 2 + i) 
        
#print(qc)
        
qc.measure_all()

qc.append(Modular_exponentiation(qubit_size, N, a), [i for i in range(6*qubit_size+3)])

qc.measure_all()

print(qc)

backend = Aer.get_backend('aer_simulator') 
job = execute(qc, backend, shots=1, memory=True)
readings = job.result().get_memory()[0].split(' ')
second_set = readings[0]
first_set = readings[1]
print(first_set)
print(second_set)

first_set_exp = first_set[(5 * qubit_size) + 3 : (6 * qubit_size) + 3]
first_set_x = first_set[(4 * qubit_size) + 2 : (5 * qubit_size) + 2]
first_set_a = first_set[(3 * qubit_size) + 2 : (4 * qubit_size) + 2]
first_set_b = first_set[(2 * qubit_size) + 1 : (3 * qubit_size) + 2]
first_set_N = first_set[qubit_size + 1 : (2 * qubit_size) + 1]

second_set_exp = second_set[(5 * qubit_size) + 3 : (6 * qubit_size) + 3]
second_set_x = second_set[(4 * qubit_size) + 2 : (5 * qubit_size) + 2]
second_set_a = second_set[(3 * qubit_size) + 2 : (4 * qubit_size) + 2]
second_set_b = second_set[(2 * qubit_size) + 1 : (3 * qubit_size) + 2]
second_set_N = second_set[qubit_size + 1 : (2 * qubit_size) + 1]

print("A = " + str(a))
print("exps 1 -> 2 = " + str(first_set_exp) + " -> " + str(second_set_exp))
print("xs 1 -> 2 = " + str(first_set_x) + " -> " + str(second_set_x))
print("as 1 -> 2 = " + str(first_set_a) + " -> " + str(second_set_a))
print("bs 1 -> 2 = " + str(first_set_b) + " -> " + str(second_set_b))
print("Ns 1 -> 2 = " + str(first_set_N) + " -> " + str(second_set_N))

#print('carry =  ' + str(second_set[0:size]))


def run_shors(qubit_size, N, a, run_amt):
    qc = QuantumCircuit(6*qubit_size + 3, qubit_size)

    # N = 21,20,19,18
    binN = bin(N)[2:].zfill(qubit_size)[::-1]

    for i in range(qubit_size):
        qc.h(i) 
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.x(5*qubit_size + 2 + i) 

    qc.append(Modular_exponentiation(qubit_size, N, a), [i for i in range(6*qubit_size+3)])

    qc.append(qft_dagger(qubit_size), [i for i in range(qubit_size)])
    
    qc.measure([i for i in range(qubit_size)],[i for i in range(qubit_size)])

    print(qc)
    
    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=run_amt, memory=True)
    return(job.result().get_memory())

print(run_shors(8,14,9,20))

%timeit run_shors(10,26,7,20)

# for loop at the bootom from qiskit textbook for postpro
# https://qiskit.org/textbook/ch-algorithms/shor.html#:~:text=if%20phase%20!%3D%200,guess)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20factor_found%20%3D%20True

def shors_with_postpro(N, a, tries):
    factors = []
    
    if not a < N:
        raise ValueError("a(%i) must be less than N(%i)"%(a,N))
    
    common_denom = gcd(a,N)
    
    if not common_denom == 1:
        raise ValueError("a(%i) and N(%i) must be coprime but they have a gcd of (%i)"%(a, N, common_denom))
        
    length_message = len(bin(N)[2:])*2 
        
    if length_message > 8:
        print(f"The length of the number you are trying to factor is large ({length_message}), it may fail or take a long time")
        
    res = run_shors(len(bin(N)[2:])*2, N, a, tries)
    
    for result in res:
        phase = (int(result,2)/(2**length_message))
        frac = Fraction(phase).limit_denominator(N)
        r = frac.denominator
        if phase != 0:
            guesses = [gcd(a**(r//2)-1, N), gcd(a**(r//2)+1, N)]
            for guess in guesses:
                if guess not in [1,N] and (N % guess) == 0: 
                    factors += [guess]
            return factors
        
shors_with_postpro(30, 7, 20)

def twos_comp_to_int(binNum):
    binNum = int(binNum,2) ^ ((2**len(binNum))-1)
    binNum = bin(binNum + 1)
    return -int(binNum,2)

@st.composite
def draw_pair_of_ints(draw): 
    size = draw(st.integers(min_value=1, max_value=6)) 
    a = draw(st.integers(min_value=0, max_value=(2**size)-1))
    b = draw(st.integers(min_value=0, max_value=(2**size)-1))
    return(a, b, size)

@st.composite
def draw_shors_variables(draw): 
    size = draw(st.integers(min_value=2, max_value=5)) 
    mod = draw(st.integers(min_value=3, max_value=(2**size)-1))
    a = draw(st.integers(min_value=2, max_value=mod-1))
    return(a, mod, size)

@st.composite
def draw_pair_of_ints_and_mod(draw): 
    size = draw(st.integers(min_value=2, max_value=3)) 
    mod = draw(st.integers(min_value=2, max_value=(2**size)-1))
    a = draw(st.integers(min_value=0, max_value=mod-1))
    b = draw(st.integers(min_value=0, max_value=mod-1))
    return(a, b, size, mod)

@given(draw_shors_variables())
@settings(deadline=None)
def test_shors_with_postpro(vals):
    a = vals[0]
    N = vals[1]
    size = vals[2]
    
    for i in range(size):
        try:
            pow(a**(2**i)%N, -1, N)
        except ValueError:
            assert(not(math.gcd(a,N) == 1))
    
    if math.gcd(a,N) == 1:
        factors = shors_with_postpro(N, a, 20)

        for factor in factors:
            print(f"N {N}, factor {factor}")
            assert(N % factor == 0)
    
@given(draw_pair_of_ints_and_mod())
@settings(deadline=None)
def test_modular_exponentiation(vals):
    a = vals[0]
    x = vals[1]
    size = vals[2]
    N = vals[3]
    
    for i in range(size):
        try:
            pow(a**(2**i)%N, -1, N)
        except ValueError:
            assert(not(math.gcd(a,N) == 1))
            
    if math.gcd(a,N) == 1:
        qc = QuantumCircuit(6*size + 3)

        binX = bin(x)[2:].zfill(size)[::-1]
        binN = bin(N)[2:].zfill(size)[::-1]

        print("size " + str(size))
        print("a " + str(a))
        print("x " + str(x))
        print("N " + str(N))
        print("\n")

        # x = 3,2,1,0
        for i in range(size):
            if binX[i] == "1":
                qc.x(i)

        for i in range(size):
            if binN[i] == "1":
                qc.x(5*size + 2 + i) 


        #print(qc)

        qc.measure_all()

        qc.append(Modular_exponentiation(size, N, a), [i for i in range(6*size+3)])

        qc.measure_all()

        backend = Aer.get_backend('aer_simulator') 
        job = execute(qc, backend, shots=1, memory=True)
        readings = job.result().get_memory()[0].split(' ')
        second_set = readings[0]
        first_set = readings[1]
        print(first_set)
        print(second_set)

        first_set_exp = first_set[(5 * size) + 3 : (6 * size) + 3]
        first_set_x = first_set[(4 * size) + 2 : (5 * size) + 2]
        first_set_a = first_set[(3 * size) + 2 : (4 * size) + 2]
        first_set_b = first_set[(2 * size) + 1 : (3 * size) + 2]
        first_set_N = first_set[size + 1 : (2 * size) + 1]

        second_set_exp = second_set[(5 * size) + 3 : (6 * size) + 3]
        second_set_x = second_set[(4 * size) + 2 : (5 * size) + 2]
        second_set_a = second_set[(3 * size) + 2 : (4 * size) + 2]
        second_set_b = second_set[(2 * size) + 1 : (3 * size) + 2]
        second_set_N = second_set[size + 1 : (2 * size) + 1]

        print("A = " + str(a))
        print("exps 1 -> 2 = " + str(first_set_exp) + " -> " + str(second_set_exp))
        print("xs 1 -> 2 = " + str(first_set_x) + " -> " + str(second_set_x))
        print("as 1 -> 2 = " + str(first_set_a) + " -> " + str(second_set_a))
        print("bs 1 -> 2 = " + str(first_set_b) + " -> " + str(second_set_b))
        print("Ns 1 -> 2 = " + str(first_set_N) + " -> " + str(second_set_N))
        print(a)
        print(x)
        print(N)
        print('\n\n')
        assert((a ** x) % N == int(second_set_x, 2))

@given(draw_pair_of_ints_and_mod(), st.booleans())
@settings(deadline=None)
def test_modular_multiplier_then_inverse_returns_original_values(vals, control):
    a = vals[0]
    x = vals[1]
    size = vals[2]
    N = vals[3]
    
    qc = QuantumCircuit(5*size + 3)

    binX = bin(x)[2:].zfill(size)[::-1]
    binN = bin(N)[2:].zfill(size)[::-1]
             
    print("size " + str(size))
    print("a " + str(a))
    print("x " + str(x))
    print("N " + str(N))
    print("control " + str(control))
    print("\n")
    
    #control qubit
    if control == True:
        qc.x(0)
    
    # x = 3,2,1,0
    for i in range(size):
        if binX[i] == "1":
            qc.x(i+1)
        
    for i in range(size):
        if binN[i] == "1":
            qc.x(4*size + 2 + i) 


    print(qc)

    qc.measure_all()

    qc.append(Modular_multiplier(size, N, a), [i for i in range(5*size+3)])
    qc.append(Modular_multiplier(size, N, a).inverse(), [i for i in range(5*size+3)])

    qc.measure_all()

    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=1, memory=True)
    readings = job.result().get_memory()
    second_set = readings[0][0:5*size+3]
    first_set = readings[0][(5*size)+4:(10*size)+7]
    print("readings " + str(readings))
    print("first set " + str(first_set))
    print("second set " + str(second_set))

    assert(first_set == second_set)


@given(draw_pair_of_ints_and_mod(), st.booleans())
@settings(deadline=None)
def test_modular_multiplier(vals, control):
    a = vals[0]
    x = vals[1]
    size = vals[2]
    N = vals[3]
    
    qc = QuantumCircuit(5*size + 3)

    binX = bin(x)[2:].zfill(size)[::-1]
    binN = bin(N)[2:].zfill(size)[::-1]
    
    #control qubit
    if control == True:
        qc.x(0)
    
    # x = 3,2,1,0
    for i in range(size):
        if binX[i] == "1":
            qc.x(i+1)
        
    for i in range(size):
        if binN[i] == "1":
            qc.x(4*size + 2 + i) 
            
    qc.measure_all()

    qc.append(Modular_multiplier(size, N, a), [i for i in range(5*size+3)])

    qc.measure_all()

    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=1, memory=True)
    readings = job.result().get_memory()
    second_set = readings[0][0:5*size+3]
    first_set = readings[0][(5*size)+3:(10*size)+7]
    print("readings " + str(readings))
    print("first set " + str(first_set))
    print("x " + str(first_set[4*size + 3: 5*size + 3]))
    x = int(first_set[4*size + 3: 5*size + 3], 2)
    print("b " + str(first_set[2*size + 1: 3*size + 2]))
    b = int(first_set[2*size + 1: 3*size + 2], 2)
    print("bout " + str(second_set[2*size + 1: 3*size + 2]))
    bout = int(second_set[2*size + 1: 3*size + 2], 2)
    
    print("control " + str(control))
    print("N " + str(N))
    print('b     = ' + str(b))
    print('a     =  ' + str(a))    
    print('x     =  ' + str(x))
    print('b out = ' + str(bout))
    print('\n\n')
    if (control == False):
        assert(bout == x)
    else: 
        assert((a * x) % N == bout)
    #assert(carry == 0)

@given(draw_pair_of_ints_and_mod())
@settings(deadline=None)
def test_modular_adder(vals):
    a = vals[0]
    b = vals[1]
    size = vals[2]
    N = vals[3]
    
    qc = QuantumCircuit(4*size + 2)

    binA = bin(a)[2:].zfill(size)[::-1]
    binB = bin(b)[2:].zfill(size)[::-1]
    binN = bin(N)[2:].zfill(size)[::-1]
             
    print("size " + str(size))
    print("a " + str(a))
    print("b " + str(b))
    print("N " + str(N))
    print("\n")
    
    
    # a = 3,2,1,0
    for i in range(size):
        if binA[i] == "1":
            qc.x(i)
        

    # b = 8,7,6,5,4
    for i in range(size):
        if binB[i] == "1":
            qc.x(i+size)
             
    # carry = 12,11,10,9

    # modulus N = 16,15,14,13
    for i in range(size):
        if binN[i] == "1":
            qc.x(3*size + 1 + i) 

    # temporary carry = 17 

    #print(qc)

    qc.measure_all()

    qc.append(Modular_adder(size,N), [i for i in range(4*size+2)])

    qc.measure_all()

    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=1, memory=True)
    readings = job.result().get_memory()
    second_set = readings[0][0:4*size+2]
    first_set = readings[0][(4*size)+3:(8*size)+6]
    print(readings)    
    a = int(first_set[3*size + 2: 4*size + 2], 2)
    b = int(first_set[2*size + 1: 3*size + 2], 2)
    bout = int(second_set[2*size + 1: 3*size + 2], 2)
             
    print('a     =  ' + str(a))
    print('b     = ' + str(b))
    print('b out = ' + str(bout))
    print('\n\n')
    assert((a+b)%N == bout)
    #assert(carry == 0)
             
@given(draw_pair_of_ints())
@settings(deadline=None)
def test_add_using_VBE_Adder(vals):
    a = vals[0]
    b = vals[1]
    size = vals[2]
    
    qc = QuantumCircuit(3*size + 1)

    binA = bin(a)[2:].zfill(size)[::-1]
    binB = bin(b)[2:].zfill(size)[::-1]
    
    print("size " + str(size))
    print("a " + str(a))
    print(binA)
    print("b " + str(b))
    print(binB)
    
    
    # a = 3,2,1,0
    for i in range(size):
        if binA[i] == "1":
            qc.x(i)
        

    # b = 8,7,6,5,4
    for i in range(size):
        if binB[i] == "1":
            qc.x(i+size)
        
    # carry = 12,11,10,9 

    qc.measure_all()

    qc.append(VBE_Adder(size), [i for i in range(3*size+1)])
    
    qc.measure_all()

    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=1, memory=True)
    readings = job.result().get_memory()
    second_set = readings[0][0:3*size+1]
    first_set = readings[0][(3*size)+2:(6*size)+4]
    print(readings)

    a = int(first_set[2*size + 1: 3*size + 1],2)
    b = int(first_set[size: 2*size + 1],2)
    bout = int(second_set[size: 2*size + 1],2)
    carry = int(second_set[0: size],2)
    
    print('a     =  ' + str(a))
    print('b     = ' + str(b))
    print('b out = ' + str(bout))
    print('\n\n')
    assert(a+b == bout)
    assert(carry == 0)
    
@given(draw_pair_of_ints())
@settings(deadline=None)
def test_subtract_using_VBE_Adder_inverse(vals):
    a = vals[0]
    b = vals[1]
    size = vals[2]
    
    qc = QuantumCircuit(3*size + 1)

    binA = bin(a)[2:].zfill(size)[::-1]
    binB = bin(b)[2:].zfill(size)[::-1]
    
    print("size " + str(size))
    print("a " + str(a))
    print(binA)
    print("b " + str(b))
    print(binB)
    
    
    # a = 3,2,1,0
    for i in range(size):
        if binA[i] == "1":
            qc.x(i)
        

    # b = 8,7,6,5,4
    for i in range(size):
        if binB[i] == "1":
            qc.x(i+size)
        
    # carry = 12,11,10,9 

    qc.measure_all()

    qc.append(VBE_Adder(size).inverse(), [i for i in range(3*size+1)])
    
    qc.measure_all()

    backend = Aer.get_backend('aer_simulator') 
    job = execute(qc, backend, shots=1, memory=True)
    readings = job.result().get_memory()
    second_set = readings[0][0:3*size+1]
    first_set = readings[0][(3*size)+2:(6*size)+4]
    print(readings)

    a = int(first_set[2*size + 1: 3*size + 1],2)
    b = int(first_set[size: 2*size + 1],2)
    if (b < a):
        bout = twos_comp_to_int(second_set[size: 2*size + 1])
    else:
        bout = int(second_set[size: 2*size + 1],2)
    carry = int(second_set[0: size],2)
    
    print('a     =  ' + str(a))
    print('b     = ' + str(b))
    print('b out = ' + str(bout))
    print('\n\n')
    assert(b-a == bout)
    assert(carry == 0)

if __name__ == '__main__':
    #test_shors_with_postpro()
    #test_modular_exponentiation()
    #test_modular_multiplier_then_inverse_returns_original_values()
    test_modular_multiplier()
    #test_modular_adder()
    #test_add_using_VBE_Adder()
    #test_subtract_using_VBE_Adder_inverse()

