diff --git a/.github/workflows/cpp_build_with_cmake.yml b/.github/workflows/cpp_build_with_cmake.yml
index d9670914..1b06cbfb 100644
--- a/.github/workflows/cpp_build_with_cmake.yml
+++ b/.github/workflows/cpp_build_with_cmake.yml
@@ -26,4 +26,4 @@ jobs:
       run: cmake --build build
 
     - name: unit test
-      run: ./build/bin/utest
+      run: ./build/bin/utest --gtest_death_test_style=threadsafe
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index e1b11480..8b4bee70 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -18,6 +18,10 @@ target_link_libraries(noisy_circuit_test.exe iqs)
 add_executable(test_of_custom_gates.exe test_of_custom_gates.cpp)
 target_link_libraries(test_of_custom_gates.exe iqs)
 
+add_executable(test_of_ncu_gates.exe test_of_ncu_gates.cpp)
+target_link_libraries(test_of_ncu_gates.exe iqs)
+
+
 ################################################################################
 
 set_target_properties( benchgates.exe
@@ -26,6 +30,7 @@ set_target_properties( benchgates.exe
                        heisenberg_dynamics.exe
                        noisy_circuit_test.exe
                        test_of_custom_gates.exe
+                       test_of_ncu_gates.exe
     PROPERTIES
     RUNTIME_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/bin"
 )
diff --git a/examples/test_of_ncu_gates.cpp b/examples/test_of_ncu_gates.cpp
new file mode 100644
index 00000000..796e107f
--- /dev/null
+++ b/examples/test_of_ncu_gates.cpp
@@ -0,0 +1,105 @@
+/**
+ * @file test_of_ncu_gates.cpp
+ * @author Lee J. O'Riordan (lee.oriordan@ichec.ie)
+ * @author Myles Doyle (myles.doyle@ichec.ie)
+ * @brief Application to show functionality of ApplyNCU gate call.
+ * @version 0.2
+ * @date 2020-06-12
+ */
+
+#include "../include/qureg.hpp"
+
+/////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////////////////
+
+int main(int argc, char **argv)
+{
+    unsigned myrank=0, nprocs=1;
+    qhipster::mpi::Environment env(argc, argv);
+    myrank = env.GetStateRank();
+    nprocs = qhipster::mpi::Environment::GetStateSize();
+    if (env.IsUsefulRank() == false) return 0;
+    int num_threads = 1;
+#ifdef _OPENMP
+#pragma omp parallel
+    {
+        num_threads = omp_get_num_threads();
+    }
+#endif
+
+    // number of qubits in compute register.
+    std::size_t num_qubits_compute = 4;
+    if(argc != 2){
+        fprintf(stderr, "usage: %s <num_qubits_compute> \n", argv[0]);
+        exit(1);
+    }
+    else{
+        num_qubits_compute = atoi(argv[1]);
+    }
+
+    // pauliX gate that will be applied in NCU
+    ComplexDP zero = {0.,0.};
+    ComplexDP one = {1.,0.};
+    TM2x2<ComplexDP> pauliX{{zero,one,one,zero}};
+
+    // Setup vector to store compute and auxiliary quantum register indices.
+    // |compute reg>|auxiliary reg>
+    // Set number of auxiliary qubits to equal number of controlqubits
+    // if there are more than 4 control qubits, else auxiliary register
+    // will be empty as it will not be used in the NCU routine for
+    // optimisations.
+    std::size_t num_qubits_control = num_qubits_compute - 1;
+    std::size_t num_qubits_auxiliary = (num_qubits_compute - 2) * (num_qubits_control > 4);
+    std::vector<std::size_t> reg_compute(num_qubits_compute);
+    std::vector<std::size_t> reg_auxiliary(num_qubits_auxiliary);
+
+    // Set qubit indices of registers
+    {
+        std::size_t qubit_index = 0;
+        for(std::size_t i = 0; i < num_qubits_compute; i++){
+            reg_compute[i] = qubit_index;
+            qubit_index++;
+        }
+        for(std::size_t i = 0; i < num_qubits_auxiliary; i++){
+            reg_auxiliary[i] = qubit_index;
+            qubit_index++;
+        }
+    }
+
+    // Set qubit indices for qubits acting as control
+    std::vector<std::size_t> control_ids(num_qubits_control);
+
+    // Set vector containing indices of the qubits acting as
+    // control for the NCU gate.
+    for(std::size_t i = 0; i < num_qubits_control; i++){
+        control_ids[i] = reg_compute[i];
+    }
+
+    // Set index of target qubit
+    std::size_t target_id = num_qubits_compute - 1;
+
+    QubitRegister<ComplexDP> psi(num_qubits_compute + num_qubits_auxiliary);
+    psi.Initialize("base", 0);
+
+    {
+        psi.EnableStatistics();
+   
+        // Apply a Hadamard gate to first num_qubits_compute-1
+        // qubits in the compute register.
+        for(std::size_t qubit_id = 0; qubit_id < num_qubits_compute-1; qubit_id++){
+            psi.ApplyHadamard(reg_compute[qubit_id]);
+        }
+
+        psi.Print("Before NCU");
+        psi.GetStatistics();
+
+        // Apply NCU
+        psi.ApplyNCU(pauliX, control_ids, reg_auxiliary, target_id);
+
+        // Observe only state with the first num_qubits_compute-1 
+        // qubits in the compute register set to 1 executes PauliX
+        // on the target qubit.
+        psi.Print("After NCU");
+    }
+}
diff --git a/include/GateCache.hpp b/include/GateCache.hpp
new file mode 100644
index 00000000..92040261
--- /dev/null
+++ b/include/GateCache.hpp
@@ -0,0 +1,147 @@
+/**
+ * @file Gates.hpp
+ * @author Lee J. O'Riordan (lee.oriordan@ichec.ie)
+ * @author Myles Doyle (myles.doyle@ichec.ie)
+ * @brief Gates wrapper/storage class. Adapted from QNLP.
+ * @version 0.2
+ * @date 2020-06-01
+ */
+
+#ifndef GATECACHE
+#define GATECACHE
+
+#include <unordered_map>
+
+#include <complex>
+#include <vector>
+#include <qureg.hpp>
+
+#include <mat_ops.hpp>
+
+#include "qureg.hpp"
+
+/**
+ * @brief Class to cache intermediate matrix values used within other parts of the computation. 
+ * Heavily depended upon by NCU to store sqrt matrix values following Barenco et al. (1995) decomposition.
+ * 
+ * @tparam SimulatorType The simulator type with SimulatorGeneral as base class
+ */
+template <class Type>
+class GateCache {
+    private:
+    //using GateType = TM2x2<Type>;
+    std::size_t cache_depth;
+
+    public:
+    GateCache() : cache_depth(0) { };
+
+    GateCache(std::size_t default_depth) : cache_depth(default_depth) {
+        initCache(cache_depth);
+    }
+
+    ~GateCache(){ clearCache(); }
+
+    //Take the 2x2 matrix type from the template SimulatorType
+    using GateType = TM2x2<Type>;
+
+    //Maintain a map for each gate label (X,Y,Z, etc.), and use vectors to store sqrt (indexed by 1/2^(i), and pairing matrix and adjoint)
+    std::unordered_map<std::string, std::vector< std::pair<GateType, GateType> > > gateCacheMap;
+    
+    void clearCache(){
+        gateCacheMap.clear();
+        cache_depth = 0;
+    }
+
+    /**
+     * @brief Initialise the gate cache with PauliX,Y,Z and H up to a given sqrt depth
+     * 
+     * @param sqrt_depth The depth to which calculate sqrt matrices and their respective adjoints
+     */
+    void initCache(const std::size_t sqrt_depth){
+        // If we do not have a sufficient circuit depth, clear and rebuild up to given depth.
+
+        if(cache_depth < sqrt_depth ){
+            gateCacheMap.clear();
+            cache_depth = 0;
+        }
+
+        if (gateCacheMap.empty()){
+            gateCacheMap["X"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_x(), adjointMatrix( construct_pauli_x() ) ) };
+            gateCacheMap["Y"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_y(), adjointMatrix( construct_pauli_y() ) ) };
+            gateCacheMap["Z"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_pauli_z(), adjointMatrix( construct_pauli_z() ) ) };
+            gateCacheMap["H"] = std::vector< std::pair<GateType, GateType> > { std::make_pair( construct_hadamard(), adjointMatrix( construct_hadamard() ) ) };
+
+            for( std::size_t depth = 1; depth <= sqrt_depth; depth++ ){
+                for( auto& kv : gateCacheMap ){
+                    kv.second.reserve(sqrt_depth + 1);
+                    auto m = matrixSqrt<GateType>(kv.second[depth-1].first);
+                    kv.second.emplace(kv.second.begin() + depth, std::make_pair( m, adjointMatrix( m ) ) );
+                }
+            }
+            cache_depth = sqrt_depth;
+        }
+    }
+
+    /**
+     * @brief Adds new gate to the cache up to a given sqrt depth
+     * 
+     * @param gateLabel Label of gate to index into map
+     * @param gate Gate matrix
+     * @param max_depth Depth of calculations for sqrt and associate adjoints
+     */
+    void addToCache(const std::string gateLabel, const GateType& gate, std::size_t max_depth){
+        if(max_depth <= cache_depth && gateCacheMap.find(gateLabel) != gateCacheMap.end() ){
+            return;
+        }
+        else if(max_depth > cache_depth){
+            initCache(max_depth);
+        }
+
+        std::vector< std::pair<GateType, GateType> > v;
+
+        v.reserve(max_depth + 1);
+        v.push_back(std::make_pair( gate, adjointMatrix( gate ) ) );
+        for( std::size_t depth = 1; depth <= max_depth; depth++ ){
+            auto m = matrixSqrt<GateType>( v[depth-1].first );
+            v.emplace(v.begin() + depth, std::make_pair( m, adjointMatrix( m ) ) );
+        }
+        gateCacheMap.emplace(std::make_pair(gateLabel, v) );
+    }
+
+    constexpr GateType construct_pauli_x(){
+        GateType px;
+        px(0, 0) = Type(0., 0.);
+        px(0, 1) = Type(1., 0.);
+        px(1, 0) = Type(1., 0.);
+        px(1, 1) = Type(0., 0.);
+        return px;
+    }
+    constexpr GateType construct_pauli_y(){
+        GateType py;
+        py(0, 0) = Type(0., 0.);
+        py(0, 1) = Type(0., -1.);
+        py(1, 0) = Type(0., 1.);
+        py(1, 1) = Type(0., 0.);
+        return py;
+    }
+    constexpr GateType construct_pauli_z(){
+        GateType pz;
+        pz(0, 0) = Type(1., 0.);
+        pz(0, 1) = Type(0., 0.);
+        pz(1, 0) = Type(0., 0.);
+        pz(1, 1) = Type(-1., 0.);
+        return pz;
+    }
+    constexpr GateType construct_hadamard(){
+        GateType h;
+        auto f = Type(1. / std::sqrt(2.), 0.0);
+        h(0, 0) = h(0, 1) = h(1, 0) = f;
+        h(1, 1) = -f;
+        return h;
+    }
+};
+
+template class GateCache<ComplexDP>;
+template class GateCache<ComplexSP>;
+
+#endif
\ No newline at end of file
diff --git a/include/mat_ops.hpp b/include/mat_ops.hpp
new file mode 100644
index 00000000..4937e49e
--- /dev/null
+++ b/include/mat_ops.hpp
@@ -0,0 +1,65 @@
+
+/**
+ * @file mat_ops.hpp
+ * @author Lee J. O'Riordan (lee.oriordan@ichec.ie)
+ * @brief Templated methods to manipulate small matrices. Adapted from QNLP.
+ * @version 0.2
+ * @date 2020-06-01
+ */
+
+#ifndef MAT_OPS
+#define MAT_OPS
+
+#include <complex>
+#include <vector>
+#include <iostream>
+
+/**
+ * @brief Calculates the unitary matrix square root (U == VV, where V is returned)
+ * 
+ * @tparam Type ComplexDP or ComplexSP
+ * @param U Unitary matrix to be rooted
+ * @return Matrix V such that VV == U
+ */
+template <class Mat2x2Type>            
+const Mat2x2Type matrixSqrt(const Mat2x2Type& U){
+    Mat2x2Type V(U);
+    std::complex<double> delta = U(0,0)*U(1,1) - U(0,1)*U(1,0);
+    std::complex<double> tau = U(0,0) + U(1,1);
+    std::complex<double> s = sqrt(delta);
+    std::complex<double> t = sqrt(tau + 2.0*s);
+
+    //must be a way to vectorise these; TinyMatrix have a scale/shift option?
+    V(0,0) += s;
+    V(1,1) += s;
+    std::complex<double> scale_factor(1.,0.);
+    scale_factor /= t;
+    V(0,0) *= scale_factor; //(std::complex<double>(1.,0.)/t);
+    V(0,1) *= scale_factor; //(1/t);
+    V(1,0) *= scale_factor; //(1/t);
+    V(1,1) *= scale_factor; //(1/t);
+
+    return V;
+}
+
+/**
+ * @brief Function to calculate the adjoint of an input matrix
+ * 
+ * @tparam Type ComplexDP or ComplexSP
+ * @param U Unitary matrix to be adjointed
+ * @return qhipster::TinyMatrix<Type, 2, 2, 32> U^{\dagger}
+ */
+template <class Mat2x2Type>            
+Mat2x2Type adjointMatrix(const Mat2x2Type& U){
+    Mat2x2Type Uadjoint(U);
+    std::complex<double> tmp;
+    tmp = Uadjoint(0,1);
+    Uadjoint(0,1) = Uadjoint(1,0);
+    Uadjoint(1,0) = tmp;
+    Uadjoint(0,0) = std::conj(Uadjoint(0,0));
+    Uadjoint(0,1) = std::conj(Uadjoint(0,1));
+    Uadjoint(1,0) = std::conj(Uadjoint(1,0));
+    Uadjoint(1,1) = std::conj(Uadjoint(1,1));
+    return Uadjoint;
+}
+#endif
\ No newline at end of file
diff --git a/include/ncu.hpp b/include/ncu.hpp
new file mode 100644
index 00000000..efacdd5d
--- /dev/null
+++ b/include/ncu.hpp
@@ -0,0 +1,182 @@
+/**
+ * @file ncu.hpp
+ * @author Lee J. O'Riordan (lee.oriordan@ichec.ie)
+ * @author Myles Doyle (myles.doyle@ichec.ie)
+ * @brief Functions for applying n-qubit controlled U (unitary) gates. Adapted from QNLP.
+ * @version 0.2
+ * @date 2020-06-01
+ */
+
+#ifndef NCU_H
+#define NCU_H
+
+#include <complex>
+#include <vector>
+#include <iostream>
+
+#include <GateCache.hpp>
+#include <mat_ops.hpp>
+
+#include <qureg.hpp>
+
+namespace ncu {
+
+/**
+ * @brief Class definition for applying n-qubit controlled unitary operations.
+ * 
+ * @tparam SimulatorType Class Simulator Type
+ */
+template <class Type>
+class NCU {
+    private:
+    using Matrix2x2Type = TM2x2<Type>;
+    GateCache<Type> gate_cache;
+
+    protected:
+    static std::size_t num_gate_ops;
+
+    public:
+    /**
+     * @brief Construct a new NCU object
+     * 
+     */
+    NCU() {
+        gate_cache = GateCache<Type>();
+    };
+
+    /**
+     * @brief Destroy the NCU object
+     * 
+     */
+    ~NCU(){
+        clearMaps();
+        gate_cache.clearCache();
+    };
+
+    /**
+     * @brief Add the PauliX and the given unitary U to the maps
+     * 
+     * @param U 
+     */
+    void initialiseMaps( std::size_t num_ctrl_lines){
+        gate_cache.initCache( num_ctrl_lines );
+    }
+
+    /**
+     * @brief Add the given unitary matrix to the maps up to the required depth
+     * 
+     * @param U 
+     */
+    void addToMaps( std::string U_label, const Matrix2x2Type& U, std::size_t num_ctrl_lines){
+        gate_cache.addToCache( U_label, U, num_ctrl_lines);
+    }
+
+    /**
+     * @brief Get the Map of cached gates. Keys are strings, and values are vectors of paired (gate, gate adjoint) types where the index give the value of (gate)^(1/2^i)
+     * 
+     * @return GateCache<SimulatorType> type 
+     */
+    GateCache<Type>& getGateCache(){
+        return gate_cache;
+    }
+
+    /**
+     * @brief Clears the maps of stored sqrt matrices
+     * 
+     */
+    void clearMaps(){
+        gate_cache.clearCache();
+    }
+
+    /**
+     * @brief Decompose n-qubit controlled op into 1 and 2 qubit gates. Control indices can be in any specified location. Ensure the gate cache has been populated with the appropriate gate type before running. This avoids O(n) checking of the container at each call for the associated gates.
+     * 
+     * @tparam Type ComplexDP or ComplexSP 
+     * @param qReg Qubit register
+     * @param ctrlIndices Vector of indices for control lines
+     * @param qTarget Target qubit for the unitary matrix U
+     * @param U Unitary matrix, U
+     * @param depth Depth of recursion.
+     */
+    void applyNQubitControl(QubitRegister<Type>& qReg, 
+            const std::vector<std::size_t> ctrlIndices,
+            const std::vector<std::size_t> auxIndices,
+            const unsigned int qTarget,
+            const std::string gateLabel,
+            const Matrix2x2Type& U,
+            const std::size_t depth
+    ){
+        int local_depth = depth + 1;
+
+        //Determine the range over which the qubits exist; consider as a count of the control ops, hence +1 since extremeties included
+        std::size_t cOps = ctrlIndices.size();
+
+        // Assuming given a set of auxiliary qubits, utilise the qubits for better depth optimisation.
+        if( (cOps >= 5) && ( auxIndices.size() >= cOps-2 ) && (gateLabel == "X") && (depth == 0) ){ //161 -> 60 2-qubit gate calls
+            qReg.ApplyToffoli( ctrlIndices.back(), *(auxIndices.begin() + ctrlIndices.size() - 3), qTarget);
+
+            for (std::size_t i = ctrlIndices.size()-2; i >= 2; i--){
+                qReg.ApplyToffoli( *(ctrlIndices.begin()+i), *(auxIndices.begin() + (i-2)), *(auxIndices.begin() + (i-1)));
+            }
+            qReg.ApplyToffoli( *(ctrlIndices.begin()), *(ctrlIndices.begin()+1), *(auxIndices.begin()) );
+            
+            for (std::size_t i = 2; i <= ctrlIndices.size()-2; i++){
+                qReg.ApplyToffoli( *(ctrlIndices.begin()+i), *(auxIndices.begin()+(i-2)), *(auxIndices.begin()+(i-1)));
+            }
+            qReg.ApplyToffoli( ctrlIndices.back(), *(auxIndices.begin() + ctrlIndices.size() - 3), qTarget);
+
+            for (std::size_t i = ctrlIndices.size()-2; i >= 2; i--){
+                qReg.ApplyToffoli( *(ctrlIndices.begin()+i), *(auxIndices.begin() + (i-2)), *(auxIndices.begin() + (i-1)));
+            }
+            qReg.ApplyToffoli( *(ctrlIndices.begin()), *(ctrlIndices.begin()+1), *(auxIndices.begin()) );
+            for (std::size_t i = 2; i <= ctrlIndices.size()-2; i++){
+                qReg.ApplyToffoli( *(ctrlIndices.begin()+i), *(auxIndices.begin()+(i-2)), *(auxIndices.begin()+(i-1)));
+            }
+        }
+        // Optimisation for replacing 17 2-qubit with 13 2-qubit gate calls
+        else if(cOps == 3){ 
+            //Apply the 13 2-qubit gate calls
+            qReg.ApplyControlled1QubitGate( ctrlIndices[0], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].first );
+            qReg.ApplyCPauliX( ctrlIndices[0], ctrlIndices[1]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[1], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].second );
+            qReg.ApplyCPauliX( ctrlIndices[0], ctrlIndices[1]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[1], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].first );
+            qReg.ApplyCPauliX( ctrlIndices[1], ctrlIndices[2]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[2], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].second );
+            qReg.ApplyCPauliX( ctrlIndices[0], ctrlIndices[2]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[2], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].first );
+            qReg.ApplyCPauliX( ctrlIndices[1], ctrlIndices[2]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[2], qTarget, gate_cache.gateCacheMap[gateLabel][local_depth+1].second );
+            qReg.ApplyCPauliX( ctrlIndices[0], ctrlIndices[2]);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices[2], qTarget,  gate_cache.gateCacheMap[gateLabel][local_depth+1].first );
+        }
+        // Default quadratic decomposition from Barenco et al (1995)
+        else if (cOps >= 2 && cOps !=3){
+            std::vector<std::size_t> subCtrlIndices(ctrlIndices.begin(), ctrlIndices.end()-1);
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices.back(), qTarget, gate_cache.gateCacheMap[gateLabel][local_depth].first );
+
+            applyNQubitControl(qReg, subCtrlIndices, auxIndices, ctrlIndices.back(), "X",  gate_cache.gateCacheMap["X"][0].first, 0 );
+
+            qReg.ApplyControlled1QubitGate( ctrlIndices.back(), qTarget, gate_cache.gateCacheMap[gateLabel][local_depth].second );
+
+            applyNQubitControl(qReg, subCtrlIndices, auxIndices, ctrlIndices.back(), "X", gate_cache.gateCacheMap["X"][0].first, 0 );
+
+            applyNQubitControl(qReg, subCtrlIndices, auxIndices, qTarget, gateLabel, gate_cache.gateCacheMap[gateLabel][local_depth+1].first, local_depth );
+        }
+
+        // If the number of control qubits is less than 2, assume we have decomposed sufficiently
+        else{
+            qReg.ApplyControlled1QubitGate( ctrlIndices[0], qTarget, gate_cache.gateCacheMap[gateLabel][depth].first );
+        }
+    }
+};
+
+};
+#endif
\ No newline at end of file
diff --git a/include/qureg.hpp b/include/qureg.hpp
index 2901dad5..8b8ae94e 100644
--- a/include/qureg.hpp
+++ b/include/qureg.hpp
@@ -292,7 +292,10 @@ class QubitRegister
   void ApplyCHadamard(unsigned const control_qubit, unsigned const target_qubit);
 
   void ApplyCPhaseRotation(unsigned const qubit, unsigned const qubit2, BaseType theta);
