# checking the version of PYTHON; only support > 3.5
import sys
if sys.version_info < (3,5):
    raise Exception('Please use Python version 3.5 or greater.')
    
# importing QISKit
from qiskit import QuantumProgram
#import Qconfig  

# import basic plotting tools
from qiskit.tools.visualization import plot_histogram
from qiskit.tools.qcvv.tomography import marginal_counts

# create a QuantumProgram object instance.
q_program = QuantumProgram()
#q_program.set_api(Qconfig.APItoken, Qconfig.config["url"]) # set the APIToken and API url

# backend
#backend = 'ibmqx2'  
backend = 'local_qasm_simulator'

# quantum plain addition algorithm for 1-qubit numbers    
def addition_1bit(circuit, q):
    circuit.h(q[2])
    circuit.cx(q[1], q[2])
    circuit.tdg(q[2])
    circuit.cx(q[0], q[2])
    circuit.t(q[2])
    circuit.cx(q[1], q[2])
    circuit.tdg(q[2])
    circuit.cx(q[0], q[2])
    circuit.t(q[2])
    circuit.h(q[2])
    circuit.t(q[1])
    circuit.cx(q[0], q[1])
    circuit.t(q[0])
    circuit.tdg(q[1])
    
# n-qubit number input state
def number_state(circuit, q, a, b):
    if a == 1:
        circuit.x(q[0]) # q[0] contains the value of a
    if b == 1:
        circuit.x(q[1]) # q[1] contain the value of b

# we define the values (0 or 1)
a = 1                 
b = 1

# one single quantum register which contains 'a' (1 qubit) and 'b' (2 qubits)
q = q_program.create_quantum_register("q", 3) # 3 qubits
# clasical register
c = q_program.create_classical_register("cr", 2) # 2 bits

# quantum circuit involving the quantum register and the classical register
add1bit_circuit = q_program.create_circuit("add", [q] ,[c])

# create the state containing a and b
number_state(add1bit_circuit, q, a, b)

# addition
addition_1bit(add1bit_circuit, q)

# measurements to see the result, which has been written in b (q[1]q[2])
add1bit_circuit.measure(q[1], c[0])
add1bit_circuit.measure(q[2], c[1])

# compile and execute the quantum program in the backend
result = q_program.execute(["add"], backend=backend, shots=1024, max_credits=3, wait=10, timeout=999999)

# show the results
print(result)
print(result.get_data("add"))
counts = marginal_counts(result.get_counts("add"), [0, 1])
plot_histogram(counts)
print("Backend:", backend)
print("Highest probability outcome: {}".format(int(max(counts, key = lambda x: counts[x]).replace(" ", ""), 2)))

# checking the version of PYTHON; only support > 3.5
import sys
if sys.version_info < (3,5):
    raise Exception('Please use Python version 3.5 or greater.')
    
# importing QISKit
from qiskit import QuantumProgram

# import basic plotting tools
from qiskit.tools.visualization import plot_histogram

# create a QuantumProgram object instance.
q_program = QuantumProgram()

# backend
backend = 'local_qasm_simulator'

def carry(circuit, q0, q1, q2, q3):
    "carry module"
    circuit.ccx(q1, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q0, q2, q3)
    
def carry_inv(circuit, q0, q1, q2, q3):
    "carry module but running backwards"
    circuit.ccx(q0, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q1, q2, q3)
    
def summation(circuit, q0, q1, q2):
    "summation module"
    circuit.cx(q1, q2)
    circuit.cx(q0, q2)
    
# quantum plain addition algorithm for n-qubit numbers    
def addition_nbit(circuit, qa, qb, qcar, n):
    if n == 1:
        circuit.ccx(qa[0], qb[0], qb[1])
        circuit.cx(qa[0], qb[0])
    else:
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
        carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        circuit.cx(qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            summation(circuit, qcar[i-1], qa[i], qb[i])
            carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1])
        summation(circuit, qcar[0], qa[1], qb[1])
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])

