qstabilizer.hpp

Go to the documentation of this file.

//
// (C) Daniel Strano and the Qrack contributors 2017-2023. All rights reserved.
//
// Adapted from:
//
// CHP: CNOT-Hadamard-Phase
// Stabilizer Quantum Computer Simulator
// by Scott Aaronson
// Last modified June 30, 2004
//
// Thanks to Simon Anders and Andrew Cross for bugfixes
//
// https://www.scottaaronson.com/chp/
//
// Daniel Strano and the Qrack contributers appreciate Scott Aaronson's open sharing of the CHP code, and we hope that
// vm6502q/qrack is one satisfactory framework by which CHP could be adapted to enter the C++ STL. Our project
// philosophy aims to raise the floor of decentralized quantum computing technology access across all modern platforms,
// for all people, not commercialization.
//
// 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 "qinterface.hpp"
 
namespace Qrack {
 
struct AmplitudeEntry {
    bitCapInt permutation;
    complex amplitude;
 
    AmplitudeEntry(const bitCapInt& p, const complex& a)
        : permutation(p)
        , amplitude(a)
    {
    }
};
 
class QStabilizer;
typedef std::shared_ptr<QStabilizer> QStabilizerPtr;
 
class QStabilizer : public QInterface {
protected:
    unsigned rawRandBools;
    unsigned rawRandBoolsRemaining;
    real1 phaseOffset;
    bitLenInt maxStateMapCacheQubitCount;
    bool isUnitarityBroken;
 
    // Phase bits: 0 for +1, 1 for i, 2 for -1, 3 for -i.  Normally either 0 or 2.
    std::vector<uint8_t> r;
    // Typedef for special type std::vector<bool> compatibility
    typedef std::vector<bool> BoolVector;
    // (2n+1)*n matrix for stabilizer/destabilizer x bits (there's one "scratch row" at the bottom)
    std::vector<BoolVector> x;
    // (2n+1)*n matrix for z bits
    std::vector<BoolVector> z;
 
    typedef std::function<void(const bitLenInt&)> StabilizerParallelFunc;
    typedef std::function<void(void)> DispatchFn;
    void Dispatch(DispatchFn fn) { fn(); }
 
    void ParFor(StabilizerParallelFunc fn, std::vector<bitLenInt> qubits);
 
    void SetPhaseOffset(real1_f phaseArg)
    {
        phaseOffset = (real1)phaseArg;
        const bool isNeg = phaseOffset < 0;
        if (isNeg) {
            phaseOffset = -phaseOffset;
        }
        phaseOffset -= (real1)(((size_t)(phaseOffset / (2 * PI_R1))) * (2 * PI_R1));
        if (phaseOffset > PI_R1) {
            phaseOffset -= 2 * PI_R1;
        }
        if (isNeg) {
            phaseOffset = -phaseOffset;
        }
    }
 
public:
    QStabilizer(bitLenInt n, bitCapInt perm = 0U, qrack_rand_gen_ptr rgp = nullptr,
        complex phaseFac = CMPLX_DEFAULT_ARG, bool doNorm = false, bool randomGlobalPhase = true, bool ignored2 = false,
        int64_t ignored3 = -1, bool useHardwareRNG = true, bool ignored4 = false, real1_f ignored5 = REAL1_EPSILON,
        std::vector<int64_t> ignored6 = {}, bitLenInt ignored7 = 0U, real1_f ignored8 = FP_NORM_EPSILON_F);
 
    ~QStabilizer() { Dump(); }
 
    QInterfacePtr Clone();
    QStabilizerPtr CloneEmpty();
 
    bool isClifford() { return true; };
    bool isClifford(bitLenInt qubit) { return true; };
 
    bitLenInt GetQubitCount() { return qubitCount; }
 
    bitCapInt GetMaxQPower() { return pow2(qubitCount); }
 
    void ResetPhaseOffset() { phaseOffset = ZERO_R1; }
    complex GetPhaseOffset() { return std::polar(ONE_R1, phaseOffset); }
 
    void SetPermutation(bitCapInt perm, complex phaseFac = CMPLX_DEFAULT_ARG);
 
    void SetRandomSeed(uint32_t seed)
    {
        if (rand_generator != NULL) {
            rand_generator->seed(seed);
        }
    }
 
