operations.cpp

/** @file
 * Unit tests of the operations module.
 *
 * @author Tyson Jones
 * 
 * @defgroup unitops Operations
 * @ingroup unittests
 */
 
#include "quest/include/quest.h"
 
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_template_test_macros.hpp>
#include <catch2/matchers/catch_matchers_string.hpp>
#include <catch2/generators/catch_generators_range.hpp>
 
#include "tests/utils/cache.hpp"
#include "tests/utils/qvector.hpp"
#include "tests/utils/qmatrix.hpp"
#include "tests/utils/compare.hpp"
#include "tests/utils/convert.hpp"
#include "tests/utils/evolve.hpp"
#include "tests/utils/linalg.hpp"
#include "tests/utils/lists.hpp"
#include "tests/utils/measure.hpp"
#include "tests/utils/macros.hpp"
#include "tests/utils/random.hpp"
 
#include <tuple>
 
using std::tuple;
 
 
 
/*
 * UTILITIES
 */
 
#define TEST_CATEGORY \
    LABEL_UNIT_TAG "[operations]"
 
 
/*
 * reference operator matrices used by testing
 */
 
namespace FixedMatrices {
 
    qmatrix H = {
        {1/std::sqrt(2),  1/std::sqrt(2)},
        {1/std::sqrt(2), -1/std::sqrt(2)}};
 
    qmatrix X = getPauliMatrix(1);
    qmatrix Y = getPauliMatrix(2);
    qmatrix Z = getPauliMatrix(3);
 
    qreal PI = 3.14159265358979323846;
    qmatrix T = {
        {1, 0},
        {0, std::exp(1_i * PI/4)}};
 
    qmatrix S = {
        {1, 0},
        {0, 1_i}};
 
    qmatrix SWAP = {
        {1, 0, 0, 0},
        {0, 0, 1, 0},
        {0, 1, 0, 0},
        {0, 0, 0, 1}};
 
    qmatrix sqrtSWAP = {
        {1, 0, 0, 0},
        {0, (1+1_i)/2, (1-1_i)/2, 0},
        {0, (1-1_i)/2, (1+1_i)/2, 0},
        {0, 0, 0, 1}};
}
 
namespace ParameterisedMatrices {
 
    auto Rx = [](qreal p) { return getExponentialOfPauliMatrix(p, FixedMatrices::X); };
    auto Ry = [](qreal p) { return getExponentialOfPauliMatrix(p, FixedMatrices::Y); };
    auto Rz = [](qreal p) { return getExponentialOfPauliMatrix(p, FixedMatrices::Z); };
 
    auto PS  = [](qreal p) { return qmatrix{{1, 0}, {0, std::exp(p*1_i)}}; };
    auto PS2 = [](qreal p) { return getControlledMatrix(PS(p), 1); };
}
 
namespace VariableSizeMatrices {
 
    auto X  = [](int n) { return getKroneckerProduct(FixedMatrices::X, n); };
    auto PF = [](int n) { return getControlledMatrix(FixedMatrices::Z, n - 1); };
}
 
namespace VariableSizeParameterisedMatrices {
 
    auto Z = [](qreal p, int n) { 
        qmatrix m = getKroneckerProduct(FixedMatrices::Z, n);
        return getExponentialOfPauliMatrix(p, m);
    };
 
    auto PS = [](qreal p, int n) {
        qmatrix m = ParameterisedMatrices::PS(p);
        return getControlledMatrix(m, n - 1);
    };
}
 
 
/*
 * execute 'function' upon each kind of cached qureg
 * (e.g. distributed, GPU-accelerated, etc) and a
 * reference state (T1 = qvector or qmatrix), when
 * both are initialised in the debug state, and
 * thereafter confirm they approximately agree. Each
 * qureg deployment is featured in a separate test
 * section, so are accounted distinctly.
 */
 
void TEST_ON_CACHED_QUREGS(quregCache quregs, auto& reference, auto& function) {
 
    for (auto& [label, qureg]: quregs) {
 
        DYNAMIC_SECTION( label ) {
 
            // no need to validate whether qureg successfully
            // enters the debug state here, because the below
            // serial setToDebugState() is guaranteed to succeed
            initDebugState(qureg);
            setToDebugState(reference);
 
            function(qureg, reference);
            REQUIRE_AGREE( qureg, reference );
        }
    }
}
 
void TEST_ON_CACHED_QUREG_AND_MATRIX(quregCache quregs, matrixCache matrices, auto apiFunc, auto refState, auto refMatr, auto refFunc) {
 
    for (auto& [labelA, qureg]: quregs) {
        for (auto& [labelB, matrix]: matrices) {
 
            // skip illegal (distributed matrix, local qureg) combo
            if (matrix.isDistributed && ! qureg.isDistributed)
                continue;
 
            DYNAMIC_SECTION( labelA + LABEL_DELIMITER + labelB ) {
 
                // set qureg and reference to debug
                initDebugState(qureg);
                setToDebugState(refState);
 
                // set API matrix to pre-initialised ref matrix
                setFullStateDiagMatr(matrix, 0, getDiagonals(refMatr));
 
                // API and reference functions should produce agreeing states
                apiFunc(qureg, matrix);
                refFunc(refState, refMatr);
                REQUIRE_AGREE( qureg, refState );
            }
        }
    }
}
 
 
/*
 * simply avoids boilerplate
 */
 
#define PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef ) \
    int numQubits = getNumCachedQubits(); \
    auto statevecQuregs = getCachedStatevecs(); \
    auto densmatrQuregs = getCachedDensmatrs(); \
    qvector statevecRef = getZeroVector(getPow2(numQubits)); \
    qmatrix densmatrRef = getZeroMatrix(getPow2(numQubits));
 
 
/*
 * Template flags for specifying what kind of additional 
 * arguments (in addition to ctrls/states/targs below) are
 * accepted by an API operation, when passing said operation
 * to automated testing facilities. This is NOT consulted
 * when generically invoking the API operation (we instead
 * used variadic templates above for that), but IS used
 * by the testing code to decide how to prepare inputs.
 * 
 * For example:
 * - applyHadamard:         none
 * - applyRotateX:          scalar
 * - applyRotateAroundAxis: axisrots
 * - applyDiagMatr1:        diagmatr
 * - applyDiagMatrPower:    diagpower
 * - applyCompMatr:         compmatr
 * - applyPauliStr:         paulistr
 * - applyPauliGadgt:       pauligad
*/
 
enum ArgsFlag { none, scalar, axisrots, diagmatr, diagpower, compmatr, paulistr, pauligad };
 
 
/*
 * Template flags for specifying how many control and
 * target qubits are accepted by an API operation.
 * Value 'anystates' is reserved for control qubits,
 * indicating that ctrls must accompany ctrl-states in
 * the API signature.
 * 
 * For example:
 * - applyHadamard:                Ctrls=zero, Targs=one
 * - applyControlledSwap:          Ctrls=one,  Targs=two
 * - applyMultiControlledCompMatr: Ctrls=any,  Targs=any
 * - applyMultiStateControlledT:   Ctrls=anystates, Targs=one
 */
 
enum NumQubitsFlag { zero, one, two, any, anystates };
 
void assertNumQubitsFlagsAreValid(NumQubitsFlag ctrlsFlag, NumQubitsFlag targsFlag) {
 
    DEMAND( 
        ctrlsFlag == zero ||
        ctrlsFlag == one  ||
        ctrlsFlag == any  ||
        ctrlsFlag == anystates );
 
    DEMAND(
        targsFlag == zero ||
        targsFlag == one  ||
        targsFlag == two  ||
        targsFlag == any);
}
 
void assertNumQubitsFlagsValid(
    NumQubitsFlag ctrlsFlag, NumQubitsFlag targsFlag, 
    vector<int> ctrls, vector<int> states, vector<int> targs
) {
    assertNumQubitsFlagsAreValid(ctrlsFlag, targsFlag);
 
    // we explicitly permit targsFlag=zero while
    // targs.size() != 0, which occurs when targs
    // are supplied to the API through an alternate
    // argument (e.g. a PauliStr)
 
    if (targsFlag == one) 
        DEMAND( targs.size() == 1 );
 
    if (targsFlag == two)
        DEMAND( targs.size() == 2 );
 
    if (ctrlsFlag == zero)
        DEMAND( ctrls.size() == 0 );
 
    if (ctrlsFlag == one)
        DEMAND( ctrls.size() == 1 );
 
    if (ctrlsFlag == anystates)
        DEMAND( states.size() == ctrls.size() );
    else
        DEMAND( states.size() == 0 );
}
 
 
/*
 * extract the runtime values of the number of control
 * and target qubtits from their compile-time templates.
 * When their number is permitted to be 'any', we use
 * a generator to successively test all possible numbers.
 * As such, this function is named similar to a Catch2
 * macro so the caller recognises it is a generator.
 */
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs, ArgsFlag Args>
int GENERATE_NUM_TARGS(int numQuregQubits) {
 
    assertNumQubitsFlagsAreValid(zero, Targs);
    DEMAND( Targs != one || numQuregQubits >= 1 );
    DEMAND( Targs != two || numQuregQubits >= 2 );
    DEMAND( numQuregQubits > 0 );
 
    // single choice if #targs is compile-time set
    if constexpr (Targs == one)
        return 1;
    if constexpr (Targs == two)
        return 2;
 
    if constexpr (Targs == any) {
 
        // we can target all qubits...
        int maxNumTargs = numQuregQubits;
 
        // unless we're applying CompMatr in distributed
        // mode. Technically we can support all targets
        // if the Qureg is a density matrix or not
        // distributed, but we safely choose the min here
        if (Args == compmatr)
            maxNumTargs = numQuregQubits - getLog2(getQuESTEnv().numNodes);
 
        // we must also ensure there is space left for a forced ctrl
        if (Ctrls == one && maxNumTargs == numQuregQubits)
            maxNumTargs -= 1;
 
        return GENERATE_COPY( range(1, maxNumTargs+1) );
    }
}
 
template <NumQubitsFlag Ctrls>
int GENERATE_NUM_CTRLS(int numFreeQubits) {
 
    assertNumQubitsFlagsAreValid(Ctrls, one);
    DEMAND( Ctrls != one || numFreeQubits >= 1 );
    DEMAND( numFreeQubits >= 0 );
 
    if constexpr (Ctrls == zero)
        return 0;
    
    if constexpr (Ctrls == one)
        return 1;
 
    if constexpr (Ctrls == any || Ctrls == anystates)
        return GENERATE_COPY( range(0, numFreeQubits+1) );
}
 
 
/*
 * invoke an API operation (e.g. applyHadamard), passing
 * any elements of (ctrls,states,targs) that it accepts 
 * (as informed by the template values) along with any
 * additional arguments. We explicitly accept numCtrls
 * and numTargs, rather than inferring them from ctrls
 * and targs, such that invalid numbers (like -1) can be
 * passed by input validation testing.
 * 
 * This big, ugly bespoke function is necessary, rather
 * than a simple variadic template, because the QuEST 
 * API accepts fixed numbers of qubits as individual 
 * arguments, rather than as lists/vectors/pointers. Note
 * that our use of variadic templates (args) means we do
 * not need to include ArgsFlag as a template parameter.
 */
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs>
void invokeApiOperation(
    auto operation, Qureg qureg, 
    vector<int> ctrls, vector<int> states, int numCtrls, 
    vector<int> targs, int numTargs, 
    auto&... args
) {
    assertNumQubitsFlagsValid(Ctrls, Targs, ctrls, states, targs);
 
    if constexpr (Ctrls == zero) {
        if constexpr (Targs == zero) operation(qureg, args...);
        if constexpr (Targs == one)  operation(qureg, targs[0], args...);
        if constexpr (Targs == two)  operation(qureg, targs[0], targs[1], args...);
        if constexpr (Targs == any)  operation(qureg, targs.data(), numTargs, args...);
    }
    if constexpr (Ctrls == one) {
        if constexpr (Targs == zero) operation(qureg, ctrls[0], args...);
        if constexpr (Targs == one)  operation(qureg, ctrls[0], targs[0], args...);
        if constexpr (Targs == two)  operation(qureg, ctrls[0], targs[0], targs[1], args...);
        if constexpr (Targs == any)  operation(qureg, ctrls[0], targs.data(), numTargs, args...);
    }
    if constexpr (Ctrls == any) {
        if constexpr (Targs == zero) operation(qureg, ctrls.data(), numCtrls, args...);
        if constexpr (Targs == one)  operation(qureg, ctrls.data(), numCtrls, targs[0], args...);
        if constexpr (Targs == two)  operation(qureg, ctrls.data(), numCtrls, targs[0], targs[1], args...);
        if constexpr (Targs == any)  operation(qureg, ctrls.data(), numCtrls, targs.data(), numTargs, args...);
    }
    if constexpr (Ctrls == anystates) {
        if constexpr (Targs == zero) operation(qureg, ctrls.data(), states.data(), numCtrls, args...);
        if constexpr (Targs == one)  operation(qureg, ctrls.data(), states.data(), numCtrls, targs[0], args...);
        if constexpr (Targs == two)  operation(qureg, ctrls.data(), states.data(), numCtrls, targs[0], targs[1], args...);
        if constexpr (Targs == any)  operation(qureg, ctrls.data(), states.data(), numCtrls, targs.data(), numTargs, args...);
    }
}
 