# n-qubit number input state
def number_state(circuit, q, x, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            circuit.x(q[i])

# we define the values
a = 9                 
b = 14

# computing the number of qubits n needed
n = len("{0:b}".format(a))
n2 = len("{0:b}".format(b))
if n2 > n:
    n = n2
    
# classical register with n+1 bits.
c = q_program.create_classical_register("cr", n+1)
# quantum registers
qa = q_program.create_quantum_register("qa", n) # a qubits
qb = q_program.create_quantum_register("qb", n+1) # b qubits
# if n = 1, no need of carry register 
if n == 1:
    qcar = 0
    # quantum circuit involving the quantum registers and the classical register
    addnbit_circuit = q_program.create_circuit("add", [qa, qb],[c])
else:
    qcar = q_program.create_quantum_register("qcar", n-1) # carry qubits
    # quantum circuit involving the quantum registers and the classical register
    addnbit_circuit = q_program.create_circuit("add", [qa, qb, qcar],[c])

# create the state containing a
number_state(addnbit_circuit, qa, a, n)
# create the state containing b
number_state(addnbit_circuit, qb, b, n)

# addition
addition_nbit(addnbit_circuit, qa, qb, qcar, n)

# measurements to see the result
for i in range(n+1):
    addnbit_circuit.measure(qb[i], c[i])

# compile and execute the quantum program in the backend
result = q_program.execute(["add"], backend=backend, shots=1024, timeout=999999)

# show the results.
print(result)
print(result.get_data("add"))
counts = result.get_counts("add")
plot_histogram(counts)
print("Backend:", backend)
print("Highest probability outcome: {}".format(int(max(counts, key = lambda x: counts[x]).replace(" ", ""), 2)))

# checking the version of PYTHON; only support > 3.5
import sys
if sys.version_info < (3,5):
    raise Exception('Please use Python version 3.5 or greater.')
    
# importing QISKit
from qiskit import QuantumProgram

# import basic plotting tools
from qiskit.tools.visualization import plot_histogram

# create a QuantumProgram object instance.
q_program = QuantumProgram()

# backend
backend = 'local_qasm_simulator'

def carry(circuit, q0, q1, q2, q3):
    "carry module"
    circuit.ccx(q1, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q0, q2, q3)
    
def carry_inv(circuit, q0, q1, q2, q3):
    "carry module running backwards"
    circuit.ccx(q0, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q1, q2, q3)
    
def summation(circuit, q0, q1, q2):
    "summation module"
    circuit.cx(q1, q2)
    circuit.cx(q0, q2)
    
def summation_inv(circuit, q0, q1, q2):
    "summation module running backwards"
    circuit.cx(q0, q2)
    circuit.cx(q1, q2)

# quantum plain addition algorithm for n-qubit numbers    
def addition_nbit(circuit, qa, qb, qcar, n):
    if n == 1:
        circuit.ccx(qa[0], qb[0], qb[1])
        circuit.cx(qa[0], qb[0])
    else:
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
        carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        circuit.cx(qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            summation(circuit, qcar[i-1], qa[i], qb[i])
            carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1])
        summation(circuit, qcar[0], qa[1], qb[1])
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
    
# quantum plain substraction algorithm for n-qubit numbers
def subs_nbit(circuit, qa, qb, qcar, n):
    "same circuit as the plain addition but going backwards"
    if n == 1:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qb[1])
    else:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        summation_inv(circuit, qcar[0], qa[1], qb[1])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
            summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2])
        circuit.cx(qa[n-1], qb[n-1])
        carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        for i in range(n-2, 0, -1):
            carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i])
        circuit.cx(qa[0], qb[0])  
        circuit.ccx(qa[0], qb[0], qcar[0])

def cond_toffoli(circuit, qcond, q1, q2, q3):
    "toffoli gate conditioned by an external qubit"
    circuit.h(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.h(q3)
    circuit.t(q2)
    circuit.ccx(qcond, q1, q2)
    circuit.t(q1)
    circuit.tdg(q2)
    circuit.ccx(qcond, q1, q2)
    
def cond_carry(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module"
    cond_toffoli(circuit, qcond, q1, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q0, q2, q3)
    
def cond_carry_inv(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module running backwards"
    cond_toffoli(circuit, qcond, q0, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q1, q2, q3)
    
def cond_summation(circuit, q0, q1, q2, qcond):
    "conditional summation module"
    circuit.ccx(qcond, q1, q2)
    circuit.ccx(qcond, q0, q2)
    
def cond_summation_inv(circuit, q0, q1, q2, qcond):
    "conditional summation module running backwards"
    circuit.ccx(qcond, q0, q2)
    circuit.ccx(qcond, q1, q2)
    
# quantum conditional plain addition algorithm for n-qubit numbers    
def cond_addition_nbit(circuit, qa, qb, qcar, qcond, n):
    "plain addition algorithm conditioned by an external qubit"
    if n == 1:
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qb[1])
        circuit.ccx(qcond[0], qa[0], qb[0])
    else:
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qcar[0])
        circuit.ccx(qcond[0], qa[0], qb[0])
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond[0])
        cond_carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond[0])
        circuit.ccx(qcond[0], qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            cond_summation(circuit, qcar[i-1], qa[i], qb[i], qcond[0])
            cond_carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1], qcond[0])
        cond_summation(circuit, qcar[0], qa[1], qb[1], qcond[0])
        circuit.ccx(qcond[0], qa[0], qb[0])
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qcar[0])
        circuit.ccx(qcond[0], qa[0], qb[0])