    void SetDevice(int64_t dID) {}
 
    bool Rand()
    {
        if (hardware_rand_generator != NULL) {
            if (!rawRandBoolsRemaining) {
                rawRandBools = hardware_rand_generator->NextRaw();
                rawRandBoolsRemaining = sizeof(unsigned) * bitsInByte;
            }
            --rawRandBoolsRemaining;
 
            return (bool)((rawRandBools >> rawRandBoolsRemaining) & 1U);
        } else {
            return (bool)rand_distribution(*rand_generator);
        }
    }
 
    void Clear()
    {
        x.clear();
        z.clear();
        r.clear();
        phaseOffset = ZERO_R1;
        qubitCount = 0U;
        maxQPower = 1U;
    }
 
protected:
    void rowcopy(const bitLenInt& i, const bitLenInt& k)
    {
        if (i == k) {
            return;
        }
 
        x[i] = x[k];
        z[i] = z[k];
        r[i] = r[k];
    }
    void rowswap(const bitLenInt& i, const bitLenInt& k)
    {
        if (i == k) {
            return;
        }
 
        std::swap(x[k], x[i]);
        std::swap(z[k], z[i]);
        std::swap(r[k], r[i]);
    }
    void rowset(const bitLenInt& i, bitLenInt b)
    {
        std::fill(x[i].begin(), x[i].end(), false);
        std::fill(z[i].begin(), z[i].end(), false);
        r[i] = 0;
 
        if (b < qubitCount) {
            x[i][b] = true;
        } else {
            b -= qubitCount;
            z[i][b] = true;
        }
    }
    void rowmult(const bitLenInt& i, const bitLenInt& k)
    {
        r[i] = clifford(i, k);
        for (bitLenInt j = 0U; j < qubitCount; ++j) {
            x[i][j] = x[i][j] ^ x[k][j];
            z[i][j] = z[i][j] ^ z[k][j];
        }
    }
    uint8_t clifford(const bitLenInt& i, const bitLenInt& k);
 
    void seed(const bitLenInt& g);
 
    AmplitudeEntry getBasisAmp(const real1_f& nrm);
 
    void setBasisState(const real1_f& nrm, complex* stateVec);
 
    void setBasisState(const real1_f& nrm, QInterfacePtr eng);
 
    void setBasisState(const real1_f& nrm, std::map<bitCapInt, complex>& stateMap);
 
    void setBasisProb(const real1_f& nrm, real1* outputProbs);
 
    real1_f getExpectation(const real1_f& nrm, const std::vector<bitCapInt>& bitPowers,
        const std::vector<bitCapInt>& perms, bitCapInt offset);
 
    real1_f getExpectation(
        const real1_f& nrm, const std::vector<bitCapInt>& bitPowers, const std::vector<real1_f>& weights);
 
    void DecomposeDispose(const bitLenInt start, const bitLenInt length, QStabilizerPtr toCopy);
 
    real1_f ApproxCompareHelper(
        QStabilizerPtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON, bool isDiscrete = false);
 
public:
    bitLenInt gaussian();
 
    bitCapInt PermCount() { return pow2(gaussian()); }
 
    void SetQuantumState(const complex* inputState);
    void SetAmplitude(bitCapInt perm, complex amp)
    {
        throw std::domain_error("QStabilizer::SetAmplitude() not implemented!");
    }
 
    void SetRandGlobalPhase(bool isRand) { randGlobalPhase = isRand; }
 
    void CNOT(bitLenInt control, bitLenInt target);
    void CY(bitLenInt control, bitLenInt target);
    void CZ(bitLenInt control, bitLenInt target);
    void AntiCNOT(bitLenInt control, bitLenInt target);
    void AntiCY(bitLenInt control, bitLenInt target);
    void AntiCZ(bitLenInt control, bitLenInt target);
    using QInterface::H;
    void H(bitLenInt qubitIndex);
    using QInterface::X;
    void X(bitLenInt qubitIndex);
    void Y(bitLenInt qubitIndex);
    void Z(bitLenInt qubitIndex);
    void S(bitLenInt qubitIndex);
    void IS(bitLenInt qubitIndex);
    // Swap two bits
    void Swap(bitLenInt qubitIndex1, bitLenInt qubitIndex2);
    // Swap two bits and apply a phase factor of i if they are different
    void ISwap(bitLenInt qubitIndex1, bitLenInt qubitIndex2);
    // Swap two bits and apply a phase factor of -i if they are different
    void IISwap(bitLenInt qubitIndex1, bitLenInt qubitIndex2);
 
