import unittest

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

import matplotlib.pyplot as plt
import numpy as np
import math
import cirq
from cirq.ops import H, X, I
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_cirq(qc, qubits, n):
    for qubit in range(n//2):
        qc.append(cirq.SWAP(qubits[qubit], qubits[n-qubit-1]))
    for j in range(n):
        for m in range(j):
            qc.append((cirq.CZ**(-1/2**(j-m)))(qubits[m],qubits[j]))
        qc.append(cirq.H(qubits[j]))

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, qubits):
    #outputs a + b
    #inverse outputs b - a (but in two's complement)
    #move to carry
    #qubits = cirq.LineQubit.range(3*qubit_size + 1)
    qc = cirq.Circuit()     
    
    for i in range(qubit_size-1):
        qc.append(cirq.CCNOT(qubits[i], qubits[i+qubit_size], qubits[2+i+(2*qubit_size)]))
    
    #toffoli to second register
    qc.append(cirq.CCNOT(qubits[qubit_size-1], qubits[(2*qubit_size)-1], qubits[2*qubit_size]))
    
    #add a to b
    for i in range(qubit_size):
        qc.append(cirq.CNOT(qubits[i], qubits[i+qubit_size]))
    
    #second register carry and last b
    for i in range(qubit_size-1):
        qc.append(cirq.CCNOT(qubits[i+qubit_size], qubits[1+i+(2*qubit_size)], qubits[2+i+(2*qubit_size)]))
        
    qc.append(cirq.CCNOT(qubits[(2*qubit_size)-1], qubits[3*qubit_size], qubits[2*qubit_size]))
    
    qc.append(cirq.CNOT(qubits[3*qubit_size], qubits[(2*qubit_size)-1]))
    
    #adder overflow
    for i in range(qubit_size-1):
        qc.append(cirq.CCNOT(qubits[(2*qubit_size) - 2 - i], qubits[(3*qubit_size) - i - 1], qubits[(3*qubit_size) - i])) 
        qc.append(cirq.CNOT(qubits[(qubit_size) - 2 - i], qubits[(2*qubit_size) - 2 - i]))
        qc.append(cirq.CCNOT(qubits[(qubit_size) - 2 - i], qubits[(2*qubit_size) - 2 - i], qubits[(3*qubit_size) - i])) 
        qc.append(cirq.CNOT(qubits[(qubit_size) - 2 - i], qubits[(2*qubit_size) - 2 - i])) 
        qc.append(cirq.CNOT(qubits[(3*qubit_size) - i - 1], qubits[(2*qubit_size) - 2 - i])) 
    
    #print(qc.draw(output='text'))
    qc.name = "VBE ADDER"
    return qc    

#VBE_Adder(4)

def Modular_adder(qubit_size, N, qubits):
    binN = bin(N)[2:].zfill(qubit_size)[::-1]
        
    #move to carry
    #qubits = cirq.LineQubit.range(4*qubit_size + 2)
    qc = cirq.Circuit()     
    
    qc.append(VBE_Adder(qubit_size, qubits[0:3*qubit_size+1]))
        
    for i in range(qubit_size):
        qc.append(cirq.SWAP(qubits[i], qubits[i + 1 + 3*qubit_size]))
    
    qc.append(cirq.inverse(VBE_Adder(qubit_size, qubits[0:3*qubit_size+1])))
    
    qc.append(cirq.X(qubits[2*qubit_size]))
    qc.append(cirq.CNOT(qubits[2*qubit_size], qubits[4*qubit_size+1]))
    qc.append(cirq.X(qubits[2*qubit_size]))
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.append(cirq.CNOT(qubits[4*qubit_size + 1], qubits[i])) 
    
    qc.append(VBE_Adder(qubit_size, qubits[0:3*qubit_size+1]))
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.append(cirq.CNOT(qubits[4*qubit_size + 1], qubits[i])) 
    
    for i in range(qubit_size):
        qc.append(cirq.SWAP(qubits[i], qubits[i + 1 + 3*qubit_size]))
    
    qc.append(cirq.inverse(VBE_Adder(qubit_size, qubits[0:3*qubit_size+1])))
    
    qc.append(cirq.CNOT(qubits[2*qubit_size], qubits[4*qubit_size+1]))
    
    qc.append(VBE_Adder(qubit_size, qubits[0:3*qubit_size+1]))
    
    qc.name = "MODULAR ADDER"
    
    return qc  