# quantum conditional plain substraction algorithm for n-qubit numbers      
def cond_subs_nbit(circuit, qa, qb, qcar, qcond, n):
    "same circuit as the conditional plain addition but going backwards"
    if n == 1:
        circuit.ccx(qcond[0], qa[0], qb[0])
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qb[1])
    else:
        circuit.ccx(qcond[0], qa[0], qb[0])
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qcar[0])
        circuit.ccx(qcond[0], qa[0], qb[0])
        cond_summation_inv(circuit, qcar[0], qa[1], qb[1], qcond[0])
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond[0])
            cond_summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2], qcond[0])
        circuit.ccx(qcond[0], qa[n-1], qb[n-1])
        cond_carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond[0])
        for i in range(n-2, 0, -1):
            cond_carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i], qcond[0])
        circuit.ccx(qcond[0], qa[0], qb[0])  
        cond_toffoli(circuit, qcond[0], qa[0], qb[0], qcar[0])
        
# quantum modular addition algorithm for n-qubit numbers 
def mod_addition_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, n):
    addition_nbit(circuit, qa, qb, qcar, n)
    subs_nbit(circuit, qN, qb, qcar, n)
    circuit.x(qb[n])
    circuit.cx(qb[n], qtemp[0])
    circuit.x(qb[n])
    cond_subs_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    addition_nbit(circuit, qN, qb, qcar, n)
    cond_addition_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    subs_nbit(circuit, qa, qb, qcar, n)
    circuit.cx(qb[n], qtemp[0])
    addition_nbit(circuit, qa, qb, qcar, n)