    bool ForceM(bitLenInt t, bool result, bool doForce = true, bool doApply = true);
 
    void GetQuantumState(complex* stateVec);
 
    void GetQuantumState(QInterfacePtr eng);
 
    std::map<bitCapInt, complex> GetQuantumState();
 
    void GetProbs(real1* outputProbs);
 
    complex GetAmplitude(bitCapInt perm);
 
    std::vector<complex> GetAmplitudes(std::vector<bitCapInt> perms);
 
    AmplitudeEntry GetAnyAmplitude();
 
    AmplitudeEntry GetQubitAmplitude(bitLenInt t, bool m);
 
    real1_f ExpectationBitsFactorized(
        const std::vector<bitLenInt>& bits, const std::vector<bitCapInt>& perms, bitCapInt offset = 0U);
 
    real1_f ExpectationFloatsFactorized(const std::vector<bitLenInt>& bits, const std::vector<real1_f>& weights);
 
    real1_f ProbPermRdm(bitCapInt perm, bitLenInt ancillaeStart);
 
    real1_f ProbMask(bitCapInt mask, bitCapInt permutation);
 
    bool IsSeparableZ(const bitLenInt& target);
    bool IsSeparableX(const bitLenInt& target);
    bool IsSeparableY(const bitLenInt& target);
    uint8_t IsSeparable(const bitLenInt& target);
 
    using QInterface::Compose;
    bitLenInt Compose(QInterfacePtr toCopy) { return Compose(std::dynamic_pointer_cast<QStabilizer>(toCopy)); }
    bitLenInt Compose(QStabilizerPtr toCopy) { return Compose(toCopy, qubitCount); }
    bitLenInt Compose(QInterfacePtr toCopy, bitLenInt start)
    {
        return Compose(std::dynamic_pointer_cast<QStabilizer>(toCopy), start);
    }
    bitLenInt Compose(QStabilizerPtr toCopy, bitLenInt start);
    void Decompose(bitLenInt start, QInterfacePtr dest)
    {
        DecomposeDispose(start, dest->GetQubitCount(), std::dynamic_pointer_cast<QStabilizer>(dest));
    }
    QInterfacePtr Decompose(bitLenInt start, bitLenInt length);
    void Dispose(bitLenInt start, bitLenInt length) { DecomposeDispose(start, length, (QStabilizerPtr)NULL); }
    void Dispose(bitLenInt start, bitLenInt length, bitCapInt ignored)
    {
        DecomposeDispose(start, length, (QStabilizerPtr)NULL);
    }
    bool CanDecomposeDispose(const bitLenInt start, const bitLenInt length);
    using QInterface::Allocate;
    bitLenInt Allocate(bitLenInt start, bitLenInt length)
    {
        if (!length) {
            return start;
        }
 
        if (start > qubitCount) {
            throw std::out_of_range("QStabilizer::Allocate() cannot start past end of register!");
        }
 
        if (!qubitCount) {
            SetQubitCount(length);
            SetPermutation(0U);
            return 0U;
        }
 
        QStabilizerPtr nQubits = std::make_shared<QStabilizer>(length, 0U, rand_generator, CMPLX_DEFAULT_ARG, false,
            randGlobalPhase, false, -1, hardware_rand_generator != NULL);
        return Compose(nQubits, start);
    }
 
    void NormalizeState(
        real1_f nrm = REAL1_DEFAULT_ARG, real1_f norm_thresh = REAL1_DEFAULT_ARG, real1_f phaseArg = ZERO_R1_F)
    {
        if (!randGlobalPhase) {
            SetPhaseOffset(phaseOffset + (real1)phaseArg);
        }
    }
    void UpdateRunningNorm(real1_f norm_thresh = REAL1_DEFAULT_ARG)
    {
        // Intentionally left blank
    }
 
