diff --git a/docs/intro_tutorial.rst b/docs/intro_tutorial.rst
index 002dd2d2d..ed0f02710 100644
--- a/docs/intro_tutorial.rst
+++ b/docs/intro_tutorial.rst
@@ -272,22 +272,15 @@ fidelity function between quantum states that happen to be identical.
 
     >>> from toqito.states import bell
     >>> from toqito.state_metrics import fidelity
+    >>> import numpy as np
     >>>
     >>> # Define two identical density operators.
     >>> rho = bell(0)*bell(0).conj().T
     >>> sigma = bell(0)*bell(0).conj().T
     >>> 
     >>> # Calculate the fidelity between `rho` and `sigma`
-    >>> '%.2f' % fidelity(rho, sigma)
-    '1.00'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(fidelity(rho, sigma), decimals=2)
+    1.0
 
 There are a number of other metrics one can compute on two density matrices
 including the trace norm, trace distance. These and others are also available
diff --git a/docs/tutorials.extended_nonlocal_games.rst b/docs/tutorials.extended_nonlocal_games.rst
index df73834eb..580d93fcd 100644
--- a/docs/tutorials.extended_nonlocal_games.rst
+++ b/docs/tutorials.extended_nonlocal_games.rst
@@ -316,21 +316,14 @@ This can be verified in :code:`toqito` as follows.
 
     >>> # Calculate the unentangled value of the BB84 extended nonlocal game.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define an ExtendedNonlocalGame object based on the BB84 game.
     >>> bb84 = ExtendedNonlocalGame(bb84_prob_mat, bb84_pred_mat)
     >>> 
     >>> # The unentangled value is cos(pi/8)**2 \approx 0.85356
-    >>> '%.2f' % bb84.unentangled_value() 
-    '0.85'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(bb84.unentangled_value(), decimals=2)
+    0.85
 
 The BB84 game also exhibits strong parallel repetition. We can specify how many
 parallel repetitions for :code:`toqito` to run. The example below provides an
@@ -340,21 +333,14 @@ example of two parallel repetitions for the BB84 game.
 
     >>> # The unentangled value of BB84 under parallel repetition.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define the bb84 game for two parallel repetitions.
     >>> bb84_2_reps = ExtendedNonlocalGame(bb84_prob_mat, bb84_pred_mat, 2)
     >>> 
     >>> # The unentangled value for two parallel repetitions is cos(pi/8)**4 \approx 0.72855
-    >>> '%.2f' % bb84_2_reps.unentangled_value() 
-    '0.73'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(bb84_2_reps.unentangled_value(), decimals=2)
+    0.73
 
 It was shown in :cite:`Johnston_2016_Extended` that the BB84 game possesses the property of strong
 parallel repetition. That is,
@@ -374,21 +360,14 @@ using :code:`toqito` as well.
 
     >>> # Calculate lower bounds on the standard quantum value of the BB84 extended nonlocal game.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define an ExtendedNonlocalGame object based on the BB84 game.
     >>> bb84_lb = ExtendedNonlocalGame(bb84_prob_mat, bb84_pred_mat)
     >>> 
     >>> # The standard quantum value is cos(pi/8)**2 \approx 0.85356
-    >>> '%.2f' % bb84_lb.quantum_value_lower_bound()
-    '0.85'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(bb84_lb.quantum_value_lower_bound(), decimals=2)
+    0.85
 
 From :cite:`Johnston_2016_Extended`, it is known that :math:`\omega(G_{BB84}) =
 \omega^*(G_{BB84})`, however, if we did not know this beforehand, we could
@@ -409,20 +388,14 @@ Using :code:`toqito`, we can see that :math:`\omega_{ns}(G) = \cos^2(\pi/8)`.
 
     >>> # Calculate the non-signaling value of the BB84 extended nonlocal game.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define an ExtendedNonlocalGame object based on the BB84 game.
     >>> bb84 = ExtendedNonlocalGame(bb84_prob_mat, bb84_pred_mat)
     >>> 
     >>> # The non-signaling value is cos(pi/8)**2 \approx 0.85356
-    >>> '%.2f' % bb84.nonsignaling_value() 
-    '0.85'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(bb84.nonsignaling_value(), decimals=2)
+    0.85
 
 So we have the relationship that
 
