diff --git a/cz_table.py b/cz_table.py
index 74a0cc1..0135e51 100644
--- a/cz_table.py
+++ b/cz_table.py
@@ -1,32 +1,42 @@
 import qi
 import numpy as np
 import tables
-import tqdm
+from tqdm import tqdm
+import itertools as it
 
-# TODO: ensure that Constraint 1 is met. i.e. 
-# if C1 is in Z, choose C1' such that it is in Z
+def get_krontable():
+    table = np.zeros((24,24,4,4), dtype=complex)
+    for i, j in it.product(range(24), range(24)):
+        u1 = tables.unitaries[i]
+        u2 = tables.unitaries[j]
+        table[i, j, :, :] = np.kron(u1, u2)
+    return table
 
-table1 = []
-table2 = []
+def find(bond, c1, c2, z, krontable):
+    # Figure out the target state
+    state = qi.bond if bond else qi.nobond
+    target = qi.cz * krontable[c1, c2] * state
 
-bond = qi.cz * np.kron(qi.plus, qi.plus)
-no_bond = np.kron(qi.plus, qi.plus)
+    # Choose the sets to search over
+    s1 = z if c1 in z else xrange(24)
+    s2 = z if c2 in z else xrange(24)
 
-def find(thing, table):
-    for index, trial in enumerate(table):
-        for qq in range(16):
-            if np.allclose(thing, np.exp(2j*np.pi*qq/16.) * trial):
-                yield index
+    for bond, c1p, c2p in it.product([0,1], s1, s2):
+        state = qi.bond if bond else qi.nobond
+        trial = krontable[c1p, c2p] * state
+        for phase in range(8):
+            if np.allclose(target, np.exp(1j * phase * np.pi / 4.) * trial):
+                return bond, c1p, c2p
 
-for state in bond, no_bond:
-    for a in tables.unitaries:
-        for b in tables.unitaries:
-            state = np.kron(a, b) * state
-            table1.append(state)
-            table2.append(qi.cz*state)
+    raise IndexError
 
-for index, thing in enumerate(table2):
-    print "{} -> {}".format(index, list(find(thing, table1)))
 
+z = [tables.find(u, tables.unitaries) for u in qi.id, qi.px, qi.pz, qi.ph, qi.ph.H]
+krontable = get_krontable()
 
+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(bond, c1, c2, z, krontable)
+
+np.save("cz_table.npy", cz_table)
 
diff --git a/qi.py b/qi.py
index 35ce5a1..fb35bf6 100644
--- a/qi.py
+++ b/qi.py
@@ -30,3 +30,6 @@ paulis = (px, py, pz)
 common_us = id, px, py, pz, ha, ph, sqz, msqz, sqy, msqy, sqx, msqx
 names = "identity", "px", "py", "pz", "hadamard", "phase", "sqz", "msqz", "sqy", "msqy", "sqx", "msqx"
 by_name = dict(zip(names, common_us))
+
+bond = cz * np.kron(plus, plus)
+nobond = np.kron(plus, plus)
diff --git a/tables.py b/tables.py
index 826a32c..a956f25 100644
--- a/tables.py
+++ b/tables.py
@@ -13,18 +13,18 @@ import cPickle
 import qi
 
 # TODO: make this more efficient / shorter
-def find_up_to_phase(u):
+def find(needle, haystack):
     """ Find the index of a given u within a list of unitaries, up to a global phase  """
-    for i, t in enumerate(unitaries):
+    for i, t in enumerate(haystack):
         for phase in range(8):
-            if np.allclose(t, np.exp(1j * phase * np.pi / 4.) * u):
-                return i, phase
+            if np.allclose(t, np.exp(1j * phase * np.pi / 4.) * needle):
+                return i
     raise IndexError
 
 
 def compose_u(decomposition):
     """ Get the unitary representation of a particular decomposition """
-    us = (elements[c] for c in decomposition)
+    us = ({"x": qi.sqx, "z": qi.msqz}[c] for c in decomposition)
     return np.matrix(reduce(np.dot, us), dtype=complex)
 
 
@@ -38,14 +38,15 @@ def construct_tables(filename="tables.cache"):
     if os.path.exists(filename):
         return cPickle.load(open(filename, "r"))
 
-    by_name = {name: find_up_to_phase(u)[0] for name, u in qi.by_name.items()}
+    by_name = {name: find(u, unitaries) for name, u in qi.by_name.items()}
     get_name = {v:k for k, v in by_name.items()}
-    conjugation_table = [find_up_to_phase(u.H)[0]
+    conjugation_table = [find(u.H, unitaries)
                          for i, u in enumerate(unitaries)]
-    times_table = [[find_up_to_phase(u * v)[0] for v in unitaries]
+    times_table = [[find(u * v, unitaries) for v in unitaries]
                    for u in tqdm(unitaries)]
     cz_table = None
     output = by_name, get_name, conjugation_table, times_table, cz_table
+
     with open(filename, "w") as f:
         cPickle.dump(output, f)
     return output
@@ -53,7 +54,5 @@ def construct_tables(filename="tables.cache"):
 decompositions = ("xxxx", "xx", "zzxx", "zz", "zxx", "z", "zzz", "xxz",
      "xzx", "xzxxx", "xzzzx", "xxxzx", "xzz", "zzx", "xxx", "x",
      "zzzx", "xxzx", "zx", "zxxx", "xxxz", "xzzz", "xz", "xzxx")
-elements = {"x": qi.sqx, "z": qi.msqz}
 unitaries = [compose_u(d) for d in decompositions]
 by_name, get_name, conjugation_table, times_table, cz_table = construct_tables()
-
diff --git a/tests/test_clifford.py b/tests/test_clifford.py
index 5df3ec5..d440eb4 100644
--- a/tests/test_clifford.py
+++ b/tests/test_clifford.py
@@ -14,11 +14,11 @@ def identify_pauli(m):
                 return sign, pauli_label
 
 
-def _test_find_up_to_phase():
+def _test_find():
     """ Test that slightly suspicious function """
-    assert lc.find_up_to_phase(id) == (0, 0)
-    assert lc.find_up_to_phase(px) == (1, 0)
-    assert lc.find_up_to_phase(exp(1j*pi/4.)*ha) == (4, 7)
+    assert lc.find(id, lc.unitaries) == 0
+    assert lc.find(px, lc.unitaries) == 1
+    assert lc.find(exp(1j*pi/4.)*ha, lc.unitaries) == 4
 
 def get_action(u):
     """ What does this unitary operator do to the Paulis? """
@@ -38,14 +38,14 @@ def test_we_have_24_matrices():
 def test_we_have_all_useful_gates():
     """ Check that all the interesting gates are included up to a global phase """
     for name, u in qi.by_name.items():
-        lc.find_up_to_phase(u)
+        lc.find(u, lc.unitaries)
 
 
 def _test_group():
     """ Test we are really in a group """
     matches = set()
     for a, b in tqdm(it.combinations(lc.unitaries, 2), "Testing this is a group"):
-        i, phase = lc.find_up_to_phase(a*b)
+        i, phase = lc.find(a*b, lc.unitaries)
         matches.add(i)
     assert len(matches)==24