    real1_f SumSqrDiff(QInterfacePtr toCompare)
    {
        return ApproxCompareHelper(std::dynamic_pointer_cast<QStabilizer>(toCompare));
    }
    bool ApproxCompare(QInterfacePtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON)
    {
        return ApproxCompare(std::dynamic_pointer_cast<QStabilizer>(toCompare), error_tol);
    }
    bool ApproxCompare(QStabilizerPtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON)
    {
        return error_tol >= ApproxCompareHelper(toCompare, error_tol, true);
    }
    bool GlobalPhaseCompare(QInterfacePtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON)
    {
        return GlobalPhaseCompare(std::dynamic_pointer_cast<QStabilizer>(toCompare), error_tol);
    }
    bool GlobalPhaseCompare(QStabilizerPtr toCompare, real1_f error_tol = TRYDECOMPOSE_EPSILON)
    {
        const AmplitudeEntry thisAmpEntry = GetAnyAmplitude();
        real1 argDiff = (real1)abs(
            (std::arg(thisAmpEntry.amplitude) - std::arg(toCompare->GetAmplitude(thisAmpEntry.permutation))) /
            (2 * PI_R1));
        argDiff -= (real1)(size_t)argDiff;
        if (argDiff > (ONE_R1 / 2)) {
            argDiff -= ONE_R1;
        }
        if (FP_NORM_EPSILON >= abs(argDiff)) {
            return false;
        }
        return error_tol >= ApproxCompareHelper(toCompare, error_tol, true);
    }
 
    real1_f Prob(bitLenInt qubit);
 
    void Mtrx(const complex* mtrx, bitLenInt target);
    void Phase(complex topLeft, complex bottomRight, bitLenInt target);
    void Invert(complex topRight, complex bottomLeft, bitLenInt target);
    void MCPhase(const std::vector<bitLenInt>& controls, complex topLeft, complex bottomRight, bitLenInt target);
    void MACPhase(const std::vector<bitLenInt>& controls, complex topLeft, complex bottomRight, bitLenInt target);
    void MCInvert(const std::vector<bitLenInt>& controls, complex topRight, complex bottomLeft, bitLenInt target);
    void MACInvert(const std::vector<bitLenInt>& controls, complex topRight, complex bottomLeft, bitLenInt target);
    void MCMtrx(const std::vector<bitLenInt>& controls, const complex* mtrx, bitLenInt target)
    {
        if (IS_NORM_0(mtrx[1U]) && IS_NORM_0(mtrx[2U])) {
            MCPhase(controls, mtrx[0U], mtrx[3U], target);
            return;
        }
 
        if (IS_NORM_0(mtrx[0U]) && IS_NORM_0(mtrx[3U])) {
            MCInvert(controls, mtrx[1U], mtrx[2U], target);
            return;
        }
 
        throw std::domain_error("QStabilizer::MCMtrx() not implemented for non-Clifford/Pauli cases!");
    }
    void MACMtrx(const std::vector<bitLenInt>& controls, const complex* mtrx, bitLenInt target)
    {
        if (IS_NORM_0(mtrx[1U]) && IS_NORM_0(mtrx[2U])) {
            MACPhase(controls, mtrx[0U], mtrx[3U], target);
            return;
        }
 
        if (IS_NORM_0(mtrx[0U]) && IS_NORM_0(mtrx[3U])) {
            MACInvert(controls, mtrx[1U], mtrx[2U], target);
            return;
        }
 
        throw std::domain_error("QStabilizer::MACMtrx() not implemented for non-Clifford/Pauli cases!");
    }
    void FSim(real1_f theta, real1_f phi, bitLenInt qubit1, bitLenInt qubit2);
 
    bool TrySeparate(const std::vector<bitLenInt>& qubits, real1_f ignored);
    bool TrySeparate(bitLenInt qubit) { return CanDecomposeDispose(qubit, 1U); }
    bool TrySeparate(bitLenInt qubit1, bitLenInt qubit2)
    {
        if (qubit2 < qubit1) {
            std::swap(qubit1, qubit2);
        }
 
        Swap(qubit1, 0U);
        Swap(qubit2, 1U);
 
        const bool toRet = CanDecomposeDispose(0U, 2U);
 
        Swap(qubit2, 1U);
        Swap(qubit1, 0U);
 
        return toRet;
    }
 
    friend std::ostream& operator<<(std::ostream& os, const QStabilizerPtr s);
    friend std::istream& operator>>(std::istream& is, const QStabilizerPtr s);
};
} // namespace Qrack