@@ -437,20 +410,14 @@ can observe this by the following snippet.
 
     >>> # The non-signaling value of BB84 under parallel repetition.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define the bb84 game for two parallel repetitions.
     >>> bb84_2_reps = ExtendedNonlocalGame(bb84_prob_mat, bb84_pred_mat, 2)
     >>> 
     >>> # The non-signaling value for two parallel repetitions is cos(pi/8)**4 \approx 0.73825
-    >>> '%.2f' % bb84_2_reps.nonsignaling_value() 
-    '0.74'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(bb84_2_reps.nonsignaling_value(), decimals=2)
+    0.74
 
 Note that :math:`0.73825 \geq \cos(\pi/8)^4 \approx 0.72855` and therefore we
 have that
@@ -554,20 +521,14 @@ the unentangled value of :math:`G_{CHSH}`.
 
     >>> # Calculate the unentangled value of the CHSH extended nonlocal game
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define an ExtendedNonlocalGame object based on the CHSH game.
     >>> chsh = ExtendedNonlocalGame(chsh_prob_mat, chsh_pred_mat)
     >>> 
     >>> # The unentangled value is 3/4 = 0.75
-    >>> '%.2f' % chsh.unentangled_value() 
-    '0.75'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(chsh.unentangled_value(), decimals=2)
+    0.75
 
 We can also run multiple repetitions of :math:`G_{CHSH}`.
 
@@ -575,21 +536,14 @@ We can also run multiple repetitions of :math:`G_{CHSH}`.
 
     >>> # The unentangled value of CHSH under parallel repetition.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define the CHSH game for two parallel repetitions.
     >>> chsh_2_reps = ExtendedNonlocalGame(chsh_prob_mat, chsh_pred_mat, 2)
     >>> 
     >>> # The unentangled value for two parallel repetitions is (3/4)**2 \approx 0.5625
-    >>> '%.2f' % chsh_2_reps.unentangled_value() 
-    '0.56'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(chsh_2_reps.unentangled_value(), decimals=2)
+    0.56
 
 Note that strong parallel repetition holds as
 
@@ -606,21 +560,14 @@ non-signaling value.
 
     >>> # Calculate the non-signaling value of the CHSH extended nonlocal game.
     >>> from toqito.nonlocal_games.extended_nonlocal_game import ExtendedNonlocalGame
+    >>> import numpy as np
     >>> 
     >>> # Define an ExtendedNonlocalGame object based on the CHSH game.
     >>> chsh = ExtendedNonlocalGame(chsh_prob_mat, chsh_pred_mat)
     >>> 
     >>> # The non-signaling value is 3/4 = 0.75
-    >>> '%.2f' % chsh.nonsignaling_value() 
-    '0.75'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
+    >>> np.around(chsh.nonsignaling_value(), decimals=2)
+    0.75
 
 As we know that :math:`\omega(G_{CHSH}) = \omega_{ns}(G_{CHSH}) = 3/4` and that
 
@@ -739,18 +686,11 @@ Now that we have encoded :math:`G_{MUB}`, we can calculate the unentangled value
 
 .. code-block:: python
 
+    >>> import numpy as np
     >>> g_mub = ExtendedNonlocalGame(prob_mat, pred_mat)
     >>> unent_val = g_mub.unentangled_value()
-    >>> '%.2f' % unent_val
-    '0.65'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(unent_val, decimals=2)
+    0.65
 
 That is, we have that 
 
@@ -763,18 +703,11 @@ obtain.
 
 .. code-block:: python
 
+    >>> import numpy as np
     >>> g_mub = ExtendedNonlocalGame(prob_mat, pred_mat)
     >>> q_val = g_mub.quantum_value_lower_bound()
-    >>> '%.2f' % q_val
-    '0.66'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(q_val, decimals=2)
+    0.66
 
 Note that as we are calculating a lower bound, it is possible that a value this
 high will not be obtained, or in other words, the algorithm can get stuck in a
diff --git a/docs/tutorials.nonlocal_games.rst b/docs/tutorials.nonlocal_games.rst
index 6431e5fe0..d8e15ec58 100644
--- a/docs/tutorials.nonlocal_games.rst
+++ b/docs/tutorials.nonlocal_games.rst
@@ -343,17 +343,11 @@ use :code:`toqito` to determine the lower bound on the quantum value.
 
 .. code-block:: python
 
+    >>> import numpy as np
     >>> from toqito.nonlocal_games.nonlocal_game import NonlocalGame
     >>> chsh = NonlocalGame(prob_mat, pred_mat)
