diff --git a/lib_q_computer_math.py b/lib_q_computer_math.py
index bb7eca4..67a565e 100644
--- a/lib_q_computer_math.py
+++ b/lib_q_computer_math.py
@@ -12,6 +12,7 @@ from load_test import sizeof_fmt
 ListOrNdarray = Union[list, np.ndarray]
 
 REPR_TENSOR_OP = "⊗"
+REPR_GREEK_BETA = "β"
 REPR_GREEK_PSI = "ψ"
 REPR_GREEK_PHI = "φ"
 REPR_MATH_KET = "⟩"
@@ -104,6 +105,10 @@ class Matrix(object):
     def complex_conjugate(self):
         return Matrix(self._complex_conjugate())
 
+    @property
+    def human_m(self):
+        return humanize(self.m)
+
     def __repr__(self):
         return str(self.m)
 
@@ -168,6 +173,15 @@ class State(Vector):
         assert 0 <= phi <= 2 * np.pi
         return theta, phi
 
+    def rotate_x(self, theta):
+        return self | Rx(theta)
+
+    def rotate_y(self, theta):
+        return self | Ry(theta)
+
+    def rotate_z(self, theta):
+        return self | Rz(theta)
+
     def __repr__(self):
         if not self.name:
             if np.count_nonzero(self.m == 1):
@@ -275,10 +289,10 @@ def s(q):
 
 
 def humanize_num(fl, tolerance=1e-3):
-    if np.abs(fl) < tolerance:
-        return 0
     if np.abs(fl.imag) < tolerance:
         fl = fl.real
+    if np.abs(fl.real) < tolerance:
+        fl = 0 + 1j * fl.imag
     try:
         return sp.nsimplify(fl, [sp.pi], tolerance, full=True)
     except:
@@ -293,7 +307,10 @@ def pp(m, tolerance=1e-3):
 def humanize(m):
     rv = []
     for element in m:
-        rv.append(humanize(element))
+        if type(element) in [np.ndarray, list]:
+            rv.append(humanize(element))
+        else:
+            rv.append(humanize_num(element))
     return rv
 
 
@@ -459,6 +476,7 @@ _1 = State([[0],
             [1]],
            name='1')
 
+# |+> and |-> states
 _p = State([[1 / np.sqrt(2)],
             [1 / np.sqrt(2)]],
            name='+')
@@ -467,7 +485,32 @@ _m = State([[1 / np.sqrt(2)],
             [- 1 / np.sqrt(2)]],
            name='-')
 
-well_known_states = [_p, _m]
+# 4 Bell states
+b_00 = State([[1 / np.sqrt(2)],
+              [0],
+              [0],
+              [1 / np.sqrt(2)]],
+             name='B00')
+
+b_01 = State([[0],
+              [1 / np.sqrt(2)],
+              [1 / np.sqrt(2)],
+              [0]],
+             name='B01')
+
+b_10 = State([[1 / np.sqrt(2)],
+              [0],
+              [0],
+              [-1 / np.sqrt(2)]],
+             name='B10')
+
+b_11 = State([[0],
+              [1 / np.sqrt(2)],
+              [-1 / np.sqrt(2)],
+              [0]],
+             name='B11')
+
+well_known_states = [_p, _m, b_00, b_01, b_10, b_11]
 
 _ = I = UnitaryOperator([[1, 0],
                          [0, 1]],
@@ -485,6 +528,35 @@ Z = UnitaryOperator([[1, 0],
                      [0, -1]],
                     name="Z")
 
+
+# TODO: These are not the rotations that are specified commonly e.g. in
+#       https://www.quantum-inspire.com/kbase/rotation-operators/
+#       http://www.vcpc.univie.ac.at/~ian/hotlist/qc/talks/bloch-sphere-rotations.pdf
+#       and elsewhere.
+#       --------------
+#       Multiply the bellow ones by the following to get the commonly given
+#       X = [[1, -1j][-1j, 1]]
+#       Y = [[1, -1j][-1j, 1]]
+#       Z = [[1j, 0][1j, 0]]
+
+def Rx(theta):
+    return UnitaryOperator([[np.cos(theta / 2), np.sin(theta / 2)],
+                            [np.sin(theta / 2), np.cos(theta / 2)]],
+                           name="Rx")
+
+
+def Ry(theta):
+    return UnitaryOperator([[np.cos(theta / 2), -1j * np.sin(theta / 2)],
+                            [1j * np.sin(theta / 2), np.cos(theta / 2)]],
+                           name="Ry")
+
+
+def Rz(theta):
+    return UnitaryOperator([[1j * np.power(np.e, -1j * theta / 2), 0],
+                            [0, 1j * np.power(np.e, 1j * theta / 2)]],
+                           name="Rz")
+
+
 H = UnitaryOperator([[1 / np.sqrt(2), 1 / np.sqrt(2)],
                      [1 / np.sqrt(2), -1 / np.sqrt(2)], ],
                     name="H")
@@ -525,7 +597,7 @@ def assert_not_raises(exception, msg, callable, *args, **kwargs):
 def test_unitary_hermitian():
     # Unitary is UU+ = I; Hermitian is U = U+
     # Matrixes could be either, neither or both
-    # Quantum operators are described by unitary transformations
+    # Quantum operators are described *only* by unitary transformations
     # TODO: What are Hermitians?
     h_not_u = [
         [1, 0],
@@ -556,6 +628,12 @@ def test_unitary_hermitian():
     assert_raises(TypeError, "Not a Unitary matrix", UnitaryMatrix, not_u_not_h)
 
 
+def testRotPauli():
+    assert Rx(np.pi) == X
+    assert Ry(np.pi) == Y
+    assert Rz(np.pi) == Z
+
+
 def test():
     # Test properties of Hilbert vector space
     # The four postulates of Quantum Mechanics
@@ -600,6 +678,9 @@ def test():
     assert Y | _1 == State([[-1j],
                             [0]])
 
+    # Test Pauli rotation gates
+    testRotPauli()
+
     # III: Measurement | A quantum measurement is described by an orthonormal basis |e_j>
     #                    for state space. If the initial state of the system is |ψ>
     #                    then we get outcome j with probability pr(j) = |<e_j|ψ>|^2