#Modular_adder(4,11,qubits)

# 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, qubits):        
    qc = cirq.Circuit()
    
    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.append(cirq.CCNOT(qubits[0], qubits[i+1], qubits[qubit_size + 1 + j]))

        qc.append(Modular_adder(qubit_size, N, qubits[qubit_size + 1:(5*qubit_size+3)]))
             
        #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.append(cirq.CCNOT(qubits[0], qubits[i+1], qubits[2*qubit_size - j])) 
                
    #set x if control is false
    qc.append(cirq.X(qubits[0]))
    for i in range(qubit_size):
        qc.append(cirq.CCNOT(qubits[0], qubits[i+1], qubits[2*qubit_size + 1 + i])) 
    qc.append(cirq.X(qubits[0]))
        
    qc.name = "MODULAR MULTIPLIER " + str(a) + "x mod " + str(N)  
    
    #print(qc.inverse().draw(output='text'))
    return qc  

#Modular_multiplier(4, 15, 13)

size = 2
N = 2
a = 1

qubits = cirq.LineQubit.range(5*size + 3)
qc = cirq.Circuit()

#control qubit
qc.append(cirq.X(qubits[0]))

# x = 4,3,2,1
#qc.append(cirq.X(qubits[4]))
#qc.append(cirq.X(qubits[3]))
#qc.append(cirq.X(qubits[2]))
#qc.append(cirq.X(qubits[1]))

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

# temporary carry = 17 

#print(qc)

qc.append(cirq.measure(*qubits[0:5*size+3], key='m1'))

qc.append(Modular_multiplier(size, N, a, qubits[0:(5*size+3)]))

#print(qc)

qc.append(cirq.measure(*qubits[0:5*size+3], key='m2'))