-    >>> '%.2f' % chsh.quantum_value_lower_bound()
-    '0.85'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(chsh.quantum_value_lower_bound(), decimals=2)
+    0.85
 
 In this case, we can see that the quantum value of the CHSH game is in fact
 attained as :math:`\cos^2(\pi/8) \approx 0.85355`.
@@ -419,17 +413,10 @@ value and calculate lower bounds on the quantum value of the FFL game.
     ...                     pred_mat[a_alice, b_bob, x_alice, y_bob] = 1
     >>> # Define the FFL game object.
     >>> ffl = NonlocalGame(prob_mat, pred_mat)
-    >>> '%.2f' % ffl.classical_value()
-    '0.67'
-    >>> '%.2f' % ffl.quantum_value_lower_bound()
-    '0.22'
-
-        .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(ffl.classical_value(), decimals=2)
+    0.67
+    >>> np.around(ffl.quantum_value_lower_bound(), decimals=2)
+    0.22
 
 In this case, we obtained the correct quantum value of :math:`2/3`, however,
 the lower bound technique is not guaranteed to converge to the true quantum
diff --git a/docs/tutorials.state_distinguishability.rst b/docs/tutorials.state_distinguishability.rst
index 6b62c775b..4bf611fa9 100644
--- a/docs/tutorials.state_distinguishability.rst
+++ b/docs/tutorials.state_distinguishability.rst
@@ -114,6 +114,7 @@ Using :code:`toqito`, we can calculate this probability directly as follows:
 
 .. code-block:: python
 
+    >>> import numpy as np
     >>> from toqito.states import basis
     >>> from toqito.state_opt import state_distinguishability
     >>> 
@@ -132,15 +133,8 @@ Using :code:`toqito`, we can calculate this probability directly as follows:
     >>>
     >>> # Calculate the probability with which Bob can 
     >>> # distinguish the state he is provided.
-    >>> '%.2f' % state_distinguishability(states, probs)[0]
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(state_distinguishability(states, probs)[0], decimals=2)
+    1.0
 
 Specifying similar state distinguishability problems can be done so using this
 general pattern.
@@ -235,15 +229,8 @@ via PPT measurements in the following manner.
     >>>
     >>> states = [rho_1, rho_2, rho_3, rho_4]
     >>> probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
-    >>> '%.2f' % ppt_distinguishability(vectors=states, probs=probs, dimensions=[2, 2, 2, 2], subsystems=[0, 2])[0]
-    '0.87'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(ppt_distinguishability(vectors=states, probs=probs, dimensions=[2, 2, 2, 2], subsystems=[0, 2])[0], decimals=2)
+    0.87
 
 Probability of distinguishing a state via separable measurements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/docs/tutorials.xor_quantum_value.rst b/docs/tutorials.xor_quantum_value.rst
index 53e4017ca..3959e91c8 100644
--- a/docs/tutorials.xor_quantum_value.rst
+++ b/docs/tutorials.xor_quantum_value.rst
@@ -288,15 +288,9 @@ follows:
 
 .. code-block:: python
 
-    >>> '%.2f' % chsh.quantum_value()
-    '0.85'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> import numpy as np
+    >>> np.around(chsh.quantum_value(), decimals=2)
+    0.85
 
 For reference, the complete code to calculate both the classical and quantum
 values of the CHSH game is provided below.
@@ -312,16 +306,8 @@ values of the CHSH game is provided below.
     >>> chsh = XORGame(prob_mat, pred_mat)
     >>> chsh.classical_value()
     0.75
-    >>> '%.2f' % chsh.quantum_value()
-    '0.85'
-
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(chsh.quantum_value(), decimals=2)
+    0.85
 
 The odd cycle game
 ------------------
@@ -363,17 +349,10 @@ the classical and quantum values of this game.
     >>>
     >>> # Compute the classical and quantum values.
     >>> odd_cycle = XORGame(prob_mat, pred_mat)
-    >>> '%.2f' % odd_cycle.classical_value()
-    '0.90'
-    >>> '%.2f' % odd_cycle.quantum_value()
-    '0.98'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(odd_cycle.classical_value(), decimals=2)
+    0.9
+    >>> np.around(odd_cycle.quantum_value(), decimals=2)
+    0.98
 
 Note that the odd cycle game is another example of an XOR game where the
 players are able to win with a strictly higher probability if they adopt a