# n-qubit number input state
def number_state(circuit, q, x, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            circuit.x(q[i])

# we define the values
a = 2
b = 3                     
N = 3

# computing number of qubits n needed
n = len("{0:b}".format(a))
n2 = len("{0:b}".format(b))
n3 = len("{0:b}".format(N))
if n2 > n:
    n = n2
if n3 > n:
    n = n3
    
# classical register with n+1 bits.
c = q_program.create_classical_register("cr", n+1)
# quantum registers
qa = q_program.create_quantum_register("qa", n) # a qubits
qb = q_program.create_quantum_register("qb", n+1) # b qubits
qN = q_program.create_quantum_register("qN", n+1) # N qubits
qNtemp = q_program.create_quantum_register("qNtemp", n) # temporary N qubits
qtemp = q_program.create_quantum_register("qtemp", 1) # temporary qubit
# if n = 1, no need of carry register 
if n == 1:
    qcar = 0
    # quantum circuit involving the quantum registers and the classical register
    mod_add_circuit = q_program.create_circuit("mod_add", [qa, qb, qN, qNtemp, qtemp],[c])
else:
    qcar = q_program.create_quantum_register("qcar", n-1) # carry qubits
    # quantum circuit involving the quantum registers and the classical register
    mod_add_circuit = q_program.create_circuit("mod_add", [qa, qb, qN, qcar, qNtemp, qtemp],[c])

# create the state containing 'a'
number_state(mod_add_circuit, qa, a, n)
# create the state containing 'b'
number_state(mod_add_circuit, qb, b, n)
# create the state containing 'N'
number_state(mod_add_circuit, qN, N, n)
# create the temporary state containing 'N'
number_state(mod_add_circuit, qNtemp, N, n)

# modular addition
mod_addition_nbit(mod_add_circuit, qa, qb, qN, qNtemp, qcar, qtemp, n)

# measurements to see the result
for i in range(n+1):
    mod_add_circuit.measure(qb[i], c[i])

# compile and execute the quantum program in the backend
result = q_program.execute(["mod_add"], backend=backend, shots=1024, max_credits=3, wait=10, timeout=999999)

# show the results.
print(result)
print(result.get_data("mod_add"))
counts = result.get_counts("mod_add")
plot_histogram(counts)
print("Backend:", backend)
print("Highest probability outcome: {}".format(int(max(counts, key = lambda x: counts[x]).replace(" ", ""), 2)))

# checking the version of PYTHON; only support > 3.5
import sys
if sys.version_info < (3,5):
    raise Exception('Please use Python version 3.5 or greater.')
    
# importing QISKit
from qiskit import QuantumProgram

# import basic plotting tools
from qiskit.tools.visualization import plot_histogram
 
# create a QuantumProgram object instance.
q_program = QuantumProgram()

# backend
backend = 'local_qasm_simulator'

def carry(circuit, q0, q1, q2, q3):
    "carry module"
    circuit.ccx(q1, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q0, q2, q3)
    
def carry_inv(circuit, q0, q1, q2, q3):
    "carry module running backwards"
    circuit.ccx(q0, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q1, q2, q3)
    
def summation(circuit, q0, q1, q2):
    "summation module"
    circuit.cx(q1, q2)
    circuit.cx(q0, q2)
    
def summation_inv(circuit, q0, q1, q2):
    "summation module running backwards"
    circuit.cx(q0, q2)
    circuit.cx(q1, q2)

# quantum plain addition algorithm for n-qubit numbers    
def addition_nbit(circuit, qa, qb, qcar, n):
    if n == 1:
        circuit.ccx(qa[0], qb[0], qb[1])
        circuit.cx(qa[0], qb[0])
    else:
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
        carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        circuit.cx(qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            summation(circuit, qcar[i-1], qa[i], qb[i])
            carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1])
        summation(circuit, qcar[0], qa[1], qb[1])
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
    
# quantum plain substraction algorithm for n-qubit numbers
def subs_nbit(circuit, qa, qb, qcar, n):
    "same circuit as the addition but going backwards"
    if n == 1:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qb[1])
    else:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        summation_inv(circuit, qcar[0], qa[1], qb[1])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
            summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2])
        circuit.cx(qa[n-1], qb[n-1])
        carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        for i in range(n-2, 0, -1):
            carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i])
        circuit.cx(qa[0], qb[0])  
        circuit.ccx(qa[0], qb[0], qcar[0])

def cond_toffoli(circuit, qcond, q1, q2, q3):
    "toffoli gate conditioned by an external qubit"
    circuit.h(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.h(q3)
    circuit.t(q2)
    circuit.ccx(qcond, q1, q2)
    circuit.t(q1)
    circuit.tdg(q2)
    circuit.ccx(qcond, q1, q2)
    
def cond_carry(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module"
    cond_toffoli(circuit, qcond, q1, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q0, q2, q3)
    
def cond_carry_inv(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module running backwards"
    cond_toffoli(circuit, qcond, q0, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q1, q2, q3)
    
def cond_summation(circuit, q0, q1, q2, qcond):
    "conditional summation module"
    circuit.ccx(qcond, q1, q2)
    circuit.ccx(qcond, q0, q2)
    
def cond_summation_inv(circuit, q0, q1, q2, qcond):
    "conditional summation module running backwards"
    circuit.ccx(qcond, q0, q2)
    circuit.ccx(qcond, q1, q2)
    
# quantum conditional plain addition algorithm for n-qubit numbers    
def cond_addition_nbit(circuit, qa, qb, qcar, qcond, n):
    if n == 1:
        cond_toffoli(circuit, qcond, qa[0], qb[0], qb[1])
        circuit.ccx(qcond, qa[0], qb[0])
    else:
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond)
        cond_carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond)
        circuit.ccx(qcond, qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            cond_summation(circuit, qcar[i-1], qa[i], qb[i], qcond)
            cond_carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1], qcond)
        cond_summation(circuit, qcar[0], qa[1], qb[1], qcond)
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])

