diff --git a/QCSim/MPSSimulatorTests.cpp b/QCSim/MPSSimulatorTests.cpp
index 434e15c..915cde7 100644
--- a/QCSim/MPSSimulatorTests.cpp
+++ b/QCSim/MPSSimulatorTests.cpp
@@ -8,8 +8,11 @@
 #include "QubitRegister.h"
 #include "MPSSimulator.h"
 
+#include <vector>
+
 #define _USE_MATH_DEFINES
 #include <math.h>
+#include <future>
 
 
 void FillOneQubitGates(std::vector<std::shared_ptr<QC::Gates::QuantumGateWithOp<>>>& gates)
@@ -527,29 +530,58 @@ bool TestMappedMeasurementsWithOneAndTwoQubitGatesCircuits()
 				std::unordered_map<std::vector<bool>, int> measurementsRegMap;
 				std::unordered_map<std::vector<bool>, int> measurementsMPSMap;
 
-				for (int t = 0; t < nrMeasurements; ++t)
-				{
-					QC::TensorNetworks::MPSSimulator mps(nrQubits);
-					QC::QubitRegister reg(nrQubits);
+				size_t remainingCounts = nrMeasurements;
+				size_t nrThreads = QC::QubitRegisterCalculator<>::GetNumberOfThreads();
+				nrThreads = std::min(nrThreads, std::max<size_t>(remainingCounts, 1ULL));
 
-					for (const auto& gate : circuit)
-					{
-						mps.ApplyGate(*gate);
-						reg.ApplyGate(*gate);
-					}
+				std::vector<std::future<void>> tasks(nrThreads);
 
-					std::vector<bool> measurementsReg(nrQubits);
-					std::vector<bool> measurementsMPS(nrQubits);
+				const size_t cntPerThread = static_cast<size_t>(ceil(static_cast<double>(remainingCounts) / nrThreads));
+				
+				std::mutex resultsMutex;
 
-					for (int q = 0; q < nrQubits; ++q)
-					{
-						measurementsReg[q] = reg.MeasureQubit(q);
-						measurementsMPS[q] = mps.MeasureQubit(q);
-					}
-					++measurementsRegMap[measurementsReg];
-					++measurementsMPSMap[measurementsMPS];
+				for (size_t i = 0; i < nrThreads; ++i)
+				{
+					const size_t curCnt = std::min(cntPerThread, remainingCounts);
+					remainingCounts -= curCnt;
+
+					tasks[i] = std::async(std::launch::async, [&circuit, &measurementsRegMap, &measurementsMPSMap, curCnt, nrQubits, &resultsMutex]()
+						{
+							QC::TensorNetworks::MPSSimulator mps(nrQubits);
+							QC::QubitRegister reg(nrQubits);
+
+							std::vector<bool> measurementsReg(nrQubits);
+							std::vector<bool> measurementsMPS(nrQubits);
+
+							for (size_t i = 0; i < curCnt; ++i)
+							{
+								for (const auto& gate : circuit)
+								{
+									mps.ApplyGate(*gate);
+									reg.ApplyGate(*gate);
+								}
+
+								for (int q = 0; q < nrQubits; ++q)
+								{
+									measurementsReg[q] = reg.MeasureQubit(q);
+									measurementsMPS[q] = mps.MeasureQubit(q);
+								}
+
+								{
+									const std::lock_guard lock(resultsMutex);
+									++measurementsRegMap[measurementsReg];
+									++measurementsMPSMap[measurementsMPS];
+								}
+
+								reg.setToBasisState(0);
+								mps.setToBasisState(0);
+							}
+						});
 				}
 
+				for (size_t i = 0; i < nrThreads; ++i)
+					tasks[i].get();
+
 				std::cout << ".";
 
 				for (const auto& [key, value] : measurementsRegMap)
diff --git a/QCSim/QubitRegister.h b/QCSim/QubitRegister.h
index 642b052..e5b044a 100644
--- a/QCSim/QubitRegister.h
+++ b/QCSim/QubitRegister.h
@@ -227,10 +227,8 @@ namespace QC {
 		void ApplyGate(const GateClass& gate, size_t qubit, size_t controllingQubit1 = 0, size_t controllingQubit2 = 0)
 		{
 			const size_t gateQubits = gate.getQubitsNumber();
-			if (qubit >= NrQubits ||
-				(gateQubits > 1 && controllingQubit1 >= NrQubits) ||
-				(gateQubits > 2 && controllingQubit2 >= NrQubits))
-				throw std::invalid_argument("Qubit number is too large");
+
+			CheckQubits(gate, qubit, controllingQubit1, controllingQubit2, gateQubits);
 
 #define OPTIMIZED_TENSOR_PRODUCT 1
 #ifdef OPTIMIZED_TENSOR_PRODUCT
@@ -468,6 +466,11 @@ namespace QC {
 		}
 
 	protected:
+		inline void CheckQubits(const GateClass& gate, size_t qubit, size_t controllingQubit1, size_t controllingQubit2, size_t gateQubits) const
+		{
+			if (NrQubits == 0) throw std::invalid_argument("Qubit number is zero");
+		}
+
 		// the following ones should be used for 'repeated measurements' that avoid reexecuting the circuit each time
 		size_t MeasureNoCollapse()
 		{