// overload to avoid passing numCtrls and numTargs
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs>
void invokeApiOperation(auto operation, Qureg qureg, vector<int> ctrls, vector<int> states, vector<int> targs, auto&... args) {
    invokeApiOperation<Ctrls,Targs>(operation, qureg, ctrls, states, ctrls.size(), targs, targs.size(), args...);
}
 
 
/*
 * prepare an API matrix (e.g. CompMatr1), as per
 * the given template parameters, with random
 * elements. This is used for testing API functions
 * which accept matrices.
 */
 
template <NumQubitsFlag Targs, ArgsFlag Args>
auto getRandomApiMatrix(int numTargs) {
 
    DEMAND(
        Args == diagmatr ||
        Args == diagpower ||
        Args == compmatr );
    
    qmatrix qm = (Args == compmatr)?
        getRandomUnitary(numTargs) : 
        getRandomDiagonalUnitary(numTargs);
 
    if constexpr (Args == compmatr && Targs == one) 
        return getCompMatr1(qm);
 
    if constexpr (Args == compmatr && Targs == two)
        return getCompMatr2(qm);
 
    if constexpr (Args == compmatr && Targs == any) {
        CompMatr cm = createCompMatr(numTargs); // must be freed
        setCompMatr(cm, qm);
        return cm;
    }
 
    qvector dv = getDiagonals(qm);
    constexpr bool diag = (Args == diagmatr || Args == diagpower);
 
    if constexpr (diag && Targs == one)
        return getDiagMatr1(dv);
 
    if constexpr (diag && Targs == two)
        return getDiagMatr2(dv);
 
    if constexpr (diag && Targs == any) {
        DiagMatr dm = createDiagMatr(numTargs); // must be freed
        setDiagMatr(dm, dv);
        return dm;
    }
}
 
 
/*
 * chooses random values for the remaining arguments
 * (after controls/states/targets) to API operations,
 * with types informed by the template parameters.
 * For example:
 * - applyRotateX() accepts a scalar
 * - applyCompMatr() accepts a CompMatr
 * - applyDiagMatrPower() accepts a DiagMatr and qcomp
 */
 
template <NumQubitsFlag Targs, ArgsFlag Args>
auto getRandomRemainingArgs(vector<int> targs) {
 
    if constexpr (Args == none)
        return tuple{ };
 
    if constexpr (Args == scalar) {
        qreal angle = getRandomPhase();
        return tuple{ angle };
    }
 
    if constexpr (Args == axisrots) {
        qreal angle = getRandomPhase();
        qreal x = getRandomReal(-1, 1);
        qreal y = getRandomReal(-1, 1);
        qreal z = getRandomReal(-1, 1);
        return tuple{ angle, x, y, z };
    }
 
    if constexpr (Args == compmatr || Args == diagmatr) {
        auto matrix = getRandomApiMatrix<Targs,Args>(targs.size()); // allocates heap mem
        return tuple{ matrix };
    }
 
    if constexpr (Args == diagpower) {
        DiagMatr matrix = getRandomApiMatrix<Targs,Args>(targs.size()); // allocates heap mem
        qcomp exponent = qcomp(getRandomReal(-3, 3), 0); // real for unitarity
        return tuple{ matrix, exponent };
    }
 
    if constexpr (Args == paulistr) {
        PauliStr str = getRandomPauliStr(targs);
        return tuple{ str };
    }
 
    if constexpr (Args == pauligad) {
        PauliStr str = getRandomPauliStr(targs);
        qreal angle = getRandomPhase();
        return tuple{ str, angle };
    }
}
 
 
template <NumQubitsFlag Targs, ArgsFlag Args>
void freeRemainingArgs(auto args) {
 
    if constexpr (Targs == any && Args == compmatr)
        destroyCompMatr(std::get<0>(args));
 
    if constexpr (Targs == any && Args == diagmatr)
        destroyDiagMatr(std::get<0>(args));
 
    if constexpr (Targs == any && Args == diagpower)
        destroyDiagMatr(std::get<0>(args));
}
 
 
/*
 * unpack the given reference operator matrix (a qmatrix) 
 * for an API operation, which is passed to testOperation(),
 * and which will be effected upon the reference state (a 
 * qvector or qmatrix). The type/form of matrixRefGen depends 
 * on the type of API operation, indicated by template parameter.
 */
 