# quantum conditional plain substraction algorithm for n-qubit numbers      
def cond_subs_nbit(circuit, qa, qb, qcar, qcond, n):
    "same circuit as the conditional plain addition but going backwards"
    if n == 1:
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qb[1])
    else:
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])
        cond_summation_inv(circuit, qcar[0], qa[1], qb[1], qcond)
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond)
            cond_summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2], qcond)
        circuit.ccx(qcond, qa[n-1], qb[n-1])
        cond_carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond)
        for i in range(n-2, 0, -1):
            cond_carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i], qcond)
        circuit.ccx(qcond, qa[0], qb[0])  
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        
# quantum modular addition algorithm for n-qubit numbers 
def mod_addition_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, n):
    addition_nbit(circuit, qa, qb, qcar, n)
    subs_nbit(circuit, qN, qb, qcar, n)
    circuit.x(qb[n])
    circuit.cx(qb[n], qtemp)
    circuit.x(qb[n])
    cond_subs_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    addition_nbit(circuit, qN, qb, qcar, n)
    cond_addition_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    subs_nbit(circuit, qa, qb, qcar, n)
    circuit.cx(qb[n], qtemp)
    addition_nbit(circuit, qa, qb, qcar, n)
    
# quantum controlled modular multiplication algorithm for n-qubit numbers
def cont_mod_mult_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, qX, qext, N, n):
    for i in range(n):
        for j in range(n):
            classical_mod = (2**(i+j))%N
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)
            mod_addition_nbit(circuit, qtempst, qb, qN, qNtemp, qcar, qtemp, n)
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)
    circuit.x(qext)
    cond_addition_nbit(circuit, qX, qb, qcar, qext, n)
    circuit.x(qext)

# n-qubit number input state
def number_state(circuit, q, x, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            circuit.x(q[i])

# n-qubit number input state, controlled by 2 control qubits
def cond_number_state(circuit, q, x, ext, control1, control2, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            cond_toffoli(circuit, ext, control1, control2, q[i])

# we define the values
a = 1
x = 1
N = 1


# computing number of qubits n needed
n = len("{0:b}".format(a))
n2 = len("{0:b}".format(x))
n3 = len("{0:b}".format(N))
if n2 > n:
    n = n2
if n3 > n:
    n = n3

# classical register with n+1 bits.
c = q_program.create_classical_register("cr", n+1)
# quantum registers
qa = q_program.create_quantum_register("qa", n) # a qubits
qb = q_program.create_quantum_register("qb", n+1) # result register
qN = q_program.create_quantum_register("qN", n+1) # N qubits
qNtemp = q_program.create_quantum_register("qNtemp", n) # temporary N qubits
qtemp = q_program.create_quantum_register("qtemp", 1) # temporary qubit
qtempst = q_program.create_quantum_register("qtempst", n) # temporary register
qX = q_program.create_quantum_register("qX", n) # x register
qext = q_program.create_quantum_register("qext", 1)
# if n = 1, no need of carry register 
if n == 1:
    qcar = 0
    # quantum circuit involving the quantum registers and the classical register
    mod_mult_circuit = q_program.create_circuit("mod_mult", [qa, qb, qN, qNtemp, qtemp, qtempst, qX, qext],[c])
else:
    qcar = q_program.create_quantum_register("qcar", n-1) # carry qubits
    # quantum circuit involving the quantum register and the classical register
    mod_mult_circuit = q_program.create_circuit("mod_mult", [qa, qb, qN, qcar, qNtemp, qtemp, qtempst, qX, qext],[c])

# create the state containing 'a'
number_state(mod_mult_circuit, qa, a, n)
# create the state containing 'b'
number_state(mod_mult_circuit, qX, x, n)
# create the state containing 'N'
number_state(mod_mult_circuit, qN, N, n+1)
# create a temporary state containing 'N'
number_state(mod_mult_circuit, qNtemp, N, n)

mod_mult_circuit.x(qext[0]) # we set the control qubit to |1>

# controlled modular multiplication
cont_mod_mult_nbit(mod_mult_circuit, qa, qb, qN, qNtemp, qcar, qtemp[0], qtempst, qX, qext[0], N, n)

# measurements to see the result
for i in range(n+1):
    mod_mult_circuit.measure(qb[i], c[i])

# compile and execute the quantum program in the backend
result = q_program.execute(["mod_mult"], backend=backend, shots=1024, max_credits=3, wait=10, timeout=999999)

# show the results.
print(result)
print(result.get_data("mod_mult"))
counts = result.get_counts("mod_mult")
plot_histogram(counts)
print("Backend:", backend)
print("Highest probability outcome: {}".format(int(max(counts, key = lambda x: counts[x]).replace(" ", ""), 2)))

# checking the version of PYTHON; only support > 3.5
import sys
if sys.version_info < (3,5):
    raise Exception('Please use Python version 3.5 or greater.')
    
# importing QISKit
from qiskit import QuantumProgram

# import basic plotting tools
from qiskit.tools.visualization import plot_histogram
 
# create a QuantumProgram object instance.
q_program = QuantumProgram()

# backend
backend = 'local_qasm_simulator'

def carry(circuit, q0, q1, q2, q3):
    "carry module"
    circuit.ccx(q1, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q0, q2, q3)
    
def carry_inv(circuit, q0, q1, q2, q3):
    "carry module running backwards"
    circuit.ccx(q0, q2, q3)
    circuit.cx(q1, q2)
    circuit.ccx(q1, q2, q3)
    
def summation(circuit, q0, q1, q2):
    "summation module"
    circuit.cx(q1, q2)
    circuit.cx(q0, q2)
    
def summation_inv(circuit, q0, q1, q2):
    "summation module running backwards"
    circuit.cx(q0, q2)
    circuit.cx(q1, q2)

# quantum plain addition algorithm for n-qubit numbers    
def addition_nbit(circuit, qa, qb, qcar, n):
    if n == 1:
        circuit.ccx(qa[0], qb[0], qb[1])
        circuit.cx(qa[0], qb[0])
    else:
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
        carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        circuit.cx(qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            summation(circuit, qcar[i-1], qa[i], qb[i])
            carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1])
        summation(circuit, qcar[0], qa[1], qb[1])
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
    
# quantum plain substraction algorithm for n-qubit numbers 
def subs_nbit(circuit, qa, qb, qcar, n):
    "same as the plain addition but running backwards"
    if n == 1:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qb[1])
    else:
        circuit.cx(qa[0], qb[0])
        circuit.ccx(qa[0], qb[0], qcar[0])
        circuit.cx(qa[0], qb[0])
        summation_inv(circuit, qcar[0], qa[1], qb[1])
        for i in range(n-2):
            carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1])
            summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2])
        circuit.cx(qa[n-1], qb[n-1])
        carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n])
        for i in range(n-2, 0, -1):
            carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i])
        circuit.cx(qa[0], qb[0])  
        circuit.ccx(qa[0], qb[0], qcar[0])
        
