diff --git a/src/openfermion/transforms/opconversions/bksf.py b/src/openfermion/transforms/opconversions/bksf.py
index 35d4ecce6..968c1b7d4 100644
--- a/src/openfermion/transforms/opconversions/bksf.py
+++ b/src/openfermion/transforms/opconversions/bksf.py
@@ -125,8 +125,7 @@ def bravyi_kitaev_fast_interaction_op(iop: InteractionOperator) -> QubitOperator
                             qubit_operator += transformed_term
                             continue
                         elif p != r and q < p:
-                            # TODO: remove pragma if reachable continue
-                            continue  # pragma: no cover
+                            continue
 
                     # Handle the two-body terms.
                     transformed_term = _two_body(edge_matrix_indices, p, q, r, s)
@@ -169,8 +168,7 @@ def bravyi_kitaev_fast_edge_matrix(
 
                     # Skip zero terms.
                     if (not coefficient2) or (p == q) or (r == s):
-                        # TODO: remove pragma if this is a reachable continue
-                        continue  # pragma: no cover
+                        continue
 
                     # Identify and skip one of the complex conjugates.
                     if [p, q, r, s] != [s, r, q, p]:
diff --git a/src/openfermion/transforms/opconversions/bksf_test.py b/src/openfermion/transforms/opconversions/bksf_test.py
index dcae7f600..d5ac7baa0 100644
--- a/src/openfermion/transforms/opconversions/bksf_test.py
+++ b/src/openfermion/transforms/opconversions/bksf_test.py
@@ -21,6 +21,7 @@
 from openfermion.linalg import eigenspectrum, get_sparse_operator
 from openfermion.ops.operators import FermionOperator, QubitOperator
 from openfermion.ops.representations import InteractionOperator
+from openfermion.transforms import bravyi_kitaev_fast_interaction_op
 from openfermion.transforms.opconversions import bksf, get_fermion_operator, normal_ordered
 from openfermion.transforms.opconversions.jordan_wigner import jordan_wigner, jordan_wigner_one_body
 from openfermion.utils import count_qubits
@@ -258,3 +259,27 @@ def test_bravyi_kitaev_fast_number_excitation_operator(self):
                 evensector_n += 1
         self.assertEqual(evensector_H, 2 ** (n_qubits - 1))
         self.assertEqual(evensector_n, 2 ** (n_qubits - 1))
+
+    def test_bravyi_kitaev_fast_interaction_op_coverage(self):
+        n_qubits = 4
+        constant = 0.0
+        one_body_tensor = numpy.zeros((n_qubits, n_qubits))
+        two_body_tensor = numpy.zeros((n_qubits, n_qubits, n_qubits, n_qubits))
+
+        # Set up a Hermitian term with p != r and q < p
+        # For p = 2, q = 1, r = 3, s = 0: p != r (2 != 3) and q < p (1 < 2)
+        two_body_tensor[2, 1, 3, 0] = 1.0
+        two_body_tensor[1, 2, 3, 0] = -1.0
+        two_body_tensor[2, 1, 0, 3] = -1.0
+        two_body_tensor[1, 2, 0, 3] = 1.0
+
+        # Hermitian conjugate terms
+        two_body_tensor[3, 0, 2, 1] = 1.0
+        two_body_tensor[3, 0, 1, 2] = -1.0
+        two_body_tensor[0, 3, 2, 1] = -1.0
+        two_body_tensor[0, 3, 1, 2] = 1.0
+
+        iop = InteractionOperator(constant, one_body_tensor, two_body_tensor)
+        qubit_op = bravyi_kitaev_fast_interaction_op(iop)
+
+        self.assertIsNotNone(qubit_op)