diff --git a/clifford.py b/clifford.py
index 0630b9f..2548963 100644
--- a/clifford.py
+++ b/clifford.py
@@ -11,97 +11,23 @@ Following the prescription of Anders (thesis pg. 26):
 > indicated by the column header and the opposite sign otherwise.
 """
 
-# TODO:
-# - check that we re-generate the table
-# - do conjugation
-# - do times table
-# - write tests
-
 from numpy import *
 
-
-def identify_pauli(m):
-    """ Given a signed Pauli matrix, name it. """
-    for sign in (+1, -1):
-        for pauli_label, pauli in zip("xyz", paulis):
-            if allclose(sign * pauli, m):
-                return sign, pauli_label
-
-
-def get_action(u):
-    """ What does this unitary operator do to the Paulis? """
-    return [identify_pauli(u * p * u.H) for p in paulis]
-
-
-def format_action(action):
-    """ Format an action as a string """
-    return "".join("{}{}".format("+" if s >= 0 else "-", p) for s, p in action)
-
-
-# Some two-qubit matrices
 i = matrix(eye(2, dtype=complex))
 px = matrix([[0, 1], [1, 0]], dtype=complex)
 py = matrix([[0, -1j], [1j, 0]], dtype=complex)
 pz = matrix([[1, 0], [0, -1]], dtype=complex)
 h = matrix([[1, 1], [1, -1]], dtype=complex) / sqrt(2)
 p = matrix([[1, 0], [0, 1j]], dtype=complex)
-paulis = (px, py, pz)
-# Some two-qubit matrices
-i = matrix(eye(2, dtype=complex))
-h = matrix([[1, 1], [1, -1]], dtype=complex) / sqrt(2)
-p = matrix([[1, 0], [0, 1j]], dtype=complex)
-
-# Basic single-qubit gates
-s_gates = (("i", i), ("p", p), ("pp", p * p), ("ppp", p * p * p))
-c_gates = [("i", i), ("h", h), ("hp", h * p), ("hpp", h * p * p),
-           ("hppp", h * p * p * p), ("hpph", h * p * p * h)]
 
-# Build the table of VOPs according to Anders (verbatim from thesis)
-table = (("a", "xyz", +1), ("b", "yxz", -1), ("c", "zyx", -1),
-        ("d", "xzy", -1), ("e", "yzx", +1), ("f", "zxy", +1))
+permutations = (i, h, p, h*p, h*p*h, h*p*h*p)
+signs = (i, px, py, pz)
+unitaries = [p*s for p in permutations for s in signs]
 
-# Build a big ol lookup table
-vop_names = []
-vop_actions = []
-vop_gates = [None] * 24
-vop_unitaries = [None] * 24
 
-for label, permutation, sign in table:
-    for column, operator in zip("ixyz", "i" + permutation):
-        effect = [((sign if (p == column or column == "i") else -sign), p)
-                  for p in permutation]
-        vop_names.append(column + label)  # think we can dump "operator"
-        vop_actions.append(format_action(effect))
-
-for s_name, s_gate in s_gates:
-    for c_name, c_gate in c_gates:
-        u = s_gate * c_gate
-        action = format_action(get_action(u))
-        index = vop_actions.index(action)
-        vop_gates[index] = s_name + c_name
-        vop_unitaries[index] = u
-
-# Add some more useful lookups
-vop_by_name = {n: {"name":n, "index": i, "action": a, "gates": g, "unitary": u}
-               for n, i, a, g, u in zip(vop_names, xrange(24), vop_actions, vop_gates, vop_unitaries)}
-vop_by_action = {a: {"name": n, "index": i, "action":a, "gates": g, "unitary": u}
-                 for n, i, a, g, u in zip(vop_names, xrange(24), vop_actions, vop_gates, vop_unitaries)}
-
-names, unitaries = [], []
-for c_name, c_gate in c_gates:
-    for s_name, s_gate in s_gates:
-        names.append(s_name+c_name)
-        unitaries.append(s_gate * c_gate)
-        print s_gate * c_gate.round(2)
-        print 
-
-i = matrix(eye(2, dtype=complex))
-px = matrix([[0, 1], [1, 0]], dtype=complex)
-py = matrix([[0, -1j], [1j, 0]], dtype=complex)
-pz = matrix([[1, 0], [0, -1]], dtype=complex)
-h = matrix([[1, 1], [1, -1]], dtype=complex) / sqrt(2)
-p = matrix([[1, 0], [0, 1j]], dtype=complex)
-
-#for m in i, px, py, pz:
-    #print any([allclose(x, m) for x in unitaries])
+# TODO:
+# - check that we re-generate the table
+# - do conjugation
+# - do times table
+# - write tests
 
diff --git a/tests/test_against_anders_thesis.py b/tests/test_against_anders_thesis.py
index 4bb8a99..ad5a103 100644
--- a/tests/test_against_anders_thesis.py
+++ b/tests/test_against_anders_thesis.py
@@ -36,13 +36,3 @@ anders = [
     ir2 * matrix([[1,-1],[-1j,-1j]], dtype=complex),
 ]
 
-
-def test_everything():
-    for i, (a, b) in enumerate(zip(lc.vop_actions, anders)):
-        a2 = lc.format_action(lc.get_action(b))
-        if i %4==0:
-            print
-        print "({} {})".format(a, a2),
-        #if not any([allclose(a, x) for x in anders]):
-            #print lc.vop_gates[i], "is not in {anders}"
-
diff --git a/tests/test_clifford.py b/tests/test_clifford.py
index 3030dbf..fb360f9 100644
--- a/tests/test_clifford.py
+++ b/tests/test_clifford.py
@@ -1,48 +1,56 @@
 import clifford as lc
 from numpy import *
-
-def test_identify_pauli():
-    assert lc.identify_pauli(lc.px) == (1, "x")
-    assert lc.identify_pauli(-lc.px) == (-1, "x")
-    assert lc.identify_pauli(-lc.pz) == (-1, "z")
-
-#def test_against_anders_table():
-    #assert allclose(lc.vop_unitaries[0], lc.i)
-    #assert allclose(lc.vop_unitaries[10], lc.h)
-
-    #yb = matrix([[1,0],[0,1j]])
-    #assert allclose(lc.vop_unitaries[5], yb)
-
-    #xb = matrix([[1,0],[0,-1j]])
-    #assert allclose(lc.vop_unitaries[6], xb)
-
-    #ye = matrix([[1,-1j],[-1,-1j]])/sqrt(2)
-    #print lc.vop_unitaries[17]
-    #print ye
-    #assert allclose(lc.vop_unitaries[17], ye)
-
-#def test_some_anders():
-    #u = matrix([[1,0],[0,1j]])
-    #print u
-    #print lc.format_action(lc.get_action(u))
-    #print lc.vop_by_name["xb"]
-
-    #u = matrix([[1,0],[0,0-1j]])
-    #print u
-    #print lc.format_action(lc.get_action(u))
-    #print lc.vop_by_name["yb"]
-
-
-#def _test_anders_problem():
-    #bi = lc.vop_by_name["bi"]
-    #print bi["name"]
-    #print bi["action"]
-    #print bi["unitary"]
-
-    #u = exp(-1j*pi/4)*matrix([[0,1],[1j,0]])
-    #print u
-    #print lc.format_action(lc.get_action(u))
-    #print lc.format_action(lc.identify_pauli(u*p*u.H) for p in lc.paulis)
-    #u = lc.vop_unitaries[4]
-    #print lc.format_action(lc.identify_pauli(u*p*u.H) for p in lc.paulis)
+from scipy.linalg import sqrtm
+
+sqy = sqrtm(1j*lc.py)
+msqy = sqrtm(-1j*lc.py)
+sqz = sqrtm(1j*lc.pz)
+msqz = sqrtm(-1j*lc.pz)
+sqx = sqrtm(1j*lc.px)
+msqx = sqrtm(-1j*lc.px)
+paulis = (lc.px, lc.py, lc.pz)
+
+def identify_pauli(m):
+    """ Given a signed Pauli matrix, name it. """
+    for sign in (+1, -1):
+        for pauli_label, pauli in zip("xyz", paulis):
+            if allclose(sign * pauli, m):
+                return sign, pauli_label
+
+def get_action(u):
+    """ What does this unitary operator do to the Paulis? """
+    return [identify_pauli(u * p * u.H) for p in paulis]
+
+def format_action(action):
+    return "".join("{}{}".format("+" if s>=0 else "-", p) for s, p in action)
+
+
+def test_we_have_all_useful_gates():
+    """ Check that all the interesting gates are included up to a global phase """
+    names = "i", "px", "py", "pz", "h", "p"
+    unitaries = lc.i, lc.px, lc.py, lc.pz, lc.h, lc.p
+    for name, unitary in zip(names, unitaries):
+        foundit = False
+        for i, clifford in enumerate(lc.unitaries):
+            if allclose(clifford, unitary): 
+                foundit = True
+                print "{}\t=\tlc.unitaries[{}]".format(name, i)
+        assert foundit
+
+    names = "sqrt(ix)", "sqrt(-ix)", "sqrt(iy)", "sqrt(-iy)", "sqrt(iz)", "sqrt(-iz)", 
+    unitaries = sqz, msqz, sqy, msqy, sqx, msqx
+    for name, unitary in zip(names, unitaries):
+        foundit = False
+        for phase in range(8):
+            for i, clifford in enumerate(lc.unitaries):
+                if allclose(exp(1j*phase*pi/4.)*clifford, unitary): 
+                    foundit = True
+                    print "{}\t=\texp({} . i . pi/4).lc.unitaries[{}]".format(name, phase, i)
+        assert foundit
+
+
+def test_we_have_24_matrices():
+    """ Check that we have 24 unique actions on the Bloch sphere """
+    actions = set(tuple(get_action(u)) for u in lc.unitaries)
+    assert len(set(actions)) == 24