diff --git a/abp/qi.py b/abp/qi.py
index fc5da2c..1276d58 100644
--- a/abp/qi.py
+++ b/abp/qi.py
@@ -80,7 +80,7 @@ class CircuitModel(object):
         output = np.zeros((self.d, 1), dtype=complex)
         for i, v in enumerate(self.state):
             q = i & where > 0
-            output[i] += v*u[q, q]
+            output[i] += v*u[q, q] # TODO this is probably wrong
             output[i ^ where] += v*u[not q, q]
         self.state = output
 
diff --git a/tests/test_against_circuit_model.py b/tests/test_against_circuit_model.py
index 0290796..9760abf 100644
--- a/tests/test_against_circuit_model.py
+++ b/tests/test_against_circuit_model.py
@@ -4,9 +4,9 @@ from abp import clifford
 import numpy as np
 import random
 
-def multi_qubit_test():
-    """ A multi qubit test """
-    n = 6
+def test_hadamard_only_multiqubit():
+    """ A multi qubit test with Hadamards only"""
+    n = 4
     g = GraphState(range(n))
     c = CircuitModel(n)
 
@@ -27,3 +27,28 @@ def multi_qubit_test():
     assert np.allclose(s1, s2)
 
 
+def test_all_multiqubit():
+    """ A multi qubit test with arbitrary local rotations """
+    n = 4
+    g = GraphState(range(n))
+    c = CircuitModel(n)
+
+    for i in range(10):
+        i = np.random.randint(0, n-1)
+        j = np.random.randint(0, 24)
+        print i, j
+        g.act_local_rotation(i, j)
+        c.act_local_rotation(i, clifford.unitaries[j])
+
+    assert np.allclose(g.to_state_vector().state, c.state)
+
+    #for i in range(100):
+        #a, b = np.random.randint(0, n-1, 2)
+        #if a != b:
+            #g.act_cz(a, b)
+            #c.act_cz(a, b)
+
+    #s1 = clifford.normalize_global_phase(g.to_state_vector().state)
+    #s2 = clifford.normalize_global_phase(c.state)
+    #assert np.allclose(s1, s2)
+