diff --git a/toqito/channel_metrics/channel_fidelity.py b/toqito/channel_metrics/channel_fidelity.py
index 622690e14..f94b90f6e 100644
--- a/toqito/channel_metrics/channel_fidelity.py
+++ b/toqito/channel_metrics/channel_fidelity.py
@@ -39,43 +39,28 @@ def channel_fidelity(choi_1: np.ndarray, choi_2: np.ndarray) -> float:
     :math:`1`.
 
 
+    >>> import numpy as np
     >>> from toqito.channels import dephasing
     >>> from toqito.channel_metrics import channel_fidelity
     >>>
     >>> # The Choi matrices of dimension-4 for the dephasing channel
     >>> choi_1 = dephasing(4)
     >>> choi_2 = dephasing(4)
-    >>> '%.2f' % channel_fidelity(choi_1, choi_2)
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(channel_fidelity(choi_1, choi_2), decimals=2)
+    1.0
 
     We can also compute the channel fidelity between two different channels. For example, we can
     compute the channel fidelity between the dephasing and depolarizing channels.
 
+    >>> import numpy as np
     >>> from toqito.channels import dephasing, depolarizing
     >>> from toqito.channel_metrics import channel_fidelity
     >>>
     >>> # The Choi matrices of dimension-4 for the dephasing and depolarizing channels
     >>> choi_1 = dephasing(4)
     >>> choi_2 = depolarizing(4)
-    >>> '%.2f' % channel_fidelity(choi_1, choi_2)
-    '0.50'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
+    >>> np.around(channel_fidelity(choi_1, choi_2), decimals=2)
+    0.5
 
     References
     ==========
diff --git a/toqito/channel_metrics/diamond_norm.py b/toqito/channel_metrics/diamond_norm.py
index a4fe906bf..b8549a398 100644
--- a/toqito/channel_metrics/diamond_norm.py
+++ b/toqito/channel_metrics/diamond_norm.py
@@ -26,18 +26,8 @@ def diamond_norm(choi_1: np.ndarray, choi_2: np.ndarray) -> float:
     >>> from toqito.channel_metrics import diamond_norm
     >>> choi_depolarizing = depolarizing(dim=2, param_p=0.2)
     >>> choi_identity = np.identity(2**2)
-    >>> dn = diamond_norm(choi_depolarizing, choi_identity)
-    >>> print("Diamond norm between depolarizing and identity channels: ", '%.2f' % dn)
-    Diamond norm between depolarizing and identity channels:  -0.00
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
+    >>> np.around(diamond_norm(choi_depolarizing, choi_identity), decimals=2)
+    -0.0
 
     Similarly, we can compute the diamond norm between the dephasing channel (with parameter 0.3) and the identity
     channel:
@@ -47,17 +37,8 @@ def diamond_norm(choi_1: np.ndarray, choi_2: np.ndarray) -> float:
     >>> from toqito.channel_metrics import diamond_norm
     >>> choi_dephasing = dephasing(dim=2)
     >>> choi_identity = np.identity(2**2)
-    >>> dn = diamond_norm(choi_dephasing, choi_identity)
-    >>> print("Diamond norm between dephasing and identity channels: ", '%.2f' % dn)
-    Diamond norm between dephasing and identity channels:  -0.00
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(diamond_norm(choi_dephasing, choi_identity), decimals=2)
+    -0.0
 
     References
     ==========
