diff --git a/make_tables.py b/make_tables.py
index 07558cd..2c3a06f 100644
--- a/make_tables.py
+++ b/make_tables.py
@@ -25,20 +25,19 @@ def find_clifford(needle, haystack):
     raise IndexError
 
 
-def find_cz(bond, c1, c2, z, krontable):
+def find_cz(bond, c1, c2, commuters, state_table):
     """ Find the output of a CZ operation """
     # Figure out the target state
     state = qi.bond if bond else qi.nobond
-    target = qi.cz.dot(krontable[c1, c2].dot(state))
+    target = qi.cz.dot(state_table[bond, c1, c2])
 
     # Choose the sets to search over
-    s1 = z if c1 in z else xrange(24)
-    s2 = z if c2 in z else xrange(24)
+    s1 = commuters if c1 in commuters else xrange(24)
+    s2 = commuters if c2 in commuters else xrange(24)
 
     # Find a match
     for bond, c1p, c2p in it.product([0, 1], s1, s2):
-        state = qi.bond if bond else qi.nobond
-        trial = krontable[c1p, c2p].dot(state)
+        trial = state_table[bond, c1p, c2p]
         for phase in range(8):
             if np.allclose(target, np.exp(1j * phase * np.pi / 4.) * trial):
                 return bond, c1p, c2p
@@ -74,23 +73,26 @@ def get_times_table(unitaries):
                    for u in tqdm(unitaries)])
 
 
-def get_krontable():
-    """ Cache a table of Kronecker products to speed up a little bit """
-    krontable = np.zeros((24, 24, 4, 4), dtype=complex)
-    for i, j in it.product(range(24), range(24)):
-        krontable[i, j,:,:] = np.kron(unitaries[i], unitaries[j])
-    return krontable
+def get_state_table():
+    """ Cache a table of state to speed up a little bit """
+    state_table = np.zeros((2, 24, 24, 4), dtype=complex)
+    for bond, i, j in it.product([0, 1], range(24), range(24)):
+        state = qi.bond if bond else qi.nobond
+        kp = np.kron(unitaries[i], unitaries[j])
+        state_table[bond, i, j,:] = np.dot(kp, state).T
+    return state_table
+
 
 
 def get_cz_table(unitaries):
     """ Compute the lookup table for the CZ (A&B eq. 9) """
-    z = (qi.id, qi.px, qi.pz, qi.ph, hermitian_conjugate(qi.ph))
-    z = [find_clifford(u, unitaries) for u in z]
-    krontable = get_krontable()
+    commuters = (qi.id, qi.px, qi.pz, qi.ph, hermitian_conjugate(qi.ph))
+    commuters = [find_clifford(u, unitaries) for u in commuters]
+    state_table = get_state_table()
 
     cz_table = np.zeros((2, 24, 24, 3))
     for bond, c1, c2 in tqdm(list(it.product([0, 1], range(24), range(24)))):
-        cz_table[bond, c1, c2] = find_cz(bond, c1, c2, z, krontable)
+        cz_table[bond, c1, c2] = find_cz(bond, c1, c2, commuters, state_table)
     return cz_table
 
 if __name__ == '__main__':
diff --git a/qi.py b/qi.py
index c57c963..5ce96d7 100644
--- a/qi.py
+++ b/qi.py
@@ -9,28 +9,28 @@ import numpy as np
 from scipy.linalg import sqrtm
 
 # Operators
-id = np.matrix(np.eye(2, dtype=complex))
-px = np.matrix([[0, 1], [1, 0]], dtype=complex)
-py = np.matrix([[0, -1j], [1j, 0]], dtype=complex)
-pz = np.matrix([[1, 0], [0, -1]], dtype=complex)
-ha = np.matrix([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
-ph = np.matrix([[1, 0], [0, 1j]], dtype=complex)
+id = np.array(np.eye(2, dtype=complex))
+px = np.array([[0, 1], [1, 0]], dtype=complex)
+py = np.array([[0, -1j], [1j, 0]], dtype=complex)
+pz = np.array([[1, 0], [0, -1]], dtype=complex)
+ha = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
+ph = np.array([[1, 0], [0, 1j]], dtype=complex)
 
 sqy = sqrtm(1j * py)
-msqy = np.matrix(sqrtm(-1j * py))
-sqz = np.matrix(sqrtm(1j * pz))
-msqz = np.matrix(sqrtm(-1j * pz))
-sqx = np.matrix(sqrtm(1j * px))
-msqx = np.matrix(sqrtm(-1j * px))
+msqy = np.array(sqrtm(-1j * py))
+sqz = np.array(sqrtm(1j * pz))
+msqz = np.array(sqrtm(-1j * pz))
+sqx = np.array(sqrtm(1j * px))
+msqx = np.array(sqrtm(-1j * px))
 paulis = (px, py, pz)
 
 # CZ gate
-cz = np.matrix(np.eye(4), dtype=complex)
+cz = np.array(np.eye(4), dtype=complex)
 cz[3,3]=-1
 
 # States
-plus = np.matrix([[1],[1]], dtype=complex) / np.sqrt(2)
-bond = cz * np.kron(plus, plus)
+plus = np.array([[1],[1]], dtype=complex) / np.sqrt(2)
+bond = cz.dot(np.kron(plus, plus))
 nobond = np.kron(plus, plus)
 
 # Labelling stuff