simulator = cirq.Simulator()
results = simulator.run(qc , repetitions =1)
first_set = np.array2string(results.measurements['m1'][0], separator='')[1:-1][::-1]
second_set = np.array2string(results.measurements['m2'][0], separator='')[1:-1][::-1]

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, qubits):            
    qc = cirq.Circuit()
    
    #need to set x to 1
    qc.append(cirq.X(qubits[qubit_size+1]))
        
    for i in range(qubit_size):
        qc.append(cirq.CNOT(qubits[i],qubits[qubit_size]))
        qc.append(Modular_multiplier(qubit_size, N, a**(2**i)%N, qubits[qubit_size:(6*qubit_size+3)]))
        for j in range(qubit_size):
            qc.append(cirq.SWAP(qubits[qubit_size+1+j], qubits[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(cirq.inverse(Modular_multiplier(qubit_size, N, pow(a**(2**i)%N, -1, N), qubits[qubit_size:(6*qubit_size+3)])))
        qc.append(cirq.CNOT(qubits[i],qubits[qubit_size]))
                            
    qc.name = "MODULAR EXPONENTIATION"    
        
    #print(qc.draw(output='text'))
    return qc  

#qubits = cirq.LineQubit.range(6*3 + 3)
#Modular_exponentiation(3, 7, 3,qubits)

# 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):
    qubits = cirq.LineQubit.range(6*qubit_size + 3)
    qc = cirq.Circuit()
    
    # N = 21,20,19,18
    binN = bin(N)[2:].zfill(qubit_size)[::-1]

    for i in range(qubit_size):
        qc.append(cirq.H(qubits[i])) 
    
    for i in range(qubit_size):
        if binN[i] == "1":
            qc.append(cirq.X(qubits[5*qubit_size + 2 + i])) 

    qc.append(Modular_exponentiation(qubit_size, N, a, qubits[0:(6*qubit_size+3)]))

    qft_dagger_cirq(qc, qubits[:qubit_size], qubit_size)
    
    qc.append(cirq.measure(*qubits[:qubit_size]))
    
    #backend = Aer.get_backend('aer_simulator') 
    #job = execute(qc, backend, shots=run_amt, memory=True)
    # return(job.result().get_memory())
    
    
    simulator = cirq.Simulator()
    results = simulator.run(qc , repetitions =run_amt)
    #readings = np.array2string(results.measurements['m'][0], separator='')[1:-1][::-1]
    return  results
    
#print(run_shors(3,7,3,1))

#%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
        print(r)
        #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(7, 3, 4)

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=4)) 
    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:
#         qubits = cirq.LineQubit.range(6*size + 3)
#         qc = cirq.Circuit()

#         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]
    
#     qubits = cirq.LineQubit.range(5*size + 3)
#     qc = cirq.Circuit()

#     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]

    print(a)
    print(x)
    print(size)
    print(N)
    
    qubits = cirq.LineQubit.range(5*size + 3)
    qc = cirq.Circuit()

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

    qc.append(cirq.measure(*qubits[0:5*size + 3], key='m1'))

    qc.append(Modular_multiplier(size, N, a, qubits[0:(5*size+3)]))

    qc.append(cirq.measure(*qubits[0:5*size + 3], key='m2'))

    simulator = cirq.Simulator()
    results = simulator.run(qc , repetitions =1)
    first_set = np.array2string(results.measurements['m1'][0], separator='')[1:-1][::-1]
    second_set = np.array2string(results.measurements['m2'][0], separator='')[1:-1][::-1]
    
    print("f " + first_set)
    print('s ' + second_set)
    
    x = int(first_set[4*size + 2: 5*size + 2], 2)
    print(first_set[4*size + 2: 5*size + 2])
    b = int(first_set[2*size + 1: 3*size + 2], 2)
    bout = int(second_set[2*size + 1: 3*size + 2], 2)
       
    
    print("control " + str(control))
    print("size  = " + str(size))
    print('b     = ' + str(b))
    print('a     = ' + str(a))    
    print('x     = ' + str(x))
    print("N     = " + str(N))
    print('b out = ' + str(bout))
    print('\n\n')
    if (control == False):
        assert(bout == x)
    else: 
        assert((a * x) % N == bout)

@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]
    
    qubits = cirq.LineQubit.range(4*size + 2)
    qc = cirq.Circuit()

    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.append(cirq.X(qubits[i]))
        

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

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

    # temporary carry = 17 

    #print(qc)

    qc.append(cirq.measure(*qubits[0:4*size+2], key='m1'))
    
    qc.append(Modular_adder(size,N, qubits[0:4*size+2]))

    qc.append(cirq.measure(*qubits[0:4*size+2], key='m2'))
    
    simulator = cirq.Simulator()
    results = simulator.run(qc , repetitions =1)
    first_set = np.array2string(results.measurements['m1'][0], separator='')[1:-1][::-1]
    second_set = np.array2string(results.measurements['m2'][0], separator='')[1:-1][::-1]
    
    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)
             
@given(draw_pair_of_ints())
@settings(deadline=None)
def test_add_using_VBE_Adder(vals):
    a = vals[0]
    b = vals[1]
    size = vals[2]
    
    qubits = cirq.LineQubit.range(3*size + 1)
    qc = cirq.Circuit()

    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.append(cirq.X(qubits[i]))
        

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

    qc.append(cirq.measure(*qubits[0:3*size+1], key='m1'))

    qc.append(VBE_Adder(size, qubits[0:3*size+1]))
    
    qc.append(cirq.measure(*qubits[0:3*size+1], key='m2'))

    simulator = cirq.Simulator()
    results = simulator.run(qc , repetitions =1)
    first_set = np.array2string(results.measurements['m1'][0], separator='')[1:-1][::-1]
    second_set = np.array2string(results.measurements['m2'][0], separator='')[1:-1][::-1]

    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]
    
    qubits = cirq.LineQubit.range(3*size + 1)
    qc = cirq.Circuit()

    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.append(cirq.X(qubits[i]))
        

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

    qc.append(cirq.measure(*qubits[0:3*size+1], key='m1'))

    qc.append(cirq.inverse(VBE_Adder(size,qubits[0:3*size+1])))
    
    qc.append(cirq.measure(*qubits[0:3*size+1], key='m2'))

    simulator = cirq.Simulator()
    results = simulator.run(qc , repetitions =1)
    first_set = np.array2string(results.measurements['m1'][0], separator='')[1:-1][::-1]
    second_set = np.array2string(results.measurements['m2'][0], separator='')[1:-1][::-1]
        
    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()