template <NumQubitsFlag Targs, ArgsFlag Args>
qmatrix getReferenceMatrix(auto matrixRefGen, vector<int> targs, auto additionalArgs) {
 
    if constexpr (Args == none && Targs != any)
        return matrixRefGen;
 
    if constexpr (Args == none && Targs == any)
        return matrixRefGen(targs.size());
 
    if constexpr (Args == scalar && Targs != any) {
        qreal angle = std::get<0>(additionalArgs);
        return matrixRefGen(angle);
    }
 
    if constexpr (Args == scalar && Targs == any) {
        qreal angle = std::get<0>(additionalArgs);
        return matrixRefGen(angle, targs.size());
    }
 
    if constexpr (Args == axisrots) {
        qreal angle = std::get<0>(additionalArgs);
        qreal x = std::get<1>(additionalArgs);
        qreal y = std::get<2>(additionalArgs);
        qreal z = std::get<3>(additionalArgs);
        return getExponentialOfNormalisedPauliVector(angle, x, y, z);
    }
 
    if constexpr (Args == compmatr || Args == diagmatr) {
        auto apiMatrix = std::get<0>(additionalArgs);
        return getMatrix(apiMatrix);
    }
 
    if constexpr (Args == diagpower) {
        auto apiMatrix = std::get<0>(additionalArgs);
        qmatrix diag = getMatrix(apiMatrix);
        qcomp power = std::get<1>(additionalArgs);
        return getPowerOfDiagonalMatrix(diag, power);
    }
 
    if constexpr (Args == paulistr) {
        PauliStr str = std::get<0>(additionalArgs);
        return getMatrix(str, targs);
    }
 
    if constexpr (Args == pauligad) {
        PauliStr str = std::get<0>(additionalArgs);
        qreal angle  = std::get<1>(additionalArgs);
        qmatrix matr = getMatrix(str, targs);
        return getExponentialOfPauliMatrix(angle, matr);
    }
}
 
 
/*
 * display all/only relevant inputs given to an 
 * API operation when its subsequent test fails.
 * This is like a customisation of CAPTURE, although
 * we must use UNSCOPED_INFO (and ergo re-implement
 * some printers) because our branching makes scopes
 * which end CAPTURE's lifetime.
 */
 
 
// @todo surely this should live somewhere else,
// and/or re-use printer_ functions as much as possible
 
std::string toString(vector<int> list) {
 
    std::string out = "{ ";
    for (int& e : list)
        out += std::to_string(e) + " ";
    out += "}";
    return out;
}
 
extern int paulis_getPauliAt(PauliStr str, int ind);
 
std::string toString(PauliStr str, vector<int> targs) {
 
    std::string labels = "IXYZ";
 
    // ugly but adequate - like me (call me)
    std::string out = "";
    for (int i=targs.size()-1; i>=0; i--)
        out += labels[paulis_getPauliAt(str, targs[i])];
    
    return out;
}
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs, ArgsFlag Args>
void CAPTURE_RELEVANT( vector<int> ctrls, vector<int> states, vector<int> targs, auto& args ) {
 
    // display ctrls
    if (Ctrls == one)
        UNSCOPED_INFO( "control := " << ctrls[0] );
    if (Ctrls == any || Ctrls == anystates )
        UNSCOPED_INFO( "controls := " << toString(ctrls) );
 
    // display states
    if (Ctrls == anystates)
        UNSCOPED_INFO( "states := " << toString(states) );
 
    // display targs
    if (Targs == one)
        UNSCOPED_INFO( "target := " << targs[0] );
 
    if (Targs == two) {
        UNSCOPED_INFO( "target A := " << targs[0] );
        UNSCOPED_INFO( "target B := " << targs[1] );
    }
    if (Targs == any)
        UNSCOPED_INFO( "targets := " << toString(targs) );
 
    // display rotation angle
    if constexpr (Args == scalar)
        UNSCOPED_INFO( "angle := " << std::get<0>(args) );
 
    // display rotation angle and axis
    if constexpr (Args == axisrots) {
        UNSCOPED_INFO( "angle := " << std::get<0>(args) );
        UNSCOPED_INFO( "x := " << std::get<1>(args) );
        UNSCOPED_INFO( "y := " << std::get<2>(args) );
        UNSCOPED_INFO( "z := " << std::get<3>(args) );
    }
 
    // note but don't display API matrices
    if constexpr (Args == compmatr || Args == diagmatr || Args == diagpower)
        UNSCOPED_INFO( "matrix := (omitted)" );
 
    // display exponent of diagonal matrix
    if constexpr (Args == diagpower) {
        qcomp p = std::get<1>(args);
        UNSCOPED_INFO( "exponent := " << std::real(p) << " + (" << std::imag(p) << ")i" );
    }
 
    // display PauliStr
    if constexpr (Args == paulistr || Args == pauligad)
        UNSCOPED_INFO( "paulis := " << toString(std::get<0>(args), targs) );
 
    // display PauliStr angle
    if constexpr (Args == pauligad)
        UNSCOPED_INFO( "angle := " << std::get<1>(args) );
}
 
 
/*
 * perform a unit test on an API operation. The
 * template parameters are compile-time clues
 * about what inputs to prepare and pass to the
 * operation, and how its reference matrix (arg
 * matrixRefGen) is formatted.
 */
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs, ArgsFlag Args>
void testOperation(auto operation, auto matrixRefGen, bool multiplyOnly) {
 
    assertNumQubitsFlagsAreValid(Ctrls, Targs);
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        // try all possible number of ctrls and targs
        int numTargs = GENERATE_NUM_TARGS<Ctrls,Targs,Args>(numQubits);
        int numCtrls = GENERATE_NUM_CTRLS<Ctrls>(numQubits - numTargs);
        
        // either try all possible ctrls and targs, or randomise them
        auto ctrlsAndTargs = GENERATE_CTRLS_AND_TARGS( numQubits, numCtrls, numTargs );
        vector<int> ctrls = std::get<0>(ctrlsAndTargs);
        vector<int> targs = std::get<1>(ctrlsAndTargs);
 
        // randomise control states (if operation accepts them)
        vector<int> states = getRandomOutcomes(numCtrls * (Ctrls == anystates));
 
        // randomise remaining operation parameters
        auto primaryArgs = tuple{ ctrls, states, targs };
        auto furtherArgs = getRandomRemainingArgs<Targs,Args>(targs); // may allocate heap memory
 
        // obtain the reference matrix for this operation 
        qmatrix matrixRef = getReferenceMatrix<Targs,Args>(matrixRefGen, targs, furtherArgs);
 
        // PauliStr arg replaces target qubit list in API operations
        constexpr NumQubitsFlag RevTargs = (Args==paulistr||Args==pauligad)? zero : Targs;
 
        // disabling unitary-check validation for compmatr, since it's hard to
        // generate numerically-precise random unitaries upon many qubits, or
        // upon few qubits are single-precision. So we disable completely until
        // we re-implement 'input validation' checks which force us to fix thresholds
        (Args == compmatr)?
            setValidationEpsilon(0):
            setValidationEpsilonToDefault();
 
        // prepare test function which will receive both statevectors and density matrices
        auto testFunc = [&](Qureg qureg, auto& stateRef) -> void { 
            
            // invoke API operation, passing all args (unpacking variadic)
            auto apiFunc = [](auto&&... args) { return invokeApiOperation<Ctrls,RevTargs>(args...); };
            auto allArgs = std::tuple_cat(tuple{operation, qureg}, primaryArgs, furtherArgs);
            std::apply(apiFunc, allArgs);
 
            // update reference state
            (multiplyOnly)?
                multiplyReferenceOperator(stateRef, ctrls, states, targs, matrixRef):
                applyReferenceOperator(stateRef, ctrls, states, targs, matrixRef);
        };
 
        // report operation's input parameters if any subsequent test fails
        CAPTURE_RELEVANT<Ctrls,Targs,Args>( ctrls, states, targs, furtherArgs );
 
        // test API operation on all available deployment combinations (e.g. OMP, MPI, MPI+GPU, etc)
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
 
        // free any heap-alloated API matrices and restore epsilon
        freeRemainingArgs<Targs,Args>(furtherArgs);
        setValidationEpsilonToDefault();
    }
 
    // @todo
    // input validation!
}
 