diff --git a/toqito/channel_metrics/fidelity_of_separability.py b/toqito/channel_metrics/fidelity_of_separability.py
index c979921b8..4535bc8c7 100644
--- a/toqito/channel_metrics/fidelity_of_separability.py
+++ b/toqito/channel_metrics/fidelity_of_separability.py
@@ -81,23 +81,14 @@ def fidelity_of_separability(
     .. math::
         \rho_{AB} = |000 \rangle \langle 000|
 
+    >>> import numpy as np
     >>> from toqito.state_metrics import fidelity_of_separability
     >>> from toqito.matrix_ops import tensor
     >>> from toqito.states import basis
     >>> state = tensor(basis(2, 0), basis(2, 0))
     >>> rho = state @ state.conj().T
-    >>> expected_value = fidelity_of_separability(rho, [2, 2])
-    >>> '%.2f' % expected_value
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
+    >>> np.around(fidelity_of_separability(rho, [2, 2]), decimals=2)
+    1.0
 
     References
     ==========
diff --git a/toqito/matrix_props/kp_norm.py b/toqito/matrix_props/kp_norm.py
index 9465caf1d..ab9c77d74 100644
--- a/toqito/matrix_props/kp_norm.py
+++ b/toqito/matrix_props/kp_norm.py
@@ -16,8 +16,8 @@ def kp_norm(mat: np.ndarray, k: int, p: int) -> float:
     >>> import numpy as np
     >>> from toqito.matrix_props import kp_norm
     >>> from toqito.states import bell
-    >>> '%.2f' % kp_norm(bell(0), 1, np.inf)
-    '1.00'
+    >>> np.around(kp_norm(bell(0), 1, np.inf), decimals=2)
+    1.0
 
 
     To compute the k-largest singular values of a matrix:
@@ -25,15 +25,8 @@ def kp_norm(mat: np.ndarray, k: int, p: int) -> float:
     >>> import numpy as np
     >>> from toqito.matrix_props import kp_norm
     >>> from toqito.rand import random_unitary
-    >>> '%.2f' % kp_norm(random_unitary(5), 5, 2)
-    '2.24'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(kp_norm(random_unitary(5), 5, 2), decimals=2)
+    2.24
 
     :param mat: 2D numpy ndarray
     :param k: The number of singular values to take.
diff --git a/toqito/nonlocal_games/quantum_hedging.py b/toqito/nonlocal_games/quantum_hedging.py
index 252d45fa2..eb2ed3252 100644
--- a/toqito/nonlocal_games/quantum_hedging.py
+++ b/toqito/nonlocal_games/quantum_hedging.py
@@ -31,6 +31,7 @@ class QuantumHedging:
     As was illustrated in :cite:`Molina_2012_Hedging`, the hedging value of the above scenario is
     :math:`\cos(\pi/8)^2 \approx 0.8536`
 
+    >>> import numpy as np
     >>> from numpy import kron, cos, sin, pi, sqrt, isclose
     >>> from toqito.states import basis
     >>> from toqito.nonlocal_games.quantum_hedging import QuantumHedging
@@ -52,16 +53,8 @@ class QuantumHedging:
     >>> molina_watrous = QuantumHedging(q_0, 1)
     >>>
     >>> # cos(pi/8)**2 \approx 0.8536
-    >>> '%.2f' % molina_watrous.max_prob_outcome_a_primal()
-    '0.85'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(molina_watrous.max_prob_outcome_a_primal(), decimals=2)
+    0.85
 
     References
     ==========
diff --git a/toqito/nonlocal_games/xor_game.py b/toqito/nonlocal_games/xor_game.py
index d9410855f..a0d9ffc62 100644
--- a/toqito/nonlocal_games/xor_game.py
+++ b/toqito/nonlocal_games/xor_game.py
@@ -71,20 +71,12 @@ class XORGame:
     >>> import numpy as np
     >>> from toqito.nonlocal_games.xor_game import XORGame
     >>> chsh = XORGame(prob_mat, pred_mat)
-    >>> '%.2f' % chsh.quantum_value()
-    '0.85'
+    >>> np.around(chsh.quantum_value(), decimals=2)
+    0.85
     >>>
     >>> chsh.classical_value()
     0.75
 
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
     The odd cycle game
 
     The odd cycle game is another XOR game :cite:`Cleve_2010_Consequences`. For this game, we can
@@ -115,18 +107,10 @@ class XORGame:
     >>> import numpy as np
     >>> from toqito.nonlocal_games.xor_game import XORGame
     >>> odd_cycle = XORGame(prob_mat, pred_mat)
-    >>> '%.2f' % odd_cycle.quantum_value()
-    '0.98'
-    >>> '%.1f' % odd_cycle.classical_value()
-    '0.9'
-
-    .. note::
-        You do not need to use `'%.2f' %` or `'%.1f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(odd_cycle.quantum_value(), decimals=2)
+    0.98
+    >>> np.around(odd_cycle.classical_value(), decimals=1)
+    0.9
 
     References
     ==========
diff --git a/toqito/state_metrics/fidelity_of_separability.py b/toqito/state_metrics/fidelity_of_separability.py
index 9dba7624e..caf1602f1 100644
--- a/toqito/state_metrics/fidelity_of_separability.py
+++ b/toqito/state_metrics/fidelity_of_separability.py
@@ -87,22 +87,14 @@ def fidelity_of_separability(
         \rho_{AB} = |00 \rangle \langle 00|
 
 
+    >>> import numpy as np
     >>> from toqito.state_metrics import fidelity_of_separability
     >>> from toqito.matrix_ops import tensor
     >>> from toqito.states import basis
     >>> state = tensor(basis(2, 0), basis(2, 0))
     >>> rho = state @ state.conj().T
-    >>> expected_value = fidelity_of_separability(rho, [2, 2])
-    >>> '%.2f' % expected_value
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(fidelity_of_separability(rho, [2, 2]), decimals=2)
+    1.0
 
     References
     ==========
diff --git a/toqito/state_metrics/hilbert_schmidt.py b/toqito/state_metrics/hilbert_schmidt.py
index 87f2bfdc8..8acce7aff 100644
--- a/toqito/state_metrics/hilbert_schmidt.py
+++ b/toqito/state_metrics/hilbert_schmidt.py
@@ -20,20 +20,13 @@ def hilbert_schmidt(rho: np.ndarray, sigma: np.ndarray) -> float:
     One may consider taking the Hilbert-Schmidt distance between two Bell states. In :code:`toqito`, one may accomplish
     this as
 
+    >>> import numpy as np
     >>> from toqito.states import bell
     >>> from toqito.state_metrics import hilbert_schmidt
     >>> rho = bell(0) * bell(0).conj().T
     >>> sigma = bell(3) * bell(3).conj().T
-    >>> '%.2f' % hilbert_schmidt(rho, sigma)
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(hilbert_schmidt(rho, sigma), decimals=2)
+    1.0
 
     References
     ==========
diff --git a/toqito/state_metrics/hilbert_schmidt_inner_product.py b/toqito/state_metrics/hilbert_schmidt_inner_product.py
index a68bb7d78..7b25910ee 100644
--- a/toqito/state_metrics/hilbert_schmidt_inner_product.py
+++ b/toqito/state_metrics/hilbert_schmidt_inner_product.py
@@ -21,19 +21,12 @@ def hilbert_schmidt_inner_product(a_mat: np.ndarray, b_mat: np.ndarray) -> compl
 
     One may consider taking the Hilbert-Schmidt distance between two Hadamard matrices.
 
+    >>> import numpy as np
     >>> from toqito.matrices import hadamard
     >>> from toqito.state_metrics import hilbert_schmidt_inner_product
     >>> h = hadamard(1)
-    >>> '%.2f' % hilbert_schmidt_inner_product(h, h)
-    '2.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(hilbert_schmidt_inner_product(h, h), decimals=2)
+    2.0
 
     References
     ==========
diff --git a/toqito/state_metrics/matsumoto_fidelity.py b/toqito/state_metrics/matsumoto_fidelity.py
index d7767849a..7641b8396 100644
--- a/toqito/state_metrics/matsumoto_fidelity.py
+++ b/toqito/state_metrics/matsumoto_fidelity.py
@@ -59,16 +59,8 @@ def matsumoto_fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:
     ...      [1, 0, 0, 1]]
     ... )
     >>> sigma = rho
-    >>> '%.2f' % matsumoto_fidelity(rho, sigma)
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(matsumoto_fidelity(rho, sigma), decimals=2)
+    1.0
 
     References
     ==========
diff --git a/toqito/state_opt/optimal_clone.py b/toqito/state_opt/optimal_clone.py
index b086af2df..e9048191e 100644
--- a/toqito/state_opt/optimal_clone.py
+++ b/toqito/state_opt/optimal_clone.py
@@ -83,16 +83,8 @@ def optimal_clone(
     >>>
     >>> states = [e_0, e_1, e_p, e_m]
     >>> probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
-    >>> '%.2f' % optimal_clone(states, probs)
-    '0.75'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(optimal_clone(states, probs), decimals=2)
+    0.75
 
     References
     ==========
diff --git a/toqito/state_opt/ppt_distinguishability.py b/toqito/state_opt/ppt_distinguishability.py
index ce3bab406..35f514447 100644
--- a/toqito/state_opt/ppt_distinguishability.py
+++ b/toqito/state_opt/ppt_distinguishability.py
@@ -91,14 +91,6 @@ def ppt_distinguishability(
     >>> '%.3f' % opt_val
     '0.875'
 
-    .. note::
-        You do not need to use `'%.3f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
     References
     ==========
     .. bibliography::
diff --git a/toqito/state_opt/state_distinguishability.py b/toqito/state_opt/state_distinguishability.py
index 0f7d51c67..be5b4ea53 100644
--- a/toqito/state_opt/state_distinguishability.py
+++ b/toqito/state_opt/state_distinguishability.py
@@ -78,21 +78,14 @@ def state_distinguishability(
 
     Minimal-error state distinguishability for the Bell states (which are perfectly distinguishable).
 
+    >>> import numpy as np
     >>> from toqito.states import bell
     >>> from toqito.state_opt import state_distinguishability
     >>> states = [bell(0), bell(1), bell(2), bell(3)]
     >>> probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
     >>> res, _ = state_distinguishability(vectors=states, probs=probs, primal_dual="dual")
-    >>> '%.2f' % res
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(res, decimals=2)
+    1.0
 
     Note that if we are just interested in obtaining the optimal value, it is computationally less intensive to compute
     the dual problem over the primal problem. However, the primal problem does allow us to extract the explicit
@@ -117,8 +110,8 @@ def state_distinguishability(
     >>> states = [np.array([[1.], [0.]]), np.array([[1.],[1.]]) / np.sqrt(2)]
     >>> probs = [1 / 2, 1 / 2]
     >>> res, _ = state_distinguishability(vectors=states, probs=probs, primal_dual="primal", strategy="unambiguous")
-    >>> '%.2f' % res
-    '0.29'
+    >>> np.around(res, decimals=2)
+    0.29
 
     References
     ==========
diff --git a/toqito/state_opt/state_exclusion.py b/toqito/state_opt/state_exclusion.py
index 0138902a9..18d4091c4 100644
--- a/toqito/state_opt/state_exclusion.py
+++ b/toqito/state_opt/state_exclusion.py
@@ -113,8 +113,8 @@ def state_exclusion(
     >>> vectors = [bell(0), bell(1)]
     >>> probs = [1/2, 1/2]
     >>>
-    >>> '%.2f' % state_exclusion(vectors, probs)[0]
-    '0.00'
+    >>> np.around(state_exclusion(vectors, probs)[0], decimals=2)
+    0.0
 
     Unambiguous state exclusion for unbiased states.
 
@@ -122,16 +122,8 @@ def state_exclusion(
     >>> import numpy as np
     >>> states = [np.array([[1.], [0.]]), np.array([[1.],[1.]]) / np.sqrt(2)]
     >>> res, _ = state_exclusion(states, primal_dual="primal", strategy="unambiguous", abs_ipm_opt_tol=1e-7)
-    >>> '%.2f' % res
-    '0.71'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> np.around(res, decimals=2)
+    0.71
 
     .. note::
         If you encounter a `ZeroDivisionError` or an `ArithmeticError` when using cvxopt as a solver (which is the
diff --git a/toqito/state_opt/symmetric_extension_hierarchy.py b/toqito/state_opt/symmetric_extension_hierarchy.py
index 66298dc70..46f3a42f6 100644
--- a/toqito/state_opt/symmetric_extension_hierarchy.py
+++ b/toqito/state_opt/symmetric_extension_hierarchy.py
@@ -108,30 +108,22 @@ def symmetric_extension_hierarchy(
     >>>
     >>> # Calculate the first level of the symmetric extension hierarchy. This
     >>> # is simply the value of optimally distinguishing via PPT measurements.
-    >>> '%.2f' % symmetric_extension_hierarchy(states=states, probs=None, level=1)
-    '0.99'
+    >>> np.around(symmetric_extension_hierarchy(states=states, probs=None, level=1), decimals=2)
+    0.99
     >>>
     >>> # Calculating the second value gets closer to the separable value.
-    >>> '%.2f' % symmetric_extension_hierarchy(states=states, probs=None, level=2)
-    '0.96'
+    >>> np.around(symmetric_extension_hierarchy(states=states, probs=None, level=2), decimals=2)
+    0.96
     >>>
     >>> # As proven in :cite:`Cosentino_2015_QuantumState`, the true separable value of distinguishing the
     >>> # three Bell states is:
-    >>> '%.2f' % (1/3 * (2 + np.sqrt(1 - eps**2)))
-    '0.96'
+    >>> np.around(1/3 * (2 + np.sqrt(1 - eps**2)), decimals=2)
+    0.96
     >>>
     >>> # Computing further levels of the hierarchy would eventually converge to
     >>> # this value, however, the higher the level, the more computationally
     >>> # demanding the SDP becomes.
 
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
     References
     ==========
     .. bibliography::
diff --git a/toqito/state_props/concurrence.py b/toqito/state_props/concurrence.py
index 54fa2a37e..58fb03aab 100644
--- a/toqito/state_props/concurrence.py
+++ b/toqito/state_props/concurrence.py
@@ -38,22 +38,14 @@ def concurrence(rho: np.ndarray) -> float:
     The following example calculates this quantity using the :code:`toqito` package.
 
     >>> import numpy as np
-    >>> from toqito.states import basis
+    >>> from toqito.matrices import standard_basis
     >>> from toqito.state_props import concurrence
-    >>> e_0, e_1 = basis(2, 0), basis(2, 1)
+    >>> e_0, e_1 = standard_basis(2)
     >>> e_00, e_11 = np.kron(e_0, e_0), np.kron(e_1, e_1)
     >>> u_vec = 1 / np.sqrt(2) * (e_00 + e_11)
-    >>> rho = u_vec * u_vec.conj().T
-    >>> '%.2f' % concurrence(rho)
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> rho = u_vec @ u_vec.conj().T
+    >>> np.around(concurrence(rho), decimals=3)
+    1.0
 
     Consider the concurrence of the following product state
 
@@ -67,7 +59,7 @@ def concurrence(rho: np.ndarray) -> float:
     >>> from toqito.state_props import concurrence
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
     >>> v_vec = np.kron(e_0, e_1)
-    >>> sigma = v_vec * v_vec.conj().T
+    >>> sigma = v_vec @ v_vec.conj().T
     >>> concurrence(sigma)
     0
 
diff --git a/toqito/state_props/entanglement_of_formation.py b/toqito/state_props/entanglement_of_formation.py
index d02ccea5e..83e42fe7c 100644
--- a/toqito/state_props/entanglement_of_formation.py
+++ b/toqito/state_props/entanglement_of_formation.py
@@ -35,21 +35,14 @@ def entanglement_of_formation(rho: np.ndarray, dim: list[int] | int = None) -> f
 
     The entanglement-of-formation of :math:`\rho` is equal to 1.
 
+    >>> import numpy as np
     >>> from toqito.state_props import entanglement_of_formation
     >>> from toqito.states import bell
     >>>
     >>> u_vec = bell(0)
-    >>> rho = u_vec * u_vec.conj().T
-    >>> '%.2f' % entanglement_of_formation(rho)
-    '1.00'
-
-    .. note::
-        You do not need to use `'%.2f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
+    >>> rho = u_vec @ u_vec.conj().T
+    >>> np.around(entanglement_of_formation(rho), decimals=3)
+    1.0
 
     References
     ==========
diff --git a/toqito/state_props/l1_norm_coherence.py b/toqito/state_props/l1_norm_coherence.py
index bd127fe27..bb7520c26 100644
--- a/toqito/state_props/l1_norm_coherence.py
+++ b/toqito/state_props/l1_norm_coherence.py
@@ -39,14 +39,6 @@ def l1_norm_coherence(rho: np.ndarray) -> float:
     >>> '%.1f' % l1_norm_coherence(v)
     '2.0'
 
-    .. note::
-        You do not need to use `'%.1f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
     References
     ==========
     .. bibliography::
diff --git a/toqito/state_props/purity.py b/toqito/state_props/purity.py
index 915fe7d40..2f5b0ba86 100644
--- a/toqito/state_props/purity.py
+++ b/toqito/state_props/purity.py
@@ -44,14 +44,6 @@ def purity(rho: np.ndarray) -> float:
     >>> '%.4f' % purity(rho)
     '0.2653'
 
-    .. note::
-        You do not need to use `'%.4f' %` when you use this function.
-
-        We use this to format our output such that `doctest` compares the calculated output to the
-        expected output upto two decimal points only. The accuracy of the solvers can calculate the
-        `float` output to a certain amount of precision such that the value deviates after a few digits
-        of accuracy.
-
     References
     ==========
     .. bibliography::
diff --git a/toqito/states/bb84.py b/toqito/states/bb84.py
index 68946d9d2..fc0d36bda 100644
--- a/toqito/states/bb84.py
+++ b/toqito/states/bb84.py
@@ -29,13 +29,6 @@ def bb84() -> np.ndarray:
            [0.70710678]]), array([[ 0.70710678],
            [-0.70710678]])]]
 
-
-    .. note::
-        We use `#doctest` to verify our examples work as expected.
-
-        `# doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE` is used here to make sure doctest matches the
-        calculated output to be close to the expected output instead of making sure both match exactly.
-
     References
     ==========
     .. bibliography::