diff --git a/include/__factory.hpp b/include/__factory.hpp
index 51f08b7..ea9e255 100644
--- a/include/__factory.hpp
+++ b/include/__factory.hpp
@@ -75,7 +75,7 @@ auto rand(std::vector<T> &blob, std::array<size_t, n_dim> extents) {
   for (auto ei : extents) {
     total_size *= ei;
   }
-  blob.reserve(total_size);
+  blob = std::vector<T>(total_size);
 
   // Fill with random numbers
   std::random_device rd;
diff --git a/tests/test_pauli_op.cpp b/tests/test_pauli_op.cpp
index b6e8981..db718a7 100644
--- a/tests/test_pauli_op.cpp
+++ b/tests/test_pauli_op.cpp
@@ -156,8 +156,6 @@ TEST_CASE("test apply multistate multistring") {
   // Check
   for (size_t t = 0; t < n_states; ++t) {
     for (size_t i = 0; i < dims; ++i) {
-      fmt::println("new_states(i, t): {}, expected_span(i, t): {}",
-                   new_states(i, t), expected_span(i, t));
       CHECK(abs(new_states(i, t) - expected_span(i, t)) < 1e-6);
     }
   }
@@ -175,7 +173,7 @@ TEST_CASE("test apply multistate multistring identity") {
 
   // Set up random states
   size_t const n_states = 10;
-  std::vector<std::complex<double>> states_raw;
+  std::vector<std::complex<double>> states_raw(dims * n_states, 1);
   std::mdspan states =
       fast_pauli::rand<std::complex<double>, 2>(states_raw, {dims, n_states});
 
diff --git a/tests/test_summed_pauli_op.cpp b/tests/test_summed_pauli_op.cpp
index 5c48c37..e5aeb5c 100644
--- a/tests/test_summed_pauli_op.cpp
+++ b/tests/test_summed_pauli_op.cpp
@@ -1,6 +1,8 @@
+#include "__pauli.hpp"
+#include "__pauli_string.hpp"
+#include <cstddef>
 #define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
 
-// #include "__summed_pauli_op.hpp"
 #include <doctest/doctest.h>
 #include <fmt/core.h>
 #include <fmt/ranges.h>
@@ -10,6 +12,72 @@
 using namespace std::literals;
 using namespace fast_pauli;
 
+/**
+ * @brief Helper function to check the apply function on a set of states for a
+ * given set of pauli strings and coefficients.
+ *
+ * @param pauli_strings
+ * @param coeff
+ */
+void __check_apply(
+    std::vector<PauliString> &pauli_strings,
+    std::mdspan<std::complex<double>, std::dextents<size_t, 2>> coeff,
+    size_t const n_states = 10) {
+
+  SummedPauliOp<double> summed_op{pauli_strings, coeff};
+
+  // Setup states
+  size_t const dim = summed_op.n_dimensions();
+  size_t const n_ops = summed_op.n_operators();
+
+  std::vector<std::complex<double>> states_raw;
+  std::mdspan states =
+      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
+
+  std::vector<std::complex<double>> new_states_raw;
+  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
+      new_states_raw, {dim, n_states});
+
+  // Init data (aka data)
+  std::vector<double> data_raw;
+  std::mdspan data = fast_pauli::rand<double, 2>(data_raw, {n_ops, n_states});
+
+  // Apply the summed operator
+  summed_op.apply(new_states, states, data);
+
+  // Check the check
+  std::vector<std::complex<double>> expected_raw;
+  std::mdspan expected =
+      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
+
+  for (size_t j = 0; j < summed_op.n_pauli_strings(); ++j) {
+    auto ps = summed_op.pauli_strings[j];
+
+    PauliOp<double> pop{{std::complex<double>(1.)}, {ps}};
+    std::vector<std::complex<double>> tmp_raw;
+    std::mdspan tmp =
+        fast_pauli::zeros<std::complex<double>, 2>(tmp_raw, {dim, n_states});
+
+    pop.apply(tmp, states);
+
+    // Manually calculate the sum over the different pauli operators
+    // This is specific to the coefficients we've chosen above
+    for (size_t i = 0; i < dim; ++i) {
+      for (size_t t = 0; t < n_states; ++t) {
+        for (size_t k = 0; k < n_ops; ++k) {
+          expected(i, t) += data(k, t) * summed_op.coeffs(j, k) * tmp(i, t);
+        }
+      }
+    }
+  }
+
+  for (size_t t = 0; t < n_states; ++t) {
+    for (size_t i = 0; i < dim; ++i) {
+      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
+    }
+  }
+}
+
 TEST_CASE("ctors") {
   {
     // Default
@@ -57,148 +125,38 @@ TEST_CASE("accessors") {
 TEST_CASE("apply 1 operator 1 PauliString") {
   fmt::print("\n\napply 1 operator 1 PauliString\n");
   // Setup operator
-  SummedPauliOp<double> summed_op{{"XYZ"}, {std::complex<double>(1.)}};
-
-  // Setup states
-  size_t const n_states = 1;
-  size_t const dim = summed_op.n_dimensions();
-  size_t const n_ops = summed_op.n_operators();
 
-  std::vector<std::complex<double>> states_raw;
-  std::mdspan states =
-      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
+  std::vector<PauliString> pauli_strings = {"XYZ"};
+  std::vector<std::complex<double>> coeff_raw = {1i};
+  std::mdspan<std::complex<double>, std::dextents<size_t, 2>> coeff(
+      coeff_raw.data(), 1, 1);
 
-  std::vector<std::complex<double>> new_states_raw;
-  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
-      new_states_raw, {dim, n_states});
-
-  // Init data
-  std::vector<double> data_raw(n_ops * n_states, 1);
-  std::mdspan<double, std::dextents<size_t, 2>> data(data_raw.data(), n_ops,
-                                                     n_states);
-  summed_op.apply(new_states, states, data);
-
-  // Calculated the expected answer and check
-  std::vector<std::complex<double>> expected_raw;
-  std::mdspan expected =
-      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
-
-  PauliOp<double> pop{{std::complex<double>(1.)}, {"XYZ"}};
-  pop.apply(expected, states);
-
-  for (size_t t = 0; t < n_states; ++t) {
-    for (size_t i = 0; i < dim; ++i) {
-      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
-    }
-  }
+  __check_apply(pauli_strings, coeff);
 }
 
 TEST_CASE("apply 2 operators 1 PauliString") {
   fmt::print("\n\napply 2 operators 1 PauliString\n");
 
   // Setup operator
-  SummedPauliOp<double> summed_op{{"XYZ"}, {std::complex<double>(1.), 1i}};
-
-  // Setup states
-  size_t const n_states = 1;
-  size_t const dim = summed_op.n_dimensions();
-  size_t const n_ops = summed_op.n_operators();
-
-  std::vector<std::complex<double>> states_raw;
-  std::mdspan states =
-      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
-
-  std::vector<std::complex<double>> new_states_raw;
-  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
-      new_states_raw, {dim, n_states});
-
-  // Init data (aka data)
-  std::vector<double> data_raw(n_ops * n_states, 1);
-  std::mdspan<double, std::dextents<size_t, 2>> data(data_raw.data(), n_ops,
-                                                     n_states);
-  summed_op.apply(new_states, states, data);
-
-  // Calculated the expected answer and check
-  std::vector<std::complex<double>> expected_raw;
-  std::mdspan expected =
-      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
-
-  // Manually calculate the sum over the different pauli operators
-  for (size_t i = 0; i < dim; ++i) {
-    for (size_t t = 0; t < n_states; ++t) {
-      expected(i, t) += 1i * expected(i, t);
-    }
-  }
-
-  for (size_t t = 0; t < n_states; ++t) {
-    for (size_t i = 0; i < dim; ++i) {
-      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
-    }
-  }
+  std::vector<PauliString> pauli_strings = {"XYZ"};
+  std::vector<std::complex<double>> coeff_raw = {1i, 1};
+  std::mdspan<std::complex<double>, std::dextents<size_t, 2>> coeff(
+      coeff_raw.data(), 1, 2);
+  __check_apply(pauli_strings, coeff);
 }
 
 TEST_CASE("apply 2 operators 2 PauliString") {
   fmt::print("\n\napply 2 operators 2 PauliString\n");
   // Setup operator
-  SummedPauliOp<double> summed_op{{"XYZ", "YYZ"}, {1i, 1, 0.5i, -0.99}};
-
-  // Setup states
-  size_t const n_states = 10;
-  size_t const dim = summed_op.n_dimensions();
-  size_t const n_ops = summed_op.n_operators();
-
-  std::vector<std::complex<double>> states_raw;
-  std::mdspan states =
-      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
-
-  std::vector<std::complex<double>> new_states_raw;
-  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
-      new_states_raw, {dim, n_states});
-
-  // Init data
-  std::vector<double> data_raw(n_ops * n_states, 1);
-  std::mdspan<double, std::dextents<size_t, 2>> data(data_raw.data(), n_ops,
-                                                     n_states);
-  summed_op.apply(new_states, states, data);
-
-  // Calculated expected answer
-  std::vector<std::complex<double>> expected_raw;
-  std::mdspan expected =
-      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
-
-  for (size_t j = 0; j < summed_op.n_pauli_strings(); ++j) {
-
-    auto ps = summed_op.pauli_strings[j];
-
-    PauliOp<double> pop{{std::complex<double>(1.)}, {ps}};
-    std::vector<std::complex<double>> tmp_raw;
-    std::mdspan tmp =
-        fast_pauli::zeros<std::complex<double>, 2>(tmp_raw, {dim, n_states});
-
-    pop.apply(tmp, states);
-
-    // Manually calculate the sum over the different pauli operators
-    // This is specific to the coefficients we've chosen above
-    for (size_t i = 0; i < dim; ++i) {
-      for (size_t t = 0; t < n_states; ++t) {
-        // expected(i, t) += 1i * tmp(i, t);
-        for (size_t k = 0; k < n_ops; ++k) {
-          expected(i, t) += data(k, t) * summed_op.coeffs(j, k) * tmp(i, t);
-        }
-        // expected(i, t) += std::complex<double>(2) * tmp(i, t);
-      }
-    }
-  }
-
-  for (size_t t = 0; t < n_states; ++t) {
-    for (size_t i = 0; i < dim; ++i) {
-      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
-    }
-  }
+  std::vector<std::complex<double>> coeff_raw = {1i, 1, 0.5i, -0.99};
+  std::mdspan<std::complex<double>, std::dextents<size_t, 2>> coeff(
+      coeff_raw.data(), 2, 2);
+  std::vector<PauliString> pauli_strings = {"XYZ", "YYZ"};
+  __check_apply(pauli_strings, coeff);
 }
 
 TEST_CASE("apply many operators many PauliString") {
-  fmt::print("\n\napply many operators many PauliString\n");
+  fmt::print("\n\napply many operators mansize_t Ranky PauliString\n");
   // Setup operator
   std::vector<PauliString> pauli_strings{"XIXXX", "IIXII", "ZYYZI",
                                          "ZYIIZ", "YXZZY", "IZYII"};
@@ -207,59 +165,7 @@ TEST_CASE("apply many operators many PauliString") {
   std::mdspan coeff = fast_pauli::rand<std::complex<double>, 2>(
       coeff_raw, {pauli_strings.size(), 100});
 
-  SummedPauliOp<double> summed_op{pauli_strings, coeff};
-
-  // Setup states
-  size_t const n_states = 10;
-  size_t const dim = summed_op.n_dimensions();
-  size_t const n_ops = summed_op.n_operators();
-
-  std::vector<std::complex<double>> states_raw;
-  std::mdspan states =
-      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
-
-  std::vector<std::complex<double>> new_states_raw;
-  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
-      new_states_raw, {dim, n_states});
-
-  // Init data (aka data)
-  std::vector<double> data_raw;
-  std::mdspan data = fast_pauli::rand<double, 2>(data_raw, {n_ops, n_states});
-
-  // Apply the summed operator
-  summed_op.apply(new_states, states, data);
-
-  // Check the check
-  std::vector<std::complex<double>> expected_raw;
-  std::mdspan expected =
-      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
-
-  for (size_t j = 0; j < summed_op.n_pauli_strings(); ++j) {
-    auto ps = summed_op.pauli_strings[j];
-
-    PauliOp<double> pop{{std::complex<double>(1.)}, {ps}};
-    std::vector<std::complex<double>> tmp_raw;
-    std::mdspan tmp =
-        fast_pauli::zeros<std::complex<double>, 2>(tmp_raw, {dim, n_states});
-
-    pop.apply(tmp, states);
-
-    // Manually calculate the sum over the different pauli operators
-    // This is specific to the coefficients we've chosen above
-    for (size_t i = 0; i < dim; ++i) {
-      for (size_t t = 0; t < n_states; ++t) {
-        for (size_t k = 0; k < n_ops; ++k) {
-          expected(i, t) += data(k, t) * summed_op.coeffs(j, k) * tmp(i, t);
-        }
-      }
-    }
-  }
-
-  for (size_t t = 0; t < n_states; ++t) {
-    for (size_t i = 0; i < dim; ++i) {
-      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
-    }
-  }
+  __check_apply(pauli_strings, coeff);
 }
 
 TEST_CASE("apply many operators many PauliString") {
@@ -272,58 +178,5 @@ TEST_CASE("apply many operators many PauliString") {
   std::mdspan coeff = fast_pauli::rand<std::complex<double>, 2>(
       coeff_raw, {pauli_strings.size(), 100});
 
-  SummedPauliOp<double> summed_op{pauli_strings, coeff};
-
-  // Setup states
-  size_t const n_states = 10;
-  size_t const dim = summed_op.n_dimensions();
-  size_t const n_ops = summed_op.n_operators();
-
-  std::vector<std::complex<double>> states_raw;
-  std::mdspan states =
-      fast_pauli::rand<std::complex<double>, 2>(states_raw, {dim, n_states});
-
-  std::vector<std::complex<double>> new_states_raw;
-  std::mdspan new_states = fast_pauli::zeros<std::complex<double>, 2>(
-      new_states_raw, {dim, n_states});
-
-  // Init data
-  std::vector<double> data_raw;
-  std::mdspan data = fast_pauli::rand<double, 2>(data_raw, {n_ops, n_states});
-
-  // Apply the summed operator
-  summed_op.apply_parallel(new_states, states, data);
-
-  // Calculate the expected answer
-  std::vector<std::complex<double>> expected_raw;
-  std::mdspan expected =
-      fast_pauli::zeros<std::complex<double>, 2>(expected_raw, {dim, n_states});
-
-  for (size_t j = 0; j < summed_op.n_pauli_strings(); ++j) {
-    auto ps = summed_op.pauli_strings[j];
-
-    PauliOp<double> pop{{std::complex<double>(1.)}, {ps}};
-    std::vector<std::complex<double>> tmp_raw;
-    std::mdspan tmp =
-        fast_pauli::zeros<std::complex<double>, 2>(tmp_raw, {dim, n_states});
-
-    pop.apply(tmp, states);
-
-    // Manually calculate the sum over the different pauli operators
-    // This is specific to the coefficients we've chosen above
-    for (size_t i = 0; i < dim; ++i) {
-      for (size_t t = 0; t < n_states; ++t) {
-        // expected(i, t) += 1i * tmp(i, t);
-        for (size_t k = 0; k < n_ops; ++k) {
-          expected(i, t) += data(k, t) * summed_op.coeffs(j, k) * tmp(i, t);
-        }
-      }
-    }
-  }
-
-  for (size_t t = 0; t < n_states; ++t) {
-    for (size_t i = 0; i < dim; ++i) {
-      CHECK(abs(new_states(i, t) - expected(i, t)) < 1e-6);
-    }
-  }
+  __check_apply(pauli_strings, coeff);
 }
\ No newline at end of file