%matplotlib inline
import numpy as np
import numpy.random as rnd
import matplotlib.pyplot as plt, mpld3
import operator
plt.rcParams['figure.figsize'] = 12, 8
mpld3.enable_notebook()
def kron(*Xs):
Y = np.eye(1)
for X in Xs:
Y = np.kron(Y, X)
return Y
def prod(iterable):
return reduce(operator.mul, iterable, 1)
def randpsi(d):
psi = np.random.normal(size=d) + 1j * np.random.normal(size=d)
return psi / np.linalg.norm(psi)
class QC(object):
def __init__(self, num_qubits):
"""starts out in state |0....0>"""
self.num_qubits = num_qubits
self.psi = np.zeros(2**num_qubits)
self.psi[0] = 1
def _tensor(self, indices, Xs):
"""compute tensor product of operators Xs at indices"""
assert len(indices) == len(set(indices)), 'no repetitions'
assert all(i < self.num_qubits for i in indices), 'index out of bounds'
X = np.eye(1)
for i in range(self.num_qubits):
if i in indices:
X = np.kron(X, Xs[indices.index(i)])
else:
X = np.kron(X, np.eye(2))
return X
def rdm(self, *indices):
"""reduced density matrix"""
indices = list(indices)
compl = [i for i in range(self.num_qubits) if i not in indices]
A = self.psi.reshape([2] * self.num_qubits).transpose(indices + compl).reshape([2**len(indices), 2**len(compl)])
return np.dot(A, A.T.conj())
def H(self, idx):
"""Hadamard gate"""
H = [[1, 1], [1, -1]] / np.sqrt(2)
U = self._tensor([idx], [H])
self.psi = np.dot(U, self.psi)
def NOT(self, idx):
"""NOT gate"""
NOT = [[0, 1], [1, 0]]
U = self._tensor([idx], [NOT])
self.psi = np.dot(U, self.psi)
def Z(self, idx):
"""Pauli Z-operator"""
Z = [[1, 0], [0, -1]]
U = self._tensor([idx], [Z])
self.psi = np.dot(U, self.psi)
def C(self, U, idx, control_idx):
"""controlled unitary"""
proj0 = [[1, 0], [0, 0]]
proj1 = [[0, 0], [0, 1]]
U = self._tensor([idx, control_idx], [np.eye(2), proj0]) + self._tensor([idx, control_idx], [U, proj1])
self.psi = np.dot(U, self.psi)
def CNOT(self, idx, control_idx):
"""controller NOT"""
NOT = [[0, 1], [1, 0]]
self.C(NOT, idx, control_idx)
def measure(self, idx):
"""measure a qubit"""
proj0 = [[1, 0], [0, 0]]
psi0 = np.dot(self._tensor([idx], [proj0]), self.psi)
prob0 = np.linalg.norm(psi0)**2
if np.random.random() < prob0:
self.psi = psi0 / np.sqrt(prob0)
return 0
else:
self.psi = (self.psi - psi0) / np.sqrt(1 - prob0)
return 1
KET_EPR = np.array([1, 0, 0, 1]) / np.sqrt(2)
def teleport(qc, ALICE_1, ALICE_2, BOB):
# teleportation protocol
assert np.allclose(qc.rdm(ALICE_2, BOB), np.outer(KET_EPR, KET_EPR))
qc.CNOT(ALICE_2, control_idx=ALICE_1)
qc.H(ALICE_1)
bits = [qc.measure(ALICE_1), qc.measure(ALICE_2)] # ...these are sent to bob...
if bits[1]:
qc.NOT(BOB)
if bits[0]:
qc.Z(BOB)
return qc.rdm(BOB)
# simple test
qc = QC(num_qubits=3)
# prepare initial state of Alice' first qubit
qc.H(0)
expected = qc.rdm(0)
# prepare maximally entangled state between Alice' other qubit and Bob's qubit
qc.H(1)
qc.CNOT(2, control_idx=1)
got = teleport(qc, 0, 1, 2)
assert np.allclose(expected, got)
Teleport random states:
%%time
qc = QC(num_qubits=3)
for _ in range(100):
# initialize QC with random state on first qubit and EPR pair between the other two
qc.psi = np.kron(randpsi(2), KET_EPR)
expected = qc.rdm(0)
# teleport
teleport(qc, 0, 1, 2)
# verify
got = qc.rdm(2)
assert np.allclose(expected, got)
Note that the "entanglement" (in fact, the entire reduced state) is preserved by the teleportation procedure.
%%time
qc = QC(num_qubits=4)
for _ in range(100):
# initialize QC with random state on first two qubits and EPR pair on the other two
qc.psi = np.kron(randpsi(4), KET_EPR)
expected = qc.rdm(0, 1)
# teleport
teleport(qc, 1, 2, 3)
# verify
got = qc.rdm(0, 3)
assert np.allclose(expected, got)
Thus given an EPR pair between Alice and Bob and another one between Bob and Charly we can get an EPR pair between Alice and Charly!
qc = QC(num_qubits=4)
ALICE, BOB_1, BOB_2, CHARLY = range(4)
# initialize with two EPR pairs between Alice and Bob1 & between Bob2 and Charly, respectively
qc.psi = np.kron(KET_EPR, KET_EPR)
# teleport Bob1 via Bob2--Charly
teleport(qc, BOB_1, BOB_2, CHARLY)
# Alice--Charly are now entangled (and can do stuff like QKD, etc.)
assert np.allclose(qc.rdm(0, 3), np.outer(KET_EPR, KET_EPR))
Truth tables of all constant and balanced functions of $n$ bits:
from itertools import combinations
def constant_fns(n):
yield [0] * (2**n)
yield [1] * (2**n)
def balanced_fns(n):
inputs = range(2**n)
for s in combinations(inputs, 2**(n-1)):
yield [int(i in s) for i in inputs]
def rand_balanced_fn(n):
inputs = range(2**n)
s = np.random.choice(inputs, 2**(n-1), replace=False)
return [int(i in s) for i in inputs]
print list(constant_fns(2))
print list(balanced_fns(2))
print rand_balanced_fn(2)
Unitary oracle for a boolean function $f$:
def U_oracle(n, f):
U = np.zeros(shape=(2**(n+1), 2**(n+1)))
for x in range(2**n):
for y in range(2):
z = (y + f[x]) % 2
U[x*2 + z, x*2 + y] = 1
return U
U_oracle(1, [1, 0])
Deutsch--Josza algorithm:
def deutsch_josza(n, f):
# construct unitary oracle
U = U_oracle(n, f)
# initialize in state |0>^{\otimes n} \otimes |1>
qc = QC(n+1)
qc.NOT(n)
# hadamards on all qubits
for i in range(n+1):
qc.H(i)
# oracle
qc.psi = np.dot(U, qc.psi)
# hadamards on all qubits (except for last one, which does not matter anymore)
for i in range(n):
qc.H(i)
# measure...
bits = [qc.measure(i) for i in range(n)]
is_const = all(b == 0 for b in bits)
return is_const
assert deutsch_josza(2, [0, 0, 0, 0]) == True
assert deutsch_josza(2, [0, 1, 0, 1]) == False
Systematic Test:
%%time
n = 3
for f in constant_fns(n):
assert deutsch_josza(n, f) == True
for f in balanced_fns(n):
assert deutsch_josza(n, f) == False
n = 5
for _ in range(100):
f = rand_balanced_fn(n)
assert deutsch_josza(n, f) == False