diff --git a/fast_pauli/cpp/examples/04_get_sparse_repr.cpp b/fast_pauli/cpp/examples/04_get_sparse_repr.cpp
new file mode 100644
index 0000000..225b32d
--- /dev/null
+++ b/fast_pauli/cpp/examples/04_get_sparse_repr.cpp
@@ -0,0 +1,17 @@
+#include <string>
+
+#include "fast_pauli.hpp"
+
+namespace fp = fast_pauli;
+
+int main()
+{
+    size_t const n_qubits = 24;
+    std::string pauli_string(n_qubits, 'X');
+    fp::PauliString ps(pauli_string);
+    std::vector<size_t> k;
+    std::vector<std::complex<double>> m;
+    std::tie(k, m) = get_sparse_repr<double>(ps.paulis);
+
+    return 0;
+}
\ No newline at end of file
diff --git a/fast_pauli/cpp/include/__pauli_string.hpp b/fast_pauli/cpp/include/__pauli_string.hpp
index 7bcf03f..cbef31a 100644
--- a/fast_pauli/cpp/include/__pauli_string.hpp
+++ b/fast_pauli/cpp/include/__pauli_string.hpp
@@ -64,7 +64,7 @@ std::tuple<std::vector<size_t>, std::vector<std::complex<T>>> get_sparse_repr(st
 
     // Helper function that let's us know if a pauli matrix has diagonal (or
     // conversely off-diagonal) elements
-    auto diag = [](Pauli const &p) {
+    auto diag = [](Pauli const &p) -> size_t {
         if (p.code == 0 || p.code == 3)
         {
             return 0UL;
@@ -98,18 +98,13 @@ std::tuple<std::vector<size_t>, std::vector<std::complex<T>>> get_sparse_repr(st
     }
     m[0] = initial_value();
 
-    // Populate the rest of the values in a recursive-like manner
     for (size_t l = 0; l < n; ++l)
     {
         Pauli const &po = ps[l];
 
-        T eps = 1.0;
-        if (po.code == 2 || po.code == 3)
-        {
-            eps = -1;
-        }
+        T const eps = (po.code == 2 || po.code == 3) ? -1.0 : 1.0;
 
-        T sign = diag(po) ? -1.0 : 1.0;
+        int sign = diag(po) ? -1 : 1;
 
         auto const lower_bound = 1UL << l;
         for (size_t li = lower_bound; li < (lower_bound << 1); li++)
diff --git a/fast_pauli/cpp/include/__summed_pauli_op.hpp b/fast_pauli/cpp/include/__summed_pauli_op.hpp
index 856a9c0..c7d088b 100644
--- a/fast_pauli/cpp/include/__summed_pauli_op.hpp
+++ b/fast_pauli/cpp/include/__summed_pauli_op.hpp
@@ -539,13 +539,16 @@ template <std::floating_point T> struct SummedPauliOp
      */
     void to_tensor(Tensor<3> A_k_out) const
     {
-        for (size_t i = 0; i < pauli_strings.size(); ++i)
+#pragma omp parallel for schedule(static)
+        for (size_t k = 0; k < n_operators(); ++k)
         {
-            PauliString const &ps = pauli_strings[i];
-            auto [cols, vals] = get_sparse_repr<T>(ps.paulis);
-
-            for (size_t k = 0; k < n_operators(); ++k)
+            for (size_t i = 0; i < pauli_strings.size(); ++i)
             {
+                PauliString const &ps = pauli_strings[i];
+                std::vector<size_t> cols;
+                std::vector<std::complex<T>> vals;
+                std::tie(cols, vals) = get_sparse_repr<T>(ps.paulis);
+
                 for (size_t j = 0; j < dim(); ++j)
                 {
                     A_k_out(k, j, cols[j]) += coeffs(i, k) * vals[j];
diff --git a/tests/fast_pauli/test_summed_pauli_op.py b/tests/fast_pauli/test_summed_pauli_op.py
index 7c81a38..0f96181 100644
--- a/tests/fast_pauli/test_summed_pauli_op.py
+++ b/tests/fast_pauli/test_summed_pauli_op.py
@@ -101,12 +101,11 @@ def test_apply(
     "summed_pauli_op", [fp.SummedPauliOp], ids=resolve_parameter_repr
 )
 @pytest.mark.parametrize(
-    "n_states,n_operators,n_qubits",
-    [(s, o, q) for s in [1, 10, 1000] for o in [1, 10, 100] for q in [1, 2, 6]],
+    "n_operators,n_qubits",
+    [(o, q) for o in [1, 10, 100] for q in [1, 2, 6]],
 )
 def test_to_tensor(
     summed_pauli_op: type[fp.SummedPauliOp],
-    n_states: int,
     n_operators: int,
     n_qubits: int,
 ) -> None: