Functions

Base.:*(b::Bra{N}, p::Union{Pauli{N}, PauliBasis{N}})

Apply a Pauli to a bra from the right: ⟨b|P = c ⟨b'| with b' = b ⊻ x. Returns a (coefficient, Bra) tuple — a Pauli connects each basis bra to exactly one basis bra, so no sum type is needed.

Base.:*(p::Union{Pauli{N}, PauliBasis{N}}, k::Ket{N})

Apply a Pauli to a computational-basis ket: P|k⟩ = c |k'⟩ with k' = k ⊻ x. Returns a (coefficient, Ket) tuple — a Pauli maps each basis state to exactly one basis state, so no sum type is needed.

Base.:*(p1::Pauli{N}, p2::Pauli{N}) where {N}

Multiply two Pauli's together

Base.:*(d1::Union{Dyad{N}, DyadBasis{N}}, d2::Union{Dyad{N}, DyadBasis{N}})

Dyad product: |i⟩⟨j| ⋅ |k⟩⟨l| = δ_jk |i⟩⟨l|, with coefficients multiplied through.

Base.:*(d1::DyadSum{N,T}, d2::DyadSum{N,T}) where {N,T}

Multiply two DyadSums.

Base.:*(ps1::PauliSum{N}, ps2::PauliSum{N}) where {N}

Multiply two PauliSums.

Base.:*(O::PauliSum{N,T}, k::Ket{N})

Apply a sum of Paulis to a basis state, returning a KetSum with one entry per distinct x bitstring in O.

Base.:*(ps1::PauliSum{N, T}, ps2::Adjoint{<:Any, PauliSum{N, T}}) where {N, T}

Multiply two PauliSums.

Base.:*(ps1::Adjoint{<:Any, PauliSum{N, T}}, ps2::PauliSum{N, T}) where {N, T}

Multiply two PauliSums.

Base.:*(ps1::Adjoint{<:Any, PauliSum{N, T}}, ps2::Adjoint{<:Any, PauliSum{N, T}}) where {N, T}

Multiply two PauliSums.

Base.:+(p::Bra{N}, q::Bra{N}) where N

Add two Ket's together, return a KetSum

Base.:+(p::Ket{N}, q::Ket{N}) where N

Add two Ket's together, return a KetSum

Base.:+(p::DyadBasis{N}, q::DyadBasis{N}) where N

Add two Dyad's together, return a DyadSum

Base.:+(p::Pauli{N}, q::Pauli{N}) where N

Add two Pauli's together, return a PauliSum

Base.:+(ps1::DyadSum{N}, ps2::DyadSum{N}) where {N}

Add two DyadSums.

Base.:-(ps1::PauliSum, ps2::PauliSum)

Subtract two PauliSums.

Base.:-(ps1::PauliSum, ps2::PauliSum)

Subtract two PauliSums.

Base.:-(p::Pauli{N}) where {N}

Negate the scalar of p.

Base.adjoint(p::Pauli)


Pauli{N}(s,z,x)  =  s ⋅ z₁...|x₁... 
                 =  s ⋅ i^-θs ⋅ P₁⊗...⊗Pₙ
                 =  coeff ⋅ P₁⊗...⊗Pₙ

Since the PauliBasis is Hermitian, we have that Pauli' = coeff' ⋅ P₁⊗...⊗Pₙ

Base.isapprox(ps1::PauliSum{N}, ps2::PauliSum{N}; atol=1e-14, rtol=0) where {N}

Compare two PauliSums for approximate equality based on the norm of their difference.

rand(Pauli{N})

Generate a random Pauli{N} with random z, x bitstrings and a random complex scalar.

Base.show(io::IO, P::Union{Ket, Bra})

Print the ket/bra string representation (e.g., |010> or <010|).

Base.show(io::IO, v::KetSum{N,T}) where {N,T}

Print each ket and its coefficient.

Base.string(p::PauliBasis{N}) where N

Return a string representation. Y sites are displayed as Y.

Base.string(p::Pauli{N}) where N

Display, y = iY

Base.sum!(p1::PauliSum{N}, p2::PauliSum{N}) where {N}

Add two PauliSums.

Base.sum!(v1::SparsePauliVector, v2::SparsePauliVector)

In-place accumulate v1 += v2 — the flat-storage mergewith!(+, ...). v2 is already sorted, so this is a straight merge (no sort); allocation-free once v1's buffers cover the union.

LinearAlgebra.diag(ps::PauliSum{N,T}) where {N,T}

Return a new PauliSum containing only the diagonal terms (those with x == 0, i.e., only Z and I factors).

LinearAlgebra.dot(v1::KetSum{N,T}, v2::KetSum{N,TT}) where {N,T,TT}

Compute the inner product v1' * v2, iterating over the shorter dictionary.

ishermitian(p::Pauli; thresh=1e-16)

Return true if the coefficient of p is real (within thresh).

LinearAlgebra.norm(ks::KetSum{N,T}, p::Real=2) where {N,T}

Compute the p-norm of the coefficient vector of a KetSum.

LinearAlgebra.norm(ps::PauliSum{N,T}, p::Real=2) where {N,T}

Compute the p-norm of the coefficient vector of a PauliSum.

• p=2 (default): Frobenius norm sqrt(sum(|c_i|^2))
• p=1: L1 norm sum(|c_i|)
• p=Inf: max norm max(|c_i|)

S_gate(p::Union{PauliSum{N}, KetSum{N}}, q) where N

Apply S gate (π/2 Z-rotation) on qubit q.

T_gate(p::Union{PauliSum{N}, KetSum{N}}, q) where N

Apply T gate (π/4 Z-rotation) on qubit q.

X_gate(p::KetSum{N}, q) where N

Apply Pauli X gate on qubit q (Schrödinger picture).

X_gate(p::PauliSum{N}, q) where N

Apply Pauli X gate on qubit q (Heisenberg picture).

X_gate_to_paulis(N, q)

Return (generators, angles) for an X gate on qubit q.

Y_gate(p::KetSum{N}, q) where N

Apply Pauli Y gate on qubit q (Schrödinger picture).

Y_gate(p::PauliSum{N}, q) where N

Apply Pauli Y gate on qubit q (Heisenberg picture).

Z_gate(p::KetSum{N}, q) where N

Apply Pauli Z gate on qubit q (Schrödinger picture).

Z_gate(p::PauliSum{N}, q) where N

Apply Pauli Z gate on qubit q (Heisenberg picture).

Z_gate_to_paulis(N, q)

Return (generators, angles) for a Z gate on qubit q.

Gather the pending appends az/ax/ac[1:an] into the workspace as (z, x, c) triples. Returns the count. Allocation-free.

_merge_spv!(v, m, f) -> (n_in, n_out)

Two-pointer merge of the sorted live buffer with m sorted workspace triples into scratch, summing coefficients of equal keys (equal-key runs in the appends are possible across rotations of one window), dropping outputs the strict filter rejects, then swapping scratch and live (pointer swap). Restores the sorted, duplicate-free live invariant and resets the append cursor. Allocation-free in steady state; grows the live buffers (chunked doubling) if the merged population exceeds capacity — a boundary allocation, never a hot-loop one.

_rotate_range!(z, x, c, lo, hi, gz, gx, n_g, cosθ, sinθ,
               dz, dx, dc, cur, cap_end, f) -> (cur, created, overflowed)

Sweep terms lo:hi of one buffer under the rotation exp(iθ/2 G)·O·exp(-iθ/2 G) (the evolve!(O, G, θ) convention — cos(θ)/sin(θ), half-angle absorbed): commuting terms untouched; anticommuting terms cos-scaled in place, with the sin branch (bits G ⊻ P, sign i·i^k = ±1 computed purely from bits — the fused-phase identity from commutator.jl) appended at cur in the destination arrays unless the local filter drops it. n_g is count_ones(gz & gx), precomputed once per rotation.

Zero-allocation hot path. overflowed only fires if the driver's capacity precheck was skipped or wrong; the driver treats it as an error.

In-place quicksort (median-of-3, insertion sort below 24, recurse-smaller / iterate-larger) of tuples by their first two fields — the x-major key for merge-workspace triples. Hand-rolled because Base's default QuickSort allocates scratch; this is the swap point for a future radix sort. Allocation-free.

_unshuffle_ws!(v, m, gz, gx)

Sort v.ws[1:m] by the x-major key, exploiting rotation structure: sin-branch appends are generated by scanning the sorted live buffer and XORing each key with the generator mask (gz, gx), so the gathered triples are already sorted with respect to key ⊻ mask. One streaming block-swap pass per set mask bit (highest combined bit first) restores natural order — weight(G) linear passes instead of an O(m log m) comparison sort. Ping-pongs between v.ws and the lazily grown v.ws2; the result always ends in v.ws. Falls back to nothing-to-do for the identity mask.

_word_type(N)

Packed bit-word type for an N-qubit register: UInt64 when the bits fit (halves bit bandwidth), UInt128 up to 128 qubits.

Beyond 128 qubits the intended extension is BitIntegers.jl (UInt256, UInt512, ... are primitive Unsigned types supporting ⊻, &, count_ones, < natively), so every kernel is written generically over W<:Unsigned with standard operators — adding larger lattices later means extending this one function (plus a wider PauliBasis for conversions).

anticommutator(A::PauliSum{N,T}, B::PauliSum{N,T}) where {N,T}

Compute the anticommutator {A, B} = AB + BA using an optimized single-pass algorithm.

For Pauli basis elements, the anticommutator contribution is (phase + conj(phase)) * c_i * c_j * P_k = 2 * Re(phase) * c_i * c_j * P_k. Anti-commuting pairs (where the phase is purely imaginary) are skipped entirely.

Uses the same bitstring-only phase computation as commutator.

bit_flip_channel!(O::PauliSum{N}, p; qubits=nothing) -> O

Apply (in place) the i.i.d. single-qubit bit-flip channel with probability p ∈ [0,1]: ρ ↦ (1-p)ρ + p·XρX. Suppresses Y and Z Pauli letters by (1 - 2p) per qubit.

bit_phase_flip_channel!(O::PauliSum{N}, p; qubits=nothing) -> O

Apply (in place) the i.i.d. single-qubit bit-phase-flip channel with probability p ∈ [0,1]: ρ ↦ (1-p)ρ + p·YρY. Suppresses X and Z Pauli letters by (1 - 2p) per qubit.

boson_to_paulis(nqubits::Integer; verbose=0)

Binary (truncated) encoding of a bosonic raising operator b† into nqubits qubits, giving access to a Fock space of dimension d = 2^nqubits:

b† = ∑_{n=0}^{d-2} √(n+1) |n+1⟩⟨n|

The bosonic occupation number n is encoded in the qubit register with qubit 1 as the most-significant bit (i.e. n = q₁·2^{nq-1} + q₂·2^{nq-2} + … + q_{nq}·2^0). Adjoint the result to obtain the lowering operator b = (b†)†.

The encoding is the same one used in the original spin–boson workflow on the spinboson branch (formerly boson_binary_transformation).

Returns a PauliSum{nqubits, ComplexF64}.

check_spv(v::SparsePauliVector)

Verify the structural invariants: consistent buffer lengths, and strictly increasing (sorted, duplicate-free) live keys. Errors on violation. Testing / debugging utility.

clip!(ps::DyadSum; thresh=1e-16)

clip!(ks::KetSum{N}; thresh=1e-16) where {N}

clip!(ps::PauliSum{N}; thresh=1e-16)

cnot(p::KetSum{N}, c::Int, t::Int) where N

Apply a CNOT gate (control c, target t) via Pauli rotations (Schrödinger picture).

cnot(p::PauliSum{N}, c::Int, t::Int) where N

Apply a CNOT gate (control c, target t) via Pauli rotations (Heisenberg picture).

cnot_to_paulis(N, c::Int, t::Int)

Return (generators, angles) for a CNOT gate (control c, target t) in an N-qubit system.

coeff(p::Pauli)

Return the coefficient from the product of the scalar times the inverse symplectic_phase

coeff_clip!(ks::KetSum{N}, thresh::Real) where {N}

Remove Ket terms with |coefficient| <= thresh.

coeff_clip!(ps::PauliSum{N}, thresh::Real)

Remove Pauli terms with |coefficient| <= thresh.

coeff_clip!(ps::DyadSum{N,T}, thresh::Real) where {N,T}

Remove Dyad terms with |coefficient| <= thresh.

commutator!(out::SparsePauliVector, A::SparsePauliVector, B::SparsePauliVector)

Compute [A, B] = AB - BA into out (contents overwritten), using out.ws as the pairing buffer. Allocation-free: errors if out's workspace or live capacity cannot hold the result — size out with SparsePauliVector(N, T; capacity=...) or use commutator(A, B), which grows automatically. Requires complex T (commutator coefficients are imaginary; store -i[A,B] yourself if you need a real-typed engine).

commutator(A::PauliSum{N,T}, B::PauliSum{N,T}) where {N,T}

Compute the commutator [A, B] = AB - BA using an optimized single-pass algorithm.

For Pauli basis elements, P_i * P_j = phase * P_k and P_j * P_i = conj(phase) * P_k, so the commutator contribution is (phase - conj(phase)) * c_i * c_j * P_k = 2i * Im(phase) * c_i * c_j * P_k. Commuting pairs (where the phase is real) are skipped entirely.

Optimization

Rather than constructing intermediate Pauli objects and doing complex arithmetic, the phase exponent k (where phase = i^k) is computed directly from bitstring operations. For non-commuting pairs, k is always odd (phase is ±i), so the imaginary part is ±1 and can be determined by a single integer comparison.

commute(p1::PauliBasis, p2::PauliBasis)

Return true if the two Pauli strings commute, via the symplectic parity test popcount(x₁ & z₂) ≡ popcount(z₁ & x₂) (mod 2) — no product is formed.

covariance(A::PauliSum{N}, B::PauliSum{N}, ψ::Ket{N}) where N

Compute the covariance of observables A and B in state ψ: <A†B> - <A†><B>.

dephasing_channel!(O::PauliSum{N}, p; qubits=nothing) -> O

Apply (in place) the i.i.d. single-qubit dephasing (phase-flip) channel with probability p ∈ [0,1]: ρ ↦ (1-p)ρ + p·ZρZ. Suppresses X and Y Pauli letters by (1 - 2p) per qubit.

Aliased as phase_flip_channel!.

depolarizing_channel!(O::PauliSum{N}, p; qubits=nothing) -> O

Apply (in place) the i.i.d. single-qubit depolarizing channel with parameter p ∈ [0,1], in the Heisenberg picture. Convention: pI = 1-p, pX = pY = pZ = p/3. Equivalent to scaling each Pauli P by (1 - 4p/3)^w, where w is the number of non-identity qubits of P within qubits.

Weight-decay equivalence

Applied to all qubits, this realizes the CPTP map O ↦ Σ_P exp(-γΔt·w(P)) ⟨P,O⟩ P / d for any γ,Δt ≥ 0 by choosing p = (3/4)·(1 - exp(-γΔt)) (see depolarizing_p_for_weight_decay). A thresholded weight damper of the form "only damp when w > lmax" is not CPTP in general and is intentionally not provided here.

depolarizing_p_for_weight_decay(rate, dt) -> Float64

Return the depolarizing parameter p such that the i.i.d. depolarizing channel depolarizing_channel!(O, p) (applied to all qubits) scales each Pauli term P by exp(-rate*dt · w(P)), where w(P) is the number of non-identity qubits of P.

This is the CPTP realization of exponential-in-weight damping. A thresholded weight damper (damp only when w(P) > lmax) is not CPTP in general — its implied multi-qubit Pauli probabilities have negative entries — and is intentionally not provided in this module.

evolve!(K::KetSum{N, ComplexF64}, G::PauliBasis{N}, θ::Real)

In-place Schrödinger-picture evolution: K → exp(-iθ/2 G) K

Modifies K in place. Element type must be ComplexF64 so that the imaginary contribution from the sine branch can be stored back in K.

evolve!(O::PauliSum{N, T}, G::PauliBasis{N}, θ::Real)

In-place Heisenberg-picture evolution: O(θ) = exp(iθ/2 G) O exp(-iθ/2 G)

evolve!(O::SparsePauliVector{N}, generators::Vector{PauliBasis{N}}, angles;
        window=1, truncation=NoTruncation(), local_truncation=NoTruncation(),
        correction=NoCorrection(), counters=nothing)

Heisenberg-picture sequence evolution on flat storage, applying generators left to right (same convention as evolve(::PauliSum, generators, angles); sequences from trotterize/qdrift plug in directly).

Rotations append sin branches under the loose local_truncation (applied per term at append time — weight cutoffs are exact there, coefficient cutoffs act on unmerged duplicates); every window rotations the appends are sort-merged into the live buffer under the strict truncation, and correction accumulates the truncation losses.

At window = 1 (the default) this reproduces the Dict path evolve!(O, g, θ); truncate!(O, truncation, correction) per rotation exactly, for every truncation strategy. window > 1 trades truncation cadence for speed: deduplication and truncation happen once per window.

If a rotation's worst-case appends cannot fit the append buffer, an early merge is triggered (harmless: it only changes truncation cadence), growing the buffer at the boundary if the population genuinely needs more room. The steady-state hot path allocates zero bytes — pass a WindowCounters to verify (counters.allocd).

evolve!(O::SparsePauliVector{N}, G::PauliBasis{N}, θ::Real)

In-place Heisenberg-picture evolution O ← exp(iθ/2 G) O exp(-iθ/2 G), the flat-storage analogue of evolve!(::PauliSum, G, θ) (identical semantics: sin branches are deduplicated into the sum immediately). Zero-allocation once buffers have warmed up.

evolve(O::SparsePauliVector{N}, generators, angles; kwargs...)

Non-mutating sequence evolution (see evolve!).

evolve(O::SparsePauliVector{N}, G::PauliBasis{N}, θ::Real)

Non-mutating single-rotation evolution (see evolve!).

evolve(K::KetSum{N,T}, generators::Vector{PauliBasis{N}}, angles::Vector{<:Real})

Schrödinger-picture sequence evolution: applies generators left to right, producing

Uₙ ⋯ U₂ U₁ |K⟩

where Uₖ = exp(-iθₖ/2 Gₖ). The effective unitary is Uₙ⋯U₁ (reversed order).

Note: sequences from trotterize and qdrift use the Heisenberg convention. To use them with KetSum, reverse the sequence:

evolve(K, reverse(generators), reverse(angles))

evolve(K::KetSum{N, T}, G::PauliBasis{N}, θ::Real)

Schrödinger-picture evolution: K(θ) = exp(-iθ/2 G) K

Applies the unitary exp(-iθ/2 G) to a KetSum state vector.

evolve(O::PauliSum{N,T}, generators::Vector{PauliBasis{N}}, angles::Vector{<:Real};
       truncation::TruncationStrategy=NoTruncation(),
       correction::CorrectionAccumulator=NoCorrection())

Heisenberg-picture sequence evolution: applies generators left to right, producing

U_n† ⋯ U_1† O U_1 ⋯ U_n

where Uₖ = exp(-iθₖ/2 Gₖ). The effective right-side unitary is U₁U₂⋯Uₙ.

Sequences from trotterize and qdrift are designed for this convention.

evolve(O::PauliSum{N, T}, G::PauliBasis{N}, θ::Real)

Heisenberg-picture evolution: O(θ) = exp(iθ/2 G) O exp(-iθ/2 G)

Commuting terms pass through unchanged. Non-commuting terms branch: O(θ) = cos(θ)·O - i·sin(θ)·G·O

expectation_value(O::AnyPauliSum, v::KetSum)

⟨v|O|v⟩ for a linear combination of basis states, including all cross terms ⟨k₂|P|k₁⟩. Works for both PauliSum and SparsePauliVector.

expectation_value(p::Union{PauliBasis{N}, Pauli{N}}, ket::Ket{N})

⟨k|P|k⟩ for a computational-basis ket: nonzero only for diagonal Paulis (x == 0), where it equals (-1)^popcount(z & k) ⋅ coeff(p).

expectation_value(p::Union{PauliBasis{N}, Pauli{N}}, d::Union{Dyad{N}, DyadBasis{N}})

tr(P ⋅ |ket⟩⟨bra|) = ⟨bra|P|ket⟩: nonzero only when the Pauli's x string connects the two basis states (ket ⊻ bra == x). Note the dyad's own coefficient is not included — DyadSum methods carry coefficients in the dictionary values.

expectation_value(p::PauliSum{N,T}, d::Union{Ket{N}, Dyad{N}, DyadBasis{N}})

Expectation value of a sum of Paulis against a basis state or dyad — the coefficient-weighted sum of the per-term expectation values.

expectation_value(v::SparsePauliVector{N}, ψ::Ket{N})

⟨ψ|v|ψ⟩ for a computational-basis ket. Includes any pending (unmerged) appends — expectation is linear, so the pre-merge state evaluates exactly. Allocation-free.

find_top_k(O::PauliSum{N,T}, k::Int) where {N,T}

Return the k terms with largest absolute coefficients, sorted by decreasing |c|. Returns a Vector{Pair{PauliBasis{N}, T}}. Efficient for k << length(O).

get_majorana_weight_counts(O::PauliSum{N}) where N

Return a vector of length 2N+1 where entry i contains the number of terms with Majorana weight i-1.

get_majorana_weight_probs(O::PauliSum{N}) where N

Return a vector of length 2N+1 where entry i contains the sum of |c|² for terms with Majorana weight i-1.

get_on_bits(x::Integer)

Return the (1-based) positions of the set bits of x, in increasing order. Used to translate z/x bitstrings into site indices, e.g. for display.

get_weight_counts(O::PauliSum{N}) where N

Return a vector of length N+1 where entry i contains the number of terms with Pauli weight i-1.

get_weight_probs(O::PauliSum{N}) where N

Return a vector of length N+1 where entry i contains the sum of |c|² for terms with Pauli weight i-1.

hadamard(p::KetSum{N}, q::Int) where N

Apply a Hadamard gate on qubit q via Pauli rotations (Schrödinger picture).

hadamard(p::PauliSum{N}, q::Int) where N

Apply a Hadamard gate on qubit q via Pauli rotations (Heisenberg picture).

hadamard_to_paulis(N, q::Int)

Return (generators, angles) for a Hadamard gate on qubit q in an N-qubit system.

inner_product(O1::PauliSum{N,T}, O2::PauliSum{N,T}) where {N,T}

Evaluate the Liouville space inner product: tr(O1'*O2)

inner_product(v1::SparsePauliVector, v2::SparsePauliVector)

Liouville inner product tr(v1†·v2) over shared basis terms (two-pointer coefficient dot product; both inputs sorted). Allocation-free.

jordan_wigner(f::Integer, N::Integer)

Jordan–Wigner image of the fermionic creation operator $\hat{a}_f^\dagger$ acting on N spin orbitals:

a†_f = ½(X_f - i Y_f) ⊗ Z_{f-1} ⋯ Z_1

Arguments

• f: orbital index in 1..N for the creation operator.
• N: total number of spin orbitals.

Returns a PauliSum{N, ComplexF64}.

largest(ps::PauliSum{N,T}) where {N,T}

Return the term with the largest absolute coefficient as a single-term PauliSum.

largest(v::SparsePauliVector)

Return the term with the largest absolute coefficient as a single-term SparsePauliVector (parity with largest(::PauliSum)).

largest_diag(ps::PauliSum{N,T}) where {N,T}

Return the PauliBasis => coefficient pair for the diagonal term (x == 0) with the largest absolute coefficient.

largest_diag(v::SparsePauliVector)

Return the PauliBasis => coefficient pair for the diagonal term (x == 0) with the largest absolute coefficient (parity with largest_diag(::PauliSum)).

majorana_weight(p::Union{PauliBasis{N}, Pauli{N}}) where N

Compute the Majorana weight of a Pauli string. The Majorana weight counts the number of Majorana operators needed to represent the Pauli string in the Jordan-Wigner encoding.

Algorithm

Uses a branchless O(1) bitwise algorithm. The original per-bit scan tracks a control flag that starts true (at the MSB) and flips at each X/Y site. This is equivalent to a suffix parity of the X-bitstring, computed via parallel prefix XOR (7 steps for Int128). The control mask then selects which Z-only and I-only positions contribute to the weight.

majorana_weight_clip!(ps::PauliSum{N}, max_weight::Int) where {N}

Remove terms with Majorana weight above max_weight.

matrix_element(b::Bra{N}, p::PauliBasis{N}, k::Ket{N})

⟨b|P|k⟩: nonzero only when b.v == k.v ⊻ p.x, i.e. each Pauli connects a ket to exactly one bra. Also accepts AnyPauliSum and KetSum arguments.

merge_pending!(v::SparsePauliVector, f::MergeFilter=NOFILTER)

Sort-merge the pending appends into the sorted live buffer under filter f, restoring the public merged-state invariant (v.an == 0). Allocation-free in steady state.

merge_pending!(v, f, gz, gx)

merge_pending! for appends produced by a single rotation under the generator mask (gz, gx): the gathered triples are then sorted with respect to key ⊻ mask, so _unshuffle_ws! restores natural order in weight(G) linear passes instead of a comparison sort.

offdiag(ps::PauliSum{N,T}) where {N,T}

Return a new PauliSum containing only the off-diagonal terms (those with x != 0).

osum(p1::Pauli{N}, p2::Pauli{M}) where {N,M}

Returns the direct sum of two Paulis

osum(p1::PauliSum{N,T}, p2::PauliSum{M,T}) where {N,M,T}

Direct sum of two PauliSums: p1 ⊕ p2 = p1 ⊗ I_M + I_N ⊗ p2, returning a PauliSum{N+M, T}.

otimes(p1::Pauli{N}, p2::Pauli{M}) where {N,M}

Tensor product of two Paulis, returning a Pauli{N+M}.

otimes(p1::DyadSum{N,T}, p2::DyadSum{M,T}) where {N,M,T}

Tensor product of two DyadSums, returning a DyadSum{N+M}.

otimes(p1::KetSum{N,T}, p2::KetSum{M,T}) where {N,M,T}

Tensor product of two KetSums, returning a KetSum{N+M}.

otimes(p1::PauliSum{N,T}, p2::PauliSum{M,T}) where {N,M,T}

Tensor product of two PauliSums, returning a PauliSum{N+M}.

pauli_channel!(O::PauliSum{N}, pX, pY, pZ; qubits=nothing) -> O

Apply (in place) the Heisenberg-picture single-qubit Pauli channel E(ρ) = pI·ρ + pX·XρX + pY·YρY + pZ·ZρZ independently to each qubit in qubits (default: all 1:N), where pI = 1 - pX - pY - pZ.

The channel is diagonal in the Pauli basis: each PauliBasis term Q is scaled by ∏_{i ∈ qubits} λ_{Q_i}, with

λ_I = 1
λ_X = 1 - 2(pY + pZ)
λ_Y = 1 - 2(pX + pZ)
λ_Z = 1 - 2(pX + pY)

qubits accepts nothing, an Integer, or any iterable of integer qubit indices (e.g. Vector, Tuple, UnitRange).

Throws ArgumentError if any of pX,pY,pZ is negative or if their sum exceeds 1.

promote_to_sum(x)

Wrap a single basis element (Pauli, PauliBasis, Dyad, DyadBasis, Ket) in its corresponding one-term sum type. Mixed products between single operators and sums promote the single operand first, so only sum × sum methods need bespoke implementations.

qdrift(H::PauliSum{N,T}, dt::Real; n_samples::Int=1, rng::AbstractRNG=Random.default_rng())

Decompose time evolution using the QDrift protocol (Campbell, 2019).

Randomly samples Pauli terms from H with probability proportional to |c_k|, producing a sequence of rotations that approximates exp(-i dt H) in expectation.

Each sampled term Pk is rotated by angle θ = -2 dt λ sign(ck), where λ = Σ|c_k| is the 1-norm of the coefficients.

Returns (generators::Vector{PauliBasis{N}}, angles::Vector{Float64}).

Arguments

• H::PauliSum{N,T}: Hamiltonian
• dt::Real: time step
• n_samples::Int=1: number of random samples (more = better approximation)
• rng::AbstractRNG: random number generator

scale!(v1::KetSum{N,T}, a::Number) where {N,T}

Scale all coefficients in v1 by a in-place.

stochastic_clip!(ps::PauliSum{N,T}, ε::Real; rng=Random.default_rng())

Unbiased stochastic compression (Russian Roulette) of a PauliSum.

For each term (basis, c):

• If |c| >= ε: keep unchanged
• If |c| < ε: with probability |c|/ε, promote to ε·sign(c); otherwise delete

This is unbiased: E[c̃] = c for every term.

stochastic_clip!(v::SparsePauliVector, ε::Real; rng=Random.default_rng())

Unbiased stochastic compression (Russian Roulette), the SparsePauliVector analogue of stochastic_clip!(::PauliSum, ...). In-place, order-preserving, allocation-free (modulo RNG state).

symplectic_phase(p::Union{Pauli{N}, PauliBasis{N}})

The power of i needed to recover the Hermitian Pauli string from the bare ZX bitstring form: P = i^θs ⋅ (z|x), with θs = (-n_Y) mod 4 where n_Y is the number of Y sites (count_ones(z & x)). Arises because each Y site is stored as ZX = iY.

trotterize(H::PauliSum{N,T}, dt::Real; n_trotter::Int=1, order::Int=1)

Decompose the time evolution operator exp(-i dt H) into a sequence of Pauli rotations using the Lie-Trotter-Suzuki product formula.

Returns (generators::Vector{PauliBasis{N}}, angles::Vector{Float64}) suitable for passing to evolve.

For a Hamiltonian H = Σk ck Pk, first-order Trotter gives: exp(-i dt H) ≈ Πk exp(-i dt ck Pk)

Each factor exp(-i dt ck Pk) = exp(i θ/2 Pk) with θ = -2 dt Re(ck) in the Heisenberg picture convention used by evolve.

Arguments

• H::PauliSum{N,T}: Hamiltonian (should be Hermitian, so coefficients are real)
• dt::Real: time step
• n_trotter::Int=1: number of Trotter steps (dt is divided by n_trotter)
• order::Int=1: Trotter order (1 = first-order, 2 = second-order symmetric)

truncate!(O::PauliSum, strategy::TruncationStrategy,
          corr::CorrectionAccumulator=NoCorrection())

Apply strategy to truncate O in-place. If a CorrectionAccumulator is provided, measure quantities before and after truncation and accumulate the differences.

Users can define new strategies by subtyping TruncationStrategy and implementing _apply!(O, s). New correction types are defined by subtyping CorrectionAccumulator and implementing _measure(O, corr) and _accumulate!(corr, before, after).

variance(O::PauliSum{N}, ψ::Ket{N}) where N

Compute the variance of observable O in state ψ: <O²> - <O>².

weight(p::PauliBasis)

Number of non-identity single-qubit Pauli factors.

weight_clip!(ps::PauliSum{N}, max_weight::Int)

Remove terms with Pauli weight above max_weight.

weight_damped_clip!(ps::PauliSum{N}, alpha::Real, thresh::Real)

Remove terms with |coefficient|·exp(-alpha·weight) <= thresh. At alpha = 0 this reduces to coeff_clip!(ps, thresh).

x_weight(p::PauliBasis)

Number of qubits acted on by an off-diagonal (X or Y) single-qubit factor. Both X and Y set the x-bit, so this counts the off-diagonal support of the Pauli string (Z and I factors do not contribute).

x_weight_clip!(ps::PauliSum{N}, max_weight::Int)

Remove terms with X-weight (number of X/Y factors) above max_weight.

x_weight_damped_clip!(ps::PauliSum{N}, alpha::Real, thresh::Real)

Remove terms with |coefficient|·exp(-alpha·xweight) <= thresh. At alpha = 0 this reduces to `coeffclip!(ps, thresh)`.