qinterface.hpp

Go to the documentation of this file.

//
// (C) Daniel Strano and the Qrack contributors 2017-2023. All rights reserved.
//
// This is a multithreaded, universal quantum register simulation, allowing
// (nonphysical) register cloning and direct measurement of probability and
// phase, to leverage what advantages classical emulation of qubits can have.
//
// Licensed under the GNU Lesser General Public License V3.
// See LICENSE.md in the project root or https://www.gnu.org/licenses/lgpl-3.0.en.html
// for details.
 
#pragma once
 
#include "common/parallel_for.hpp"
#include "common/pauli.hpp"
#include "common/rdrandwrapper.hpp"
#include "hamiltonian.hpp"
 
#include <map>
#include <random>
 
#if ENABLE_UINT128
#include <ostream>
#endif
 
namespace Qrack {
 
class QInterface;
typedef std::shared_ptr<QInterface> QInterfacePtr;
 
enum QInterfaceEngine {
 
    QINTERFACE_CPU = 0,
 
    QINTERFACE_OPENCL,
 
    QINTERFACE_CUDA,
 
    QINTERFACE_HYBRID,
 
    QINTERFACE_BDT,
 
    QINTERFACE_BDT_HYBRID,
 
    QINTERFACE_STABILIZER,
 
    QINTERFACE_STABILIZER_HYBRID,
 
    QINTERFACE_QPAGER,
 
    QINTERFACE_QUNIT,
 
    QINTERFACE_QUNIT_MULTI,
 
    QINTERFACE_QUNIT_CLIFFORD,
 
    QINTERFACE_TENSOR_NETWORK,
 
    QINTERFACE_NOISY,
 
#if ENABLE_OPENCL || ENABLE_CUDA
    QINTERFACE_OPTIMAL_SCHROEDINGER = QINTERFACE_QPAGER,
#if ENABLE_OPENCL
    QINTERFACE_OPTIMAL_BASE = QINTERFACE_OPENCL,
#else
    QINTERFACE_OPTIMAL_BASE = QINTERFACE_CUDA,
#endif
#else
    QINTERFACE_OPTIMAL_SCHROEDINGER = QINTERFACE_CPU,
 
    QINTERFACE_OPTIMAL_BASE = QINTERFACE_CPU,
#endif
 
    QINTERFACE_OPTIMAL = QINTERFACE_QUNIT,
 
    QINTERFACE_OPTIMAL_MULTI = QINTERFACE_QUNIT_MULTI,
 
    QINTERFACE_MAX
};
 
class QInterface : public ParallelFor {
protected:
    bool doNormalize;
    bool randGlobalPhase;
    bool useRDRAND;
    bitLenInt qubitCount;
    uint32_t randomSeed;
    real1 amplitudeFloor;
    bitCapInt maxQPower;
    qrack_rand_gen_ptr rand_generator;
    std::uniform_real_distribution<real1_s> rand_distribution;
    std::shared_ptr<RdRandom> hardware_rand_generator;
 
    // Compilers have difficulty figuring out types and overloading if the "norm" handle is passed to std::transform. If
    // you need a safe pointer to norm(), try this:
    static inline real1_f normHelper(const complex& c) { return (real1_f)norm(c); }
 
    static inline real1_f clampProb(real1_f toClamp)
    {
        if (toClamp < ZERO_R1_F) {
            toClamp = ZERO_R1_F;
        }
        if (toClamp > ONE_R1_F) {
            toClamp = ONE_R1_F;
        }
        return toClamp;
    }
 
    complex GetNonunitaryPhase()
    {
        if (randGlobalPhase) {
            real1_f angle = Rand() * 2 * (real1_f)PI_R1;
            return complex((real1)cos(angle), (real1)sin(angle));
        } else {
            return ONE_CMPLX;
        }
    }
 
    template <typename Fn> void MACWrapper(const std::vector<bitLenInt>& controls, Fn fn)
    {
        bitCapInt xMask = ZERO_BCI;
        for (const bitLenInt& control : controls) {
            bi_or_ip(&xMask, pow2(control));
        }
 
        XMask(xMask);
        fn(controls);
        XMask(xMask);
    }
 
    virtual bitCapInt SampleClone(const std::vector<bitCapInt>& qPowers)
    {
        QInterfacePtr clone = Clone();
 
        const bitCapInt rawSample = clone->MAll();
        bitCapInt sample = ZERO_BCI;
        for (size_t i = 0U; i < qPowers.size(); ++i) {
            if (bi_compare_0(rawSample & qPowers[i]) != 0) {
                bi_or_ip(&sample, pow2(i));
            }
        }
 
        return sample;
    }
 
    virtual real1_f ExpVarUnitaryAll(bool isExp, const std::vector<bitLenInt>& bits,
        const std::vector<std::shared_ptr<complex>>& basisOps, std::vector<real1_f> eigenVals = {});
    virtual real1_f ExpVarUnitaryAll(bool isExp, const std::vector<bitLenInt>& bits,
        const std::vector<real1_f>& basisOps, std::vector<real1_f> eigenVals = {});
    virtual real1_f ExpVarBitsAll(bool isExp, const std::vector<bitLenInt>& bits, const bitCapInt& offset = ZERO_BCI)
    {
        std::vector<bitCapInt> perms;
        perms.reserve(bits.size() << 1U);
        for (size_t i = 0U; i < bits.size(); ++i) {
            perms.push_back(ZERO_BCI);
            perms.push_back(pow2(i));
        }
 
        return isExp ? ExpectationBitsFactorized(bits, perms, offset) : VarianceBitsFactorized(bits, perms, offset);
    }
 
    virtual void Copy(QInterfacePtr orig)
    {
        orig->Finish();
        doNormalize = orig->doNormalize;
        randGlobalPhase = orig->randGlobalPhase;
        useRDRAND = orig->useRDRAND;
        qubitCount = orig->qubitCount;
        randomSeed = orig->randomSeed;
        amplitudeFloor = orig->amplitudeFloor;
        maxQPower = orig->maxQPower;
        rand_generator = orig->rand_generator;
        rand_distribution = orig->rand_distribution;
        hardware_rand_generator = orig->hardware_rand_generator;
    }
 
public:
    QInterface(bitLenInt n, qrack_rand_gen_ptr rgp = nullptr, bool doNorm = false, bool useHardwareRNG = true,
        bool randomGlobalPhase = true, real1_f norm_thresh = REAL1_EPSILON);
 
