diff --git a/.github/workflows/compile.yml b/.github/workflows/compile.yml
index fa7ca423..0950b7db 100644
--- a/.github/workflows/compile.yml
+++ b/.github/workflows/compile.yml
@@ -211,13 +211,28 @@ jobs:
       - name: Install ROCm
         if: ${{ matrix.hip == 'ON' }}
         run: |
+          sudo apt update
+          sudo apt upgrade
           sudo apt install "linux-headers-$(uname -r)" "linux-modules-extra-$(uname -r)"
-          sudo apt install python3-setuptools python3-wheel
-          sudo usermod -a -G render,video $USER
-          wget https://repo.radeon.com/amdgpu-install/6.3.3/ubuntu/noble/amdgpu-install_6.3.60303-1_all.deb
-          sudo apt install ./amdgpu-install_6.3.60303-1_all.deb
+          # Make the directory if it doesn't exist yet.
+          # This location is recommended by the distribution maintainers.
+          sudo mkdir --parents --mode=0755 /etc/apt/keyrings
+          # Download the key, convert the signing-key to a full
+          # keyring required by apt and store in the keyring directory
+          wget https://repo.radeon.com/rocm/rocm.gpg.key -O - | \
+          gpg --dearmor | sudo tee /etc/apt/keyrings/rocm.gpg > /dev/null
+          sudo tee /etc/apt/sources.list.d/amdgpu.list << EOF
+          deb [arch=amd64 signed-by=/etc/apt/keyrings/rocm.gpg] https://repo.radeon.com/rocm/apt/7.2 noble main
+          deb [arch=amd64 signed-by=/etc/apt/keyrings/rocm.gpg] https://repo.radeon.com/graphics/7.2/ubuntu noble main
+          EOF
+          sudo tee /etc/apt/preferences.d/rocm-pin-600 << EOF
+          Package: *
+          Pin: release o=repo.radeon.com
+          Pin-Priority: 600
+          EOF
           sudo apt update
-          sudo apt install amdgpu-dkms rocm
+          sudo apt autoremove
+          sudo apt install rocm rocm-hip-sdk rocm-hip-runtime-dev
           echo "${{ env.rocm_path }}" >> $GITHUB_PATH
 
       # invoke cmake, disabling LTO (it duplicates symbols with CUDA + MPI)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index eb826e18..5d308795 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -40,8 +40,6 @@ set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
 
 
 # GNUInstallDirs to provide sensible default install directory names
-cmake_path(SET ORG_INSTALL_PATH NORMALIZE "${CMAKE_INSTALL_PREFIX}")
-cmake_path(APPEND CMAKE_INSTALL_PREFIX "quest")
 include(GNUInstallDirs)
 include(CMakePackageConfigHelpers)
 
@@ -74,10 +72,11 @@ endif()
 
 
 # Library type