def cond_toffoli(circuit, qcond, q1, q2, q3):
    "conditional toffoli gate"
    circuit.h(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.ccx(qcond, q2, q3)
    circuit.tdg(q3)
    circuit.ccx(qcond, q1, q3)
    circuit.t(q3)
    circuit.h(q3)
    circuit.t(q2)
    circuit.ccx(qcond, q1, q2)
    circuit.t(q1)
    circuit.tdg(q2)
    circuit.ccx(qcond, q1, q2)

def cond_carry(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module"
    cond_toffoli(circuit, qcond, q1, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q0, q2, q3)
    
def cond_carry_inv(circuit, q0, q1, q2, q3, qcond):
    "conditional carry module running backwards"
    cond_toffoli(circuit, qcond, q0, q2, q3)
    circuit.ccx(qcond, q1, q2)
    cond_toffoli(circuit, qcond, q1, q2, q3)
    
def cond_summation(circuit, q0, q1, q2, qcond):
    "conditional summation"
    circuit.ccx(qcond, q1, q2)
    circuit.ccx(qcond, q0, q2)
    
def cond_summation_inv(circuit, q0, q1, q2, qcond):
    "conditional summation running backwards"
    circuit.ccx(qcond, q0, q2)
    circuit.ccx(qcond, q1, q2)
    
# quantum conditional plain addition algorithm for n-qubit numbers    
def cond_addition_nbit(circuit, qa, qb, qcar, qcond, n):
    if n == 1:
        cond_toffoli(circuit, qcond, qa[0], qb[0], qb[1])
        circuit.ccx(qcond, qa[0], qb[0])
    else:
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond)
        cond_carry(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond)
        circuit.ccx(qcond, qa[n-1], qb[n-1])
        for i in range(n-1, 1, -1):
            cond_summation(circuit, qcar[i-1], qa[i], qb[i], qcond)
            cond_carry_inv(circuit, qcar[i-2], qa[i-1], qb[i-1], qcar[i-1], qcond)
        cond_summation(circuit, qcar[0], qa[1], qb[1], qcond)
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])

