diff --git a/README.md b/README.md
index ab1c4fe..b504db2 100644
--- a/README.md
+++ b/README.md
@@ -113,6 +113,26 @@ Once installed, you can import and use the Digest library in Python:
 ['ACGTA', 'CGTAG', 'AGCTA', 'TAGCT', 'GCTGA', 'TTACA', 'TACAT', 'GTATG', 'GCAAG', 'TGATC', 'CGTAG', 'TAGTG', 'ATGCT']
 ```
 
+We have also implemented parallel execution in the python library:
+```
+>>> import timeit
+>>> timeit.timeit('window_minimizer(seq, k=5, w=11)', 
+...     setup='from digest import window_minimizer; seq = open("tests/bench/chrY.fa", "rb").read()',
+...     number=10)
+19.85955114942044
+>>> timeit.timeit('window_minimizer(seq, k=5, w=11, num_threads=4)',
+...     setup='from digest import window_minimizer; seq = open("tests/bench/chrY.fa", "rb").read()',
+...     number=10)
+10.327348310500383
+```
+
+Functions can also return numpy arrays with the `*_np` function variant:
+```
+>>> from digest import window_minimizer_np
+>>> window_minimizer_np(b'ACGTACGTAGCTAGCTAGCTAGCTGATTACATACTGTATGCAAGCTAGCTGATCGATCGTAGCTAGTGATGCTAGCTAC', k=5, w=11)
+array([ 4,  5, 16, 19, 21, 26, 27, 35, 39, 49, 57, 63, 68], dtype=uint32)
+```
+
 <!---
 # Implementation
 Supports Mod Minimizers, Window Minimizers, and Syncmers  
diff --git a/include/digest/thread_out.hpp b/include/digest/thread_out.hpp
index 7846cf3..2d2c0f2 100644
--- a/include/digest/thread_out.hpp
+++ b/include/digest/thread_out.hpp
@@ -173,7 +173,7 @@ void thread_mod(
 		ind += assigned_kmer_am;
 	}
 	for (auto &t : thread_vector) {
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 }
 
@@ -235,7 +235,7 @@ void thread_mod(
 		ind += assigned_kmer_am;
 	}
 	for (auto &t : thread_vector) {
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 }
 
@@ -305,14 +305,14 @@ void thread_wind(
 			extras--;
 		}
 
-		thread_vector.emplace_back(std::async(thread_wind_roll1<P, T>, seq, ind,
+		thread_vector.push_back(std::async(thread_wind_roll1<P, T>, seq, ind,
 											  k, large_wind_kmer_am,
 											  minimized_h, assigned_lwind_am));
 
 		ind += assigned_lwind_am;
 	}
 	for (auto &t : thread_vector) {
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 
 	// handle duplicates
@@ -365,7 +365,7 @@ void thread_wind(
 	unsigned lwinds_per_thread = num_lwinds / thread_count;
 	unsigned extras = num_lwinds % thread_count;
 	vec.reserve(thread_count);
-	std::vector<std::future<std::pair<uint32_t, uint32_t>>> thread_vector;
+	std::vector<std::future<std::vector<std::pair<uint32_t, uint32_t>>>> thread_vector;
 
 	size_t ind = start;
 	for (unsigned i = 0; i < thread_count; i++) {
@@ -377,14 +377,14 @@ void thread_wind(
 			extras--;
 		}
 
-		thread_vector.emplace_back(std::async(thread_wind_roll2<P, T>, seq, ind,
+		thread_vector.push_back(std::async(thread_wind_roll2<P, T>, seq, ind,
 											  k, large_wind_kmer_am,
 											  minimized_h, assigned_lwind_am));
 
 		ind += assigned_lwind_am;
 	}
 	for (auto &t : thread_vector) {
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 
 	// handle duplicates
@@ -473,7 +473,7 @@ void thread_sync(
 	}
 	for (auto &t : thread_vector) {
 
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 }
 
@@ -535,7 +535,7 @@ void thread_sync(
 		ind += assigned_lwind_am;
 	}
 	for (auto &t : thread_vector) {
-		vec.emplace_back(t.get());
+		vec.emplace_back(std::move(t.get()));
 	}
 }
 
diff --git a/pybind/bindings.cpp b/pybind/bindings.cpp
index 33a785c..2319c24 100644
--- a/pybind/bindings.cpp
+++ b/pybind/bindings.cpp
@@ -6,15 +6,53 @@ namespace py = pybind11;
 
 PYBIND11_MODULE(digest, m) {
 	m.doc() = "bindings for digest";
-	m.def("window_minimizer", &window_minimizer,
-		  "A function that runs window minimizer digestion", py::arg("seq"),
-		  py::arg("k") = 31, py::arg("w") = 11,
+	
+	m.def("window_minimizer", 
+		  [](std::string_view seq, unsigned k, unsigned w, unsigned num_threads, bool include_hash) {
+			  return window_minimizer_wrapper(seq.data(), seq.size(), k, w, num_threads, include_hash);
+		  },
+		  "A function that runs window minimizer digestion",
+		  py::call_guard<py::gil_scoped_release>(),
+		  py::arg("seq"), py::arg("k") = 31, py::arg("w") = 11, py::arg("num_threads") = 1,
 		  py::arg("include_hash") = false);
-	m.def("modimizer", &modimizer,
-		  "A function that runs mod-minimizer digestion", py::arg("seq"),
-		  py::arg("k") = 31, py::arg("mod") = 100,
+
+	m.def("modimizer", 
+		  [](std::string_view seq, unsigned k, uint32_t mod, unsigned num_threads, bool include_hash) {
+			  return modimizer_wrapper(seq.data(), seq.size(), k, mod, num_threads, include_hash);
+		  },
+		  "A function that runs mod-minimizer digestion",
+		  py::call_guard<py::gil_scoped_release>(),
+		  py::arg("seq"), py::arg("k") = 31, py::arg("mod") = 100, py::arg("num_threads") = 1,
+		  py::arg("include_hash") = false);
+
+	m.def("syncmer", 
+		  [](std::string_view seq, unsigned k, unsigned w, unsigned num_threads, bool include_hash) {
+			  return syncmer_wrapper(seq.data(), seq.size(), k, w, num_threads, include_hash);
+		  },
+		  "A function that runs syncmer digestion",
+		  py::call_guard<py::gil_scoped_release>(),
+		  py::arg("seq"), py::arg("k") = 31, py::arg("w") = 11, py::arg("num_threads") = 1,
+		  py::arg("include_hash") = false);
+
+	m.def("window_minimizer_np",
+		  [](std::string_view seq, unsigned k, unsigned w, unsigned num_threads, bool include_hash) {
+			  return window_minimizer_numpy(seq.data(), seq.size(), k, w, num_threads, include_hash);
+		  },
+		  "A function that runs window minimizer digestion and returns a numpy array",
+		  py::arg("seq"), py::arg("k") = 31, py::arg("w") = 11, py::arg("num_threads") = 1,
+		  py::arg("include_hash") = false);
+	m.def("modimizer_np",
+		  [](std::string_view seq, unsigned k, uint32_t mod, unsigned num_threads, bool include_hash) {
+			  return modimizer_numpy(seq.data(), seq.size(), k, mod, num_threads, include_hash);
+		  },
+		  "A function that runs mod-minimizer digestion and returns a numpy array",
+		  py::arg("seq"), py::arg("k") = 31, py::arg("mod") = 100, py::arg("num_threads") = 1,
 		  py::arg("include_hash") = false);
-	m.def("syncmer", &syncmer, "A function that runs syncmer digestion",
-		  py::arg("seq"), py::arg("k") = 31, py::arg("w") = 11,
+	m.def("syncmer_np",
+		  [](std::string_view seq, unsigned k, unsigned w, unsigned num_threads, bool include_hash) {
+			  return syncmer_numpy(seq.data(), seq.size(), k, w, num_threads, include_hash);
+		  },
+		  "A function that runs syncmer digestion and returns a numpy array",
+		  py::arg("seq"), py::arg("k") = 31, py::arg("w") = 11, py::arg("num_threads") = 1,
 		  py::arg("include_hash") = false);
 }
diff --git a/pybind/digest_utils.hpp b/pybind/digest_utils.hpp
index f2547ca..63e23d8 100644
--- a/pybind/digest_utils.hpp
+++ b/pybind/digest_utils.hpp
@@ -1,52 +1,274 @@
 #include <digest/mod_minimizer.hpp>
 #include <digest/syncmer.hpp>
 #include <digest/window_minimizer.hpp>
+#include <digest/thread_out.hpp>
 #include <variant>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/pytypes.h>
+
+namespace py = pybind11;
 
 std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
-window_minimizer(const std::string &seq, unsigned k, unsigned large_window,
-				 bool include_hash = false) {
-	digest::WindowMin<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>
-		digester(seq, k, large_window);
-	if (include_hash) {
-		std::vector<std::pair<uint32_t, uint32_t>> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
-	} else {
-		std::vector<uint32_t> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
+window_minimizer_wrapper(const char *seq, size_t len, unsigned k, unsigned large_window,
+				unsigned threads, bool include_hash = false) {
+	if (threads == 1) {
+		digest::WindowMin<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>
+			digester(seq, len, k, large_window);
+		if (include_hash) {
+			std::vector<std::pair<uint32_t, uint32_t>> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		} else {
+			std::vector<uint32_t> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		}
+	}
+	// return {};
+	else {
+		if (include_hash) {
+			std::vector<std::vector<std::pair<uint32_t, uint32_t>>> output;
+			digest::thread_out::thread_wind<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>(
+				threads, output, seq, len, k, large_window);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			std::vector<std::pair<uint32_t, uint32_t>> concatenated;
+			concatenated.reserve(total_size);
+			for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		} else {
+			std::vector<std::vector<uint32_t>> output;
+			digest::thread_out::thread_wind<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>(
+				threads, output, seq, len, k, large_window);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			
+			std::vector<uint32_t> concatenated;
+			concatenated.reserve(total_size);
+		for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		}
 	}
 }
-// std::vector<std::pair<size_t, size_t>> output;
 
 std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
-modimizer(const std::string &seq, unsigned k, uint32_t mod,
-		  bool include_hash = false) {
-	digest::ModMin<digest::BadCharPolicy::SKIPOVER> digester(seq, k, mod);
-	if (include_hash) {
-		std::vector<std::pair<uint32_t, uint32_t>> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
-	} else {
-		std::vector<uint32_t> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
+modimizer_wrapper(const char *seq, size_t len, unsigned k, uint32_t mod,
+		  unsigned threads, bool include_hash = false) {
+	if (threads == 1) {
+		digest::ModMin<digest::BadCharPolicy::SKIPOVER> digester(seq, len, k, mod);
+		if (include_hash) {
+			std::vector<std::pair<uint32_t, uint32_t>> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		} else {
+			std::vector<uint32_t> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		}
+	}
+	else {
+		if (include_hash) {
+			std::vector<std::vector<std::pair<uint32_t, uint32_t>>> output;
+			digest::thread_out::thread_mod<digest::BadCharPolicy::SKIPOVER>(
+				threads, output, seq, len, k, mod);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			std::vector<std::pair<uint32_t, uint32_t>> concatenated;
+			concatenated.reserve(total_size);
+			for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		} else {
+			std::vector<std::vector<uint32_t>> output;
+			digest::thread_out::thread_mod<digest::BadCharPolicy::SKIPOVER>(
+				threads, output, seq, len, k, mod);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			
+			std::vector<uint32_t> concatenated;
+			concatenated.reserve(total_size);
+			for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		}
 	}
 }
 
 std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
-syncmer(const std::string &seq, unsigned k, unsigned large_window,
-		bool include_hash = false) {
-	digest::Syncmer<digest::BadCharPolicy::WRITEOVER, digest::ds::Adaptive>
-		digester(seq, k, large_window);
-	if (include_hash) {
-		std::vector<std::pair<uint32_t, uint32_t>> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
-	} else {
-		std::vector<uint32_t> output;
-		digester.roll_minimizer(seq.length(), output);
-		return output;
+syncmer_wrapper(const char *seq, size_t len, unsigned k, unsigned large_window,
+		unsigned threads, bool include_hash = false) {
+	if (threads == 1) {
+		digest::Syncmer<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>
+			digester(seq, len, k, large_window);
+		if (include_hash) {
+			std::vector<std::pair<uint32_t, uint32_t>> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		} else {
+			std::vector<uint32_t> output;
+			digester.roll_minimizer(len, output);
+			return output;
+		}
+	}
+	else {
+		if (include_hash) {
+			std::vector<std::vector<std::pair<uint32_t, uint32_t>>> output;
+			digest::thread_out::thread_sync<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>(
+				threads, output, seq, len, k, large_window);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			std::vector<std::pair<uint32_t, uint32_t>> concatenated;
+			concatenated.reserve(total_size);
+			for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		} else {
+			std::vector<std::vector<uint32_t>> output;
+			digest::thread_out::thread_sync<digest::BadCharPolicy::SKIPOVER, digest::ds::Adaptive>(
+				threads, output, seq, len, k, large_window);
+			
+			size_t total_size = 0;
+			for (const auto& vec : output) {
+				total_size += vec.size();
+			}
+			
+			std::vector<uint32_t> concatenated;
+			concatenated.reserve(total_size);
+			for (const auto& vec : output) {
+				concatenated.insert(concatenated.end(), vec.begin(), vec.end());
+			}
+			return concatenated;
+		}
 	}
-}
\ No newline at end of file
+}
+
+pybind11::array window_minimizer_numpy(const char *seq, size_t len, unsigned k, unsigned large_window,
+                                       unsigned threads, bool include_hash = false) {
+	py::gil_scoped_acquire gil;									
+    std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>> result;
+    {
+        py::gil_scoped_release release;
+        result = window_minimizer_wrapper(seq, len, k, large_window, threads, include_hash);
+    }
+    if (std::holds_alternative<std::vector<uint32_t>>(result)) {
+        const auto& vec = std::get<std::vector<uint32_t>>(result);
+        return pybind11::array_t<uint32_t>(vec.size(), vec.data());
+    } else {
+        const auto& vec = std::get<std::vector<std::pair<uint32_t, uint32_t>>>(result);
+        pybind11::array_t<uint32_t> arr({static_cast<pybind11::ssize_t>(vec.size()), static_cast<pybind11::ssize_t>(2)});
+        auto buf = arr.mutable_unchecked<2>();
+        for (pybind11::ssize_t i = 0; i < static_cast<pybind11::ssize_t>(vec.size()); ++i) {
+            buf(i, 0) = vec[i].first;
+            buf(i, 1) = vec[i].second;
+        }
+        return arr;
+    }
+}
+
+pybind11::array modimizer_numpy(const char *seq, size_t len, unsigned k, uint32_t mod,
+                                unsigned threads, bool include_hash = false) {
+	py::gil_scoped_acquire gil;
+    std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>> result;
+    {
+        py::gil_scoped_release release;
+        result = modimizer_wrapper(seq, len, k, mod, threads, include_hash);
+    }
+    if (std::holds_alternative<std::vector<uint32_t>>(result)) {
+        const auto& vec = std::get<std::vector<uint32_t>>(result);
+        return pybind11::array_t<uint32_t>(vec.size(), vec.data());
+    } else {
+        const auto& vec = std::get<std::vector<std::pair<uint32_t, uint32_t>>>(result);
+        pybind11::array_t<uint32_t> arr({static_cast<pybind11::ssize_t>(vec.size()), static_cast<pybind11::ssize_t>(2)});
+        auto buf = arr.mutable_unchecked<2>();
+        for (pybind11::ssize_t i = 0; i < static_cast<pybind11::ssize_t>(vec.size()); ++i) {
+            buf(i, 0) = vec[i].first;
+            buf(i, 1) = vec[i].second;
+        }
+        return arr;
+    }
+}
+
+pybind11::array syncmer_numpy(const char *seq, size_t len, unsigned k, unsigned large_window,
+                              unsigned threads, bool include_hash = false) {
+	py::gil_scoped_acquire gil;							
+    std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>> result;
+    {
+        py::gil_scoped_release release;
+        result = syncmer_wrapper(seq, len, k, large_window, threads, include_hash);
+    }
+    if (std::holds_alternative<std::vector<uint32_t>>(result)) {
+        const auto& vec = std::get<std::vector<uint32_t>>(result);
+        return pybind11::array_t<uint32_t>(vec.size(), vec.data());
+    } else {
+        const auto& vec = std::get<std::vector<std::pair<uint32_t, uint32_t>>>(result);
+        pybind11::array_t<uint32_t> arr({static_cast<pybind11::ssize_t>(vec.size()), static_cast<pybind11::ssize_t>(2)});
+        auto buf = arr.mutable_unchecked<2>();
+        for (pybind11::ssize_t i = 0; i < static_cast<pybind11::ssize_t>(vec.size()); ++i) {
+            buf(i, 0) = vec[i].first;
+            buf(i, 1) = vec[i].second;
+        }
+        return arr;
+    }
+}
+
+// std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
+// window_minimizer(py::bytes seq, unsigned k, unsigned large_window,
+//                     unsigned threads, bool include_hash = false) {
+//     py::gil_scoped_acquire gil;
+//     py::buffer_info info(py::buffer(seq).request());
+//     const char* data = static_cast<const char*>(info.ptr);
+//     size_t len = static_cast<size_t>(info.size);
+//     {
+//         py::gil_scoped_release release;
+//         return window_minimizer_wrapper(data, len, k, large_window, threads, include_hash);
+//     }
+// }
+
+// std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
+// modimizer(py::bytes seq, unsigned k, uint32_t mod,
+//              unsigned threads, bool include_hash = false) {
+//     py::gil_scoped_acquire gil;
+//     py::buffer_info info(py::buffer(seq).request());
+//     const char* data = static_cast<const char*>(info.ptr);
+//     size_t len = static_cast<size_t>(info.size);
+//     {
+//         py::gil_scoped_release release;
+//         return modimizer_wrapper(data, len, k, mod, threads, include_hash);
+//     }
+// }
+
+// std::variant<std::vector<uint32_t>, std::vector<std::pair<uint32_t, uint32_t>>>
+// syncmer(py::bytes seq, unsigned k, unsigned large_window,
+//            unsigned threads, bool include_hash = false) {
+//     py::gil_scoped_acquire gil;
+//     py::buffer_info info(py::buffer(seq).request());
+//     const char* data = static_cast<const char*>(info.ptr);
+//     size_t len = static_cast<size_t>(info.size);
+//     {
+//         py::gil_scoped_release release;
+//         return syncmer_wrapper(data, len, k, large_window, threads, include_hash);
+//     }
+// }