    QInterface()
        : doNormalize(false)
        , randGlobalPhase(true)
        , useRDRAND(true)
        , qubitCount(0U)
        , randomSeed(0)
        , amplitudeFloor(REAL1_EPSILON)
        , maxQPower(ONE_BCI)
        , rand_distribution(0.0, 1.0)
        , hardware_rand_generator{ nullptr }
    {
        // Intentionally left blank
    }
 
    virtual ~QInterface()
    {
        // Virtual destructor for inheritance
    }
 
    void SetRandomSeed(uint32_t seed)
    {
        if (!!rand_generator) {
            rand_generator->seed(seed);
        }
    }
 
    virtual void SetQubitCount(bitLenInt qb)
    {
        qubitCount = qb;
        maxQPower = pow2(qubitCount);
    }
 
    virtual void SetConcurrency(uint32_t threadsPerEngine) { SetConcurrencyLevel(threadsPerEngine); }
 
    virtual bitLenInt GetQubitCount() { return qubitCount; }
 
    virtual bitCapInt GetMaxQPower() { return maxQPower; }
 
    virtual bitCapInt GetAmplitudeCount() { return maxQPower; }
 
    virtual bool GetIsArbitraryGlobalPhase() { return randGlobalPhase; }
 
    real1_f Rand()
    {
        if (!!hardware_rand_generator) {
            return hardware_rand_generator->Next();
        } else {
            return rand_distribution(*rand_generator);
        }
    }
 
    virtual void LossySaveStateVector(std::string f, int p = 6, int b = 4);
    virtual void LossyLoadStateVector(std::string f);
 
    virtual void SetQuantumState(const complex* inputState) = 0;
 
    virtual void GetQuantumState(complex* outputState) = 0;
 
    virtual void GetReducedDensityMatrix(const std::vector<bitLenInt>& qubits, complex* outputState);
 
    virtual void GetProbs(real1* outputProbs) = 0;
 
    virtual complex GetAmplitude(const bitCapInt& perm) = 0;
 
    virtual void SetAmplitude(const bitCapInt& perm, const complex& amp) = 0;
 
    virtual void SetPermutation(const bitCapInt& perm, const complex& phaseFac = CMPLX_DEFAULT_ARG);
 
    virtual bitLenInt Compose(QInterfacePtr toCopy) { return Compose(toCopy, qubitCount); }
    virtual bitLenInt ComposeNoClone(QInterfacePtr toCopy) { return Compose(toCopy); }
    virtual std::map<QInterfacePtr, bitLenInt> Compose(std::vector<QInterfacePtr> toCopy);
    virtual bitLenInt Compose(QInterfacePtr toCopy, bitLenInt start);
 
    virtual void Decompose(bitLenInt start, QInterfacePtr dest) = 0;
 
    virtual QInterfacePtr Decompose(bitLenInt start, bitLenInt length) = 0;
 
    virtual void Dispose(bitLenInt start, bitLenInt length) = 0;
 
    virtual void Dispose(bitLenInt start, bitLenInt length, const bitCapInt& disposedPerm) = 0;
 
    virtual bitLenInt Allocate(bitLenInt length) { return Allocate(qubitCount, length); }
 
    virtual bitLenInt Allocate(bitLenInt start, bitLenInt length) = 0;
 
    virtual void Mtrx(const complex mtrx[4U], bitLenInt qubit) = 0;
 
    virtual void MCMtrx(const std::vector<bitLenInt>& controls, const complex mtrx[4U], bitLenInt target) = 0;
 
    virtual void MACMtrx(const std::vector<bitLenInt>& controls, const complex mtrx[4U], bitLenInt target)
    {
        if (IS_NORM_0(mtrx[1U]) && IS_NORM_0(mtrx[2U])) {
            MACPhase(controls, mtrx[0U], mtrx[3U], target);
        } else if (IS_NORM_0(mtrx[0U]) && IS_NORM_0(mtrx[3U])) {
            MACInvert(controls, mtrx[1U], mtrx[2U], target);
        } else {
            MACWrapper(controls, [this, mtrx, target](const std::vector<bitLenInt>& lc) { MCMtrx(lc, mtrx, target); });
        }
    }
 
    virtual void UCMtrx(
        const std::vector<bitLenInt>& controls, const complex* mtrx, bitLenInt target, const bitCapInt& controlPerm);
 
    virtual void Phase(const complex& topLeft, const complex& bottomRight, bitLenInt qubit)
    {
        if ((randGlobalPhase || IS_NORM_0(ONE_CMPLX - topLeft)) && IS_NORM_0(topLeft - bottomRight)) {
            return;
        }
 
        const complex mtrx[4U]{ topLeft, ZERO_CMPLX, ZERO_CMPLX, bottomRight };
        Mtrx(mtrx, qubit);
    }
 
    virtual void Invert(const complex& topRight, const complex& bottomLeft, bitLenInt qubit)
    {
        const complex mtrx[4U]{ ZERO_CMPLX, topRight, bottomLeft, ZERO_CMPLX };
        Mtrx(mtrx, qubit);
    }
 
    virtual void MCPhase(
        const std::vector<bitLenInt>& controls, const complex& topLeft, const complex& bottomRight, bitLenInt target)
    {
        if (IS_NORM_0(ONE_CMPLX - topLeft) && IS_NORM_0(ONE_CMPLX - bottomRight)) {
            return;
        }
 
        const complex mtrx[4U]{ topLeft, ZERO_CMPLX, ZERO_CMPLX, bottomRight };
        MCMtrx(controls, mtrx, target);
    }
 
    virtual void MCInvert(
        const std::vector<bitLenInt>& controls, const complex& topRight, const complex& bottomLeft, bitLenInt target)
    {
        const complex mtrx[4U]{ ZERO_CMPLX, topRight, bottomLeft, ZERO_CMPLX };
        MCMtrx(controls, mtrx, target);
    }
 
    virtual void MACPhase(
        const std::vector<bitLenInt>& controls, const complex& topLeft, const complex& bottomRight, bitLenInt target)
    {
        if (IS_NORM_0(ONE_CMPLX - topLeft) && IS_NORM_0(ONE_CMPLX - bottomRight)) {
            return;
        }
 
        MACWrapper(controls, [this, topLeft, bottomRight, target](const std::vector<bitLenInt>& lc) {
            MCPhase(lc, topLeft, bottomRight, target);
        });
    }
 
    virtual void MACInvert(
        const std::vector<bitLenInt>& controls, const complex& topRight, const complex& bottomLeft, bitLenInt target)
    {
        MACWrapper(controls, [this, topRight, bottomLeft, target](const std::vector<bitLenInt>& lc) {
            MCInvert(lc, topRight, bottomLeft, target);
        });
    }
 