# quantum conditional plain substraction algorithm for n-qubit numbers      
def cond_subs_nbit(circuit, qa, qb, qcar, qcond, n):
    "same as conditional plain addition but running backwards"
    if n == 1:
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qb[1])
    else:
        circuit.ccx(qcond, qa[0], qb[0])
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        circuit.ccx(qcond, qa[0], qb[0])
        cond_summation_inv(circuit, qcar[0], qa[1], qb[1], qcond)
        for i in range(n-2):
            cond_carry(circuit, qcar[i], qa[i+1], qb[i+1], qcar[i+1], qcond)
            cond_summation_inv(circuit, qcar[i+1], qa[i+2], qb[i+2], qcond)
        circuit.ccx(qcond, qa[n-1], qb[n-1])
        cond_carry_inv(circuit, qcar[n-2], qa[n-1], qb[n-1], qb[n], qcond)
        for i in range(n-2, 0, -1):
            cond_carry_inv(circuit, qcar[i-1], qa[i], qb[i], qcar[i], qcond)
        circuit.ccx(qcond, qa[0], qb[0])  
        cond_toffoli(circuit, qcond, qa[0], qb[0], qcar[0])
        
# quantum modular addition algorithm for n-qubit numbers 
def mod_addition_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, n):
    addition_nbit(circuit, qa, qb, qcar, n)
    subs_nbit(circuit, qN, qb, qcar, n)
    circuit.x(qb[n])
    circuit.cx(qb[n], qtemp)
    circuit.x(qb[n])
    cond_subs_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    addition_nbit(circuit, qN, qb, qcar, n)
    cond_addition_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    subs_nbit(circuit, qa, qb, qcar, n)
    circuit.cx(qb[n], qtemp)
    addition_nbit(circuit, qa, qb, qcar, n)
    
# quantum modular substraction algorithm for n-qubit numbers
def mod_subs_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, n):
    "same as modular addition but running backwards"
    subs_nbit(circuit, qa, qb, qcar, n)
    circuit.cx(qb[n], qtemp)
    addition_nbit(circuit, qa, qb, qcar, n)
    cond_subs_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    subs_nbit(circuit, qN, qb, qcar, n)
    cond_addition_nbit(circuit, qNtemp, qN, qcar, qtemp, n)
    circuit.x(qb[n])
    circuit.cx(qb[n], qtemp)
    circuit.x(qb[n])
    addition_nbit(circuit, qN, qb, qcar, n)
    subs_nbit(circuit, qa, qb, qcar, n)

# quantum controlled modular multiplication algorithm for n-qubit numbers
def cont_mod_mult_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, qX, qext, N, n):
    for i in range(n):
        for j in range(n):
            classical_mod = (2**(i+j))%N
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)
            mod_addition_nbit(circuit, qtempst, qb, qN, qNtemp, qcar, qtemp, n)
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)
    circuit.x(qext)
    cond_addition_nbit(circuit, qX, qb, qcar, qext, n)
    circuit.x(qext)

def cont_inv_mod_mult_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, qX, qext, N, n):
    "same as the controlled modular multiplication but running backwards"
    circuit.x(qext)
    cond_subs_nbit(circuit, qX, qb, qcar, qext, n)
    circuit.x(qext)
    for i in range(n):
        for j in range(n):
            classical_mod = (2**(i+j))%N
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)
            mod_subs_nbit(circuit, qtempst, qb, qN, qNtemp, qcar, qtemp, n)
            cond_number_state(circuit, qtempst, classical_mod, qext, qa[i], qX[j], n)