-  
+
+  // Algorithms
+  void ApplyNCU(TM2x2<Type> gate, const std::vector<std::size_t>& ctrl_indices, const std::vector<std::size_t>& aux_indices, unsigned const target);
+
   // fusion  
   void TurnOnFusion(unsigned log2llc = 20);
   void TurnOffFusion();
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8e2ad1a3..449aa9f0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,6 +14,7 @@ set(IQS_FILES
   qureg_expectval.cpp
   qureg_fusion.cpp
   qureg_init.cpp
+  qureg_ncu.cpp
   qureg_measure.cpp
   qureg_noisysimul.cpp
   qureg_permute.cpp
@@ -27,7 +28,7 @@ set(IQS_FILES
 
 add_library(iqs SHARED ${IQS_FILES})
 
-target_include_directories(iqs PUBLIC ../include)
+target_include_directories(iqs PUBLIC ${CMAKE_SOURCE_DIR}/include)
 
 target_include_directories(iqs INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}>)
 
diff --git a/src/qureg_ncu.cpp b/src/qureg_ncu.cpp
new file mode 100644
index 00000000..e2363530
--- /dev/null
+++ b/src/qureg_ncu.cpp
@@ -0,0 +1,46 @@
+#include "qureg.hpp"
+#include "ncu.hpp"
+#include "GateCache.hpp"
+#include "mat_ops.hpp"
+#include <iostream>
+
+template <class Type>
+void QubitRegister<Type>::ApplyNCU(
+                TM2x2<Type> gate,
+                const std::vector<std::size_t>& ctrl_indices, 
+                const std::vector<std::size_t>& aux_indices, 
+                unsigned const target)
+{
+    ncu::NCU<Type> ncu;
+    ncu.initialiseMaps(ctrl_indices.size());
+
+    auto gcache = ncu.getGateCache();
+
+    // initial support for Pauli gates only
+    std::string gate_label = "";
+    if( gate == gcache.construct_pauli_x() ){
+        gate_label = "X";
+    }
+    else if( gate == gcache.construct_pauli_y() ){
+        gate_label = "Y";
+    }
+    else if( gate == gcache.construct_pauli_z() ){
+        gate_label = "Z";
+    }
+    else {
+        std::cerr << "Gate not currently support for NCU: " << gate.tostr() << std::endl;
+        std::abort();
+    }
+    ncu.applyNQubitControl(*this, ctrl_indices, aux_indices, target, gate_label, gate, 0);
+}
+
+template void QubitRegister<ComplexDP>::ApplyNCU(
+                TM2x2<ComplexDP> gate,
+                const std::vector<std::size_t>& ctrl_indices, 
+                const std::vector<std::size_t>& aux_indices, 
+                unsigned const target);
+template void QubitRegister<ComplexSP>::ApplyNCU(
+                TM2x2<ComplexSP> gate,
+                const std::vector<std::size_t>& ctrl_indices, 
+                const std::vector<std::size_t>& aux_indices, 
+                unsigned const target);
diff --git a/unit_test/include/apply_ncu_test.hpp b/unit_test/include/apply_ncu_test.hpp
new file mode 100644
index 00000000..6960431d
--- /dev/null
+++ b/unit_test/include/apply_ncu_test.hpp
@@ -0,0 +1,630 @@
+/**
+ * @file apply_ncu_test.hpp
+ * @author Lee J. O'Riordan (lee.oriordan@ichec.ie)
+ * @author Myles Doyle (myles.doyle@ichec.ie)
+ * @brief Test functionality of ApplyNCU gate call.
+ * @version 0.2
+ * @date 2020-06-12
+ */
+
+#ifndef APPLY_NCU_TEST_HPP
+#define APPLY_NCU_TEST_HPP
+
+#include "../../include/qureg.hpp"
+
+//////////////////////////////////////////////////////////////////////////////
+// Test fixture class.
+
+class ApplyNCUGateTest : public ::testing::Test
+{
+ protected:
+
+    ApplyNCUGateTest(){
+    }
+
+    // just after the 'constructor'
+    void SetUp() override{
+        // All tests are skipped if the rank is dummy.
+        if (qhipster::mpi::Environment::IsUsefulRank() == false){
+            GTEST_SKIP();
+        }
+
+        // All tests are skipped if the 4-qubit state is distributed in more than 2^3 ranks.
+        // In fact the MPI version needs to allocate half-the-local-storage for communication.
+        if (qhipster::mpi::Environment::GetStateSize() > 8){
+            GTEST_SKIP();
+        }
+    }
+
+    // Pauli gates that will be applied in NCU
+    ComplexDP zero = {0.,0.};
+    ComplexDP one_re = {1.,0.};
+    ComplexDP one_im = {0.,1.};
+    ComplexDP neg_one_re = {-1.,0.};
+    ComplexDP neg_one_im = {0.,-1.};
+    TM2x2<ComplexDP> pauliX_{{zero,one_re,one_re,zero}};
+    TM2x2<ComplexDP> pauliY_{{zero,neg_one_im,one_im,zero}};
+    TM2x2<ComplexDP> pauliZ_{{one_re,zero,zero,neg_one_re}};
+
+    double accepted_error_ = 1e-14;
+    double sqrt2_ = std::sqrt(2.);
+
+    // Standard NCU
+    const std::size_t num_qubits_ = 6;
+    const std::size_t num_qubits_control_ = 3;
+    size_t init0000_ = 0, init0001 = 8, init111000_ = 7, init111001_ = 39;
+
+    std::size_t target_index_ = num_qubits_-1;
+    std::vector<std::size_t> control_indices_ = {0,1,2};
+    std::vector<std::size_t> auxiliary_indices_empty_;
+
+    std::size_t num_states_ = pow(2, num_qubits_control_);
+    double percent_equal_dist_ = 1.0/(double) num_states_;
+    double amplitude_ = 1.0/std::sqrt(num_states_);
+
+    // NCU with optimisations
+    const std::size_t num_qubits_opt_ = 11;
+    const std::size_t num_qubits_control_opt_ = 6;
+    const std::size_t num_qubits_auxiliary_opt_ = 4;
+    size_t init11111100001_ = 63+1024, init11111100000_ = 63;
+
+    std::size_t target_index_opt_ = num_qubits_opt_ -1;
+    std::vector<std::size_t> control_indices_opt_ = {0,1,2,3,4,5};
+    std::vector<std::size_t> auxiliary_indices_opt_ = {6,7,8,9};
+
+    std::size_t num_states_opt_ = pow(2, num_qubits_control_opt_);
+    double percent_equal_dist_opt_ = 1.0/(double) num_states_opt_;
+    double amplitude_opt_ = 1.0/std::sqrt(num_states_opt_);
+};
+
+//////////////////////////////////////////////////////////////////////////////
+// For all 1-qubit gates we test the expected behavior:
+// - on the last qubit
+// - for initial states in orthogonal basis (X and Z eigensattes)
+//
+//
+//////////////////////////////////////////////////////////////////////////////
+// PauliX on |0>
+TEST_F(ApplyNCUGateTest, PauliXSingleState0){
+
+    // Apply to |111>|00>|0>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111000_);
+    psi.ApplyNCU(pauliX_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 1.,accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(),  1., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliX on |1>
+TEST_F(ApplyNCUGateTest, PauliXSingleState1){
+
+    // Apply to |111>|00>|1>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111001_);
+    psi.ApplyNCU(pauliX_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 0.,accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(),  1., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(),  0., accepted_error_);
+
+
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliY on |0>
+TEST_F(ApplyNCUGateTest, PauliYSingleState0){
+
+    // Apply to |111>|00>|0>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111000_);
+    psi.ApplyNCU(pauliY_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 1., accepted_error_);
+
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  1., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliY on |1>
+TEST_F(ApplyNCUGateTest, PauliYSingleState1){
+
+    // Apply to |111>|00>|1>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111001_);
+    psi.ApplyNCU(pauliY_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 0., accepted_error_);
+
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(),  -1., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliZ on |0>
+TEST_F(ApplyNCUGateTest, PauliZSingleState0){
+
+    // Apply to |111>|00>|0>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111000_);
+    psi.ApplyNCU(pauliZ_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 0., accepted_error_);
+
+    // Check amplitude_
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(),  1., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliZ on |1>
+TEST_F(ApplyNCUGateTest, PauliZSingleState1){
+
+    // Apply to |111>|00>|1>
+    //
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",init111001_);
+    psi.ApplyNCU(pauliZ_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure control qubits remain unchanged.
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        ASSERT_NEAR(psi.GetProbability(control_indices_[i]), 1., accepted_error_);
+    }
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 1., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(),  -1., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  0., accepted_error_);
+
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliX on superposition of states on |0>
+TEST_F(ApplyNCUGateTest, PauliXMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+
+    psi.ApplyNCU(pauliX_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), percent_equal_dist_,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(),  amplitude_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliY on superposition of states on |0>
+TEST_F(ApplyNCUGateTest, PauliYMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+
+    psi.ApplyNCU(pauliY_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), percent_equal_dist_,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(), 0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  amplitude_, accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliZ on superposition of states on |0>
+TEST_F(ApplyNCUGateTest, PauliZMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+
+    psi.ApplyNCU(pauliZ_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 0.,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(), amplitude_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(), 0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(), 0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliX on superposition of states on |1>
+TEST_F(ApplyNCUGateTest, PauliXMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_);
+
+    psi.ApplyNCU(pauliX_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(1. - psi.GetProbability(target_index_), percent_equal_dist_,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(),  amplitude_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliY on superposition of states on |1>
+TEST_F(ApplyNCUGateTest, PauliYMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_);
+
+    psi.ApplyNCU(pauliY_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(1. - psi.GetProbability(target_index_), percent_equal_dist_,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(), -amplitude_, accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(), 0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// PauliZ on superposition of states on |1>
+TEST_F(ApplyNCUGateTest, PauliZMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_; i++){
+        psi.ApplyHadamard(control_indices_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_);
+
+    psi.ApplyNCU(pauliZ_, control_indices_, auxiliary_indices_empty_, target_index_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_), 1.,accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).real(),  amplitude_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 32).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).real(), 0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111000_).imag(), 0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).real(), -amplitude_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init111001_).imag(), 0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// Test expcted Failure of NCU if arbitrary
+// U that is not a Pauli matrix is passed.
+TEST_F(ApplyNCUGateTest, FailArbitraryU)
+{
+
+	TM2x2<ComplexDP> U;
+	U(0, 0) = {0.592056606032915, 0.459533060553574}; 
+	U(0, 1) = {-0.314948020757856, -0.582328159830658};
+	U(1, 0) = {0.658235557641767, 0.070882241549507}; 
+	U(1, 1) = {0.649564427121402, 0.373855203932477};
+	// |psi> = |1110>
+	QubitRegister<ComplexDP> psi (num_qubits_,"base",init111000_);
+        ASSERT_DEATH(psi.ApplyNCU(U, control_indices_, auxiliary_indices_empty_, target_index_),"Gate not currently support for NCU:*" );
+}
+
+//////////////////////////////////////////////////////////////////////////////
+// NCU with Optimisations
+//
+// PauliX on superposition of states with NCU Optimised using
+// auxilairy register acting on |0> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliXMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+
+    psi.ApplyNCU(pauliX_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_opt_), percent_equal_dist_opt_, accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  amplitude_opt_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// PauliY on superposition of states with NCU Optimised using
+// auxilairy register acting on |0> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliYMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+
+    psi.ApplyNCU(pauliY_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    //psi.Print("Y");
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_opt_), percent_equal_dist_opt_, accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  amplitude_opt_, accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// PauliZ on superposition of states with NCU Optimised using
+// auxilairy register acting on |0> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliZMultiState0){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+
+    psi.ApplyNCU(pauliZ_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_opt_), 0., accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  amplitude_opt_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//
+/// PauliX on superposition of states with NCU Optimised using
+// auxilairy register acting on |1> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliXMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_opt_);
+
+    psi.ApplyNCU(pauliX_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(1. - psi.GetProbability(target_index_opt_), percent_equal_dist_opt_, accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  amplitude_opt_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// PauliY on superposition of states with NCU Optimised using
+// auxilairy register acting on |1> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliYMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_opt_);
+
+    psi.ApplyNCU(pauliY_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    //psi.Print("Y");
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(1. - psi.GetProbability(target_index_opt_), percent_equal_dist_opt_, accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  -amplitude_opt_, accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  0., accepted_error_);
+}
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// PauliZ on superposition of states with NCU Optimised using
+// auxilairy register acting on |1> 
+TEST_F(ApplyNCUGateTest, OptimisedNCUPauliZMultiState1){
+
+    QubitRegister<ComplexDP> psi (num_qubits_opt_,"base",0);
+    for(std::size_t i = 0; i < num_qubits_control_opt_; i++){
+        psi.ApplyHadamard(control_indices_opt_[i]);
+    }
+    // Flip target qubit to |1>
+    psi.ApplyPauliX(target_index_opt_);
+
+    psi.ApplyNCU(pauliZ_, control_indices_opt_, auxiliary_indices_opt_, target_index_opt_);
+
+    ASSERT_NEAR(psi.ComputeNorm(), 1., accepted_error_);
+
+    // Ensure target qubit was set correctly.
+    ASSERT_NEAR(psi.GetProbability(target_index_opt_), 1., accepted_error_);
+
+    // Check states without all control qubits set have 
+    // correct value.
+    for(std::size_t i = 0; i < num_states_opt_-1; i++){
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).real(),  amplitude_opt_, accepted_error_);
+        ASSERT_NEAR( psi.GetGlobalAmplitude(i + 1024).imag(),  0., accepted_error_);
+    }
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).real(),  0., accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100000_).imag(),  0., accepted_error_);
+
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).real(),  -amplitude_opt_, accepted_error_);
+    ASSERT_NEAR( psi.GetGlobalAmplitude(init11111100001_).imag(),  0., accepted_error_);
+}/////////////////////////////////////////////////////////////////////////////
+
+#endif	// header guard APPLY_NCU_TEST_HPP
diff --git a/unit_test/suite_of_tests.cpp b/unit_test/suite_of_tests.cpp
index d9d39f9d..2bc14a22 100644
--- a/unit_test/suite_of_tests.cpp
+++ b/unit_test/suite_of_tests.cpp
@@ -57,6 +57,7 @@ ASSERT_NEAR(val1.imag(),val2.imag(),error);
 // Utility methods and distributed implementation.
 #include "include/utility_methods_test.hpp"
 #include "include/chunking_communication_test.hpp"
+#include "include/apply_ncu_test.hpp"
 
 // Noisy simulations and extra features.
 #include "include/qureg_permute_test.hpp"