    virtual void UCPhase(const std::vector<bitLenInt>& controls, const complex& topLeft, const complex& bottomRight,
        bitLenInt target, const bitCapInt& perm)
    {
        if (IS_NORM_0(ONE_CMPLX - topLeft) && IS_NORM_0(ONE_CMPLX - bottomRight)) {
            return;
        }
 
        const complex mtrx[4U]{ topLeft, ZERO_CMPLX, ZERO_CMPLX, bottomRight };
        UCMtrx(controls, mtrx, target, perm);
    }
 
    virtual void UCInvert(const std::vector<bitLenInt>& controls, const complex& topRight, const complex& bottomLeft,
        bitLenInt target, const bitCapInt& perm)
    {
        const complex mtrx[4U]{ ZERO_CMPLX, topRight, bottomLeft, ZERO_CMPLX };
        UCMtrx(controls, mtrx, target, perm);
    }
 
    virtual void UniformlyControlledSingleBit(
        const std::vector<bitLenInt>& controls, bitLenInt qubit, const complex* mtrxs)
    {
        UniformlyControlledSingleBit(controls, qubit, mtrxs, std::vector<bitCapInt>(), ZERO_BCI);
    }
    virtual void UniformlyControlledSingleBit(const std::vector<bitLenInt>& controls, bitLenInt qubit,
        const complex* mtrxs, const std::vector<bitCapInt>& mtrxSkipPowers, const bitCapInt& mtrxSkipValueMask);
 
    virtual void TimeEvolve(Hamiltonian h, real1_f timeDiff);
 
    virtual void CSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void AntiCSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void CSqrtSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void AntiCSqrtSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void CISqrtSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void AntiCISqrtSwap(const std::vector<bitLenInt>& controls, bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void CCNOT(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MCInvert(controls, ONE_CMPLX, ONE_CMPLX, target);
    }
 
    virtual void AntiCCNOT(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MACInvert(controls, ONE_CMPLX, ONE_CMPLX, target);
    }
 
    virtual void CNOT(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCInvert(controls, ONE_CMPLX, ONE_CMPLX, target);
    }
 
    virtual void AntiCNOT(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MACInvert(controls, ONE_CMPLX, ONE_CMPLX, target);
    }
 
    virtual void CY(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCInvert(controls, -I_CMPLX, I_CMPLX, target);
    }
 
    virtual void AntiCY(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MACInvert(controls, -I_CMPLX, I_CMPLX, target);
    }
 
    virtual void CCY(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MCInvert(controls, -I_CMPLX, I_CMPLX, target);
    }
 
    virtual void AntiCCY(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MACInvert(controls, -I_CMPLX, I_CMPLX, target);
    }
 
    virtual void CZ(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, -ONE_CMPLX, target);
    }
 
    virtual void AntiCZ(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MACPhase(controls, ONE_CMPLX, -ONE_CMPLX, target);
    }
 
    virtual void CCZ(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MCPhase(controls, ONE_CMPLX, -ONE_CMPLX, target);
    }
 
    virtual void AntiCCZ(bitLenInt control1, bitLenInt control2, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control1, control2 };
        MACPhase(controls, ONE_CMPLX, -ONE_CMPLX, target);
    }
 
    virtual void U(bitLenInt target, real1_f theta, real1_f phi, real1_f lambda);
 
    virtual void U2(bitLenInt target, real1_f phi, real1_f lambda) { U(target, (real1_f)(M_PI / 2), phi, lambda); }
 
    virtual void IU2(bitLenInt target, real1_f phi, real1_f lambda)
    {
        U(target, (real1_f)(M_PI / 2), (real1_f)(-lambda - PI_R1), (real1_f)(-phi + PI_R1));
    }
 