-# Shared library by default
-option(BUILD_SHARED_LIBS "Build shared library. Turned ON by default." ON)
-message(STATUS "Shared library is turned ${BUILD_SHARED_LIBS}. Set BUILD_SHARED_LIBS to modify.")
-
+# Shared library by default when we are the top level cmake project
+if(PROJECT_IS_TOP_LEVEL)
+  option(BUILD_SHARED_LIBS "Build shared library. Turned ON by default." ON)
+  message(STATUS "Shared library is turned ${BUILD_SHARED_LIBS}. Set BUILD_SHARED_LIBS to modify.")
+endif ()
 
 # Library naming
 set(LIB_NAME QuEST 
@@ -185,6 +184,7 @@ option(
 )
 message(STATUS "Disabling of deprecated API warnings is turned ${DISABLE_DEPRECATION_WARNINGS}. Set DISABLE_DEPRECATION_WARNINGS to modify.")
 
+option(INSTALL_BINARIES "Whether to include example and user binaries in the install." OFF)
 
 
 # ============================
@@ -565,10 +565,12 @@ add_executable(min_example
 )
 target_link_libraries(min_example PRIVATE QuEST::QuEST)
 
-install(TARGETS min_example
-  RUNTIME
-  DESTINATION ${CMAKE_INSTALL_BINDIR}
-)
+if (INSTALL_BINARIES)
+  install(TARGETS min_example
+    RUNTIME
+    DESTINATION ${CMAKE_INSTALL_BINDIR}
+  )
+endif ()
 
 
 # all examples optionally built
@@ -577,7 +579,7 @@ if (BUILD_EXAMPLES)
 endif()
 
 
-## RATH
+## RPATH
 set(BUILD_RPATH_USE_ORIGIN ON)
 if(APPLE)
   set(_RPATH_ORIGIN "@loader_path")
@@ -623,7 +625,11 @@ if (USER_SOURCE AND OUTPUT_EXE)
 
   add_executable(${OUTPUT_EXE} ${USER_SOURCE})
   target_link_libraries(${OUTPUT_EXE} PUBLIC QuEST)
-  install(TARGETS ${OUTPUT_EXE} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+   
+  if (INSTALL_BINARIES)
+    install(TARGETS ${OUTPUT_EXE} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+  endif()
+  
   setup_quest_rpath(${OUTPUT_EXE})
 endif()
 
@@ -732,3 +738,7 @@ install(
   NAMESPACE QuEST::
   DESTINATION "${QuEST_INSTALL_CONFIGDIR}"
 )
+
+if(PROJECT_IS_TOP_LEVEL)
+  include(CPack)
+endif ()
\ No newline at end of file
diff --git a/README.md b/README.md
index 9eb85bfd..f1dd5e87 100644
--- a/README.md
+++ b/README.md
@@ -94,7 +94,7 @@ initRandomPureState(qureg);
 
 applyHadamard(qureg, 0);
 applyControlledRotateX(qureg, 0, 1, angle);
-applyFullQuantumFourierTransform(qureg);
+applyFullQuantumFourierTransform(qureg, inverse);
 
 reportQureg(qureg);
 
diff --git a/docs/cmake.md b/docs/cmake.md
index 9240b679..d3c23ee4 100644
--- a/docs/cmake.md
+++ b/docs/cmake.md
@@ -36,8 +36,9 @@ make
 | `VERBOSE_LIB_NAME` | (`OFF`), `ON` | When turned on `LIB_NAME` will be modified according to the other configuration options chosen. For example compiling QuEST with multithreading, distribution, and double precision with `VERBOSE_LIB_NAME` turned on creates `libQuEST-fp2+mt+mpi.so`. |
 | `FLOAT_PRECISION` | (`2`), `1`, `4` | Determines which floating-point precision QuEST will use: double, single, or quad. *Note: Quad precision is not supported when also compiling for GPU.* |
 | `BUILD_EXAMPLES` | (`OFF`), `ON` | Determines whether the example programs will be built alongside QuEST. Note that `min_example` is always built. |
-| `ENABLE_MULTITHREADING` | (`ON`), OFF | Determines whether QuEST will be built with support for parallelisation with OpenMP. |
-| `ENABLE_DISTRIBUTION` | (`OFF`), ON | Determines whether QuEST will be built with support for parallelisation with MPI. |
+| `INSTALL_BINARIES` | (`OFF`), `ON` | Determines whether compiled binaries such as the examples will be installed as well as the QuEST library. |
+| `ENABLE_MULTITHREADING` | (`ON`), `OFF` | Determines whether QuEST will be built with support for parallelisation with OpenMP. |
+| `ENABLE_DISTRIBUTION` | (`OFF`), `ON` | Determines whether QuEST will be built with support for parallelisation with MPI. |
 | `ENABLE_CUDA` | (`OFF`), `ON` | Determines whether QuEST will be built with support for NVIDIA GPU acceleration. If turned on, `CMAKE_CUDA_ARCHITECTURES` should probably also be set. |
 | `ENABLE_CUQUANTUM` | (`OFF`), `ON` | Determines whether QuEST will make use of the NVIDIA CuQuantum library. Cannot be turned on if `ENABLE_CUDA` is off. |
 | `ENABLE_HIP` | (`OFF`), `ON` | Determines whether QuEST will be built with support for AMD GPU acceleration. If turned on, `CMAKE_HIP_ARCHITECTURES` should probably also be set. |
diff --git a/docs/tutorial.md b/docs/tutorial.md
index a006a3fc..b3e99706 100644
--- a/docs/tutorial.md
+++ b/docs/tutorial.md
@@ -626,7 +626,8 @@ applyCompMatr1(qureg, 0, m);
 QuEST includes a few convenience functions for effecting [QFT](https://quest-kit.github.io/QuEST/group__op__qft.html) and [Trotter](https://quest-kit.github.io/QuEST/group__op__paulistrsum.html) circuits.
 
 ```cpp
-applyQuantumFourierTransform(qureg, targets, 3);
+bool inverse = false;
+applyQuantumFourierTransform(qureg, targets, 3, inverse);
 
 qreal time = .3;
 int order = 4;
@@ -865,4 +866,4 @@ This is important because it ensures:
 > [!CAUTION]
 > After calling `finalizeQuESTEnv()`, MPI will close and if being accessed directly by the user, will enter an undefined state. Subsequent calls to MPI routines may return gibberish, and distributed machines will lose their ability to communicate. It is recommended to call `finalizeQuESTEnv()` immediately before exiting.
 
-You are now a QuEST expert 🎉 though there are _many_ more functions in the [API](https://quest-kit.github.io/QuEST/group__api.html) not covered here. Go forth and simulate!
\ No newline at end of file
+You are now a QuEST expert 🎉 though there are _many_ more functions in the [API](https://quest-kit.github.io/QuEST/group__api.html) not covered here. Go forth and simulate!
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 00d4c343..afc8f85d 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -20,11 +20,13 @@ function(add_example direc in_fn)
   add_executable(${target} ${in_fn})
   target_link_libraries(${target} PUBLIC QuEST)
 
-  install(
-    TARGETS ${target}
-    RUNTIME
-    DESTINATION ${out_dir}
-  )
+  if (INSTALL_BINARIES)
+    install(
+      TARGETS ${target}
+      RUNTIME
+      DESTINATION ${out_dir}
+    )
+  endif ()
 
   set_target_properties(${target} 
     PROPERTIES 
diff --git a/examples/extended/dynamics.c b/examples/extended/dynamics.c
index 0d6763c1..03abcc16 100644
--- a/examples/extended/dynamics.c
+++ b/examples/extended/dynamics.c
@@ -153,11 +153,12 @@ int main() {
     int order = 4;
     int reps  = 5;
     int steps = 20;
+    bool randomise = true;
 
     for (int i=0; i<steps; i++) {
 
         // evolve qureg under (approx) exp(-i dt H)
-        applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps);
+        applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps, randomise);
 
         // calculate and report <O>
         qreal time = dt * (i+1);
@@ -187,7 +188,7 @@ int main() {
 
     // verify results by uninterrupted higher-order simulation to target time
     initPlusState(qureg);
-    applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps);
+    applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps, randomise);
     reportScalar("final <O>", calcExpecPauliStrSum(qureg, observ));
 
     // clean up
diff --git a/examples/extended/dynamics.cpp b/examples/extended/dynamics.cpp
index 95245fb8..63614538 100644
--- a/examples/extended/dynamics.cpp
+++ b/examples/extended/dynamics.cpp
@@ -150,11 +150,12 @@ int main() {
     int order = 4;
     int reps  = 5;
     int steps = 20;
+    bool randomise = true;
 
     for (int i=0; i<steps; i++) {
 
         // evolve qureg under (approx) exp(-i dt H)
-        applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps);
+        applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps, randomise);
 
         // calculate and report <O>
         qreal time = dt * (i+1);
@@ -181,7 +182,7 @@ int main() {
 
     // verify results by uninterrupted higher-order simulation to target time
     initPlusState(qureg);
-    applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps);
+    applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps, randomise);
     reportScalar("final <O>", calcExpecPauliStrSum(qureg, observ));
 
     // clean up
diff --git a/quest/include/deprecated.h b/quest/include/deprecated.h
index 1bc8bda5..1a63b204 100644
--- a/quest/include/deprecated.h
+++ b/quest/include/deprecated.h
@@ -21,6 +21,7 @@
 #include "quest.h"
 
 #include "stdlib.h"
+#include <stdbool.h>
 
 
 
@@ -1321,11 +1322,11 @@ static inline void _multiControlledMultiRotatePauli(Qureg qureg, int* ctrls, int
 
 #define applyFullQFT(...) \
     _WARN_FUNC_RENAMED("applyFullQFT()", "applyFullQuantumFourierTransform()") \
-    applyFullQuantumFourierTransform(__VA_ARGS__)
+    applyFullQuantumFourierTransform(__VA_ARGS__, false)
 
 #define applyQFT(...) \
     _WARN_FUNC_RENAMED("applyQFT()", "applyQuantumFourierTransform()") \
-    applyQuantumFourierTransform(__VA_ARGS__)
+    applyQuantumFourierTransform(__VA_ARGS__, false)
 
 
 
diff --git a/quest/include/multiplication.h b/quest/include/multiplication.h
index 81d27a50..6422e69f 100644
--- a/quest/include/multiplication.h
+++ b/quest/include/multiplication.h
@@ -746,7 +746,6 @@ extern "C" {
 
 
 /// @notyetdoced
-/// @notyetvalidated
 /// @see
 /// - leftapplyCompMatr1()
 /// - applyQubitProjector()
@@ -754,7 +753,6 @@ void leftapplyQubitProjector(Qureg qureg, int qubit, int outcome);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 /// @see
 /// - leftapplyCompMatr1()
 /// - applyMultiQubitProjector()
@@ -762,7 +760,6 @@ void leftapplyMultiQubitProjector(Qureg qureg, int* qubits, int* outcomes, int n
 
 
 /// @notyetdoced
-/// @notyetvalidated
 /// @see
 /// - rightapplyCompMatr1()
 /// - applyQubitProjector()
@@ -782,6 +779,25 @@ void rightapplyMultiQubitProjector(Qureg qureg, int* qubits, int* outcomes, int
 }
 #endif
 
+#ifdef __cplusplus
+
+
+/// @notyetdoced
+/// @cppvectoroverload
+/// @see leftapplyMultiQubitProjector()
+void leftapplyMultiQubitProjector(Qureg qureg, std::vector<int> qubits, std::vector<int> outcomes);
+
+
+/// @notyetdoced
+/// @cppvectoroverload
+/// @see rightapplyMultiQubitProjector()
+void rightapplyMultiQubitProjector(Qureg qureg, std::vector<int> qubits, std::vector<int> outcomes);
+
+
+#endif
+
+/** @} */
+
 
 
 /** 
diff --git a/quest/include/operations.h b/quest/include/operations.h
index b138f200..8cacdf44 100644
--- a/quest/include/operations.h
+++ b/quest/include/operations.h
@@ -25,6 +25,8 @@
 #include "quest/include/matrices.h"
 #include "quest/include/channels.h"
 
+#include <stdbool.h>
+
 #ifdef __cplusplus
     #include <vector>
 #endif
@@ -2291,32 +2293,26 @@ extern "C" {
 
 
 /// @notyetdoced
-/// @notyetvalidated
 int applyQubitMeasurement(Qureg qureg, int target);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 int applyQubitMeasurementAndGetProb(Qureg qureg, int target, qreal* probability);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 qreal applyForcedQubitMeasurement(Qureg qureg, int target, int outcome);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 qindex applyMultiQubitMeasurement(Qureg qureg, int* qubits, int numQubits);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 qindex applyMultiQubitMeasurementAndGetProb(Qureg qureg, int* qubits, int numQubits, qreal* probability);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 qreal applyForcedMultiQubitMeasurement(Qureg qureg, int* qubits, int* outcomes, int numQubits);
 
 
@@ -2328,16 +2324,18 @@ qreal applyForcedMultiQubitMeasurement(Qureg qureg, int* qubits, int* outcomes,
 #ifdef __cplusplus
 
 
-/// @notyettested
-/// @notyetvalidated
+/// @notyetdoced
+/// @cppvectoroverload
+/// @see applyMultiQubitMeasurement()
+qindex applyMultiQubitMeasurement(Qureg qureg, std::vector<int> qubits);
+
+
 /// @notyetdoced
 /// @cppvectoroverload
 /// @see applyMultiQubitMeasurementAndGetProb()
 qindex applyMultiQubitMeasurementAndGetProb(Qureg qureg, std::vector<int> qubits, qreal* probability);
 
 
-/// @notyettested
-/// @notyetvalidated
 /// @notyetdoced
 /// @cppvectoroverload
 /// @see applyForcedMultiQubitMeasurement()
@@ -2363,12 +2361,10 @@ extern "C" {
 
 
 /// @notyetdoced
-/// @notyetvalidated
 void applyQubitProjector(Qureg qureg, int target, int outcome);
 
 
 /// @notyetdoced
-/// @notyetvalidated
 void applyMultiQubitProjector(Qureg qureg, int* qubits, int* outcomes, int numQubits);
 
 
@@ -2380,8 +2376,6 @@ void applyMultiQubitProjector(Qureg qureg, int* qubits, int* outcomes, int numQu
 #ifdef __cplusplus
 
 
-/// @notyettested
-/// @notyetvalidated
 /// @notyetdoced
 /// @cppvectoroverload
 /// @see applyMultiQubitProjector()
@@ -2406,14 +2400,67 @@ extern "C" {
 #endif
 
 
-/// @notyetdoced
-/// @notyetvalidated
-void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets);
+/** 
+ * Applies the Quantum Fourier Transform upon the specified @p targets of @p qureg.
+ * Alternatively, applies the Inverse Quantum Fourier Transform according to @p inverse.
+ * 
+ * @formulae
+ * 
+ * Letting @f$ N @f$ = @p numTargets, the @f$ N @f$ qubit Quantum Fourier Transform maps each
+ * computational basis state of the targeted qubits, @f$ \ket{j} @f$, according to
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ * Similarly the Inverse Quantum Fourier Transform maps each basis state like
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{-2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ *
+ * @param[in,out] qureg      the state to modify.
+ * @param[in]     targets    the indices of the target qubits.
+ * @param[in]     numTargets the length of list @p targets
+ * @param[in]     inverse    whether to apply the inverse QFT or forward QFT
+ * @throws @validationerror
+ * - if @p qureg is uninitialised.
+*  - if @p targets are invalid qubit indices.
+*  - if @p targets are not unique.
+ * - if @p numTargets < 1.
+ * @see
+ * - applyFullQuantumFourierTransform()
+ * @author Vasco Ferreira
+ */
+void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets, bool inverse);
 
 
-/// @notyetdoced
-/// @notyetvalidated
-void applyFullQuantumFourierTransform(Qureg qureg);
+/** 
+ * Applies the Quantum Fourier Transform upon all qubits in @p qureg. Alternatively,
+ * applies the Inverse Quantum Fourier Transform according to @p inverse.
+ * 
+ * @formulae
+ * 
+ * The Quantum Fourier Transform maps each computational basis state @f$ \ket{j} @f$
+ * in an @f$ N @f$ qubit @p qureg according to
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ * Similarly the Inverse Quantum Fourier Transform maps each basis state like
+ * @f[ 
+        \ket{j} \rightarrow \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} e^{-2 \pi i j k / 2^N} \ket{k}.
+ * @f]
+ *
+ * @equivalences
+ *
+ * - This function wraps applyQuantumFourierTransform(), passing all qubits in the @p qureg as targets.
+ *
+ * @param[in,out] qureg      the state to modify.
+ * @param[in]     inverse    whether to apply the inverse QFT or forward QFT
+ * @throws @validationerror
+ * - if @p qureg is uninitialised.
+ * @see
+ * - applyQuantumFourierTransform()
+ * @author Vasco Ferreira
+ */
+void applyFullQuantumFourierTransform(Qureg qureg, bool inverse);
 
 
 // end de-mangler
@@ -2424,12 +2471,10 @@ void applyFullQuantumFourierTransform(Qureg qureg);
 #ifdef __cplusplus
 
 
-/// @notyettested
-/// @notyetvalidated
 /// @notyetdoced
 /// @cppvectoroverload
 /// @see applyQuantumFourierTransform()
-void applyQuantumFourierTransform(Qureg qureg, std::vector<int> targets);
+void applyQuantumFourierTransform(Qureg qureg, std::vector<int> targets, bool inverse);
 
 
 #endif // __cplusplus
diff --git a/quest/include/paulis.h b/quest/include/paulis.h
index 1d08169f..e4645923 100644
--- a/quest/include/paulis.h
+++ b/quest/include/paulis.h
@@ -99,6 +99,9 @@ typedef struct {
  * 
  * @defgroup paulis_reporters Reporters
  * @brief Functions for printing Pauli data structures.
+ *
+ * @defgroup paulis_setters Setters
+ * @brief Functions for overwriting the elements of Pauli data structures.
  */
 
 
@@ -420,6 +423,91 @@ extern "C" {
 
 
 
+/*
+ * SETTERS
+ */
+
+// enable invocation by both C and C++ binaries
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+    /** @ingroup paulis_setters
+     *
+     * Reorders the terms within a @p sum of weighted Pauli strings to sort Pauli
+     * strings into lexicographic (dictionary) ordering.
+     *
+     * @formulae
+     * Let @f$ H = @f$ @p sum, which can be represented as
+     * @f[
+          H = \sum\limits_j c_j \, \hat{\sigma}_j
+     * @f]
+     * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$.
+     *
+     * This function constructs and applies the permutation @f$ \pi @f$ to @f$ H @f$
+     * @f[
+          H = \sum\limits_j c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
+     * @f]
+     * such that
+     * @f[
+     *   \hat{\sigma}_{\pi(i)} <_{lex} \hat{\sigma}_{\pi(j)} \  \forall \  \pi(i) < \pi(j).
+     * @f]
+     *
+     *
+     * @param[in,out] sum  a weighted sum of Pauli strings to reorder.
+     *
+     * @throws @validationerror
+     * - if @p sum is not initialised.
+     *
+     * @see
+     * - sortPauliStrSumMagnitude()
+     * @author Vasco Ferreira
+     */
+    void sortPauliStrSumLexicographic(PauliStrSum sum);
+
+
+    /** @ingroup paulis_setters
+     *
+     * Reorders the terms within a @p sum of weighted Pauli strings to sort Pauli
+     * strings into decreasing magnitude weights.
+     *
+     * @formulae
+     * Let @f$ H = @f$ @p sum, represented as the weighted sum
+     * @f[
+          H = \sum\limits_j c_j \, \hat{\sigma}_j
+     * @f]
+     * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$.
+     *
+     * This function constructs and applies the permutation @f$ \pi @f$ to @f$ H @f$
+     * @f[
+          H = \sum\limits_j c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
+     * @f]
+     * such that
+     * @f[
+     *    |c_{\pi(i)}| > |c_{\pi(j)}| \, \forall \,  \pi(i) < \pi(j).
+     * @f]
+     *
+     * @param[in,out] sum  a weighted sum of Pauli strings to reorder.
+     *
+     * @throws @validationerror
+     * - if @p sum is not initialised.
+     *
+     * @see
+     * - sortPauliStrSumLexicographic()
+     *
+     * @author Vasco Ferreira
+     */
+    void sortPauliStrSumMagnitude(PauliStrSum sum);
+
+
+// end de-mangler
+#ifdef __cplusplus
+}
+#endif
+
+
+
 #endif // PAULIS_H
 
 /** @} */ // (end file-wide doxygen defgroup)
diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h
index 667eedd9..9e7cad25 100644
--- a/quest/include/trotterisation.h
+++ b/quest/include/trotterisation.h
@@ -10,9 +10,12 @@
  * @{
  */
 
+
 #ifndef TROTTERISATION_H
 #define TROTTERISATION_H
 
+#include <stdbool.h>
+
 #include "quest/include/qureg.h"
 #include "quest/include/paulis.h"
 #include "quest/include/matrices.h"
@@ -40,9 +43,12 @@ extern "C" {
  * 
  * Effects an approximation to the exponential of @p sum, weighted by @p angle times @f$ i @f$, upon @p qureg,
  * via the symmetrized Trotter-Suzuki decomposition (<a href="https://arxiv.org/abs/math-ph/0506007">arXiv</a>).
+ * 
  * Increasing @p reps (the number of Trotter repetitions) or @p order (an even, positive integer or one) 
  * improves the accuracy of the approximation by reducing the "Trotter error" due to non-commuting 
- * terms of @p sum, though increases the runtime linearly and exponentially respectively.
+ * terms of @p sum, though increases the runtime linearly and exponentially respectively. 
+ * Using @p permutePaulis the ordering of terms in the sum can also be randomised, which generally 
+ * improves the accuracy of the approximation for low order decompositions (<a href="https://arxiv.org/abs/1805.08385">arXiv</a>).
  * 
  * @formulae 
  * 
@@ -103,6 +109,15 @@ extern "C" {
  * > These formulations are taken from 'Finding Exponential Product Formulas
  * > of Higher Orders', Naomichi Hatano and Masuo Suzuki (2005) (<a href="https://arxiv.org/abs/math-ph/0506007">arXiv</a>).
  * 
+ * When @p permutePaulis=true the terms of @p sum are effected in a random order at each repetition. That is, each repetition of the Trotter-Suzuki decomposition is evaluated with the sum
+ * @f[
+      \hat{H} = \sum\limits_j^T c_{\pi(j)} \, \hat{\sigma}_{\pi(j)}
+ * @f]
+ * where @f$ \pi @f$  is a randomly selected permutation.
+ *
+ * @important
+ *   Using @p permutePaulis=true will cause @p sum to be mutated by the Trotterisation.
+ *
  * @equivalences
  * 
  * - By passing @f$ \theta = - \Delta t / \hbar @f$, this function approximates unitary time evolution of a closed 
@@ -137,11 +152,12 @@ extern "C" {
  * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly
  *   when all PauliStr in @p sum = @f$ \hat{H} @f$ commute, or @p reps @f$ \rightarrow \infty @f$.
  * 
- * @param[in,out] qureg  the state to modify.
- * @param[in]     sum    a weighted sum of Pauli strings to approximately exponentiate.
- * @param[in]     angle  the prefactor of @p sum times @f$ i @f$ in the exponent.
- * @param[in]     order  the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
- * @param[in]     reps   the number of Trotter repetitions.
+ * @param[in,out] qureg                the state to modify.
+ * @param[in]     sum                  a weighted sum of Pauli strings to approximately exponentiate.
+ * @param[in]     angle                the prefactor of @p sum times @f$ i @f$ in the exponent.
+ * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
+ * @param[in]     reps                 the number of Trotter repetitions.
+ * @param[in]     permutePaulis        whether to randomly reorder Pauli strings at each repetition.
  * 
  * @throws @validationerror
  * - if @p qureg or @p sum are uninitialised.
@@ -156,32 +172,45 @@ extern "C" {
  *  - applyTrotterizedUnitaryTimeEvolution()
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (randomisation)
  */
-void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 /// @notyetdoced
 /// @notyettested
+///
+/// @author Tyson Jones
+/// @author Vasco Ferreira (randomisation)
+///
 /// @see
 ///  - applyTrotterizedPauliStrSumGadget()
 ///  - applyControlledCompMatr1()
-void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 /// @notyetdoced
 /// @notyettested
+///
+/// @author Tyson Jones
+/// @author Vasco Ferreira (randomisation)
+///
 /// @see
 ///  - applyTrotterizedPauliStrSumGadget()
 ///  - applyMultiControlledCompMatr1()
-void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 /// @notyetdoced
 /// @notyettested
+///
+/// @author Tyson Jones
+/// @author Vasco Ferreira (randomisation)
+///
 /// @see
 ///  - applyTrotterizedPauliStrSumGadget()
 ///  - applyMultiStateControlledCompMatr1()
-void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 /** @notyettested
@@ -218,11 +247,12 @@ void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* con
  * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly
  *   when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. 
  * 
- * @param[in,out] qureg  the state to modify.
- * @param[in]     sum    a weighted sum of Pauli strings to approximately exponentiate.
- * @param[in]     angle  an effective prefactor of @p sum in the exponent.
- * @param[in]     order  the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
- * @param[in]     reps   the number of Trotter repetitions.
+ * @param[in,out] qureg                the state to modify.
+ * @param[in]     sum                  a weighted sum of Pauli strings to approximately exponentiate.
+ * @param[in]     angle                an effective prefactor of @p sum in the exponent.
+ * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
+ * @param[in]     reps                 the number of Trotter repetitions.
+ * @param[in]     permutePaulis        whether to randomly reorder Pauli strings at each repetition.
  * 
  * @throws @validationerror
  * - if @p qureg or @p sum are uninitialised.
@@ -231,8 +261,9 @@ void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* con
  * - if @p reps is not a positive integer.
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (randomisation)
  */
-void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps);
+void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps, bool permutePaulis);
 
 
 // end de-mangler
@@ -247,16 +278,24 @@ void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, q
 /// @notyetvalidated
 /// @notyetdoced
 /// @cppvectoroverload
+///
+/// @author Tyson Jones
+/// @author Vasco Ferreira (randomisation)
+///
 /// @see applyTrotterizedMultiControlledPauliStrSumGadget()
-void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, std::vector<int> controls, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, std::vector<int> controls, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 /// @notyettested
 /// @notyetvalidated
 /// @notyetdoced
 /// @cppvectoroverload
+///
+/// @author Tyson Jones
+/// @author Vasco Ferreira (randomisation)
+///
 /// @see applyTrotterizedMultiStateControlledPauliStrSumGadget()
-void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, std::vector<int> controls, std::vector<int> states, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, std::vector<int> controls, std::vector<int> states, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
 
 #endif // __cplusplus
@@ -277,9 +316,7 @@ extern "C" {
 #endif
 
 
-/** @notyettested
- * 
- * Unitarily time evolves @p qureg for the duration @p time under the time-independent Hamiltonian @p hamil, 
+/** Unitarily time evolves @p qureg for the duration @p time under the time-independent Hamiltonian @p hamil, 
  * as approximated by symmetrized Trotterisation of the specified @p order and number of cycles @p reps. 
  * 
  * @formulae 
@@ -353,11 +390,12 @@ extern "C" {
  *  - applyTrotterizedNoisyTimeEvolution()
  *  - applyTrotterizedNonUnitaryPauliStrSumGadget()
  * 
- * @param[in,out] qureg  the state to modify.
- * @param[in]     hamil  the Hamiltonian as a a weighted sum of Pauli strings.
- * @param[in]     time   the duration over which to simulate evolution.
- * @param[in]     order  the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
- * @param[in]     reps   the number of Trotter repetitions.
+ * @param[in,out] qureg                the state to modify.
+ * @param[in]     hamil                the Hamiltonian as a a weighted sum of Pauli strings.
+ * @param[in]     time                 the duration over which to simulate evolution.
+ * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
+ * @param[in]     reps                 the number of Trotter repetitions.
+ * @param[in]     permutePaulis        whether to randomly reorder Pauli strings at each repetition.
  * 
  * @throws @validationerror
  * - if @p qureg or @p hamil are uninitialised.
@@ -367,13 +405,12 @@ extern "C" {
  * - if @p reps is not a positive integer.
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (randomisation)
  */
-void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps);
+void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps, bool permutePaulis);
 
 
-/** @notyettested
- * 
- * Simulates imaginary-time evolution of @p qureg for the duration @p tau under the time-independent 
+/** Simulates imaginary-time evolution of @p qureg for the duration @p tau under the time-independent 
  * Hamiltonian @p hamil, as approximated by symmetrized Trotterisation of the specified @p order and
  * number of cycles @p reps. 
  * 
@@ -484,11 +521,12 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal
  *  - applyTrotterizedUnitaryTimeEvolution()
  *  - applyTrotterizedNonUnitaryPauliStrSumGadget()
  * 
- * @param[in,out] qureg  the state to modify.
- * @param[in]     hamil  the Hamiltonian as a a weighted sum of Pauli strings.
- * @param[in]     tau    the duration over which to simulate imaginary-time evolution.
- * @param[in]     order  the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
- * @param[in]     reps   the number of Trotter repetitions.
+ * @param[in,out] qureg                the state to modify.
+ * @param[in]     hamil                the Hamiltonian as a a weighted sum of Pauli strings.
+ * @param[in]     tau                  the duration over which to simulate imaginary-time evolution.
+ * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
+ * @param[in]     reps                 the number of Trotter repetitions.
+ * @param[in]     permutePaulis        whether to randomly reorder Pauli strings at each repetition.
  * 
  * @throws @validationerror
  * - if @p qureg or @p hamil are uninitialised.
@@ -498,8 +536,9 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal
  * - if @p reps is not a positive integer.
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (randomisation)
  */
-void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps);
+void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps, bool permutePaulis);
 
 
 /** @notyettested
@@ -622,14 +661,15 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea
  *  - applyTrotterizedUnitaryTimeEvolution()
  *  - applyTrotterizedImaginaryTimeEvolution()
  * 
- * @param[in,out] qureg     the density-matrix state to evolve and modify.
- * @param[in]     hamil     the Hamiltonian of the qubit system (excludes any environment).
- * @param[in]     damps     the damping rates of each jump operator in @p jumps.
- * @param[in]     jumps     the jump operators specified as PauliStrSum.
- * @param[in]     numJumps  the length of list @p jumps (and @p damps).
- * @param[in]     time      the duration through which to evolve the state.
- * @param[in]     order     the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
- * @param[in]     reps      the number of Trotter repetitions.
+ * @param[in,out] qureg                the density-matrix state to evolve and modify.
+ * @param[in]     hamil                the Hamiltonian of the qubit system (excludes any environment).
+ * @param[in]     damps                the damping rates of each jump operator in @p jumps.
+ * @param[in]     jumps                the jump operators specified as PauliStrSum.
+ * @param[in]     numJumps             the length of list @p jumps (and @p damps).
+ * @param[in]     time                 the duration through which to evolve the state.
+ * @param[in]     order                the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...).
+ * @param[in]     reps                 the number of Trotter repetitions.
+ * @param[in]     permutePaulis        whether to randomly reorder Pauli strings at each repetition.
  * 
  * @throws @validationerror
  * - if @p qureg, @p hamil or any element of @p jumps are uninitialised.
@@ -645,8 +685,9 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea
  * - if @p reps is not a positive integer.
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (randomisation)
  */
-void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps);
+void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps, bool permutePaulis);
 
 
 // end de-mangler
diff --git a/quest/src/api/operations.cpp b/quest/src/api/operations.cpp
index e458605a..88141e23 100644
--- a/quest/src/api/operations.cpp
+++ b/quest/src/api/operations.cpp
@@ -1717,6 +1717,11 @@ qreal applyForcedMultiQubitMeasurement(Qureg qureg, int* qubits, int* outcomes,
 
 } // end de-mangler
 
+qindex applyMultiQubitMeasurement(Qureg qureg, vector<int> qubits) {
+
+    return applyMultiQubitMeasurement(qureg, qubits.data(), qubits.size());
+}
+
 qindex applyMultiQubitMeasurementAndGetProb(Qureg qureg, vector<int> qubits, qreal* probability) {
 
     return applyMultiQubitMeasurementAndGetProb(qureg, qubits.data(), qubits.size(), probability);
@@ -1736,7 +1741,7 @@ qreal applyForcedMultiQubitMeasurement(Qureg qureg, vector<int> qubits, vector<i
 
 extern "C" {
 
-void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets) {
+void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets, bool inverse) {
     validate_quregFields(qureg, __func__);
     validate_targets(qureg, targets, numTargets, __func__);
 
@@ -1744,20 +1749,37 @@ void applyQuantumFourierTransform(Qureg qureg, int* targets, int numTargets) {
     /// change this placeholder implementation to the bespoke, optimised routine,
     /// wherein each contiguous controlled-phase gate is merged
 
-    for (int n=numTargets-1; n>=0; n--) {
-        applyHadamard(qureg, targets[n]);
-        for (int m=0; m<n; m++) {
-            qreal arg = const_PI / powerOf2(m+1);
-            applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+    int mid = numTargets/2; // floors
+
+    // don't think it is worth overcomplicating this to avoid slight repetition
+    if (inverse) {
+        for (int n=0; n<mid; n++)
+            applySwap(qureg, targets[n], targets[numTargets-1-n]);
+
+        for (int n=0; n<numTargets; n++) {
+            for (int m=0; m<n; m++) {
+                qreal arg = - const_PI / powerOf2(m+1);
+                applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+            }
+            applyHadamard(qureg, targets[n]);
         }
-    }
 
-    int mid = numTargets/2; // floors
-    for (int n=0; n<mid; n++)
-        applySwap(qureg, targets[n], targets[numTargets-1-n]);
+    } else {
+        for (int n=numTargets-1; n>=0; n--) {
+            applyHadamard(qureg, targets[n]);
+
+            for (int m=0; m<n; m++) {
+                qreal arg = const_PI / powerOf2(m+1);
+                applyTwoQubitPhaseShift(qureg, targets[n], targets[n-m-1], arg);
+            }
+        }
+
+        for (int n=0; n<mid; n++)
+            applySwap(qureg, targets[n], targets[numTargets-1-n]);
+    }
 }
 
-void applyFullQuantumFourierTransform(Qureg qureg) {
+void applyFullQuantumFourierTransform(Qureg qureg, bool inverse) {
     validate_quregFields(qureg, __func__);
 
     // tiny; no need to validate alloc
@@ -1765,12 +1787,12 @@ void applyFullQuantumFourierTransform(Qureg qureg) {
     for (size_t i=0; i<targets.size(); i++)
         targets[i] = i;
 
-    applyQuantumFourierTransform(qureg, targets.data(), targets.size());
+    applyQuantumFourierTransform(qureg, targets.data(), targets.size(), inverse);
 }
 
 } // end de-mangler
 
-void applyQuantumFourierTransform(Qureg qureg, vector<int> targets) {
+void applyQuantumFourierTransform(Qureg qureg, vector<int> targets, bool inverse) {
 
-    applyQuantumFourierTransform(qureg, targets.data(), targets.size());
+    applyQuantumFourierTransform(qureg, targets.data(), targets.size(), inverse);
 }
diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp
index 5a2a122c..e7f85a88 100644
--- a/quest/src/api/paulis.cpp
+++ b/quest/src/api/paulis.cpp
@@ -291,3 +291,33 @@ extern "C" void reportPauliStrSum(PauliStrSum sum) {
     // exclude mandatory newline above
     print_oneFewerNewlines();
 }
+
+
+
+/*
+ * SORTING
+ */
+
+
+extern "C" void sortPauliStrSumLexicographic(PauliStrSum sum) {
+    validate_pauliStrSumFields(sum, __func__);
+
+    auto lexSort = [&](qindex i, qindex j) {
+        PauliStr strI = sum.strings[i];
+        PauliStr strJ = sum.strings[j];
+        return std::tie(strI.highPaulis, strI.lowPaulis) < std::tie(strJ.highPaulis, strJ.lowPaulis);
+    };
+
+    paulis_sortTermsViaComparator(sum, lexSort);
+}
+
+
+extern "C" void sortPauliStrSumMagnitude(PauliStrSum sum) {
+    validate_pauliStrSumFields(sum, __func__);
+
+    auto magSort = [&](qindex i, qindex j) {
+        return std::norm(sum.coeffs[i]) > std::norm(sum.coeffs[j]);
+    };
+
+    paulis_sortTermsViaComparator(sum, magSort);
+}
diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp
index af25b961..700077aa 100644
--- a/quest/src/api/trotterisation.cpp
+++ b/quest/src/api/trotterisation.cpp
@@ -15,6 +15,7 @@
 #include "quest/src/core/localiser.hpp"
 #include "quest/src/core/paulilogic.hpp"
 #include "quest/src/core/errors.hpp"
+#include "quest/src/core/randomiser.hpp"
 
 #include <vector>
 
@@ -84,7 +85,7 @@ void internal_applyHigherOrderTrotterRepetition(
 
 void internal_applyAllTrotterRepetitions(
     Qureg qureg, int* controls, int* states, int numControls, 
-    PauliStrSum sum, qcomp angle, int order, int reps, bool onlyLeftApply
+    PauliStrSum sum, qcomp angle, int order, int reps, bool onlyLeftApply, bool permutePaulis
 ) {
     // exp(i angle sum) = identity when angle=0
     if (angle == qcomp(0,0))
@@ -98,9 +99,12 @@ void internal_applyAllTrotterRepetitions(
     qcomp arg = angle / reps;
 
     // perform carefully-ordered sequence of gadgets
-    for (int r=0; r<reps; r++)
+    for (int r=0; r<reps; r++){
+        if (permutePaulis)
+            rand_permutePauliStrSum(sum);
         internal_applyHigherOrderTrotterRepetition(
             qureg, ketCtrlsVec, braCtrlsVec, statesVec, sum, arg, order, onlyLeftApply);
+    }
 
     /// @todo
     /// the accuracy of Trotterisation is greatly improved by randomisation
@@ -147,7 +151,7 @@ qindex internal_getNumTotalSuperPropagatorTerms(PauliStrSum hamil, PauliStrSum*
 
 extern "C" {
 
-void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps) {
+void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps, bool permutePaulis) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumTargets(sum, qureg, __func__);
@@ -156,10 +160,10 @@ void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, q
 
     // |psi> -> U |psi>, rho -> U rho U^dagger
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
-void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumTargets(sum, qureg, __func__);
@@ -167,10 +171,13 @@ void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle
     validate_trotterParams(qureg, order, reps, __func__);
 
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
-void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedControlledPauliStrSumGadget(
+        Qureg qureg, int control, PauliStrSum sum,
+        qreal angle, int order, int reps, bool permutePaulis
+) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumIsHermitian(sum, __func__);
@@ -178,10 +185,13 @@ void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, Pauli
     validate_trotterParams(qureg, order, reps, __func__);
     
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
-void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedMultiControlledPauliStrSumGadget(
+        Qureg qureg, int* controls, int numControls, PauliStrSum sum,
+        qreal angle, int order, int reps, bool permutePaulis
+) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumIsHermitian(sum, __func__);
@@ -189,10 +199,13 @@ void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls
     validate_trotterParams(qureg, order, reps, __func__);
 
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
-void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedMultiStateControlledPauliStrSumGadget(
+        Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum,
+        qreal angle, int order, int reps, bool permutePaulis
+) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumIsHermitian(sum, __func__);
@@ -201,20 +214,26 @@ void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* con
     validate_trotterParams(qureg, order, reps, __func__);
 
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
 } // end de-mangler
 
-void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, vector<int> controls, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedMultiControlledPauliStrSumGadget(
+        Qureg qureg, vector<int> controls, PauliStrSum sum,
+        qreal angle, int order, int reps, bool permutePaulis
+) {
 
-    applyTrotterizedMultiControlledPauliStrSumGadget(qureg, controls.data(), controls.size(), sum, angle, order, reps);
+    applyTrotterizedMultiControlledPauliStrSumGadget(qureg, controls.data(), controls.size(), sum, angle, order, reps, permutePaulis);
 }
 
-void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, vector<int> controls, vector<int> states, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyTrotterizedMultiStateControlledPauliStrSumGadget(
+        Qureg qureg, vector<int> controls, vector<int> states, PauliStrSum sum,
+        qreal angle, int order, int reps, bool permutePaulis
+) {
     validate_controlsMatchStates(controls.size(), states.size(), __func__);
 
-    applyTrotterizedMultiStateControlledPauliStrSumGadget(qureg, controls.data(), states.data(), controls.size(), sum, angle, order, reps);
+    applyTrotterizedMultiStateControlledPauliStrSumGadget(qureg, controls.data(), states.data(), controls.size(), sum, angle, order, reps, permutePaulis);
 }
 
 
@@ -225,7 +244,7 @@ void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, vector<i
 
 extern "C" {
 
-void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps) {
+void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps, bool permutePaulis) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(hamil, __func__);
     validate_pauliStrSumTargets(hamil, qureg, __func__);
@@ -235,10 +254,10 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal
     // exp(-i t H) = exp(x i H) | x=-t
     qcomp angle = - time;
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
-void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps) {
+void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps, bool permutePaulis) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(hamil, __func__);
     validate_pauliStrSumTargets(hamil, qureg, __func__);
@@ -248,7 +267,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea
     // exp(-tau H) = exp(x i H) | x=tau*i
     qcomp angle = qcomp(0, tau);
     bool onlyLeftApply = false;
-    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
 } // end de-mangler
@@ -261,7 +280,10 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea
 
 extern "C" {
 
-void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps) {
+void applyTrotterizedNoisyTimeEvolution(
+        Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps,
+        int numJumps, qreal time, int order, int reps, bool permutePaulis
+) {
     validate_quregFields(qureg, __func__);
     validate_quregIsDensityMatrix(qureg, __func__);
     validate_pauliStrSumFields(hamil, __func__);
@@ -346,7 +368,7 @@ void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* d
     // effect exp(t S) = exp(x i S) | x=-i*time, left-multiplying only
     qcomp angle = qcomp(0, -time);
     bool onlyLeftApply = true;
-    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, superSum, angle, order, reps, onlyLeftApply);
+    internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, superSum, angle, order, reps, onlyLeftApply, permutePaulis);
 }
 
 } // end de-mangler
diff --git a/quest/src/core/bitwise.hpp b/quest/src/core/bitwise.hpp
index 5479c4ee..4d455c2d 100644
--- a/quest/src/core/bitwise.hpp
+++ b/quest/src/core/bitwise.hpp
@@ -28,7 +28,9 @@
  */
 
 
-#define QINDEX_ONE 1ULL
+// alternatives to type-unsafe literals 
+constexpr qindex QINDEX_ZERO = 0; // used by gpu_thrust.cuh
+constexpr qindex QINDEX_ONE  = 1; // used only here
 
 
 INLINE qindex powerOf2(int exponent) {
diff --git a/quest/src/core/paulilogic.cpp b/quest/src/core/paulilogic.cpp
index 1a90dce2..aa7a06ca 100644
--- a/quest/src/core/paulilogic.cpp
+++ b/quest/src/core/paulilogic.cpp
@@ -13,9 +13,12 @@
 #include "quest/src/core/bitwise.hpp"
 #include "quest/src/core/errors.hpp"
 
+#include <numeric>
 #include <utility>
 #include <vector>
 #include <array>
+#include <functional>
+#include <algorithm>
 
 using std::vector;
 
@@ -307,6 +310,37 @@ qindex paulis_getTargetBitMask(PauliStrSum sum) {
 }
 
 
+void paulis_applyPermutationToTerms(PauliStrSum sum, vector<qindex> scatterPermutation) {
+    // permutation passed by value since we modify it
+
+    // scatterPermutation[i] = destination index for element originally at i
+    for (qindex i = 0; i < sum.numTerms; i++) {
+        while (scatterPermutation[i] != i) {
+            qindex j = scatterPermutation[i];
+            std::swap(sum.strings[i], sum.strings[j]);
+            std::swap(sum.coeffs[i], sum.coeffs[j]);
+            std::swap(scatterPermutation[i], scatterPermutation[j]);
+        }
+    }
+}
+
+
+void paulis_sortTermsViaComparator(PauliStrSum sum, std::function<bool(qindex, qindex)> comparator) {
+
+    // TODO: below is an unguarded vector alloc, forgiven since a subsequent
+    // change (giving PauliStrSum an 'ordering' list) supersedes it
+
+    // gatherPermutation[j] = source index of element placed at j
+    vector<qindex> gatherPermutation(sum.numTerms);
+    std::iota(gatherPermutation.begin(), gatherPermutation.end(), 0);
+    std::stable_sort(gatherPermutation.begin(), gatherPermutation.end(), comparator);
+
+    // invert permutation and apply
+    vector<qindex> scatterPermutation = util_getInversePermutation(gatherPermutation);
+    paulis_applyPermutationToTerms(sum, scatterPermutation);
+}
+
+
 void paulis_setPauliStrSumToScaledTensorProdOfConjWithSelf(PauliStrSum out, qreal factor, PauliStrSum in, int numQubits) {
 
     // sets out = factor * conj(in) (x) in, where in has dim of numQubits
diff --git a/quest/src/core/paulilogic.hpp b/quest/src/core/paulilogic.hpp
index f3748b5e..52ffe05f 100644
--- a/quest/src/core/paulilogic.hpp
+++ b/quest/src/core/paulilogic.hpp
@@ -15,6 +15,7 @@
 #include <utility>
 #include <vector>
 #include <array>
+#include <functional>
 
 using std::vector;
 
@@ -76,6 +77,10 @@ int paulis_getIndOfLefmostNonIdentityPauli(PauliStrSum sum);
 
 qindex paulis_getTargetBitMask(PauliStrSum sum);
 
+void paulis_applyPermutationToTerms(PauliStrSum sum, vector<qindex> permutation);
+
+void paulis_sortTermsViaComparator(PauliStrSum sum, std::function<bool(qindex, qindex)> comparator);
+
 
 // below are used exclusively by Trotterisation
 
@@ -88,4 +93,4 @@ void paulis_setPauliStrSumToScaledProdOfAdjointWithSelf(PauliStrSum out, qreal f
 void paulis_setPauliStrSumToShiftedConj(PauliStrSum out, PauliStrSum in, int numQubits);
 
 
-#endif // PAULILOGIC_HPP
\ No newline at end of file
+#endif // PAULILOGIC_HPP
diff --git a/quest/src/core/randomiser.cpp b/quest/src/core/randomiser.cpp
index 45668fd3..1e9b4a94 100644
--- a/quest/src/core/randomiser.cpp
+++ b/quest/src/core/randomiser.cpp
@@ -5,6 +5,7 @@
  * 
  * @author Tyson Jones
  * @author Balint Koczor (patched v3 MSVC seeding)
+ * @author Vasco Ferreira (PauliStrSum permutation)
  */
 
 #include "quest/include/types.h"
@@ -266,3 +267,22 @@ qcomp rand_getThreadPrivateRandomAmp(std::mt19937_64 &gen, std::normal_distribut
     qcomp amp = std::sqrt(prob) * std::exp(phase * 1_i);
     return amp;
 }
+
+
+
+/*
+ * PAULI STRINGS
+ */
+
+
+void rand_permutePauliStrSum(PauliStrSum &sum) {
+
+    // permute ordering of terms inplace using Fisher-Yates
+    for (qindex i = sum.numTerms - 1; i > 0; --i) {
+        std::uniform_int_distribution<qindex> distrib(0, i);
+        qindex j = distrib(mainGenerator);
+
+        std::swap(sum.coeffs[i], sum.coeffs[j]);
+        std::swap(sum.strings[i], sum.strings[j]);
+    }
+}
diff --git a/quest/src/core/randomiser.hpp b/quest/src/core/randomiser.hpp
index 179e5ccb..b867fae7 100644
--- a/quest/src/core/randomiser.hpp
+++ b/quest/src/core/randomiser.hpp
@@ -4,6 +4,7 @@
  * quantum measurements and randomly initialising Quregs.
  * 
  * @author Tyson Jones
+ * @author Vasco Ferreira (PauliStrSum permutation)
  */
 
 #ifndef RANDOMISER_HPP
@@ -61,4 +62,11 @@ qcomp rand_getThreadPrivateRandomAmp(std::mt19937_64 &gen, std::normal_distribut
 
 
 
-#endif // RANDOMISER_HPP
\ No newline at end of file
+/*
+ * PAULI STRINGS
+ */
+
+
+void rand_permutePauliStrSum(PauliStrSum &sum);
+
+#endif // RANDOMISER_HPP
diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp
index 999bd9a7..a5ca635b 100644
--- a/quest/src/core/utilities.cpp
+++ b/quest/src/core/utilities.cpp
@@ -368,7 +368,7 @@ bool util_willSumOverflow(vector<qindex> terms) {
 
 
 /*
- * VECTOR REDUCTION
+ * LIST PROCESSING
  */
 
 qreal util_getSum(vector<qreal> list) {
@@ -387,6 +387,20 @@ qreal util_getSum(vector<qreal> list) {
     return sum;
 }
 
+vector<qindex> util_getInversePermutation(vector<qindex> permutation) {
+    
+    // TODO: below is an unguarded vector alloc, forgiven since a subsequent
+    // change (giving PauliStrSum an 'ordering' list) supersedes it
+
+    qindex numTerms = permutation.size();
+    vector<qindex> out(numTerms);
+
+    for (qindex i = 0; i < numTerms; i++)
+        out[permutation[i]] = i;
+
+    return out;
+}
+
 
 
 /*
diff --git a/quest/src/core/utilities.hpp b/quest/src/core/utilities.hpp
index 4b7fb5db..f2d7087e 100644
--- a/quest/src/core/utilities.hpp
+++ b/quest/src/core/utilities.hpp
@@ -240,11 +240,13 @@ qcomp* util_getGpuMemPtr(T matr) {
 
 
 /*
- * VECTOR REDUCTION
+ * LIST PROCESSING
  */
 
 qreal util_getSum(vector<qreal> list);
 
+vector<qindex> util_getInversePermutation(vector<qindex> permutation);
+
 
 
 /*
@@ -430,4 +432,4 @@ void util_tryAllocMatrix(vector<vector<qcomp>> &vec, qindex numRows, qindex numC
 
 
 
-#endif // UTILITIES_HPP
\ No newline at end of file
+#endif // UTILITIES_HPP
diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp
index 2639b46f..3ac48505 100644
--- a/quest/src/core/validation.cpp
+++ b/quest/src/core/validation.cpp
@@ -3765,7 +3765,7 @@ void validate_measurementOutcomesMatchTargets(int numQubits, int numOutcomes, co
         {"${NUM_QUBITS}",    numQubits},
         {"${NUM_OUTCOMES}",  numOutcomes}};
 
-    assertThat(numQubits == numOutcomes, report::MEASUREMENT_OUTCOMES_MISMATCH_NUM_TARGETS, caller);
+    assertThat(numQubits == numOutcomes, report::MEASUREMENT_OUTCOMES_MISMATCH_NUM_TARGETS, vars, caller);
 }
 
 
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index 9f8d8f1a..3b653dfd 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/quest/src/gpu/gpu_thrust.cuh
@@ -175,21 +175,21 @@ auto getEndPtr(FullStateDiagMatr matr) {
  */
 
 
-struct functor_getAmpConj : public thrust::unary_function<cu_qcomp,cu_qcomp> {
+struct functor_getAmpConj {
 
     __host__ __device__ cu_qcomp operator()(cu_qcomp amp) {
         return getCompConj(amp);
     }
 };
 
-struct functor_getAmpNorm : public thrust::unary_function<cu_qcomp,qreal> {
+struct functor_getAmpNorm {
 
     __host__ __device__ qreal operator()(cu_qcomp amp) {
         return getCompNorm(amp);
     }
 };
 
-struct functor_getAmpReal : public thrust::unary_function<cu_qcomp,qreal> {
+struct functor_getAmpReal {
 
     __host__ __device__ qreal operator()(cu_qcomp amp) {
         return getCompReal(amp);
@@ -197,14 +197,14 @@ struct functor_getAmpReal : public thrust::unary_function<cu_qcomp,qreal> {
 };
 
 
-struct functor_getAmpConjProd : public thrust::binary_function<cu_qcomp,cu_qcomp,cu_qcomp> {
+struct functor_getAmpConjProd {
 
     __host__ __device__ cu_qcomp operator()(cu_qcomp braAmp, cu_qcomp ketAmp) { 
         return getCompConj(braAmp) * ketAmp;
     }
 };
 
-struct functor_getNormOfAmpDif : public thrust::binary_function<cu_qcomp,cu_qcomp,qreal> {
+struct functor_getNormOfAmpDif {
 
     __host__ __device__ qreal operator()(cu_qcomp amp1, cu_qcomp amp2) { 
         return getCompNorm(amp1 - amp2);
@@ -212,7 +212,7 @@ struct functor_getNormOfAmpDif : public thrust::binary_function<cu_qcomp,cu_qcom
 };
 
 
-struct functor_getExpecStateVecZTerm : public thrust::binary_function<qindex,cu_qcomp,qreal> {
+struct functor_getExpecStateVecZTerm {
 
     // this functor computes a single term from the sum
     // in the expectation value of Z of a statevector
@@ -229,7 +229,7 @@ struct functor_getExpecStateVecZTerm : public thrust::binary_function<qindex,cu_
 };
 
 
-struct functor_getExpecDensMatrZTerm : public thrust::unary_function<qindex,cu_qcomp> {
+struct functor_getExpecDensMatrZTerm {
 
     // this functor computes a single term from the sum
     // in the expectation value of Z of a density matrix
@@ -252,7 +252,7 @@ struct functor_getExpecDensMatrZTerm : public thrust::unary_function<qindex,cu_q
 };
 
 
-struct functor_getExpecStateVecPauliTerm : public thrust::unary_function<qindex,cu_qcomp> {
+struct functor_getExpecStateVecPauliTerm {
 
     // this functor computes a single term from the sum in the
     // expectation value of a Pauli str (which necessarily contains 
@@ -276,7 +276,7 @@ struct functor_getExpecStateVecPauliTerm : public thrust::unary_function<qindex,
 };
 
 
-struct functor_getExpecDensMatrPauliTerm : public thrust::unary_function<qindex,cu_qcomp> {
+struct functor_getExpecDensMatrPauliTerm {
 
     // this functor computes a single term from the sum in the
     // expectation value of a Pauli str (which necessarily contains 
@@ -304,7 +304,7 @@ struct functor_getExpecDensMatrPauliTerm : public thrust::unary_function<qindex,
 
 
 template <bool HasPower, bool UseRealPow>
-struct functor_getExpecDensMatrDiagMatrTerm : public thrust::unary_function<qindex,cu_qcomp> {
+struct functor_getExpecDensMatrDiagMatrTerm {
 
     // this functor computes a single term from the sum in the expectation 
     // value of a FullStateDiagMatr upon a density matrix
@@ -372,7 +372,7 @@ struct functor_setAmpToPauliStrSumElem {
 };
 
 
-struct functor_mixAmps : public thrust::binary_function<cu_qcomp,cu_qcomp,cu_qcomp> {
+struct functor_mixAmps {
 
     // this functor linearly combines the given pair
     // of amplitudes, weighted by the fixed qreals,
@@ -389,7 +389,7 @@ struct functor_mixAmps : public thrust::binary_function<cu_qcomp,cu_qcomp,cu_qco
 
 
 template <bool HasPower, bool UseRealPow, bool Norm>
-struct functor_multiplyElemPowerWithAmpOrNorm : public thrust::binary_function<cu_qcomp,cu_qcomp,cu_qcomp> {
+struct functor_multiplyElemPowerWithAmpOrNorm {
 
     // this functor multiplies a diagonal matrix element 
     // raised to a power (templated to optimise away the 
@@ -415,7 +415,7 @@ struct functor_multiplyElemPowerWithAmpOrNorm : public thrust::binary_function<c
 };
 
 
-struct functor_getDiagInd : public thrust::unary_function<qindex,qindex> {
+struct functor_getDiagInd {
 
     // this functor accepts the index of a statevector 
     // basis-state and produces the index of a density 
@@ -435,7 +435,7 @@ struct functor_getDiagInd : public thrust::unary_function<qindex,qindex> {
 
 
 template <int NumBits>
-struct functor_insertBits : public thrust::unary_function<qindex,qindex> {
+struct functor_insertBits {
 
     // this functor inserts bits into a qindex value, and
     // is used to enumerate specific basis-state indices
@@ -463,8 +463,7 @@ struct functor_insertBits : public thrust::unary_function<qindex,qindex> {
 
 
 template <bool Conj>
-struct functor_getFidelityTerm : public thrust::unary_function<qindex,cu_qcomp>
-{
+struct functor_getFidelityTerm {
     int rank, numQubits;
     qindex logNumAmpsPerNode, numAmpsPerCol;
     cu_qcomp *rho, *psi;
@@ -503,7 +502,7 @@ struct functor_getFidelityTerm : public thrust::unary_function<qindex,cu_qcomp>
 
 
 template <int NumTargets>
-struct functor_projectStateVec : public thrust::binary_function<qindex,cu_qcomp,cu_qcomp> {
+struct functor_projectStateVec {
 
     // this functor multiplies an amp with zero or a 
     // renormalisation codfficient, depending on whether
@@ -540,7 +539,7 @@ struct functor_projectStateVec : public thrust::binary_function<qindex,cu_qcomp,
 
 
 template <int NumTargets>
-struct functor_projectDensMatr : public thrust::binary_function<qindex,cu_qcomp,cu_qcomp> {
+struct functor_projectDensMatr {
 
     // this functor multiplies an amp with zero or a 
     // renormalisation coefficient, depending on whether
@@ -585,7 +584,7 @@ struct functor_projectDensMatr : public thrust::binary_function<qindex,cu_qcomp,
 };
 
 
-struct functor_setRandomStateVecAmp : public thrust::unary_function<qindex,cu_qcomp> {
+struct functor_setRandomStateVecAmp {
 
     // this functor generates a random, unnormalised
     // statevector amplitude which, after normalisation
@@ -666,7 +665,7 @@ void thrust_fullstatediagmatr_setElemsToPauliStrSum(FullStateDiagMatr out, Pauli
         rank, out.numElems, logNumElemsPerNode,
         in.numTerms, toCuQcomps(out.gpuElems), devCoeffsPtr, devStringsPtr);
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + out.numElemsPerNode;
     thrust::for_each(indIter, endIter, functor);
 }
@@ -709,7 +708,7 @@ void thrust_densmatr_setAmpsToPauliStrSum_sub(Qureg qureg, PauliStrSum sum) {
         qureg.rank, powerOf2(qureg.numQubits), qureg.logNumAmpsPerNode,
         sum.numTerms, toCuQcomps(qureg.gpuAmps), devCoeffsPtr, devStringsPtr);
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
     thrust::for_each(indIter, endIter, functor);
 }
@@ -772,7 +771,7 @@ qreal thrust_densmatr_calcTotalProb_sub(Qureg qureg) {
     /// reduction may be suboptimal; it may be necessary to
     /// implement a custom numerically-stable CUDA reduction.
 
-    auto rawIter = thrust::make_counting_iterator(0);
+    auto rawIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto indIter = thrust::make_transform_iterator(rawIter, functor_getDiagInd(qureg));
     auto ampIter = thrust::make_permutation_iterator(getStartPtr(qureg), indIter);
     auto probIter= thrust::make_transform_iterator(ampIter, functor_getAmpReal());
@@ -792,7 +791,7 @@ qreal thrust_statevec_calcProbOfMultiQubitOutcome_sub(Qureg qureg, vector<int> q
     auto indFunctor = functor_insertBits<NumQubits>(getPtr(sortedQubits), valueMask, qubits.size());
     auto probFunctor = functor_getAmpNorm();
 
-    auto rawIter = thrust::make_counting_iterator(0);
+    auto rawIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto indIter = thrust::make_transform_iterator(rawIter, indFunctor);
     auto ampIter = thrust::make_permutation_iterator(getStartPtr(qureg), indIter);
     auto probIter = thrust::make_transform_iterator(ampIter, probFunctor);
@@ -815,7 +814,7 @@ qreal thrust_densmatr_calcProbOfMultiQubitOutcome_sub(Qureg qureg, vector<int> q
     auto diagIndFunctor = functor_getDiagInd(qureg);
     auto probFunctor = functor_getAmpReal();
 
-    auto rawIter = thrust::make_counting_iterator(0);
+    auto rawIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto indIter = thrust::make_transform_iterator(rawIter, basisIndFunctor);
     auto diagIter= thrust::make_transform_iterator(indIter, diagIndFunctor);
     auto ampIter = thrust::make_permutation_iterator(getStartPtr(qureg), diagIter);
@@ -865,7 +864,7 @@ cu_qcomp thrust_densmatr_calcFidelityWithPureState_sub(Qureg rho, Qureg psi) {
         rho.rank, rho.numQubits, rho.logNumAmpsPerNode, 
         psi.numAmps, toCuQcomps(rho.gpuAmps), toCuQcomps(psi.gpuAmps));
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     qindex numIts = rho.numAmpsPerNode;
 
     cu_qcomp init = getCuQcomp(0, 0);
@@ -889,7 +888,7 @@ qreal thrust_statevec_calcExpecAnyTargZ_sub(Qureg qureg, vector<int> targs) {
     auto functor = functor_getExpecStateVecZTerm(mask);
 
     qreal init = 0;
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
     return thrust::inner_product(
@@ -906,7 +905,7 @@ cu_qcomp thrust_densmatr_calcExpecAnyTargZ_sub(Qureg qureg, vector<int> targs) {
     auto functor = functor_getExpecDensMatrZTerm(dim, ind, mask, toCuQcomps(qureg.gpuAmps));
 
     cu_qcomp init = getCuQcomp(0, 0);
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
     return thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<cu_qcomp>());
@@ -921,7 +920,7 @@ cu_qcomp thrust_statevec_calcExpecPauliStr_subA(Qureg qureg, vector<int> x, vect
     auto functor = functor_getExpecStateVecPauliTerm(maskXY, maskYZ, ampsPtr, ampsPtr); // amps=pairAmps
 
     cu_qcomp init = getCuQcomp(0, 0);
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
     cu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<cu_qcomp>());
@@ -939,7 +938,7 @@ cu_qcomp thrust_statevec_calcExpecPauliStr_subB(Qureg qureg, vector<int> x, vect
     auto functor = functor_getExpecStateVecPauliTerm(maskXY, maskYZ, ampsPtr, buffPtr);
 
     cu_qcomp init = getCuQcomp(0, 0);
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
     cu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<cu_qcomp>());
@@ -957,7 +956,7 @@ cu_qcomp thrust_densmatr_calcExpecPauliStr_sub(Qureg qureg, vector<int> x, vecto
     auto functor = functor_getExpecDensMatrPauliTerm(mXY, mYZ, dim, ind, toCuQcomps(qureg.gpuAmps));
 
     cu_qcomp init = getCuQcomp(0, 0);
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
     cu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<cu_qcomp>());
@@ -996,7 +995,7 @@ cu_qcomp thrust_densmatr_calcExpecFullStateDiagMatr_sub(Qureg qureg, FullStateDi
     auto functor = functor_getExpecDensMatrDiagMatrTerm<HasPower,UseRealPow>(dim, ind, ampsPtr, elemsPtr, expo);
 
     cu_qcomp init = getCuQcomp(0, 0);
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
     return thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<cu_qcomp>());
@@ -1017,7 +1016,7 @@ void thrust_statevec_multiQubitProjector_sub(Qureg qureg, vector<int> qubits, ve
     auto projFunctor = functor_projectStateVec<NumQubits>(
         getPtr(devQubits), qubits.size(), retainValue, renorm);
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
@@ -1034,7 +1033,7 @@ void thrust_densmatr_multiQubitProjector_sub(Qureg qureg, vector<int> qubits, ve
         getPtr(devQubits), qubits.size(), qureg.rank, qureg.numQubits,
         qureg.logNumAmpsPerNode, retainValue, renorm);
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
@@ -1074,7 +1073,7 @@ void thrust_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qure
     unsigned seed = rand_getThreadSharedRandomSeed(qureg.isDistributed);
     auto functor = functor_setRandomStateVecAmp(seed);
 
-    auto indIter = thrust::make_counting_iterator(0);
+    auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
diff --git a/tests/unit/initialisations.cpp b/tests/unit/initialisations.cpp
index 649b9290..ac1f1abd 100644
--- a/tests/unit/initialisations.cpp
+++ b/tests/unit/initialisations.cpp
@@ -486,8 +486,7 @@ TEST_CASE( "setQuregToWeightedSum", TEST_CATEGORY ) {
 
     SECTION( LABEL_VALIDATION ) {
 
-        // arbitrary existing qureg
-        Qureg qureg = getCachedStatevecs().begin()->second;
+        Qureg qureg = getArbitraryCachedStatevec();
 
         SECTION( "out qureg uninitialised" ) {
 
@@ -702,7 +701,7 @@ TEST_CASE( "setQuregToMixture", TEST_CATEGORY ) {
 
         SECTION( "out qureg is statevector" ) {
 
-            Qureg badQureg = getCachedStatevecs().begin()->second;
+            Qureg badQureg = getArbitraryCachedStatevec();
 
             REQUIRE_THROWS_WITH( 
                 setQuregToMixture(badQureg, nullptr, nullptr, 1), 
@@ -717,7 +716,7 @@ TEST_CASE( "setQuregToMixture", TEST_CATEGORY ) {
 
             // hide a statevector among them
             int badInd = GENERATE_COPY( range(0,numIn) );
-            inQuregs[badInd] = getCachedStatevecs().begin()->second;;
+            inQuregs[badInd] = getArbitraryCachedStatevec();
 
             REQUIRE_THROWS_WITH( 
                 setQuregToMixture(qureg, nullptr, inQuregs.data(), numIn), 
diff --git a/tests/unit/operations.cpp b/tests/unit/operations.cpp
index 05d46418..0e33220d 100644
--- a/tests/unit/operations.cpp
+++ b/tests/unit/operations.cpp
@@ -24,6 +24,7 @@
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/catch_template_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_string.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
 #include <catch2/generators/catch_generators_range.hpp>
 
 #include "tests/utils/config.hpp"
@@ -1042,8 +1043,12 @@ void testOperationValidation(auto operation) {
 
     SECTION( "targeted amps fit in node" ) {
 
-        // simplest to trigger validation using a statevector
-        qureg = getCachedStatevecs().begin()->second;
+        // use any qureg which is otherwise compatible
+        // (but beware statevecs vs density matrices permit
+        //  different num targets before validation is triggered)
+        qureg = (Apply == rightapply)?
+            getArbitraryCachedDensmatr():
+            getArbitraryCachedStatevec();
 
         // can only be validated when environment AND qureg
         // are distributed (over more than 1 node, of course)
@@ -1053,7 +1058,8 @@ void testOperationValidation(auto operation) {
         // can only be validated if forced ctrl qubits permit 
         // enough remaining targets
         int minNumCtrls = (Ctrls == one)? 1 : 0;
-        int minNumTargs = numQubits - qureg.logNumNodes + 1;
+        int numVecQubits = (qureg.isDensityMatrix)? 2*numQubits : numQubits;
+        int minNumTargs = numVecQubits - qureg.logNumNodes + 1;
         int maxNumTargs = numQubits - minNumCtrls;
         if (minNumTargs > maxNumTargs)
             return;
@@ -1132,9 +1138,8 @@ void testOperationValidation(auto operation) {
         if (Apply != rightapply)
             return;
 
-        // use any statevector
-        qureg = getCachedStatevecs().begin()->second;
-
+        // override qureg in apiyFunc with a statevector
+        qureg = getArbitraryCachedStatevec();
         REQUIRE_THROWS_WITH( apiFunc(), ContainsSubstring("Expected a density matrix") );
     }
 
@@ -1370,16 +1375,18 @@ TEST_CASE( "applyQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
         int numTargs = GENERATE_COPY( range(1,numQubits+1) );
         auto targs = GENERATE_TARGS( numQubits, numTargs );
+        bool inverse = GENERATE(false, true);
 
         CAPTURE( targs );
 
         SECTION( LABEL_STATEVEC ) { 
 
             auto testFunc = [&](Qureg qureg, qvector& ref) {
-                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
-                ref = getDisceteFourierTransform(ref, targs);
+                applyQuantumFourierTransform(qureg, targs.data(), targs.size(), inverse);
+                ref = getDiscreteFourierTransform(ref, targs, inverse);
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
         }
 
@@ -1393,20 +1400,70 @@ TEST_CASE( "applyQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
                 // overwrite the Qureg debug state set by caller to above mixture
                 setQuregToReference(qureg, getMixture(states, probs));
-                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
+                applyQuantumFourierTransform(qureg, targs.data(), targs.size(), inverse);
                 
                 ref = getZeroMatrix(ref.size());
                 for (size_t i=0; i<states.size(); i++) {
-                    qvector vec = getDisceteFourierTransform(states[i], targs);
+                    qvector vec = getDiscreteFourierTransform(states[i], targs, inverse);
                     ref += probs[i] * getOuterProduct(vec, vec);
                 }
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
         }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targs[] = {0, 1, 2};
+        int numTargs = 3;
+        bool inverse = false;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyQuantumFourierTransform(badQureg, targs, numTargs, inverse),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargs[] = {0, 1, qureg.numQubits}; // latter is too large
+            REQUIRE_THROWS_WITH(
+                applyQuantumFourierTransform(qureg, badTargs, 3, inverse),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+
+            int dupTargs[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                applyQuantumFourierTransform(qureg, dupTargs, 3, inverse),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+
+            int badNumTargs = GENERATE_COPY( -1, 0 ); 
+            REQUIRE_THROWS_WITH(
+                applyQuantumFourierTransform(qureg, targs, badNumTargs, inverse),
+                ContainsSubstring("targets") || ContainsSubstring("target qubits")
+            );
+
+            badNumTargs = qureg.numQubits+1;
+            REQUIRE_THROWS_WITH(
+                applyQuantumFourierTransform(qureg, targs, badNumTargs, inverse),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+    }
 }
 
 
@@ -1417,14 +1474,16 @@ TEST_CASE( "applyFullQuantumFourierTransform", TEST_CATEGORY_OPS ) {
     SECTION( LABEL_CORRECTNESS ) {
 
         GENERATE( range(0,10) );
+        bool inverse = GENERATE(false, true);
 
         SECTION( LABEL_STATEVEC ) { 
 
             auto testFunc = [&](Qureg qureg, qvector& ref) {
-                applyFullQuantumFourierTransform(qureg);
-                ref = getDisceteFourierTransform(ref);
+                applyFullQuantumFourierTransform(qureg, inverse);
+                ref = getDiscreteFourierTransform(ref, inverse);
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
         }
 
@@ -1438,20 +1497,35 @@ TEST_CASE( "applyFullQuantumFourierTransform", TEST_CATEGORY_OPS ) {
 
                 // overwrite the Qureg debug state set by caller to above mixture
                 setQuregToReference(qureg, getMixture(states, probs));
-                applyFullQuantumFourierTransform(qureg);
+                applyFullQuantumFourierTransform(qureg, inverse);
                 
                 ref = getZeroMatrix(ref.size());
                 for (size_t i=0; i<states.size(); i++) {
-                    qvector vec = getDisceteFourierTransform(states[i]);
+                    qvector vec = getDiscreteFourierTransform(states[i], inverse);
                     ref += probs[i] * getOuterProduct(vec, vec);
                 }
             };
 
+            CAPTURE(inverse);
             TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
         }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        bool inverse = false;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyFullQuantumFourierTransform(badQureg, inverse),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+    }
 }
 
 
@@ -1477,7 +1551,40 @@ TEST_CASE( "applyQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyQubitProjector(badQureg, 0, 0),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                applyQubitProjector(qureg, badTarget, 0),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "invalid outcome" ) {
+
+            int badOutcome = GENERATE_COPY( -1, 2 );
+            REQUIRE_THROWS_WITH(
+                applyQubitProjector(qureg, 0, badOutcome),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
@@ -1503,7 +1610,75 @@ TEST_CASE( "applyMultiQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targets[] = {0, 1, 2};
+        int outcomes[] = {0, 1, 0};
+        int numTargets = 3;
+
+        SECTION( "qureg uninitialised" ) {
+            
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(badQureg, targets, outcomes, numTargets),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY( -1, qureg.numQubits) };
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, badTargets, outcomes, numTargets),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, dupTargets, outcomes, numTargets),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( 0, -1 );
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+
+        SECTION( "invalid outcomes" ) {
+
+            int badOutcomes[] = {0, 1, GENERATE( -1, 2 ) };
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        SECTION( "targets mismatch outcomes (C++ only)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitProjector(qureg, {0,1,2}, {0,1}),
+                ContainsSubstring("outcomes") && ContainsSubstring("inconsistent with the given number of qubits")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
@@ -1542,7 +1717,77 @@ TEST_CASE( "applyForcedQubitMeasurement", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+
+        // below validation tests assume qubit 0 can collapse to either outcome
+        // (which does not require normalisation; qureg can be in the debug state)
+        initDebugState(qureg);
+        REQUIRE( calcProbOfQubitOutcome(qureg, 0, 0) > getValidationEpsilon() );
+        REQUIRE( calcProbOfQubitOutcome(qureg, 0, 1) > getValidationEpsilon() );
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyForcedQubitMeasurement(badQureg, 0, 0),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                applyForcedQubitMeasurement(qureg, badTarget, 0),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "invalid outcome" ) {
+
+            int badOutcome = GENERATE_COPY( -1, 2 );
+            REQUIRE_THROWS_WITH(
+                applyForcedQubitMeasurement(qureg, 0, badOutcome),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        SECTION( "improbable outcome" ) {
+
+            // precisely zero probability outcome
+            initZeroState(qureg);
+            int badOutcome = 1; // impossible
+            REQUIRE_THROWS_WITH(
+                applyForcedQubitMeasurement(qureg, 0, badOutcome),
+                ContainsSubstring("impossibly unlikely")
+            );
+
+            // outcome of non-zero probability smaller than epsilon
+            qreal badTheta = 1E-8;
+            applyRotateX(qureg, 0, badTheta); // causes prob(1) = sin^2(theta) ~ 2.5E-17
+            REQUIRE_THROWS_WITH(
+                applyForcedQubitMeasurement(qureg, 0, badOutcome),
+                ContainsSubstring("impossibly unlikely")
+            );
+
+            // confirm that >epsilon probability is fine
+            initZeroState(qureg);
+            qreal goodTheta = 0.1;
+            applyRotateX(qureg, 0, goodTheta);
+            REQUIRE( 
+                calcProbOfQubitOutcome(qureg, 0, badOutcome) > getValidationEpsilon() 
+            );
+            REQUIRE_NOTHROW(
+                applyForcedQubitMeasurement(qureg, 0, badOutcome)
+            );
+
+            // restore qureg state
+            initDebugState(qureg);
+        }
+    }
 }
 
 
@@ -1588,7 +1833,111 @@ TEST_CASE( "applyForcedMultiQubitMeasurement", TEST_CATEGORY_OPS ) {
         setValidationEpsilonToDefault();
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targets[] = {0, 1, 2};
+        int outcomes[] = {0, 1, 0};
+        int numTargets = 3;
+
+        // below validation tests assume the above parameters are valid (not impossibly unlikely)
+        initDebugState(qureg);
+        REQUIRE( calcProbOfMultiQubitOutcome(qureg, targets, outcomes, numTargets) > getValidationEpsilon() );
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(badQureg, targets, outcomes, numTargets),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY( -1, qureg.numQubits )};
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, badTargets, outcomes, numTargets),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+            
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, dupTargets, outcomes, numTargets),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( -1, 0 );
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+
+        SECTION( "invalid outcomes" ) {
+
+            int badOutcomes[] = {0, 1, GENERATE( -1, 2 )};
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        SECTION( "improbable outcomes" ) {
+
+            // impossible outcome
+            initZeroState(qureg);
+            int badOutcomes[] = {0, 0, 1};
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("impossibly unlikely")
+            );
+
+            // outcome of non-zero probability smaller than epsilon
+            qreal badTheta = 1E-8;
+            applyRotateX(qureg, targets[2], badTheta);
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("impossibly unlikely")
+            );
+
+            // confirm that >epsilon probability is fine
+            initZeroState(qureg);
+            qreal goodTheta = 0.1;
+            applyRotateX(qureg, targets[2], goodTheta);
+            int goodOutcomes[] = {0, 0, 1};
+            REQUIRE( 
+                calcProbOfMultiQubitOutcome(qureg, targets, goodOutcomes, numTargets) > getValidationEpsilon() 
+            );
+            REQUIRE_NOTHROW(
+                applyForcedMultiQubitMeasurement(qureg, targets, goodOutcomes, numTargets)
+            );
+
+            // restore qureg state
+            initDebugState(qureg);
+        }
+
+        SECTION( "targets mismatch outcomes (C++ only)") {
+
+            REQUIRE_THROWS_WITH(
+                applyForcedMultiQubitMeasurement(qureg, {0,1}, {0,1,1}),
+                ContainsSubstring("inconsistent")
+            );
+        }
+    }
 }
 
 
@@ -1625,7 +1974,55 @@ TEST_CASE( "applyMultiQubitMeasurement", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targets[] = {0, 1, 2};
+        int numTargets = 3;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurement(badQureg, targets, numTargets),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY(-1, qureg.numQubits)};
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurement(qureg, badTargets, numTargets),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+            
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurement(qureg, dupTargets, numTargets),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( -1, 0 );
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurement(qureg, targets, badNumTargs),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurement(qureg, targets, badNumTargs),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+    }
 }
 
 
@@ -1664,7 +2061,56 @@ TEST_CASE( "applyMultiQubitMeasurementAndGetProb", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targets[] = {0, 1, 2};
+        int numTargets = 3;
+        qreal outProb = 0;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurementAndGetProb(badQureg, targets, numTargets, nullptr),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY( -1, qureg.numQubits )};
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurementAndGetProb(qureg, badTargets, numTargets, &outProb),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+            
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurementAndGetProb(qureg, dupTargets, numTargets, &outProb),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( -1, 0 );
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurementAndGetProb(qureg, targets, badNumTargs, &outProb),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                applyMultiQubitMeasurementAndGetProb(qureg, targets, badNumTargs, &outProb),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+    }
 }
 
 
@@ -1700,7 +2146,29 @@ TEST_CASE( "applyQubitMeasurement", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyQubitMeasurement(badQureg, 0),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                applyQubitMeasurement(qureg, badTarget),
+                ContainsSubstring("target")
+            );
+        }
+    }
 }
 
 
@@ -1738,7 +2206,30 @@ TEST_CASE( "applyQubitMeasurementAndGetProb", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        qreal outProb = 0;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyQubitMeasurementAndGetProb(badQureg, 0, &outProb),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                applyQubitMeasurementAndGetProb(qureg, badTarget, &outProb),
+                ContainsSubstring("target")
+            );
+        }
+    }
 }
 
 
@@ -1854,7 +2345,23 @@ TEST_CASE( "applyNonUnitaryPauliGadget", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        PauliStr str = getPauliStr("XY", {0, 1});
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyNonUnitaryPauliGadget(badQureg, str, qcomp(0.5, 0)),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        /// @todo remaining input validation
+    }
 }
 
 
@@ -2092,7 +2599,7 @@ TEST_CASE( "rightapplyFullStateDiagMatrPower", TEST_CATEGORY_MULT LABEL_MIXED_DE
 }
 
 
-TEST_CASE( "leftapplyQubitProjector", TEST_CATEGORY_OPS ) {
+TEST_CASE( "leftapplyQubitProjector", TEST_CATEGORY_MULT ) {
 
     PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
@@ -2114,11 +2621,44 @@ TEST_CASE( "leftapplyQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                leftapplyQubitProjector(badQureg, 0, 0),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+            
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                leftapplyQubitProjector(qureg, badTarget, 0),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "invalid outcome" ) {
+
+            int badOutcome = GENERATE_COPY( -1, 2 );
+            REQUIRE_THROWS_WITH(
+                leftapplyQubitProjector(qureg, 0, badOutcome),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
-TEST_CASE( "rightapplyQubitProjector", TEST_CATEGORY_OPS ) {
+TEST_CASE( "rightapplyQubitProjector", TEST_CATEGORY_MULT ) {
 
     PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
@@ -2139,11 +2679,53 @@ TEST_CASE( "rightapplyQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedDensmatr();
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                rightapplyQubitProjector(badQureg, 0, 0),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "qureg is not density matrix" ) {
+
+            Qureg badQureg = getArbitraryCachedStatevec();
+            REQUIRE_THROWS_WITH(
+                rightapplyQubitProjector(badQureg, 0, 0),
+                ContainsSubstring("received a statevector")
+            );
+        }
+
+        SECTION( "invalid target qubit" ) {
+            
+            int badTarget = GENERATE_COPY( -1, qureg.numQubits );
+            REQUIRE_THROWS_WITH(
+                rightapplyQubitProjector(qureg, badTarget, 0),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "invalid outcome" ) {
+
+            int badOutcome = GENERATE_COPY( -1, 2 );
+            REQUIRE_THROWS_WITH(
+                rightapplyQubitProjector(qureg, 0, badOutcome),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
-TEST_CASE( "leftapplyMultiQubitProjector", TEST_CATEGORY_OPS ) {
+TEST_CASE( "leftapplyMultiQubitProjector", TEST_CATEGORY_MULT ) {
 
     PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
@@ -2165,11 +2747,79 @@ TEST_CASE( "leftapplyMultiQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        int targets[] = {0, 1, 2};
+        int outcomes[] = {0, 1, 0};
+        int numTargets = 3;
+
+        SECTION( "qureg uninitialised" ) {
+            
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(badQureg, targets, outcomes, numTargets),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY( -1, qureg.numQubits) };
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, badTargets, outcomes, numTargets),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, dupTargets, outcomes, numTargets),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( 0, -1 );
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+
+        SECTION( "invalid outcomes" ) {
+
+            int badOutcomes[] = {0, 1, GENERATE( -1, 2 ) };
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        SECTION( "targets mismatch outcomes (C++ only)" ) {
+
+            REQUIRE_THROWS_WITH(
+                leftapplyMultiQubitProjector(qureg, {0,1,2}, {0,1}),
+                ContainsSubstring("inconsistent")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
-TEST_CASE( "rightapplyMultiQubitProjector", TEST_CATEGORY_OPS ) {
+TEST_CASE( "rightapplyMultiQubitProjector", TEST_CATEGORY_MULT ) {
 
     PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
@@ -2190,7 +2840,84 @@ TEST_CASE( "rightapplyMultiQubitProjector", TEST_CATEGORY_OPS ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedDensmatr();
+        int targets[] = {0, 1, 2};
+        int outcomes[] = {0, 1, 0};
+        int numTargets = 3;
+
+        SECTION( "qureg uninitialised" ) {
+            
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(badQureg, targets, outcomes, numTargets),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "qureg is not density matrix" ) {
+
+            Qureg badQureg = getArbitraryCachedStatevec();
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(badQureg, targets, outcomes, numTargets),
+                ContainsSubstring("received a statevector")
+            );
+        }
+
+        SECTION( "invalid target qubits" ) {
+
+            int badTargets[] = {0, 1, GENERATE_COPY( -1, qureg.numQubits) };
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, badTargets, outcomes, numTargets),
+                ContainsSubstring("target")
+            );
+        }
+
+        SECTION( "duplicate target qubits" ) {
+
+            int dupTargets[] = {0, 1, 1};
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, dupTargets, outcomes, numTargets),
+                ContainsSubstring("duplicate")
+            );
+        }
+
+        SECTION( "invalid number of targets" ) {
+            
+            int badNumTargs = GENERATE( 0, -1 );
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("targets")
+            );
+
+            badNumTargs = qureg.numQubits + 1;
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, targets, outcomes, badNumTargs),
+                ContainsSubstring("exceeds the number of qubits in the Qureg")
+            );
+        }
+
+        SECTION( "invalid outcomes" ) {
+
+            int badOutcomes[] = {0, 1, GENERATE( -1, 2 ) };
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, targets, badOutcomes, numTargets),
+                ContainsSubstring("outcome")
+            );
+        }
+
+        SECTION( "targets mismatch outcomes (C++ only)" ) {
+
+            REQUIRE_THROWS_WITH(
+                rightapplyMultiQubitProjector(qureg, {0,1,2}, {0,1}),
+                ContainsSubstring("inconsistent")
+            );
+        }
+
+        // projector does NOT validate outcome probability
+    }
 }
 
 
@@ -2220,7 +2947,27 @@ TEST_CASE( "leftapplyPauliStrSum", TEST_CATEGORY_MULT LABEL_MIXED_DEPLOY_TAG ) {
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+       
+        Qureg qureg = getArbitraryCachedStatevec();
+        PauliStrSum sum = createRandomPauliStrSum(numQubits, 2);
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            Qureg workspace = createCloneQureg(qureg);
+            REQUIRE_THROWS_WITH(
+                leftapplyPauliStrSum(badQureg, sum, workspace),
+                ContainsSubstring("invalid Qureg")
+            );
+            destroyQureg(workspace);
+        }
+
+        destroyPauliStrSum(sum);
+
+        /// @todo remaining input validation
+    }
 }
 
 
@@ -2249,7 +2996,27 @@ TEST_CASE( "rightapplyPauliStrSum", TEST_CATEGORY_MULT LABEL_MIXED_DEPLOY_TAG )
         SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
     }
 
-    /// @todo input validation
+    SECTION( LABEL_VALIDATION ) {
+        
+        Qureg qureg = getArbitraryCachedStatevec();
+        PauliStrSum sum = createRandomPauliStrSum(numQubits, 2);
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            Qureg workspace = createCloneQureg(qureg);
+            REQUIRE_THROWS_WITH(
+                rightapplyPauliStrSum(badQureg, sum, workspace),
+                ContainsSubstring("invalid Qureg")
+            );
+            destroyQureg(workspace);
+        }
+
+        destroyPauliStrSum(sum);
+
+        /// @todo remaining input validation
+    }
 }
 
 
diff --git a/tests/unit/paulis.cpp b/tests/unit/paulis.cpp
index 8b6fc4f9..e3339100 100644
--- a/tests/unit/paulis.cpp
+++ b/tests/unit/paulis.cpp
@@ -587,6 +587,61 @@ TEST_CASE( "destroyPauliStrSum", TEST_CATEGORY ) {
     }
 }
 
+TEST_CASE( "sortPauliStrSumLexicographic", TEST_CATEGORY ) {
+
+    SECTION( LABEL_CORRECTNESS ) {
+
+        vector<qcomp> coeffs = {0.1_i, 2+1_i, 5, 3+4_i};
+        vector<PauliStr> strings = {
+            getPauliStr("XY", {31,32}),
+            getPauliStr("YX", {0,1}),
+            getPauliStr("II", {0,1}),
+            getPauliStr("YY", {31,32})
+        };
+
+        PauliStrSum sum = createPauliStrSum(strings, coeffs);
+        sortPauliStrSumLexicographic(sum);
+
+        REQUIRE(sum.coeffs[0] == 5+0_i);
+        REQUIRE(sum.coeffs[1] == 2+1_i);
+        REQUIRE(sum.coeffs[3] == 3+4_i);
+
+        REQUIRE(sum.strings[0].lowPaulis == 0);
+        REQUIRE(sum.strings[1].lowPaulis == 2 + 1*4);
+        REQUIRE(sum.strings[3].highPaulis == 2);
+        REQUIRE(sum.strings[3].lowPaulis == 2*std::pow(4, 31));
+
+        destroyPauliStrSum(sum);
+    }
+}
+
+TEST_CASE( "sortPauliStrSumMagnitude", TEST_CATEGORY ) {
+
+    SECTION( LABEL_CORRECTNESS ) {
+
+        vector<qcomp> coeffs = {0.1_i, 2+1_i, 5, 3+4_i};
+        vector<PauliStr> strings = {
+            getPauliStr("XY", {0,1}),
+            getPauliStr("ZX", {0,1}),
+            getPauliStr("II", {0,1}),
+            getPauliStr("YZ", {0,1})
+        };
+
+        PauliStrSum sum = createPauliStrSum(strings, coeffs);
+        sortPauliStrSumMagnitude(sum);
+
+        REQUIRE(sum.coeffs[0] == 5+0_i);
+        REQUIRE(sum.coeffs[1] == 3+4_i);
+        REQUIRE(sum.coeffs[3] == 0+0.1_i);
+
+        REQUIRE(sum.strings[0].lowPaulis == 0);
+        REQUIRE(sum.strings[1].lowPaulis == 2 + 3*4);
+        REQUIRE(sum.strings[3].lowPaulis == 1 + 2*4);
+
+        destroyPauliStrSum(sum);
+    }
+}
+
 
 /** @} (end defgroup) */
 
diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp
index cc691402..35af81ba 100644
--- a/tests/unit/trotterisation.cpp
+++ b/tests/unit/trotterisation.cpp
@@ -2,6 +2,9 @@
  * Unit tests of the trotterisation module.
  *
  * @author Tyson Jones
+ * @author Vasco Ferreira (initial Pauli permutation tests)
+ * @author Maurice Jamieson (real and imaginary time evolution tests)
+ * @author Oliver Thomson Brown (real and imaginary time evolution tests)
  * 
  * @defgroup unittrotter Trotterisation
  * @ingroup unittests
@@ -9,7 +12,22 @@
 
 #include "quest.h"
 
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/generators/catch_generators_range.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/matchers/catch_matchers_string.hpp>
 
+#include "tests/utils/macros.hpp"
+#include "tests/utils/cache.hpp"
+#include "tests/utils/compare.hpp"
+#include "tests/utils/random.hpp"
+
+#include <vector>
+#include <string>
+
+using std::vector;
+using std::string;
+using namespace Catch::Matchers;
 
 /*
  * UTILITIES
@@ -18,25 +36,625 @@
 #define TEST_CATEGORY \
     LABEL_UNIT_TAG "[trotterisation]"
 
+void TEST_ON_CACHED_QUREGS(quregCache quregs, auto& refFunc, auto& regularFunc, auto& randFunc) {
+    for (auto& [label, refQureg]: quregs) {
+
+        DYNAMIC_SECTION( label ) {
+            initDebugState(refQureg);
+
+            Qureg regularQureg = createCloneQureg(refQureg);
+            Qureg randQureg = createCloneQureg(refQureg);
+
+            refFunc(refQureg);
+            regularFunc(regularQureg);
+            randFunc(randQureg);
+
+            double regularDistance = calcDistance(regularQureg, refQureg);
+            double randDistance = calcDistance(randQureg, refQureg);
+
+            REQUIRE( randDistance < regularDistance );
+
+            destroyQureg(regularQureg);
+            destroyQureg(randQureg);
+        }
+    }
+}
+
+/*
+ * Prepare a Hamiltonian H under which dynamical
+ * evolution will be simulated via Trotterisation
+ * of unitary-time evolution operator e^(-itH).
+ * If the Hamiltonian was fixed/known in advance,
+ * we could instead use createInlinePauliStrSum()
+ *
+ * (Adapted from dynamics.cpp, @author Tyson Jones)
+ */
+PauliStrSum createHeisenbergHamiltonian(int numQubits) {
+
+    // we prepare a Heisenberg XYZ spin-ring Hamiltonian,
+    // i.e. H = -1/2 sum( Jx XX + Jy YY + Jz ZZ + h Z )
+    // upon all nearest neighbour qubits, with periodicity.
+    // The coefficients must be real for H to be Hermitian
+    // and ergo its time-evolution operator to be unitary,
+    // although they must be represented with a qcomp type.
+    vector<string> operators = {"XX", "YY", "ZZ", "Z"};
+    vector<qcomp> coefficients = {.1, .2, .3, .4}; // Jx,Jy,Jz,h
+
+    // we will populate the below vectors with 4*numQubits
+    // elements which we could pre-allocate with .reserve,
+    // but we might incur Donald Knuth's justified wrath.
+    vector<PauliStr> allStrings;
+    vector<qcomp> allCoeffs;
+    
+    // prepare all XX + YY + ZZ
+    for (int p=0; p<3; p++) {
+        for (int i=0; i<numQubits; i++) {
+
+            // A_i, A_i+1
+            vector<int> targs = {i, (i+1)%numQubits};
+            PauliStr str = getPauliStr(operators[p], targs);
+
+            allStrings.push_back(str);
+            allCoeffs.push_back(coefficients[p]);
+        }
+    }
+
+    // prepare Z
+    for (int i=0; i<numQubits; i++) {
+        allStrings.push_back(getPauliStr(operators[3], {i}));
+        allCoeffs.push_back(coefficients[3]);
+    }
+
+    // must be freed by caller
+    return createPauliStrSum(allStrings, allCoeffs);
+}
+
+/*
+ * Prepare the observable operator O under which the
+ * evolved state (under H above) will be measured.
+ * If this were one term (a single tensor product of
+ * Pauli operators), we could return instead a PauliStr
+ * but we here return an arbitrary weighted sum thereof.
+ *
+ * (Adapted from dynamics.cpp, @author Tyson Jones)
+ */
+PauliStrSum createAlternatingPauliObservable(int numQubits) {
+
+    // we prepare a weighted sum of alternating Paulis
+    // upon each qubit, i.e. 1 X0 + 2 Y1 + 3 Z2 + 1 X3 + ...
+    // where the coefficients are real such that the
+    // output observable is Hermitian.
+
+    vector<PauliStr> strings(numQubits);
+    vector<qcomp> coeffs(numQubits);
+
+    for (int i=0; i<numQubits; i++) {
+        strings[i] = getPauliStr({"XYZ"[i%3]}, {i});
+        coeffs[i] = getQcomp(i%4 + 1, 0);
+    }
+
+    // must be freed by caller
+    return createPauliStrSum(strings, coeffs);
+}
+
+/*
+ * Constructs a PauliStrSum representing a 1D Hamiltonian of the form
+ * H = - \mu \sum^{N}_{j} Z_{j} - J \sum_{<ij>}^{N} Z_{i}Z_{j}
+ * where,
+ * \mu = magField,
+ * J = interactionStrength,
+ * <ij> indicates nearest-neighbour interactions only,
+ * and boundary conditions are periodic such that site N-1 interacts with site 0.
+ *
+ * The asymmetricBias term can be used to break the symmetry of the system
+ * in order to 'choose' a preferred antiferromagnetic state, and ensure repeatable
+ * predictable outcomes.
+ * It adds a term of the form:
+ * -BZ_{0}
+ */
+PauliStrSum createIsingHamiltonian(int numQubits, qreal magField, 
+                                   qreal interactionStrength, qreal asymmetricBias) {
+    const int NTERMS = 2 * numQubits + 1;
+    
+    vector<qcomp> coeffs;
+    vector<PauliStr> pauli_terms;
+    coeffs.reserve(NTERMS);
+    pauli_terms.reserve(NTERMS);
+    
+    for (int i = 0; i < numQubits; ++i) {
+        pauli_terms.push_back(getPauliStr("Z", {i}));
+        coeffs.push_back(getQcomp(-magField, 0));
+        
+        int next = (i + 1) % numQubits;
+        pauli_terms.push_back(getPauliStr("ZZ", {i, next}));
+        coeffs.push_back(getQcomp(-interactionStrength, 0));
+    }
+    
+    pauli_terms.push_back(getPauliStr("Z", {0}));
+    coeffs.push_back(getQcomp(-asymmetricBias, 0));
+    
+    return createPauliStrSum(pauli_terms, coeffs);
+}
+
+
+/* 
+ * TESTS
+ */
 
 
 /**
  * @todo
- * UNTESTED FUNCTIONS
+ * Basic validation for randomisation, should be expanded and merged
+ * once the Trotterisation function tests have been implemented.
  */
+TEST_CASE( "randomisedTrotter", TEST_CATEGORY ) {
+
+    SECTION( LABEL_CORRECTNESS ) {
+
+        int numQubits = getNumCachedQubits();
+        int numTerms = 25;
+        int reps = 50;
+        double time = 1.0;
+
+        int refOrder = 4;
+        int order = GENERATE_COPY(1, 2);
+
+        GENERATE( range(0, 10) );
+        PauliStrSum sum = createRandomPauliStrSum(numQubits, numTerms);
+
+        auto refFunc = [&](Qureg qureg) { applyTrotterizedUnitaryTimeEvolution(qureg, sum, time, refOrder, reps, false); };
+        auto regularFunc = [&](Qureg qureg) { applyTrotterizedUnitaryTimeEvolution(qureg, sum, time, order, reps, false); };
+        auto randFunc = [&](Qureg qureg) { applyTrotterizedUnitaryTimeEvolution(qureg, sum, time, order, reps, true); };
+        
+        TEST_ON_CACHED_QUREGS(getCachedDensmatrs(), refFunc, regularFunc, randFunc);
+        TEST_ON_CACHED_QUREGS(getCachedStatevecs(), refFunc, regularFunc, randFunc);
+
+        destroyPauliStrSum(sum);
+
+    }
+}
+
+/*
+* Time evolution tests
+* @todo Add Pauli permutation variants
+*/ 
+TEST_CASE( "applyTrotterizedUnitaryTimeEvolution", TEST_CATEGORY ) { 
+
+    // BEWARE: this test creates a new Qureg below which will have
+    // deployments chosen by the auto-deployer; it is ergo unpredictable
+    // whether it will be multithreaded, GPU-accelerated or distributed.
+    // This test is ergo checking only a single, unspecified deployment,
+    // unlike other tests which check all deployments. This is tolerable
+    // since (non-randomised) Trotterisation is merely invoking routines
+    // (Pauli gadgets) already independently tested across deployments
+
+    SECTION( LABEL_CORRECTNESS ) {
+
+        int numQubits = 20;
+        Qureg qureg = createQureg(numQubits);
+        initPlusState(qureg);
+        bool permutePaulis = false;
+        
+        PauliStrSum hamil = createHeisenbergHamiltonian(numQubits);
+        PauliStrSum observ = createAlternatingPauliObservable(numQubits);
+        
+        qreal dt = 0.1;
+        int order = 4;
+        int reps = 5;
+        int steps = 10;
+        
+        // nudge the epsilon used by internal validation functions up a bit
+        // as the time evolution operation plays badly with single precision
+        // Defaults for validation epsilon are:
+        //  - 1E-5 at single precision
+        //  - 1E-12 at double precision
+        //  - 1E-13 at quad precision
+        qreal initialValidationEps = getValidationEpsilon();
+        setValidationEpsilon(2 * initialValidationEps);
+
+        /*
+        * Tolerance for floating-point comparisons
+        * Note that the underlying numerics are sensitive to the float
+        * precision AND to the number of threads. As such we set quite 
+        * large epsilon values to account for the worst-case scenario which 
+        * is single precision, single thread. The baseline for these results
+        * is double precision, multiple threads.
+        *
+        * Values (assuming default initialValidationEps) are:
+        * Single precision:
+        *   obsEps = 0.03
+        *   normEps = 0.001
+        *
+        * Double precision:
+        *   obsEps = 3E-9
+        *   normEps = 1E-10
+        *
+        * Quad precision:
+        *   obsEps = 3E-10
+        *   normEps = 1E-11
+        */
+        qreal obsEps = 3E3 * initialValidationEps;
+        qreal normEps = 100 * initialValidationEps;
+       
+        vector<qreal> refObservables = {
+            19.26827777028073,
+            20.34277275871839,
+            21.21120737889526,
+            21.86585902741717,
+            22.30371711358924,
+            22.52644660547882,
+            22.54015748825067,
+            22.35499202583118,
+            21.9845541501027,
+            21.44521638719462
+        };
+        
+        for (int i = 0; i < steps; i++) {
+            applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps, permutePaulis);
+            qreal expec = calcExpecPauliStrSum(qureg, observ);
+            
+            REQUIRE_THAT( expec, WithinAbs(refObservables[i], obsEps) );
+        }
+        
+        // Verify state remains normalized
+        REQUIRE_THAT( calcTotalProb(qureg), WithinAbs(1.0, normEps) );
+
+        // Restore validation epsilon
+        setValidationEpsilon(initialValidationEps);
+
+        destroyQureg(qureg);
+        destroyPauliStrSum(hamil);
+        destroyPauliStrSum(observ);
+    }
+
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        PauliStrSum hamil = createHeisenbergHamiltonian(qureg.numQubits);
+        bool permutePaulis = false;
 
-void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps);
+        SECTION( "qureg uninitialised" ) {
 
-void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps);
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH( 
+                applyTrotterizedUnitaryTimeEvolution(badQureg, hamil, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "pauli sum uninitialized" ) {
+
+            PauliStrSum badHamil = hamil;
+            badHamil.numTerms = 0;
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, badHamil, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("Pauli")
+            );
+        }
+
+        SECTION( "hamiltonian not hermitian" ) {
+
+            vector<PauliStr> strings;
+            vector<qcomp> coeffs;
+            strings.push_back(getPauliStr("X", {0}));
+            coeffs.push_back(getQcomp(1.0, 1.0));  
+            PauliStrSum nonHermitian = createPauliStrSum(strings, coeffs);
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, nonHermitian, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("Hermitian")
+            );
+            destroyPauliStrSum(nonHermitian);
+        }
+
+        SECTION( "pauli sum exceeds qureg qubits" ) {
+
+            PauliStrSum largeHamil = createHeisenbergHamiltonian(qureg.numQubits + 1);
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, largeHamil, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("only compatible")
+            );
+            destroyPauliStrSum(largeHamil);
+        }
+
+        SECTION( "invalid trotter order (zero)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, hamil, 0.1, 0, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter order (negative)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, hamil, 0.1, -2, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter order (odd, not 1)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, hamil, 0.1, 3, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter reps (zero)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, hamil, 0.1, 4, 0, permutePaulis),
+                ContainsSubstring("repetitions")
+            );
+        }
+
+        SECTION( "invalid trotter reps (negative)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedUnitaryTimeEvolution(qureg, hamil, 0.1, 4, -3, permutePaulis),
+                ContainsSubstring("repetitions")
+            );
+        }
+
+        destroyPauliStrSum(hamil);
+    }
+}
+
+
+TEST_CASE( "applyTrotterizedImaginaryTimeEvolution", TEST_CATEGORY ) {
+
+    // BEWARE: this test creates a new Qureg below which will have
+    // deployments chosen by the auto-deployer; it is ergo unpredictable
+    // whether it will be multithreaded, GPU-accelerated or distributed.
+    // This test is ergo checking only a single, unspecified deployment,
+    // unlike other tests which check all deployments. This is tolerable
+    // since (non-randomised) Trotterisation is merely invoking routines
+    // (Pauli gadgets) already independently tested across deployments
+
+    SECTION( LABEL_CORRECTNESS ) {
+           
+        int numQubits = 16;
+        qreal tau = 0.1;
+        int order = 6;
+        int reps = 5;
+        int steps = 10;
+        bool permutePaulis = false;
+        
+        // Tolerance for ground state amplitude
+        qreal eps = 1E-2;
+        
+        // Ground state: all qubits align down (driven by strong magnetic field)
+        {
+            Qureg qureg = createQureg(numQubits);
+            initPlusState(qureg);
+            
+            PauliStrSum ising = createIsingHamiltonian(numQubits, 10.0, 0.0, 0.0);
+            
+            for (int i = 0; i < steps; ++i) {
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis);
+                setQuregToRenormalized(qureg);
+            }
+            
+            qcomp amp = getQuregAmp(qureg, 0);
+            qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag();
+            
+            REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) );
+            
+            for (long long i = 1; i < (1LL << numQubits); i++) {
+                qcomp other_amp = getQuregAmp(qureg, i);
+                qreal other_mag = other_amp.real() * other_amp.real() + 
+                                 other_amp.imag() * other_amp.imag();
+                REQUIRE( other_mag < eps );
+            }
+            
+            destroyQureg(qureg);
+            destroyPauliStrSum(ising);
+        }
+        
+        // Ground state: all qubits align up (driven by opposite magnetic field)
+        {
+            Qureg qureg = createQureg(numQubits);
+            initPlusState(qureg);
+            
+            PauliStrSum ising = createIsingHamiltonian(numQubits, -10.0, 0.0, 0.0);
+            
+            for (int i = 0; i < steps; ++i) {
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis);
+                setQuregToRenormalized(qureg);
+            }
+            
+            long long last_state = (1LL << numQubits) - 1;
+            qcomp amp = getQuregAmp(qureg, last_state);
+            qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag();
+            
+            REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) );
+            
+            for (long long i = 0; i < (1LL << numQubits); i++) {
+                if (i == last_state) continue;
+                qcomp other_amp = getQuregAmp(qureg, i);
+                qreal other_mag = other_amp.real() * other_amp.real() + 
+                                 other_amp.imag() * other_amp.imag();
+                REQUIRE( other_mag < eps );
+            }
+            
+            destroyQureg(qureg);
+            destroyPauliStrSum(ising);
+        }
+        
+        // Ground state: all qubits align down (driven by ferromagnetic interactions and bias)
+        {
+            Qureg qureg = createQureg(numQubits);
+            initPlusState(qureg);
+            
+            PauliStrSum ising = createIsingHamiltonian(numQubits, 0.0, 10.0, 10.0);
+            
+            for (int i = 0; i < steps; ++i) {
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis);
+                setQuregToRenormalized(qureg);
+            }
+            
+            qcomp amp = getQuregAmp(qureg, 0);
+            qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag();
+            
+            REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) );
+            
+            for (long long i = 1; i < (1LL << numQubits); i++) {
+                qcomp other_amp = getQuregAmp(qureg, i);
+                qreal other_mag = other_amp.real() * other_amp.real() + 
+                                 other_amp.imag() * other_amp.imag();
+                REQUIRE( other_mag < eps );
+            }
+            
+            destroyQureg(qureg);
+            destroyPauliStrSum(ising);
+        }
+        
+        // Ground state: alternating pattern (driven by antiferromagnetic interactions)
+        {
+            Qureg qureg = createQureg(numQubits);
+            initPlusState(qureg);
+            
+            PauliStrSum ising = createIsingHamiltonian(numQubits, 0.0, -10.0, 10.0);
+            
+            for (int i = 0; i < steps; ++i) {
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis);
+                setQuregToRenormalized(qureg);
+            }
+            
+            unsigned long long idx = 0;
+            for (int i = 0; i < numQubits / 2; ++i) {
+                idx += (1ULL << (2*i + 1));
+            }
+            
+            qcomp amp = getQuregAmp(qureg, idx);
+            qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag();
+            
+            REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) );
+            
+            for (long long i = 0; i < (1LL << numQubits); i++) {
+                if (i == idx) continue;
+                qcomp other_amp = getQuregAmp(qureg, i);
+                qreal other_mag = other_amp.real() * other_amp.real() + 
+                                 other_amp.imag() * other_amp.imag();
+                REQUIRE( other_mag < eps );
+            }
+            
+            destroyQureg(qureg);
+            destroyPauliStrSum(ising);
+        }
+    }
+
+    SECTION( LABEL_VALIDATION ) {
+
+        Qureg qureg = getArbitraryCachedStatevec();
+        PauliStrSum ising = createIsingHamiltonian(qureg.numQubits, 1.0, 1.0, 0.0);
+        bool permutePaulis = false;
+
+        SECTION( "qureg uninitialised" ) {
+
+            Qureg badQureg = qureg;
+            badQureg.numQubits = -1;
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(badQureg, ising, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("invalid Qureg")
+            );
+        }
+
+        SECTION( "pauli sum uninitialized" ) {
+
+            PauliStrSum badIsing = ising;
+            badIsing.numTerms = 0;
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, badIsing, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("Pauli")
+            );
+        }
+
+        SECTION( "pauli sum exceeds qureg qubits" ) {
+
+            PauliStrSum largeIsing = createIsingHamiltonian(qureg.numQubits+1, 1.0, 1.0, 0.0);
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, largeIsing, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("only compatible")
+            );
+            destroyPauliStrSum(largeIsing);
+        }
+
+        SECTION( "hamiltonian not hermitian" ) {
+
+            vector<PauliStr> strings;
+            vector<qcomp> coeffs;
+            strings.push_back(getPauliStr("X", {0}));
+            coeffs.push_back(getQcomp(1.0, 1.0));  
+            PauliStrSum nonHermitian = createPauliStrSum(strings, coeffs);
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, nonHermitian, 0.1, 4, 5, permutePaulis),
+                ContainsSubstring("Hermitian")
+            );
+            destroyPauliStrSum(nonHermitian);
+        }
+
+        SECTION( "invalid trotter order (zero)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, 0.1, 0, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter order (negative)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, 0.1, -2, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter order (odd, not 1)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, 0.1, 3, 5, permutePaulis),
+                ContainsSubstring("order")
+            );
+        }
+
+        SECTION( "invalid trotter reps (zero)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, 0.1, 4, 0, permutePaulis),
+                ContainsSubstring("repetitions")
+            );
+        }
+
+        SECTION( "invalid trotter reps (negative)" ) {
+
+            REQUIRE_THROWS_WITH(
+                applyTrotterizedImaginaryTimeEvolution(qureg, ising, 0.1, 4, -3, permutePaulis),
+                ContainsSubstring("repetitions")
+            );
+        }
+
+        destroyPauliStrSum(ising);
+    }
+}
+
+
+/**
+ * @todo
+ * UNTESTED FUNCTIONS (NOT YET VALIDATED BY REFERENCE TESTS)
+ */
 
-void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedNonUnitaryPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps, bool permutePaulis);
 
-void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
-void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps);
+void applyTrotterizedControlledPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
-void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps);
+void applyTrotterizedMultiControlledPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
-void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps);
+void applyTrotterizedMultiStateControlledPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps, bool permutePaulis);
 
-void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps);
+void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps, bool permutePaulis);
diff --git a/tests/utils/cache.cpp b/tests/utils/cache.cpp
index ec5001bc..0cb8a603 100644
--- a/tests/utils/cache.cpp
+++ b/tests/utils/cache.cpp
@@ -168,6 +168,22 @@ quregCache getAltCachedDensmatrs() {
     return densmatrs2;
 }
 
+Qureg getArbitraryCachedStatevec() {
+
+    // must not be called pre-creation nor post-destruction
+    DEMAND( !statevecs1.empty() );
+
+    return statevecs1.begin()->second;
+}
+
+Qureg getArbitraryCachedDensmatr() {
+
+    // must not be called pre-creation nor post-destruction
+    DEMAND( !densmatrs1.empty() );
+
+    return densmatrs1.begin()->second;
+}
+
 
 
 /*
diff --git a/tests/utils/cache.hpp b/tests/utils/cache.hpp
index ec8c71c7..87d8382f 100644
--- a/tests/utils/cache.hpp
+++ b/tests/utils/cache.hpp
@@ -40,6 +40,9 @@ quregCache getCachedDensmatrs();
 quregCache getAltCachedStatevecs();
 quregCache getAltCachedDensmatrs();
 
+Qureg getArbitraryCachedStatevec();
+Qureg getArbitraryCachedDensmatr();
+
 qvector getRefStatevec();
 qmatrix getRefDensmatr();
 
diff --git a/tests/utils/linalg.cpp b/tests/utils/linalg.cpp
index c7477521..e8285cdf 100644
--- a/tests/utils/linalg.cpp
+++ b/tests/utils/linalg.cpp
@@ -157,7 +157,7 @@ qvector getNormalised(qvector vec) {
 }
 
 
-qvector getDisceteFourierTransform(qvector in) {
+qvector getDiscreteFourierTransform(qvector in, bool inverse) {
     DEMAND( in.size() > 0 );
     
     size_t dim = in.size();
@@ -166,7 +166,7 @@ qvector getDisceteFourierTransform(qvector in) {
     // PI must be accurate here
     qreal pi = 3.14159265358979323846;
     qreal a = 1 / std::sqrt(dim);
-    qreal b = 2 * pi / dim;
+    qreal b = (inverse ? -1 : 1) * 2 * pi / dim;
     
     for (size_t x=0; x<dim; x++)
         for (size_t y=0; y<dim; y++)
@@ -176,7 +176,7 @@ qvector getDisceteFourierTransform(qvector in) {
 }
 
 
-qvector getDisceteFourierTransform(qvector in, vector<int> targs) {
+qvector getDiscreteFourierTransform(qvector in, vector<int> targs, bool inverse) {
     DEMAND( in.size() > 0 );
 
     size_t dim = in.size();
@@ -185,7 +185,7 @@ qvector getDisceteFourierTransform(qvector in, vector<int> targs) {
     qindex len = getPow2(targs.size());
     qreal pi = 3.14159265358979323846;
     qreal a = 1 / std::sqrt(len);
-    qreal b = 2 * pi / len;
+    qreal b = (inverse ? -1 : 1) * 2 * pi / len;
 
     for (size_t i=0; i<dim; i++) {
         size_t x = getBitsAt(i, targs);
diff --git a/tests/utils/linalg.hpp b/tests/utils/linalg.hpp
index c761b5fe..18555bf6 100644
--- a/tests/utils/linalg.hpp
+++ b/tests/utils/linalg.hpp
@@ -32,8 +32,8 @@ qindex getPow2(int);
 qreal getSum(vector<qreal> vec);
 qcomp getSum(qvector);
 qvector getNormalised(qvector);
-qvector getDisceteFourierTransform(qvector);
-qvector getDisceteFourierTransform(qvector in, vector<int> targs);
+qvector getDiscreteFourierTransform(qvector in, bool inverse);
+qvector getDiscreteFourierTransform(qvector in, vector<int> targs, bool inverse);
 
 qcomp   getInnerProduct(qvector bra, qvector ket);
 qmatrix getOuterProduct(qvector ket, qvector bra);