template <NumQubitsFlag Ctrls, NumQubitsFlag Targs, ArgsFlag Args>
void testOperation(auto operation, auto matrixRefGen) {
 
    // default multiplyOnly=false
    testOperation<Ctrls,Targs,Args>(operation, matrixRefGen, false);
}
 
 
/*
 * perform unit tests for the four distinctly-controlled
 * variants of the given API operation (with the specified
 * function name suffix). 'numtargs' indicates the number 
 * of target qubits accepted by the operation, and 'argtype'
 * indicates the types of remaining arguments (if any exist).
 * 'matrixgen' is the matrix representation (of varying 
 * formats) of the operation, against which it will be compared.
 */
 
#define TEST_ANY_CTRL_OPERATION( namesuffix, numtargs, argtype, matrixgen ) \
    TEST_CASE( "apply" #namesuffix,                     TEST_CATEGORY ) { testOperation<zero,     numtargs,argtype>( apply ## namesuffix,                     matrixgen); } \
    TEST_CASE( "applyControlled" #namesuffix,           TEST_CATEGORY ) { testOperation<one,      numtargs,argtype>( applyControlled ## namesuffix,           matrixgen); } \
    TEST_CASE( "applyMultiControlled" #namesuffix,      TEST_CATEGORY ) { testOperation<any,      numtargs,argtype>( applyMultiControlled ## namesuffix,      matrixgen); } \
    TEST_CASE( "applyMultiStateControlled" #namesuffix, TEST_CATEGORY ) { testOperation<anystates,numtargs,argtype>( applyMultiStateControlled ## namesuffix, matrixgen); } 
 
 
 
/** 
 * TESTS
 * 
 * @ingroup unitops
 * @{
 */
 
TEST_ANY_CTRL_OPERATION( PauliStr,    any, paulistr, nullptr );
TEST_ANY_CTRL_OPERATION( PauliGadget, any, pauligad, nullptr );
TEST_ANY_CTRL_OPERATION( CompMatr1, one, compmatr, nullptr );
TEST_ANY_CTRL_OPERATION( CompMatr2, two, compmatr, nullptr );
TEST_ANY_CTRL_OPERATION( CompMatr,  any, compmatr, nullptr );
TEST_ANY_CTRL_OPERATION( DiagMatr1, one, diagmatr, nullptr );
TEST_ANY_CTRL_OPERATION( DiagMatr2, two, diagmatr, nullptr );
TEST_ANY_CTRL_OPERATION( DiagMatr,  any, diagmatr, nullptr );
TEST_ANY_CTRL_OPERATION( DiagMatrPower, any, diagpower, nullptr );
 
TEST_ANY_CTRL_OPERATION( Hadamard, one, none, FixedMatrices::H );
TEST_ANY_CTRL_OPERATION( PauliX,   one, none, FixedMatrices::X );
TEST_ANY_CTRL_OPERATION( PauliY,   one, none, FixedMatrices::Y );
TEST_ANY_CTRL_OPERATION( PauliZ,   one, none, FixedMatrices::Z );
TEST_ANY_CTRL_OPERATION( T,        one, none, FixedMatrices::T );
TEST_ANY_CTRL_OPERATION( S,        one, none, FixedMatrices::S );
TEST_ANY_CTRL_OPERATION( Swap,     two, none, FixedMatrices::SWAP );
TEST_ANY_CTRL_OPERATION( SqrtSwap, two, none, FixedMatrices::sqrtSWAP );
 
TEST_ANY_CTRL_OPERATION( RotateX, one, scalar, ParameterisedMatrices::Rx );
TEST_ANY_CTRL_OPERATION( RotateY, one, scalar, ParameterisedMatrices::Ry );
TEST_ANY_CTRL_OPERATION( RotateZ, one, scalar, ParameterisedMatrices::Rz );
 
TEST_ANY_CTRL_OPERATION( RotateAroundAxis, one, axisrots, nullptr );
 
TEST_ANY_CTRL_OPERATION( MultiQubitNot, any, none, VariableSizeMatrices::X );
 
TEST_ANY_CTRL_OPERATION( PhaseGadget, any, scalar, VariableSizeParameterisedMatrices::Z );
 
TEST_CASE( "multiplyPauliStr",     TEST_CATEGORY ) { testOperation<zero,any,paulistr>(multiplyPauliStr,    nullptr, true); }
TEST_CASE( "multiplyPauliGadget",  TEST_CATEGORY ) { testOperation<zero,any,pauligad>(multiplyPauliGadget, nullptr, true); }
 
TEST_CASE( "multiplyCompMatr1", TEST_CATEGORY ) { testOperation<zero,one,compmatr>(multiplyCompMatr1, nullptr, true); }
TEST_CASE( "multiplyCompMatr2", TEST_CATEGORY ) { testOperation<zero,two,compmatr>(multiplyCompMatr2, nullptr, true); }
TEST_CASE( "multiplyCompMatr",  TEST_CATEGORY ) { testOperation<zero,any,compmatr>(multiplyCompMatr,  nullptr, true); }
 
TEST_CASE( "multiplyDiagMatr1", TEST_CATEGORY ) { testOperation<zero,one,diagmatr>(multiplyDiagMatr1, nullptr, true); }
TEST_CASE( "multiplyDiagMatr2", TEST_CATEGORY ) { testOperation<zero,two,diagmatr>(multiplyDiagMatr2, nullptr, true); }
TEST_CASE( "multiplyDiagMatr",  TEST_CATEGORY ) { testOperation<zero,any,diagmatr>(multiplyDiagMatr,  nullptr, true); }
 
TEST_CASE( "multiplyDiagMatrPower",  TEST_CATEGORY ) { testOperation<zero,any,diagpower>(multiplyDiagMatrPower, nullptr, true); }
 
TEST_CASE( "multiplyMultiQubitNot", TEST_CATEGORY ) { testOperation<zero,any,none>(multiplyMultiQubitNot, VariableSizeMatrices::X, true); }
TEST_CASE( "multiplyPhaseGadget",   TEST_CATEGORY ) { testOperation<zero,any,scalar>(multiplyPhaseGadget, VariableSizeParameterisedMatrices::Z, true); }
 
TEST_CASE( "applyPhaseFlip",           TEST_CATEGORY ) { testOperation<zero,one,none>(applyPhaseFlip,           VariableSizeMatrices::PF(1)); }
TEST_CASE( "applyTwoQubitPhaseFlip",   TEST_CATEGORY ) { testOperation<zero,two,none>(applyTwoQubitPhaseFlip,   VariableSizeMatrices::PF(2)); }
TEST_CASE( "applyMultiQubitPhaseFlip", TEST_CATEGORY ) { testOperation<zero,any,none>(applyMultiQubitPhaseFlip, VariableSizeMatrices::PF); }
 
TEST_CASE( "applyPhaseShift",           TEST_CATEGORY ) { testOperation<zero,one,scalar>(applyPhaseShift,           ParameterisedMatrices::PS); }
TEST_CASE( "applyTwoQubitPhaseShift",   TEST_CATEGORY ) { testOperation<zero,two,scalar>(applyTwoQubitPhaseShift,   ParameterisedMatrices::PS2); }
TEST_CASE( "applyMultiQubitPhaseShift", TEST_CATEGORY ) { testOperation<zero,any,scalar>(applyMultiQubitPhaseShift, VariableSizeParameterisedMatrices::PS); }
 
 
TEST_CASE( "applyQuantumFourierTransform", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numTargs = GENERATE_COPY( range(1,numQubits+1) );
        auto targs = GENERATE_TARGS( numQubits, numTargs );
 
        CAPTURE( targs );
 
        SECTION( LABEL_STATEVEC ) { 
 
            auto testFunc = [&](Qureg qureg, qvector& ref) {
                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
                ref = getDisceteFourierTransform(ref, targs);
            };
 
            TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
        }
 
        SECTION( LABEL_DENSMATR ) { 
 
            // prepare a random mixture
            auto states = getRandomOrthonormalStateVectors(numQubits, getRandomInt(1,10));
            auto probs = getRandomProbabilities(states.size());
 
            auto testFunc = [&](Qureg qureg, qmatrix& ref) {
 
                // overwrite the Qureg debug state set by caller to above mixture
                setQuregToReference(qureg, getMixture(states, probs));
                applyQuantumFourierTransform(qureg, targs.data(), targs.size());
                
                ref = getZeroMatrix(ref.size());
                for (size_t i=0; i<states.size(); i++) {
                    qvector vec = getDisceteFourierTransform(states[i], targs);
                    ref += probs[i] * getOuterProduct(vec, vec);
                }
            };
 
            TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
        }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyFullQuantumFourierTransform", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        GENERATE( range(0,10) );
 
        SECTION( LABEL_STATEVEC ) { 
 
            auto testFunc = [&](Qureg qureg, qvector& ref) {
                applyFullQuantumFourierTransform(qureg);
                ref = getDisceteFourierTransform(ref);
            };
 
            TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc);
        }
 
        SECTION( LABEL_DENSMATR ) { 
 
            // prepare a random mixture
            auto states = getRandomOrthonormalStateVectors(numQubits, getRandomInt(1,10));
            auto probs = getRandomProbabilities(states.size());
 
            auto testFunc = [&](Qureg qureg, qmatrix& ref) {
 
                // overwrite the Qureg debug state set by caller to above mixture
                setQuregToReference(qureg, getMixture(states, probs));
                applyFullQuantumFourierTransform(qureg);
                
                ref = getZeroMatrix(ref.size());
                for (size_t i=0; i<states.size(); i++) {
                    qvector vec = getDisceteFourierTransform(states[i]);
                    ref += probs[i] * getOuterProduct(vec, vec);
                }
            };
 
            TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc);
        }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyQubitProjector", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        GENERATE( range(0,10) );
        int target = GENERATE_COPY( range(0,numQubits) );
        int outcome = GENERATE( 0, 1 );
 
        qmatrix projector = getProjector(outcome);
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
            applyQubitProjector(qureg, target, outcome);
            applyReferenceOperator(ref, {target}, projector);
        };
 
        CAPTURE( target, outcome );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyMultiQubitProjector", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numTargs = GENERATE_COPY( range(1,numQubits+1) );
        auto targets = GENERATE_TARGS( numQubits, numTargs );
        auto outcomes = getRandomOutcomes(numTargs);
 
        qmatrix projector = getProjector(targets, outcomes, numQubits);
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
            applyMultiQubitProjector(qureg, targets.data(), outcomes.data(), numTargs);
            applyReferenceOperator(ref, projector);
        };
 
        CAPTURE( targets, outcomes );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyForcedQubitMeasurement", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        GENERATE( range(0,10) );
        int target = GENERATE_COPY( range(0,numQubits) );
        int outcome = GENERATE( 0, 1 );
 
        qmatrix projector = getProjector(outcome);
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // that precludes outcomes=|0><0| due to zero-probability
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
 
            // compare the probabilities...
            qreal apiProb = applyForcedQubitMeasurement(qureg, target, outcome);
            qreal refProb = getReferenceProbability(ref, {target}, {outcome});
            REQUIRE_AGREE( apiProb, refProb );
 
            // and the post-projection states (caller calls subsequent REQUIRE_AGREE)
            applyReferenceOperator(ref, {target}, projector);
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);
        };
 
        CAPTURE( target, outcome );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyForcedMultiQubitMeasurement", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numTargs = GENERATE_COPY( range(1,numQubits+1) );
        auto targets = GENERATE_TARGS( numQubits, numTargs );
        auto outcomes = getRandomOutcomes(numTargs);
 
        qmatrix projector = getProjector(targets, outcomes, numQubits);
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // that precludes outcomes=|0><0| due to zero-probability
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
 
            // compare the probabilities...
            qreal apiProb = applyForcedMultiQubitMeasurement(qureg, targets.data(), outcomes.data(), numTargs);
            qreal refProb = getReferenceProbability(ref, targets, outcomes);
            REQUIRE_AGREE( apiProb, refProb );
 
            // and the post-measurement states (caller calls subsequent REQUIRE_AGREE)
            applyReferenceOperator(ref, projector);
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);
        };
 
        CAPTURE( targets, outcomes );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyMultiQubitMeasurement", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numTargs = GENERATE_COPY( range(1,numQubits+1) );
        auto targets = GENERATE_TARGS( numQubits, numTargs );
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // sampling requires the outcome probs are normalised
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
            
            // the output API state...
            qindex apiOut = applyMultiQubitMeasurement(qureg, targets.data(), numTargs);
 
            // informs the projector which determines the post-measurement reference
            auto apiOutBits = getBits(apiOut, numTargs);
            qmatrix projector = getProjector(targets, apiOutBits, numQubits);
            applyReferenceOperator(ref, projector);
            qreal refProb = getReferenceProbability(ref, targets, apiOutBits);
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);
        };
 
        CAPTURE( targets );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyMultiQubitMeasurementAndGetProb", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numTargs = GENERATE_COPY( range(1,numQubits+1) );
        auto targets = GENERATE_TARGS( numQubits, numTargs );
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // sampling requires the outcome probs are normalised
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
 
            // compare the measurement probability...
            qreal apiProb = -1;
            qindex apiOut = applyMultiQubitMeasurementAndGetProb(qureg, targets.data(), numTargs, &apiProb);
            auto apiOutBits = getBits(apiOut, numTargs);
            qreal refProb = getReferenceProbability(ref, targets, apiOutBits);
            REQUIRE_AGREE( apiProb, refProb );
            
            // and the post-measurement states (caller calls subsequent REQUIRE_AGREE)
            qmatrix projector = getProjector(targets, apiOutBits, numQubits);
            applyReferenceOperator(ref, projector);
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);
        };
 
        CAPTURE( targets );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyQubitMeasurement", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        GENERATE( range(0,10) );
        int target = GENERATE_COPY( range(0,numQubits) );
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // sampling requires the outcome probs are normalised
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
 
            // the output API state...
            int apiOut = applyQubitMeasurement(qureg, target);
 
            // informs the projector which determines the post-measurement reference
            qmatrix projector = getProjector(apiOut);
            applyReferenceOperator(ref, {target}, projector);
            qreal refProb = getReferenceProbability(ref, {target}, {apiOut});
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);
        };
 
        CAPTURE( target );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyQubitMeasurementAndGetProb", TEST_CATEGORY ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        GENERATE( range(0,10) );
        int target = GENERATE_COPY( range(0,numQubits) );
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // overwrite caller's setting of initDebugState, since
            // sampling requires the outcome probs are normalised
            setToRandomState(ref);
            setQuregToReference(qureg, ref);
 
            // compare the measurement probability...
            qreal apiProb = -1;
            int apiOut = applyQubitMeasurementAndGetProb(qureg, target, &apiProb);
            qreal refProb = getReferenceProbability(ref, {target}, {apiOut});
            REQUIRE_AGREE( apiProb, refProb );
            
            // and the post-measurement states (caller calls subsequent REQUIRE_AGREE)
            qmatrix projector = getProjector(apiOut);
            applyReferenceOperator(ref, {target}, projector);
            ref /= (qureg.isDensityMatrix)?
                refProb : std::sqrt(refProb);            
        };
 
        CAPTURE( target );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "multiplyFullStateDiagMatr", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 
    PREPARE_TEST( numQubits, cachedSV, cachedDM, refSV, refDM );
 
    auto cachedMatrs = getCachedFullStateDiagMatrs();
 
    SECTION( LABEL_CORRECTNESS ) {
 
        qmatrix refMatr = getRandomDiagonalMatrix(getPow2(numQubits));
        auto apiFunc = multiplyFullStateDiagMatr;
 
        GENERATE( range(0, TEST_NUM_MIXED_DEPLOYMENT_REPETITIONS) );
 
        SECTION( LABEL_STATEVEC ) {
 
            auto refFunc = [&] (qvector& state, qmatrix matr) { multiplyReferenceOperator(state, matr); };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedSV, cachedMatrs, apiFunc, refSV, refMatr, refFunc);
        }
 
        SECTION( LABEL_DENSMATR ) {
 
            auto refFunc = [&] (qmatrix& state, qmatrix matr) { multiplyReferenceOperator(state, matr); };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedDM, cachedMatrs, apiFunc, refDM, refMatr, refFunc);
        }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "multiplyFullStateDiagMatrPower", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 
    PREPARE_TEST( numQubits, cachedSV, cachedDM, refSV, refDM );
 
    auto cachedMatrs = getCachedFullStateDiagMatrs();
 
    SECTION( LABEL_CORRECTNESS ) {
 
        qmatrix refMatr = getRandomDiagonalMatrix(getPow2(numQubits));
        qcomp exponent = getRandomComplex();
 
        auto apiFunc = [&](Qureg qureg, FullStateDiagMatr matr) { 
            return multiplyFullStateDiagMatrPower(qureg, matr, exponent);
        };
 
        CAPTURE( exponent );
        
        GENERATE( range(0, TEST_NUM_MIXED_DEPLOYMENT_REPETITIONS) );
 
        SECTION( LABEL_STATEVEC ) {
 
            auto refFunc = [&] (qvector& state, qmatrix matr) { 
                matr = getPowerOfDiagonalMatrix(matr, exponent);
                multiplyReferenceOperator(state, matr);
            };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedSV, cachedMatrs, apiFunc, refSV, refMatr, refFunc);
        }
 
        SECTION( LABEL_DENSMATR ) {
 
            auto refFunc = [&] (qmatrix& state, qmatrix matr) { 
                matr = getPowerOfDiagonalMatrix(matr, exponent);
                multiplyReferenceOperator(state, matr);
            };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedDM, cachedMatrs, apiFunc, refDM, refMatr, refFunc);
        }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyFullStateDiagMatr", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 
    PREPARE_TEST( numQubits, cachedSV, cachedDM, refSV, refDM );
 
    auto cachedMatrs = getCachedFullStateDiagMatrs();
 
    SECTION( LABEL_CORRECTNESS ) {
 
        qmatrix refMatr = getRandomDiagonalUnitary(numQubits);
        auto apiFunc = applyFullStateDiagMatr;
 
        GENERATE( range(0, TEST_NUM_MIXED_DEPLOYMENT_REPETITIONS) );
 
        SECTION( LABEL_STATEVEC ) {
 
            auto refFunc = [&] (qvector& state, qmatrix matr) { applyReferenceOperator(state, matr); };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedSV, cachedMatrs, apiFunc, refSV, refMatr, refFunc);
        }
 
        SECTION( LABEL_DENSMATR ) {
 
            auto refFunc = [&] (qmatrix& state, qmatrix matr) { applyReferenceOperator(state, matr); };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedDM, cachedMatrs, apiFunc, refDM, refMatr, refFunc);
        }
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "applyFullStateDiagMatrPower", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 
    PREPARE_TEST( numQubits, cachedSV, cachedDM, refSV, refDM );
 
    auto cachedMatrs = getCachedFullStateDiagMatrs();
 
    SECTION( LABEL_CORRECTNESS ) {
 
        qmatrix refMatr = getRandomDiagonalUnitary(numQubits);
 
        // supplying a complex exponent requires disabling
        // numerical validation to relax unitarity
        bool testRealExp = GENERATE( true, false );
        qcomp exponent = (testRealExp)?
            qcomp(getRandomReal(-2, 2), 0):
            getRandomComplex();
 
        auto apiFunc = [&](Qureg qureg, FullStateDiagMatr matr) { 
            return applyFullStateDiagMatrPower(qureg, matr, exponent);
        };
 
        CAPTURE( exponent );
 
        GENERATE( range(0, TEST_NUM_MIXED_DEPLOYMENT_REPETITIONS) );
 
        if (!testRealExp)
            setValidationEpsilon(0);
 
        SECTION( LABEL_STATEVEC ) {
 
            auto refFunc = [&] (qvector& state, qmatrix matr) { 
                matr = getPowerOfDiagonalMatrix(matr, exponent);
                applyReferenceOperator(state, matr);
            };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedSV, cachedMatrs, apiFunc, refSV, refMatr, refFunc);
        }
 
        SECTION( LABEL_DENSMATR ) {
 
            auto refFunc = [&] (qmatrix& state, qmatrix matr) { 
                matr = getPowerOfDiagonalMatrix(matr, exponent);
                applyReferenceOperator(state, matr);
            };
 
            TEST_ON_CACHED_QUREG_AND_MATRIX( cachedDM, cachedMatrs, apiFunc, refDM, refMatr, refFunc);
        }
 
        setValidationEpsilonToDefault();
    }
 
    /// @todo input validation
}
 
 
TEST_CASE( "multiplyPauliStrSum", TEST_CATEGORY LABEL_MIXED_DEPLOY_TAG ) {
 
    PREPARE_TEST( numQubits, statevecQuregs, densmatrQuregs, statevecRef, densmatrRef );
 
    SECTION( LABEL_CORRECTNESS ) {
 
        int numQubits = getNumCachedQubits();
        int numTerms = GENERATE_COPY( 1, 2, 10 );
 
        PauliStrSum sum = createRandomPauliStrSum(numQubits, numTerms);
 
        auto testFunc = [&](Qureg qureg, auto& ref) {
 
            // must use (and ergo make) an identically-deployed workspace
            Qureg workspace = createCloneQureg(qureg);
            multiplyPauliStrSum(qureg, sum, workspace);
            destroyQureg(workspace);
 
            ref = getMatrix(sum, numQubits) * ref;
        };
 
        CAPTURE( numTerms );
        SECTION( LABEL_STATEVEC ) { TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, testFunc); }
        SECTION( LABEL_DENSMATR ) { TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, testFunc); }
    }
 
    /// @todo input validation
}
 
 
/** @} (end defgroup) */
 
 
 
/**
 * @todo
 * UNTESTED FUNCTIONS
 */
 
void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps);