julia> using BenchmarkTools

julia> n = 100000; a = rand(Int32, n);b = rand(Int32, n); c = rand(Int32, n);

julia> function _allthreesumtozero(a, b, c)
         @inbounds @simd for i in 1:length(a)
           iseven(a[i] + b[i] + c[i]) || return false
         end
         true
       end
_allthreesumtozero (generic function with 1 method)

julia> function _allthreesumtozero_optimized(a, b, c)
           @inbounds @simd for i in 1:length(a)
               ((a[i] + b[i] + c[i]) & 1) == 0 || return false
           end
           true
       end
_allthreesumtozero_optimized (generic function with 1 method)

julia> println(@benchmark _allthreesumtozero($a, $b, $c))
Trial(10.601 ns)

julia> println(@benchmark _allthreesumtozero_optimized($a, $b, $c))
Trial(3.505 ns)

 Test Summary:                                  | Pass  Fail  Total  Time
noisycircuits                                  |   90     1     91  7.2s
  Noisy Gates                                  |    2            2  1.5s
  Monte Carlo Purification examples            |    9            9  0.5s
  Perturbative expansion Purification examples |   11     1     12  4.3s
    Comparison to MC                           |    8     1      9  2.1s
    Symbolic                                   |    3            3  2.2s
  Measurements                                 |   53           53  0.5s
  Classical Bits                               |   15           15  0.2s

function apply!(state::GeneralizedStabilizer, gate::AbstractPauliChannel; prune_threshold::Union{Nothing, Float64}=nothing)
    if prune_threshold === nothing
        prune_threshold = 1e-14  # Default value
    end
    dict = state.destabweights
    stab = state.stab
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero) 
    for ((dᵢ,dⱼ), χ) in dict # the state
        for ((Pₗ,Pᵣ), w) in zip(gate.paulis, gate.weights) # the channel
            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ, stab)
            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ, stab)
            c = (dot(dₗˢᵗᵃᵇ, dᵢ) + dot(dᵣˢᵗᵃᵇ, dⱼ)) * 2
            dᵢ′ = dₗ .⊻ dᵢ
            dⱼ′ = dᵣ .⊻ dⱼ
            χ′ = χ * w * (-tone)^c * (im)^(-phaseₗ + phaseᵣ + 4)
            if abs(χ′) >= prune_threshold
                if haskey(newdict, (dᵢ′, dⱼ′))
                    newdict[(dᵢ′, dⱼ′)] += χ′
                else
                    newdict[(dᵢ′, dⱼ′)] = χ′
                end
            end
    end
    filter!(x -> abs(x[2]) >= prune_threshold, newdict)
    state.destabweights = newdict
    state
end

function apply!(state::GeneralizedStabilizer, gate::AbstractPauliChannel; prune_threshold::Union{Nothing, Float64}=nothing)
    if prune_threshold === nothing
        prune_threshold = 1e-14  # Default value
    end
    dict = state.destabweights
    stab = state.stab
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero)
    for ((dᵢ,dⱼ), χ) in dict # the state
        for ((Pₗ,Pᵣ), w) in zip(gate.paulis, gate.weights) # the channel
            phaseₗ, dₗ, dₗˢᵗᵃᵇ = rowdecompose(Pₗ,stab)
            phaseᵣ, dᵣ, dᵣˢᵗᵃᵇ = rowdecompose(Pᵣ,stab)
            c = (dot(dₗˢᵗᵃᵇ,dᵢ) + dot(dᵣˢᵗᵃᵇ,dⱼ))*2
            dᵢ′ = dₗ .⊻ dᵢ
            dⱼ′ = dᵣ .⊻ dⱼ
            χ′ = χ * w * (-tone)^c * (im)^(-phaseₗ+phaseᵣ+4)
           if abs(χ′) >= prune_threshold
               newdict[(dᵢ′,dⱼ′)] = get!(newdict,(dᵢ′,dⱼ′),0)+χ′
           end
    end
    filter!(x -> abs(x[2]) >= prune_threshold, newdict)
    state.destabweights = newdict
    state
end

1) Full stabilizer `T = (S, D)` where `S` is full Stabilizer and `D` is full Destabilizer.
2) `T = (S, D)` and Pauli operator `P ∈ Gₙ`, `P = ads` in O(n²) where `α ∈ {±1, ±i}`, `d ∈ D`, and `s ∈ S`.
3) `B(T ) = B(S, D)` stabilizers a unique pure stabilizer state, `|ψₛ⟩`, or `ρₛ = |ψₛ⟩⟨ψₛ| `
4) `dₖρₛdₖ = (½)ⁿ( I - sₖ) ∏ (I + sᵢ) `where `∏` is from `ᵢ ≠ ₖ`
5) `B(T) = B(S, D) = {d|ψₛ⟩  d ∈ D}` ; stabilizer basis corresponding to T.
6) Generalized Stabilizer for τ(arbitrary state) is a two-tuple` (χ, B(S, D))`
where `χ` is a density matrix and ` B(S, D) = B(T )` is the basis in which `χ` is expressed.
7) arbitrary state: `τ = ∑∑ χᵢⱼ dᵢ |ψₛ⟩⟨ψₛ| dⱼ = ∑∑ χᵢⱼ dᵢ ρₛ dⱼ`
8 Inverse sparsity: Λ(τ) = min B(S,D) Λ(χ) = min S Λ(χ) where Λ(χ) is # of non zero elements in density matrix `χ`
9) Bounds: rank(τ) ≤ Λ(τ) ≤ 4ⁿ
10) Pauli channel: `E(τ ) = ∑ϕᵢⱼ Bᵢ ρ Bⱼ†`
where `ϕᵢ` is complex number and `Bᵢ `is Pauli operators `Bᵢ ∈ Gₙ.`
11) Λ(E) ≡ Λ(ϕ) is Pauli Channel complexity.
12) Anything a S representation can do is also possible with a (χ, B(S, D)).
13) This follows from 10) that χ is all zeros except for a single 1 along on the diagonal with appropriate stabilizer basis B(S, D)

1) Clifford Gates: Suppose Clifford circuit C.
CτC† = ∑ χᵢⱼ CdᵢC† CρₛC† CdⱼC† = ∑ χᵢⱼ d'ᵢ ρₛ' dⱼ'
------------------------------------------------------------
input: τ and Clifford gate, G
output: S', D'
Procedure:
Conjugate each generator of  S: S' = ⟨Gs₁G†, Gs₂G†, ... , GsₙG†⟩
Conjugate each generator of D: D' = ⟨Gd₁G†, Gd₂G†, ... , GdₙG†⟩
------------------------------------------------------------

2) Pauli measurements
----------------------------------------------------------------------------------------------------------------
input: τ = (χ, B(S, D)) and α, b, c from rowdecompose where M = a*d_b*s_c , M ∈ Gₙ
output:  τ' = (χ, B(S', D'))
procedure:
Initialize sparse matrix χ'
if b = 0 then
    for all χᵢⱼ ≠ 0 do
        if α(−1)^(i·c) = 1 and α(−1)^(j·c) = 1 then
            χᵢⱼ' ← χᵢⱼ'  + χᵢⱼ
else
k ← vector of zeros with a single 1 in the location of b’s first 1
    for all χᵢⱼ ≠ 0 do
    x ← i
    y ← j
    q ← 1
    if i · k = 1 then
        q ← α(−1)^(i·c) q
        x ← i + b
    if j · k = 1 then  
        q ← α∗(−1)^(j·c) q
        y ← j + b
     χ_xy' ← χ_xy'  + ½qχᵢⱼ
Return χ', S', D'
Note: (S', D') can be found from (S, D) by traditional stabilizer update via Algorithm 1
----------------------------------------------------------------------------------------------------------------

3) Pauli Channels:
--------------------------------------------------------------------------------
input: τ = (χ, B(S, D)), {αₕ, bₕ, cₕ}ₕ for each Bₕ
output:  χ' s.t. τ = (χ', B(S, D))
procedure:
Initialize sparse matrix χ
    for all ϕₖₗ ≠ 0 do
        for all χᵢⱼ ≠ 0 do
            x ← i + bₖ
            y ← j + bₗ
            χ_xy' ← χ_xy'  + [αₕαₗ*(−1)^(i·cₖ+j·cₗ)]ϕₖₗχᵢⱼ
Return χ'
--------------------------------------------------------------------------------

I was wondering why the destabweights are updated following `apply!(GeneralizedStabilizer(S"-X"), pcT)`.
The answer seems to be the following (taken from the paper):

1) An algorithm that uses just a density matrix representation could also use sparsity to improve its runtime for simulation of such non-stabilizer states through stabilizer circuits.
2) However, simulating g Clifford gates will generally result in an exponential (i.e. 2 ^O(g)) decrease in the sparsity of the density matrix representation, such that, when it comes time to measure the state at the conclusion of the simulation, the calculation will be intractable.
3) Our method avoids this by updating the stabilizer basis instead of the density matrix directly.
4) The price we pay is having to interpret the concluding measurements in terms of our modified basis, but this is polynomial in n by the decomposition lemma.

function expect(p::PauliOperator, s::GeneralizedStabilizer)
    dict = s.destabweights
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    tone = one(dtype)
    newdict = typeof(dict)(tzero)
    phase, b, c = rowdecompose(p, s.stab)
    if all(x -> x == 0, b)
        for ((dᵢ, dⱼ), χ) in dict
            if ((-tone)^(dot(dᵢ, c)))*im^phase == 1 && ((-tone)^(dot(dⱼ, c)))*im^phase == 1
                newdict[(dᵢ, dⱼ)] = get(newdict, (dᵢ, dⱼ), 0) + χ
            end
        end
        s.destabweights = newdict
        return s
    else
        k = _create_k(b)
        for ((dᵢ, dⱼ), χ) in dict
            q = 1
            dᵢ′ = copy(dᵢ)
            dⱼ′ = copy(dⱼ)
            if dot(dᵢ, k) == 1
                q *= ((-tone)^(dot(dᵢ, c)))*im^phase
                dᵢ′ = b .⊻ dᵢ
            end
            if dot(dⱼ, k) == 1
                q *= (-tone)^(dot(dⱼ, c))*conj(im^phase)
                dⱼ′ = b .⊻ dⱼ
            end
            newdict[(dᵢ′, dⱼ′)] = get(newdict, (dᵢ′,dⱼ′), 0) + χ * q / 2
        end
         s.destabweights = newdict
        return s
    end
end

function _create_k(b)
    k = zeros(Int, length(b))
    idx = findfirst(x -> x == 1, b)
    if idx !== nothing
        k[idx] = true
    end
    return k
end

julia> s = S"X"; p = P"Y";

julia> gs = GeneralizedStabilizer(s)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.853553+0.0im | + _ | + _
 0.146447+0.0im | + Z | + Z

julia> expect(p, gs)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
+ X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _

function conj_destabs!(state::GeneralizedStabilizer, gate::AbstractCliffordOperator)
    dict = state.destabweights
    dtype = _dictvaltype(dict)
    tzero = zero(dtype)
    newdict = typeof(dict)(tzero)
    gate_fields = fieldnames(typeof(gate))
    g_field = nothing
    for field in gate_fields
        if getfield(gate, field) isa Vector{UInt64}
            g_field = field
            break
        end
    end
    if g_field === nothing
        error("No Vector{UInt64} field found in the gate object")
    end
    g_part = getfield(gate, g_field)
    C = Vector{Bool}(g_part)
    for ((dᵢ, dⱼ), χ) in dict
        # Compute dᵢ' and dⱼ'
        dᵢ′ = C .* dᵢ .* C
        dⱼ′ = C .* dⱼ .* C
        newdict[(dᵢ′, dⱼ′)] = χ
    end
    state.destabweights = newdict
    state
end

julia> gs = GeneralizedStabilizer(S"-Y")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.853553+0.0im | + _ | + _
 0.146447+0.0im | + Z | + Z

julia> QuantumClifford.conj_destabs!(gs, P"-X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.146447+0.0im | + Z | + Z

julia> QuantumClifford.conj_destabs!(gs, P"X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- Y
with ϕᵢⱼ | Pᵢ | Pⱼ:
 0.853553+0.0im | + _ | + _
 0.0+0.353553im | + _ | + Z
 0.0-0.353553im | + Z | + _
 0.146447+0.0im | + Z | + Z

julia> gs = GeneralizedStabilizer(S"Z")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> apply!(gs, pcT)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> QuantumClifford.conj_destabs!(gs, P"Y")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
+ Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> sm = GeneralizedStabilizer(S"-X")
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> sm ⊗ sm
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z_
+ _Z
𝒮𝓉𝒶𝒷
- X_
- _X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + __ | + __

julia> tHadamard*sm
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ X
𝒮𝓉𝒶𝒷
- Z
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> pcT⊗P"X"
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
 0.853553+0.353553im | + _X
 0.146447-0.353553im | + ZX

julia> copy(sm)
A mixture ∑ ϕᵢⱼ Pᵢ ρ Pⱼ† where ρ is
𝒟ℯ𝓈𝓉𝒶𝒷
+ Z
𝒮𝓉𝒶𝒷
- X
with ϕᵢⱼ | Pᵢ | Pⱼ:
 1.0+0.0im | + _ | + _

julia> tId1 * pcT
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
 0.853553+0.353553im | X₁ ⟼ + X
Z₁ ⟼ + Z
 0.146447-0.353553im | X₁ ⟼ - X
Z₁ ⟼ + Z

julia> tHadamard * pcT
A unitary Pauli channel P = ∑ ϕᵢ Pᵢ with the following branches:
with ϕᵢ | Pᵢ
0.853553+0.353553im | X₁ ⟼ + Z
Z₁ ⟼ + X
0.146447-0.353553im | X₁ ⟼ + Z
Z₁ ⟼ - X