diff --git a/abp/qi.py b/abp/qi.py
index 819a045..e86f544 100644
--- a/abp/qi.py
+++ b/abp/qi.py
@@ -74,6 +74,12 @@ class CircuitModel(object):
             if (i & control) and (i & target):
                 self.state[i, 0] *= -1
 
+    def act_cnot(self, control, target):
+        """ Act a CNOT """
+        self.act_hadamard(target)
+        self.act_cz(control, target)
+        self.act_hadamard(target)
+
     def act_hadamard(self, qubit):
         """ Act a hadamard somewhere """
         where = 1 << qubit
diff --git a/chp/chp.c b/chp/chp.c
index 21cd1e3..0d5c2af 100644
--- a/chp/chp.c
+++ b/chp/chp.c
@@ -631,6 +631,21 @@ static PyObject* act_cnot(PyObject* self, PyObject* args)
 }
 
 
+static char act_cz_docs[] = "act_cz(): Do a cz";
+static PyObject* act_cz(PyObject* self, PyObject* args)
+{
+    int a, b;
+    if (!PyArg_ParseTuple(args, "ii", &a, &b)) { return NULL; }
+    hadamard(q, b);
+    cnot(q, a, b);
+    hadamard(q, b);
+
+    Py_INCREF(Py_None);
+    return Py_None;
+}
+
+
+
 
 static char get_ket_docs[] = "get_ket(): Get the state vector";
 static PyObject* get_ket(PyObject* self, PyObject* noarg)
@@ -958,6 +973,7 @@ static PyMethodDef chp_funcs[] = {
     {"get_ket", (PyCFunction)get_ket, METH_NOARGS, get_ket_docs},
     {"act_hadamard", (PyCFunction)act_hadamard, METH_VARARGS, act_hadamard_docs},
     {"act_cnot", (PyCFunction)act_cnot, METH_VARARGS, act_cnot_docs},
+    {"act_cz", (PyCFunction)act_cz, METH_VARARGS, act_cz_docs},
     {"act_phase", (PyCFunction)act_phase, METH_VARARGS, act_phase_docs},
     {NULL}
 };
diff --git a/tests/test_circuit_model_against_chp.py b/tests/test_circuit_model_against_chp.py
index 854a3ba..45c8c94 100644
--- a/tests/test_circuit_model_against_chp.py
+++ b/tests/test_circuit_model_against_chp.py
@@ -2,28 +2,43 @@ import chp
 from abp import qi
 import numpy as np
 
-n = 10
+n = 5
 
 def get_chp_state():
     """ Convert CHP to CircuitModel """
     output = qi.CircuitModel(n)
     ket = chp.get_ket()
     nonzero = np.sqrt(len(ket))
+    output.state[0, 0]=0
     for key, phase in ket.items():
         output.state[key] = np.exp(1j*phase*np.pi/2)/nonzero
     return output
 
 
-def test1():
-    chp.init(5)
+def _bell_test():
+    chp.init(n)
     chp.act_hadamard(0)
     chp.act_cnot(0, 1)
 
-    yy = qi.CircuitModel(n)
-    yy.act_hadamard(0)
-    yy.act_hadamard(1)
-    yy.act_cz(0, 1)
-    yy.act_hadamard(1)
-    assert yy == get_chp_state()
+    psi = qi.CircuitModel(n)
+    psi.act_hadamard(0)
+    psi.act_cnot(0, 1)
+    assert psi == get_chp_state()
+
+def random_test():
+    chp.init(n)
+    psi = qi.CircuitModel(n)
+    for i in range(1000):
+        if np.random.rand()>.5:
+            a = np.random.randint(0, n-1)
+            chp.act_hadamard(a)
+            psi.act_hadamard(a)
+        else:
+            a, b = np.random.randint(0, n-1, 2)
+            if a!=b:
+                chp.act_cnot(a, b)
+                psi.act_cnot(a, b)
+        assert psi == get_chp_state()
+