    virtual void AI(bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void IAI(bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void CAI(bitLenInt control, bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void AntiCAI(bitLenInt control, bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void CIAI(bitLenInt control, bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void AntiCIAI(bitLenInt control, bitLenInt target, real1_f azimuth, real1_f inclination);
 
    virtual void CU(const std::vector<bitLenInt>& controls, bitLenInt target, real1_f theta, real1_f phi,
        real1_f lambda, real1_f gamma = ZERO_R1_F);
 
    virtual void AntiCU(const std::vector<bitLenInt>& controls, bitLenInt target, real1_f theta, real1_f phi,
        real1_f lambda, real1_f gamma = ZERO_R1_F);
 
    virtual void H(bitLenInt qubit)
    {
        QRACK_CONST complex C_SQRT1_2 = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex C_SQRT1_2_NEG = complex(-SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex mtrx[4]{ C_SQRT1_2, C_SQRT1_2, C_SQRT1_2, C_SQRT1_2_NEG };
        Mtrx(mtrx, qubit);
    }
 
    virtual void SqrtH(bitLenInt qubit)
    {
        QRACK_CONST complex m00 =
            complex((real1)((ONE_R1 + SQRT2_R1) / (2 * SQRT2_R1)), (real1)((-ONE_R1 + SQRT2_R1) / (2 * SQRT2_R1)));
        QRACK_CONST complex m01 = complex((real1)(SQRT1_2_R1 / 2), (real1)(-SQRT1_2_R1 / 2));
        QRACK_CONST complex m10 = m01;
        QRACK_CONST complex m11 =
            complex((real1)((-ONE_R1 + SQRT2_R1) / (2 * SQRT2_R1)), (real1)((ONE_R1 + SQRT2_R1) / (2 * SQRT2_R1)));
        QRACK_CONST complex mtrx[4]{ m00, m01, m10, m11 };
        Mtrx(mtrx, qubit);
    }
 
    virtual void SH(bitLenInt qubit)
    {
        QRACK_CONST complex C_SQRT1_2 = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex C_I_SQRT1_2 = complex(ZERO_R1, SQRT1_2_R1);
        QRACK_CONST complex C_I_SQRT1_2_NEG = complex(ZERO_R1, -SQRT1_2_R1);
        QRACK_CONST complex mtrx[4]{ C_SQRT1_2, C_SQRT1_2, C_I_SQRT1_2, C_I_SQRT1_2_NEG };
        Mtrx(mtrx, qubit);
    }
 
    virtual void HIS(bitLenInt qubit)
    {
        QRACK_CONST complex C_SQRT1_2 = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex C_I_SQRT1_2 = complex(ZERO_R1, SQRT1_2_R1);
        QRACK_CONST complex C_I_SQRT1_2_NEG = complex(ZERO_R1, -SQRT1_2_R1);
        QRACK_CONST complex mtrx[4]{ C_SQRT1_2, C_I_SQRT1_2_NEG, C_SQRT1_2, C_I_SQRT1_2 };
        Mtrx(mtrx, qubit);
    }
 
    virtual bool M(bitLenInt qubit) { return ForceM(qubit, false, false); };
 
    virtual bool ForceM(bitLenInt qubit, bool result, bool doForce = true, bool doApply = true) = 0;
 
    virtual void S(bitLenInt qubit) { Phase(ONE_CMPLX, I_CMPLX, qubit); }
 
    virtual void IS(bitLenInt qubit) { Phase(ONE_CMPLX, -I_CMPLX, qubit); }
 
    virtual void T(bitLenInt qubit) { Phase(ONE_CMPLX, complex(SQRT1_2_R1, SQRT1_2_R1), qubit); }
 
    virtual void IT(bitLenInt qubit) { Phase(ONE_CMPLX, complex(SQRT1_2_R1, -SQRT1_2_R1), qubit); }
 
    virtual void PhaseRootN(bitLenInt n, bitLenInt qubit)
    {
        if (n == 0) {
            return;
        }
 
        Phase(ONE_CMPLX, pow(-ONE_CMPLX, (real1)(ONE_R1 / pow2Ocl(n - 1U))), qubit);
    }
 
    virtual void IPhaseRootN(bitLenInt n, bitLenInt qubit)
    {
        if (n == 0) {
            return;
        }
 
        Phase(ONE_CMPLX, pow(-ONE_CMPLX, (real1)(-ONE_R1 / pow2Ocl(n - 1U))), qubit);
    }
 
    virtual void PhaseRootNMask(bitLenInt n, const bitCapInt& mask);
 
    virtual void PhaseParity(real1_f radians, const bitCapInt& mask);
 
    virtual void X(bitLenInt qubit) { Invert(ONE_CMPLX, ONE_CMPLX, qubit); }
 
    virtual void XMask(const bitCapInt& mask);
 
    virtual void Y(bitLenInt qubit) { Invert(-I_CMPLX, I_CMPLX, qubit); }
 
    virtual void YMask(const bitCapInt& mask);
 
    virtual void Z(bitLenInt qubit) { Phase(ONE_CMPLX, -ONE_CMPLX, qubit); }
 
    virtual void ZMask(const bitCapInt& mask);
 
    virtual void SqrtX(bitLenInt qubit)
    {
        QRACK_CONST complex ONE_PLUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex ONE_MINUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ ONE_PLUS_I_DIV_2, ONE_MINUS_I_DIV_2, ONE_MINUS_I_DIV_2, ONE_PLUS_I_DIV_2 };
        Mtrx(mtrx, qubit);
    }
 
    virtual void ISqrtX(bitLenInt qubit)
    {
        QRACK_CONST complex ONE_PLUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex ONE_MINUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ ONE_MINUS_I_DIV_2, ONE_PLUS_I_DIV_2, ONE_PLUS_I_DIV_2, ONE_MINUS_I_DIV_2 };
        Mtrx(mtrx, qubit);
    }
 
    virtual void SqrtY(bitLenInt qubit)
    {
        QRACK_CONST complex ONE_PLUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex ONE_PLUS_I_DIV_2_NEG = complex((real1)(-ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ ONE_PLUS_I_DIV_2, ONE_PLUS_I_DIV_2_NEG, ONE_PLUS_I_DIV_2, ONE_PLUS_I_DIV_2 };
        Mtrx(mtrx, qubit);
    }
 
    virtual void ISqrtY(bitLenInt qubit)
    {
        QRACK_CONST complex ONE_MINUS_I_DIV_2 = complex((real1)(ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex ONE_MINUS_I_DIV_2_NEG = complex((real1)(-ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ ONE_MINUS_I_DIV_2, ONE_MINUS_I_DIV_2, ONE_MINUS_I_DIV_2_NEG, ONE_MINUS_I_DIV_2 };
        Mtrx(mtrx, qubit);
    }
 
    virtual void SqrtW(bitLenInt qubit)
    {
        QRACK_CONST complex diag = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex m01 = complex((real1)(-ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex m10 = complex((real1)(ONE_R1 / 2), (real1)(-ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ diag, m01, m10, diag };
        Mtrx(mtrx, qubit);
    }
 
    virtual void ISqrtW(bitLenInt qubit)
    {
        QRACK_CONST complex diag = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex m01 = complex((real1)(ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex m10 = complex((real1)(-ONE_R1 / 2), (real1)(ONE_R1 / 2));
        QRACK_CONST complex mtrx[4]{ diag, m01, m10, diag };
        Mtrx(mtrx, qubit);
    }
 
    virtual void CH(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        QRACK_CONST complex C_SQRT1_2 = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex C_SQRT1_2_NEG = complex(-SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex mtrx[4]{ C_SQRT1_2, C_SQRT1_2, C_SQRT1_2, C_SQRT1_2_NEG };
        MCMtrx(controls, mtrx, target);
    }
 
    virtual void AntiCH(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        QRACK_CONST complex C_SQRT1_2 = complex(SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex C_SQRT1_2_NEG = complex(-SQRT1_2_R1, ZERO_R1);
        QRACK_CONST complex mtrx[4]{ C_SQRT1_2, C_SQRT1_2, C_SQRT1_2, C_SQRT1_2_NEG };
        MACMtrx(controls, mtrx, target);
    }
 
    virtual void CS(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, I_CMPLX, target);
    }
 
    virtual void AntiCS(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MACPhase(controls, ONE_CMPLX, I_CMPLX, target);
    }
 
    virtual void CIS(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, -I_CMPLX, target);
    }
 
    virtual void AntiCIS(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MACPhase(controls, ONE_CMPLX, -I_CMPLX, target);
    }
 
    virtual void CT(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, complex(SQRT1_2_R1, SQRT1_2_R1), target);
    }
 
    virtual void CIT(bitLenInt control, bitLenInt target)
    {
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, complex(SQRT1_2_R1, -SQRT1_2_R1), target);
    }
 
    virtual void CPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target)
    {
        if (n == 0) {
            return;
        }
 
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, pow(-ONE_CMPLX, (real1)(ONE_R1 / pow2Ocl(n - 1U))), target);
    }
 
    virtual void AntiCPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target)
    {
        if (n == 0) {
            return;
        }
 
        const std::vector<bitLenInt> controls{ control };
        MACPhase(controls, ONE_CMPLX, pow(-ONE_CMPLX, (real1)(ONE_R1 / pow2Ocl(n - 1U))), target);
    }
 
    virtual void CIPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target)
    {
        if (n == 0) {
            return;
        }
 
        const std::vector<bitLenInt> controls{ control };
        MCPhase(controls, ONE_CMPLX, pow(-ONE_CMPLX, (real1)(-ONE_R1 / pow2Ocl(n - 1U))), target);
    }
 
    virtual void AntiCIPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target)
    {
        if (n == 0) {
            return;
        }
 
        const std::vector<bitLenInt> controls{ control };
        MACPhase(controls, ONE_CMPLX, pow(-ONE_CMPLX, (real1)(-ONE_R1 / pow2Ocl(n - 1U))), target);
    }
 
    virtual void AND(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void OR(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void XOR(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void CLAND(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void CLOR(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void CLXOR(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void NAND(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void NOR(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void XNOR(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt outputBit);
 
    virtual void CLNAND(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void CLNOR(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void CLXNOR(bitLenInt inputQBit, bool inputClassicalBit, bitLenInt outputBit);
 
    virtual void UniformlyControlledRY(const std::vector<bitLenInt>& controls, bitLenInt qubit, const real1* angles);
 
    virtual void UniformlyControlledRZ(const std::vector<bitLenInt>& controls, bitLenInt qubit, const real1* angles);
 
    virtual void RT(real1_f radians, bitLenInt qubit);
 
    virtual void RX(real1_f radians, bitLenInt qubit);
 
    virtual void RY(real1_f radians, bitLenInt qubit);
 
    virtual void RZ(real1_f radians, bitLenInt qubit);
 
    virtual void CRZ(real1_f radians, bitLenInt control, bitLenInt target);
 
    virtual void CRY(real1_f radians, bitLenInt control, bitLenInt target);
 
#if ENABLE_ROT_API
    virtual void RTDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void RXDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void Exp(real1_f radians, bitLenInt qubit);
 
    virtual void Exp(
        const std::vector<bitLenInt>& controls, bitLenInt qubit, const complex mtrx[4U], bool antiCtrled = false);
 
    virtual void ExpDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void ExpX(real1_f radians, bitLenInt qubit);
 
    virtual void ExpXDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void ExpY(real1_f radians, bitLenInt qubit);
 
    virtual void ExpYDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void ExpZ(real1_f radians, bitLenInt qubit);
 
    virtual void ExpZDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void CRX(real1_f radians, bitLenInt control, bitLenInt target);
 
    virtual void CRXDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target);
 
    virtual void RYDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void CRYDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target);
 
    virtual void RZDyad(int numerator, int denomPower, bitLenInt qubit);
 
    virtual void CRZDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target);
 
    virtual void CRT(real1_f radians, bitLenInt control, bitLenInt target);
 
    virtual void CRTDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target);
#endif
 
    virtual void H(bitLenInt start, bitLenInt length);
 
    virtual void X(bitLenInt start, bitLenInt length) { XMask(bitRegMask(start, length)); }
 
#if ENABLE_REG_GATES
 
    virtual void U(bitLenInt start, bitLenInt length, real1_f theta, real1_f phi, real1_f lambda);
 
    virtual void U2(bitLenInt start, bitLenInt length, real1_f phi, real1_f lambda);
 
    virtual void SH(bitLenInt start, bitLenInt length);
 
    virtual void HIS(bitLenInt start, bitLenInt length);
 
    virtual void SqrtH(bitLenInt start, bitLenInt length);
 
    virtual void S(bitLenInt start, bitLenInt length);
 
    virtual void IS(bitLenInt start, bitLenInt length);
 
    virtual void T(bitLenInt start, bitLenInt length);
 
    virtual void IT(bitLenInt start, bitLenInt length);
 
    virtual void PhaseRootN(bitLenInt n, bitLenInt start, bitLenInt length);
 
    virtual void IPhaseRootN(bitLenInt n, bitLenInt start, bitLenInt length);
 
    virtual void Y(bitLenInt start, bitLenInt length);
 
    virtual void SqrtX(bitLenInt start, bitLenInt length);
 
    virtual void ISqrtX(bitLenInt start, bitLenInt length);
 
    virtual void SqrtY(bitLenInt start, bitLenInt length);
 
    virtual void ISqrtY(bitLenInt start, bitLenInt length);
 
    virtual void Z(bitLenInt start, bitLenInt length);
 
    virtual void CNOT(bitLenInt inputBits, bitLenInt targetBits, bitLenInt length);
 
    virtual void AntiCNOT(bitLenInt inputBits, bitLenInt targetBits, bitLenInt length);
 
    virtual void CCNOT(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void AntiCCNOT(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void CY(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void AntiCY(bitLenInt inputBits, bitLenInt targetBits, bitLenInt length);
 
    virtual void CCY(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void AntiCCY(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void CZ(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void AntiCZ(bitLenInt inputBits, bitLenInt targetBits, bitLenInt length);
 
    virtual void CCZ(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void AntiCCZ(bitLenInt control1, bitLenInt control2, bitLenInt target, bitLenInt length);
 
    virtual void Swap(bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void ISwap(bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void IISwap(bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void SqrtSwap(bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void ISqrtSwap(bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void FSim(real1_f theta, real1_f phi, bitLenInt start1, bitLenInt start2, bitLenInt length);
 
    virtual void AND(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLAND(bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
    virtual void OR(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLOR(bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
    virtual void XOR(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLXOR(bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
    virtual void NAND(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLNAND(
        bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
    virtual void NOR(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLNOR(bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
    virtual void XNOR(bitLenInt inputStart1, bitLenInt inputStart2, bitLenInt outputStart, bitLenInt length);
 
    virtual void CLXNOR(
        bitLenInt qInputStart, const bitCapInt& classicalInput, bitLenInt outputStart, bitLenInt length);
 
#if ENABLE_ROT_API
    virtual void RT(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void RTDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void RX(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void RXDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void CRX(real1_f radians, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CRXDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void RY(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void RYDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void CRY(real1_f radians, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CRYDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void RZ(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void RZDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void CRZ(real1_f radians, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CRZDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CRT(real1_f radians, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CRTDyad(int numerator, int denomPower, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void Exp(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void ExpDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void ExpX(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void ExpXDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void ExpY(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void ExpYDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
 
    virtual void ExpZ(real1_f radians, bitLenInt start, bitLenInt length);
 
    virtual void ExpZDyad(int numerator, int denomPower, bitLenInt start, bitLenInt length);
#endif
 
    virtual void CH(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CS(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CIS(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CT(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CIT(bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target, bitLenInt length);
 
    virtual void CIPhaseRootN(bitLenInt n, bitLenInt control, bitLenInt target, bitLenInt length);
 
#endif
 
    virtual void ROL(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void ROR(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void ASL(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void ASR(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void LSL(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void LSR(bitLenInt shift, bitLenInt start, bitLenInt length);
 
    virtual void INC(const bitCapInt& toAdd, bitLenInt start, bitLenInt length);
 
    virtual void DEC(const bitCapInt& toSub, bitLenInt start, bitLenInt length)
    {
        const bitCapInt invToSub = pow2(length) - toSub;
        INC(invToSub, start, length);
    }
 
    virtual void INCDECC(const bitCapInt& toAdd, bitLenInt start, bitLenInt length, bitLenInt carryIndex);
 
    virtual void INCC(const bitCapInt& toAdd, bitLenInt start, bitLenInt length, bitLenInt carryIndex)
    {
        const bool hasCarry = M(carryIndex);
        if (hasCarry) {
            X(carryIndex);
            INCDECC(toAdd + 1U, start, length, carryIndex);
        } else {
            INCDECC(toAdd, start, length, carryIndex);
        }
    }
 
    virtual void DECC(const bitCapInt& toSub, bitLenInt start, bitLenInt length, bitLenInt carryIndex)
    {
        const bool hasCarry = M(carryIndex);
        bitCapInt invToSub = pow2(length) - toSub;
        if (hasCarry) {
            X(carryIndex);
        } else {
            --invToSub;
        }
 
        INCDECC(invToSub, start, length, carryIndex);
    }
 
    virtual void CINC(
        const bitCapInt& toAdd, bitLenInt inOutStart, bitLenInt length, const std::vector<bitLenInt>& controls);
 
    virtual void CDEC(
        const bitCapInt& toSub, bitLenInt inOutStart, bitLenInt length, const std::vector<bitLenInt>& controls)
    {
        const bitCapInt invToSub = pow2(length) - toSub;
        CINC(invToSub, inOutStart, length, controls);
    }
 
    virtual void INCS(const bitCapInt& toAdd, bitLenInt start, bitLenInt length, bitLenInt overflowIndex)
    {
        const bitCapInt signMask = pow2(length - 1U);
        INC(signMask, start, length);
        INCDECC(toAdd & ~signMask, start, length, overflowIndex);
        if (bi_compare_0(toAdd & signMask) == 0) {
            DEC(signMask, start, length);
        }
    }
 
    virtual void DECS(const bitCapInt& toSub, bitLenInt start, bitLenInt length, bitLenInt overflowIndex)
    {
        const bitCapInt invToSub = pow2(length) - toSub;
        INCS(invToSub, start, length, overflowIndex);
    }
 
    virtual void MULModNOut(
        const bitCapInt& toMul, const bitCapInt& modN, bitLenInt inStart, bitLenInt outStart, bitLenInt length);
 
    virtual void IMULModNOut(
        const bitCapInt& toMul, const bitCapInt& modN, bitLenInt inStart, bitLenInt outStart, bitLenInt length);
 
    virtual void CMULModNOut(const bitCapInt& toMul, const bitCapInt& modN, bitLenInt inStart, bitLenInt outStart,
        bitLenInt length, const std::vector<bitLenInt>& controls);
 
    virtual void CIMULModNOut(const bitCapInt& toMul, const bitCapInt& modN, bitLenInt inStart, bitLenInt outStart,
        bitLenInt length, const std::vector<bitLenInt>& controls);
 
    virtual void FullAdd(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt carryInSumOut, bitLenInt carryOut)
    {
        // See https://quantumcomputing.stackexchange.com/questions/1654/how-do-i-add-11-using-a-quantum-computer
 
        // Assume outputBit is in 0 state.
        CCNOT(inputBit1, inputBit2, carryOut);
        CNOT(inputBit1, inputBit2);
        CCNOT(inputBit2, carryInSumOut, carryOut);
        CNOT(inputBit2, carryInSumOut);
        CNOT(inputBit1, inputBit2);
    }
 
    virtual void IFullAdd(bitLenInt inputBit1, bitLenInt inputBit2, bitLenInt carryInSumOut, bitLenInt carryOut)
    {
        // See https://quantumcomputing.stackexchange.com/questions/1654/how-do-i-add-11-using-a-quantum-computer
        // Quantum computing is reversible! Simply perform the inverse operations in reverse order!
        // (CNOT and CCNOT are self-inverse.)
 
        // Assume outputBit is in 0 state.
        CNOT(inputBit1, inputBit2);
        CNOT(inputBit2, carryInSumOut);
        CCNOT(inputBit2, carryInSumOut, carryOut);
        CNOT(inputBit1, inputBit2);
        CCNOT(inputBit1, inputBit2, carryOut);
    }
 
    virtual void CFullAdd(const std::vector<bitLenInt>& controls, bitLenInt inputBit1, bitLenInt inputBit2,
        bitLenInt carryInSumOut, bitLenInt carryOut);
 
    virtual void CIFullAdd(const std::vector<bitLenInt>& controls, bitLenInt inputBit1, bitLenInt inputBit2,
        bitLenInt carryInSumOut, bitLenInt carryOut);
 
    virtual void ADC(bitLenInt input1, bitLenInt input2, bitLenInt output, bitLenInt length, bitLenInt carry);
 
    virtual void IADC(bitLenInt input1, bitLenInt input2, bitLenInt output, bitLenInt length, bitLenInt carry);
 
    virtual void CADC(const std::vector<bitLenInt>& controls, bitLenInt input1, bitLenInt input2, bitLenInt output,
        bitLenInt length, bitLenInt carry);
 
    virtual void CIADC(const std::vector<bitLenInt>& controls, bitLenInt input1, bitLenInt input2, bitLenInt output,
        bitLenInt length, bitLenInt carry);
 
    virtual void QFT(bitLenInt start, bitLenInt length, bool trySeparate = false);
 
    virtual void QFTR(const std::vector<bitLenInt>& qubits, bool trySeparate = false);
 
    virtual void IQFT(bitLenInt start, bitLenInt length, bool trySeparate = false);
 
    virtual void IQFTR(const std::vector<bitLenInt>& qubits, bool trySeparate = false);
 
    virtual void ZeroPhaseFlip(bitLenInt start, bitLenInt length);
 
    virtual void PhaseFlip() { Phase(-ONE_CMPLX, -ONE_CMPLX, 0); }
 
    virtual void SetReg(bitLenInt start, bitLenInt length, const bitCapInt& value);
 
    virtual bitCapInt MReg(bitLenInt start, bitLenInt length) { return ForceMReg(start, length, ZERO_BCI, false); }
 
    virtual bitCapInt MAll() { return MReg(0, qubitCount); }
 
    virtual bitCapInt ForceMReg(
        bitLenInt start, bitLenInt length, const bitCapInt& result, bool doForce = true, bool doApply = true);
 
    virtual bitCapInt M(const std::vector<bitLenInt>& bits) { return ForceM(bits, std::vector<bool>()); }
 
    virtual bitCapInt ForceM(const std::vector<bitLenInt>& bits, const std::vector<bool>& values, bool doApply = true);
 
    virtual void Swap(bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void ISwap(bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void IISwap(bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void SqrtSwap(bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void ISqrtSwap(bitLenInt qubit1, bitLenInt qubit2);
 
    virtual void FSim(real1_f theta, real1_f phi, bitLenInt qubit1, bitLenInt qubit2)
    {
        const std::vector<bitLenInt> controls{ qubit1 };
        const real1 sinTheta = (real1)sin(theta);
        const complex expIPhi = exp(complex(ZERO_R1, -(real1)phi));
 
        MCPhase(controls, ONE_CMPLX, expIPhi, qubit2);
 
        if ((2 * abs(sinTheta)) <= (FP_NORM_EPSILON * PI_R1)) {
            // ISWAP power is 0.
            return;
        }
 
        const real1 sinThetaDiffNeg = ONE_R1 + sinTheta;
        if ((2 * abs(sinThetaDiffNeg)) <= (FP_NORM_EPSILON * PI_R1)) {
            return ISwap(qubit1, qubit2);
        }
 
        const real1 sinThetaDiffPos = ONE_R1 - sinTheta;
        if ((2 * abs(sinThetaDiffPos)) <= (FP_NORM_EPSILON * PI_R1)) {
            return IISwap(qubit1, qubit2);
        }
 
        // Decompose root of ISWAP gate:
 
        // RXX
        H(qubit1);
        H(qubit2);
        CNOT(qubit1, qubit2);
        RZ(2 * theta, qubit2);
        CNOT(qubit1, qubit2);
        H(qubit2);
        H(qubit1);
 
        // RYY
        S(qubit1);
        S(qubit2);
        H(qubit1);
        H(qubit2);
        CNOT(qubit1, qubit2);
        RZ(2 * theta, qubit2);
        CNOT(qubit1, qubit2);
        H(qubit2);
        H(qubit1);
        IS(qubit2);
        IS(qubit1);
    }
 
    virtual void Reverse(bitLenInt first, bitLenInt last)
    {
        while ((last > 0U) && (first < (last - 1U))) {
            --last;
            Swap(first, last);
            ++first;
        }
    }
 
    virtual real1_f Prob(bitLenInt qubit) = 0;
    virtual real1_f CProb(bitLenInt control, bitLenInt target)
    {
        AntiCNOT(control, target);
        const real1_f prob = Prob(target);
        AntiCNOT(control, target);
 
        return prob;
    }
    virtual real1_f ACProb(bitLenInt control, bitLenInt target)
    {
        CNOT(control, target);
        const real1_f prob = Prob(target);
        CNOT(control, target);
 
        return prob;
    }
 
    virtual bitCapInt HighestProbAll()
    {
        real1_f totProb = ZERO_R1_F;
        real1_f highestProb = ZERO_R1_F;
        bitCapInt bestPerm = ZERO_BCI;
        for (bitCapInt p = ZERO_BCI; p < maxQPower; ++p) {
            real1_f prob = ProbAll(p);
            if (prob > highestProb) {
                highestProb = prob;
                bestPerm = p;
            }
            totProb += prob;
            if (highestProb > (ONE_R1_F - totProb)) {
                break;
            }
        }
 
        return bestPerm;
    }
 
    virtual std::vector<bitCapInt> HighestProbAll(size_t n);
 
    virtual real1_f ProbAll(const bitCapInt& fullRegister)
    {
        return clampProb((real1_f)norm(GetAmplitude(fullRegister)));
    }
 
    virtual real1_f ProbReg(bitLenInt start, bitLenInt length, const bitCapInt& permutation);
 
    virtual real1_f ProbMask(const bitCapInt& mask, const bitCapInt& permutation);
 
    virtual void ProbMaskAll(const bitCapInt& mask, real1* probsArray);
 
    virtual void ProbBitsAll(const std::vector<bitLenInt>& bits, real1* probsArray);
 
    virtual real1_f VarianceBitsAll(const std::vector<bitLenInt>& bits, const bitCapInt& offset = ZERO_BCI)
    {
        return ExpVarBitsAll(false, bits, offset);
    }
 
    virtual real1_f VarianceBitsAllRdm(
        bool roundRz, const std::vector<bitLenInt>& bits, const bitCapInt& offset = ZERO_BCI)
    {
        return VarianceBitsAll(bits, offset);
    }
 
    virtual real1_f VariancePauliAll(std::vector<bitLenInt> bits, std::vector<Pauli> paulis);
 
    virtual real1_f VarianceUnitaryAll(
        const std::vector<bitLenInt>& bits, const std::vector<real1_f>& basisOps, std::vector<real1_f> eigenVals = {})
    {
        return ExpVarUnitaryAll(false, bits, basisOps, eigenVals);
    }
 
    virtual real1_f VarianceUnitaryAll(const std::vector<bitLenInt>& bits,
        const std::vector<std::shared_ptr<complex>>& basisOps, std::vector<real1_f> eigenVals = {})
    {
        return ExpVarUnitaryAll(false, bits, basisOps, eigenVals);
    }
 
    virtual real1_f VarianceFloatsFactorized(const std::vector<bitLenInt>& bits, const std::vector<real1_f>& weights);
 
    virtual real1_f VarianceFloatsFactorizedRdm(
        bool roundRz, const std::vector<bitLenInt>& bits, const std::vector<real1_f>& weights)
    {
        return VarianceFloatsFactorized(bits, weights);
    }
 
    virtual real1_f VarianceBitsFactorized(
        const std::vector<bitLenInt>& bits, const std::vector<bitCapInt>& perms, const bitCapInt& offset = ZERO_BCI);
 
    virtual real1_f VarianceBitsFactorizedRdm(bool roundRz, const std::vector<bitLenInt>& bits,
        const std::vector<bitCapInt>& perms, const bitCapInt& offset = ZERO_BCI)
    {
        return VarianceBitsFactorized(bits, perms, offset);
    }
 
    virtual real1_f ExpectationBitsAll(const std::vector<bitLenInt>& bits, const bitCapInt& offset = ZERO_BCI)
    {
        return ExpVarBitsAll(true, bits, offset);
    }
 
    virtual real1_f ExpectationPauliAll(std::vector<bitLenInt> bits, std::vector<Pauli> paulis);
 
    virtual real1_f ExpectationUnitaryAll(const std::vector<bitLenInt>& bits,
        const std::vector<std::shared_ptr<complex>>& basisOps, std::vector<real1_f> eigenVals = {})
    {
        return ExpVarUnitaryAll(true, bits, basisOps, eigenVals);
    }
 
    virtual real1_f ExpectationUnitaryAll(
        const std::vector<bitLenInt>& bits, const std::vector<real1_f>& basisOps, std::vector<real1_f> eigenVals = {})
    {
        return ExpVarUnitaryAll(true, bits, basisOps, eigenVals);
    }
 
    virtual real1_f ExpectationBitsFactorized(
        const std::vector<bitLenInt>& bits, const std::vector<bitCapInt>& perms, const bitCapInt& offset = ZERO_BCI);
 
    virtual real1_f ExpectationBitsFactorizedRdm(bool roundRz, const std::vector<bitLenInt>& bits,
        const std::vector<bitCapInt>& perms, const bitCapInt& offset = ZERO_BCI)
    {
        return ExpectationBitsFactorized(bits, perms, offset);
    }
 
    virtual real1_f ExpectationFloatsFactorized(
        const std::vector<bitLenInt>& bits, const std::vector<real1_f>& weights);
 
    virtual real1_f ExpectationFloatsFactorizedRdm(
        bool roundRz, const std::vector<bitLenInt>& bits, const std::vector<real1_f>& weights)
    {
        return ExpectationFloatsFactorized(bits, weights);
    }
 
    virtual real1_f ProbRdm(bitLenInt qubit) { return Prob(qubit); }
    virtual real1_f ProbAllRdm(bool roundRz, const bitCapInt& fullRegister) { return ProbAll(fullRegister); }
    virtual real1_f ProbMaskRdm(bool roundRz, const bitCapInt& mask, const bitCapInt& permutation)
    {
        return ProbMask(mask, permutation);
    }
    virtual real1_f ExpectationBitsAllRdm(
        bool roundRz, const std::vector<bitLenInt>& bits, const bitCapInt& offset = ZERO_BCI)
    {
        return ExpectationBitsAll(bits, offset);
    }
 
    virtual std::map<bitCapInt, int> MultiShotMeasureMask(const std::vector<bitCapInt>& qPowers, unsigned shots);
 
    virtual void MultiShotMeasureMask(
        const std::vector<bitCapInt>& qPowers, unsigned shots, unsigned long long* shotsArray);
 
    virtual void SetBit(bitLenInt qubit, bool value)
    {
        if (value != M(qubit)) {
            X(qubit);
        }
    }
 
    virtual bool ApproxCompare(QInterfacePtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON)
    {
        return SumSqrDiff(toCompare) <= error_tol;
    }
 
    virtual real1_f SumSqrDiff(QInterfacePtr toCompare) = 0;
 
    virtual bool TryDecompose(bitLenInt start, QInterfacePtr dest, real1_f error_tol = TRYDECOMPOSE_EPSILON);
 
    virtual void UpdateRunningNorm(real1_f norm_thresh = REAL1_DEFAULT_ARG) = 0;
 
    virtual void NormalizeState(
        real1_f nrm = REAL1_DEFAULT_ARG, real1_f norm_thresh = REAL1_DEFAULT_ARG, real1_f phaseArg = ZERO_R1_F) = 0;
 
    virtual void Finish() {}
 
    virtual bool isFinished() { return true; }
 
    virtual void Dump() {}
 
    virtual bool isBinaryDecisionTree() { return false; }
 
    virtual bool isClifford() { return false; }
 
    virtual bool isClifford(bitLenInt qubit) { return false; }
 
    virtual bool isOpenCL() { return false; }
 
    virtual bool TrySeparate(const std::vector<bitLenInt>& qubits, real1_f error_tol) { return false; }
    virtual bool TrySeparate(bitLenInt qubit) { return false; }
    virtual bool TrySeparate(bitLenInt qubit1, bitLenInt qubit2) { return false; }
    virtual bool AreFactorized(std::vector<bitLenInt> a, std::vector<bitLenInt> b, bool flushCache = false)
    {
        return false;
    }
    virtual double GetUnitaryFidelity() { return 1.0; }
    virtual void ResetUnitaryFidelity() {}
    virtual void SetSdrp(real1_f sdrp) {};
    virtual void SetNcrp(real1_f ncrp) {};
    virtual void SetReactiveSeparate(bool isAggSep) {}
    virtual bool GetReactiveSeparate() { return false; }
    virtual void SetTInjection(bool useGadget) {}
    virtual bool GetTInjection() { return false; }
    virtual void SetUseExactNearClifford(bool useExact) {}
    virtual bool GetUseExactNearClifford() { return false; }
    virtual void SetNoiseParameter(real1_f lambda) {}
    virtual real1_f GetNoiseParameter() { return ZERO_R1_F; }
    virtual void SetAceMaxQubits(bitLenInt qb) {}
    virtual void SetSparseAceMaxMb(size_t mb) {}
    virtual void SetSparseProbabilityFloor(real1_f p) {}
    virtual void SetStochastic(bool s) {}
    virtual void SetMajorQuadrant(bool q) {}
    virtual void SetMajorQuadrant(bitLenInt t, bool q) {}
    virtual void FlipQuadrant(bitLenInt t) {}
 
    virtual QInterfacePtr Clone() = 0;
 
    virtual QInterfacePtr Copy() { return Clone(); }
 
    virtual void SetDevice(int64_t dID) {};
 
    virtual int64_t GetDevice() { return -1; }
    virtual void SetDeviceList(std::vector<int64_t> dIDs) {};
 
    virtual std::vector<int64_t> GetDeviceList() { return std::vector<int64_t>(); }
 
    bitCapIntOcl GetMaxSize() { return pow2Ocl(sizeof(bitCapInt) * 8); };
 
    virtual real1_f FirstNonzeroPhase()
    {
        complex amp;
        bitCapInt perm = ZERO_BCI;
        do {
            amp = GetAmplitude(perm);
            bi_increment(&perm, 1U);
        } while ((abs(amp) <= REAL1_EPSILON) && (perm < maxQPower));
 
        return (real1_f)std::arg(amp);
    }
 
    virtual void DepolarizingChannelWeak1Qb(bitLenInt qubit, real1_f lambda);
 
};
} // namespace Qrack