# quantum modular exponentiation algorithm for n-qubit numbers
def mod_exp_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, q1, qX, N, a, n):
    for k in range(len(qX)):
        clas_value = (a**(2**(k)))%N
        if k % 2 == 0:
            number_state(circuit, qa, clas_value, n)
            cont_mod_mult_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, q1, qX[k], N, n)
            number_state(circuit, qa, clas_value, n)
            clas_value = modinv(a**(2**(k)), N)
            number_state(circuit, qa, clas_value, n)
            cont_inv_mod_mult_nbit(circuit, qa, q1, qN, qNtemp, qcar, qtemp, qtempst, qb, qX[k], N, n)
            number_state(circuit, qa, clas_value, n)
        else:
            number_state(circuit, qa, clas_value, n)
            cont_mod_mult_nbit(circuit, qa, q1, qN, qNtemp, qcar, qtemp, qtempst, qb, qX[k], N, n)
            number_state(circuit, qa, clas_value, n)
            clas_value = modinv(a**(2**(k)), N)
            number_state(circuit, qa, clas_value, n)
            cont_inv_mod_mult_nbit(circuit, qa, qb, qN, qNtemp, qcar, qtemp, qtempst, q1, qX[k], N, n)
            number_state(circuit, qa, clas_value, n)

# n-qubit number input state
def number_state(circuit, q, x, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            circuit.x(q[i])

# n-qubit number input state, controlled by 2 control qubits
def cond_number_state(circuit, q, x, ext, control1, control2, n):
    # integer to binary
    x = "{0:b}".format(x)
    x = x.zfill(n)
    # creating the state
    for i in range(n):
        if int(x[n-1-i]) == 1:
            cond_toffoli(circuit, ext, control1, control2, q[i])

# efficient algorithm for computing the modular multiplicative inverse a^-1 mod m
def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception('modular inverse does not exist')
    else:
        return x % m

# we define the values, a and N must be coprimes
a = 1
x = 1
N = 1


# computing number of qubits n needed
n = len("{0:b}".format(a))
n2 = len("{0:b}".format(x))
n3 = len("{0:b}".format(N))

if n3 > n:
    n = n3

# classical register with n+1 bits.
c = q_program.create_classical_register("cr", n+1)
# quantum registers
qa = q_program.create_quantum_register("qa", n) # a qubits
qb = q_program.create_quantum_register("qb", n+1) # initial state |0>
qN = q_program.create_quantum_register("qN", n+1) # N qubits
qNtemp = q_program.create_quantum_register("qNtemp", n) # temporary N qubits
qtemp = q_program.create_quantum_register("qtemp", 1) # temporary qubit
qtempst = q_program.create_quantum_register("qtempst", n) # temporary register
q1 = q_program.create_quantum_register("q1", n+1) # initial state |1>
qX = q_program.create_quantum_register("qX", n2) # x register
# if n = 1, no need of carry register 
if n == 1:
    qcar = 0
    # quantum circuit involving the quantum registers and the classical register
    mod_exp_circuit = q_program.create_circuit("mod_exp", [qa, qb, qN, qNtemp, qtemp, qtempst, q1, qX],[c])
else:
    qcar = q_program.create_quantum_register("qcar", n-1) # carry qubits
    # quantum circuit involving the quantum registers and the classical register
    mod_exp_circuit = q_program.create_circuit("mod_exp", [qa, qb, qN, qcar, qNtemp, qtemp, qtempst, q1, qX],[c])                       

# create the initial state |1>. If N = 1, initial state is |0>
if N != 1:
    number_state(mod_exp_circuit, q1, 1, 1)
# create the state containing 'x'
number_state(mod_exp_circuit, qX, x, n2)
# create the state containing 'N'
number_state(mod_exp_circuit, qN, N, n+1)
# create a temporary state containing 'N'
number_state(mod_exp_circuit, qNtemp, N, n)

# modular exponentiation
mod_exp_nbit(mod_exp_circuit, qa, qb, qN, qNtemp, qcar, qtemp[0], qtempst, q1, qX, N, a, n)

# measurements to see the result, the result would be in one of those registers, q1 or qb
if n2 % 2 == 0:
    for i in range(n+1):
        mod_exp_circuit.measure(q1[i], c[i])
else:
    for i in range(n+1):
        mod_exp_circuit.measure(qb[i], c[i])

# compile and execute the quantum program in the backend
result = q_program.execute(["mod_exp"], backend=backend, shots=1024, max_credits=3, wait=10, timeout=9999999)

# show the results.
print(result)
print(result.get_data("mod_exp"))
counts = result.get_counts("mod_exp")
plot_histogram(counts)
print("Backend:", backend)
print("Highest probability outcome: {}".format(int(max(counts, key = lambda x: counts[x]).replace(" ", ""), 2)))