diff --git a/thewalrus/_hafnian.py b/thewalrus/_hafnian.py
index 81ffc767..8c07810f 100644
--- a/thewalrus/_hafnian.py
+++ b/thewalrus/_hafnian.py
@@ -721,6 +721,7 @@ def hafnian(
     rtol=1e-05,
     atol=1e-08,
     num_samples=1000,
+    approx=False,
     method="glynn",
 ):  # pylint: disable=too-many-arguments
     """Returns the hafnian of a matrix.
@@ -733,11 +734,15 @@ def hafnian(
         method (string): Set this to ``"glynn"`` to use the
             glynn formula,
             or ``"inclexcl"`` to use the inclusion exclusion principle,
-            or ``"recursive"`` to use a recursive algorithm,
-            or ``"barvinok"`` to use an approximate Barvinok estimator (non-negative matrices only),
+            or ``"recursive"`` to use a recursive algorithm.
+        If ``approx`` is ``True``, one can use approximate methods:
+            (default) ``"barvinok"`` to use an approximate Barvinok estimator (non-negative matrices only),
             or ``"godsilgutman"`` to use an approximate Godsil-Gutman estimator (non-negative matrices only).
         rtol (float): the relative tolerance parameter used in ``np.allclose``
         atol (float): the absolute tolerance parameter used in ``np.allclose``
+        approx (bool): If ``True``, an approximation algorithm is used to estimate the hafnian. Note
+            that the approximation algorithm can only be applied to matrices ``A`` that only have
+            non-negative entries.
         num_samples (int): If ``method=barvinok`` or ``method=godsilgutman``, the approximation
             algorithm performs ``num_samples`` iterations for estimation of the hafnian
             of the non-negative matrix ``A``.
@@ -797,14 +802,17 @@ def hafnian(
 
         return A[0, 1] * A[2, 3] + A[0, 2] * A[1, 3] + A[0, 3] * A[1, 2]
 
-    if method in {"godsilgutman", "barvinok"}:
+    if approx:
         if np.any(np.iscomplex(A)):
             raise ValueError("Input matrix must be real")
 
         if np.any(A < 0):
             raise ValueError("Input matrix must not have negative entries")
 
-        return hafnian_approx(A, num_samples=num_samples, method=method)
+        if method == "godsilgutman":
+            return hafnian_approx(A, num_samples=num_samples, method=method)
+        else:
+            return hafnian_approx(A, num_samples=num_samples, method="barvinok")
 
     if loop:
         if method == "recursive":
diff --git a/thewalrus/tests/test_hafnian_approx.py b/thewalrus/tests/test_hafnian_approx.py
index d1abab24..4fa7a06e 100644
--- a/thewalrus/tests/test_hafnian_approx.py
+++ b/thewalrus/tests/test_hafnian_approx.py
@@ -29,7 +29,7 @@ def test_rank_one(n):
     x = np.random.rand(n)
     A = np.outer(x, x)
     exact = factorial2(n - 1) * np.prod(x)
-    approx = hafnian(A, method="barvinok", num_samples=10000)
+    approx = hafnian(A, approx=True, num_samples=10000)
     assert np.allclose(approx, exact, rtol=2e-1, atol=0)
 
 
@@ -37,7 +37,7 @@ def test_approx_complex_error():
     """Check exception raised if matrix is complex"""
     A = 1j * np.ones([6, 6])
     with pytest.raises(ValueError, match="Input matrix must be real"):
-        hafnian(A, method="barvinok")
+        hafnian(A, approx=True)
 
 
 def test_approx_negative_error():
@@ -45,13 +45,13 @@ def test_approx_negative_error():
     A = np.ones([6, 6])
     A[0, 0] = -1
     with pytest.raises(ValueError, match="Input matrix must not have negative entries"):
-        hafnian(A, method="barvinok")
+        hafnian(A, approx=True)
 
 
 @pytest.mark.parametrize("n", [6, 8])
 def test_ones_approx(n):
     """Check hafnian_approx(J_2n)=(2n)!/(n!2^n)"""
     A = np.float64(np.ones([2 * n, 2 * n]))
-    haf = hafnian(A, method="barvinok", num_samples=10000)
+    haf = hafnian(A, approx=True, num_samples=10000)
     expected = fac(2 * n) / (fac(n) * (2**n))
     assert np.abs(haf - expected) / expected < 0.2