diff --git a/CHANGELOG.md b/CHANGELOG.md
index 630597930..f63730503 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,12 @@
 
 # News
 
+## v0.8.7 - 2023-06-22
+
+- Better UX and threading support for `pftrajectories`.
+- `affectedqubits` now more consistently returns tuples instead of vectors.
+- Private `affectedbits` is now implemented.
+- Many operation constructors now throw an error if they are given negative qubit indices.
 ## v0.8.6 - 2023-06-20
 
 - Fixes to Quantikz circuit plotting of empty circuits and `PauliOperator`
diff --git a/Project.toml b/Project.toml
index e49adab05..b6aa42b17 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumClifford"
 uuid = "0525e862-1e90-11e9-3e4d-1b39d7109de1"
 authors = ["Stefan Krastanov <stefan@krastanov.org>"]
-version = "0.8.6"
+version = "0.8.7"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/docs/src/ecc_example_sim.md b/docs/src/ecc_example_sim.md
index 0db8ca546..5a9ea1e50 100644
--- a/docs/src/ecc_example_sim.md
+++ b/docs/src/ecc_example_sim.md
@@ -32,24 +32,15 @@ scirc = naive_syndrome_circuit(code)
 The most straightforward way to start sampling syndromes is to set up a table of Pauli frames.
 
 ```@example 1
-nframes = 4
-dataqubits = code_n(code)
-ancqubits = code_s(code)
-regbits = ancqubits
-frames = PauliFrame(nframes, dataqubits+ancqubits, regbits)
 circuit = [ecirc..., scirc...]
-pftrajectories(frames, circuit) # run the sims
-pfmeasurements(frames)          # extract the measurements
+nframes = 4
+frames = pftrajectories(circuit; trajectories=nframes) # run the sims
+pfmeasurements(frames)                                 # extract the measurements
 ```
 
-To avoid runtime dispatch, it is a good practice to compactify the circuit:
-
-```@example 1
-frames = PauliFrame(nframes, dataqubits+ancqubits, regbits)
-ccircuit = compactify_circuit(circuit)
-pftrajectories(frames, ccircuit)
-pfmeasurements(frames)
-```
+The [`pftrajectories`](@ref) function is multithreaded.
+If you want more low-level control over these Pauli frame simulations, check out the [`PauliFrame`](@ref) structure,
+the other methods of [`pftrajectories`](@ref), and the circuit compactifaction function [`compactify_circuit`](@ref).
 
 If you want to model Pauli errors, use:
 
@@ -64,7 +55,6 @@ fullcircuit = [ecirc..., errors..., scirc...]
 
 And running this noisy simulation:
 ```@example 1
-frames = PauliFrame(nframes, dataqubits+ancqubits, regbits)
-pftrajectories(frames, fullcircuit)
+frames = pftrajectories(fullcircuit; trajectories=nframes)
 pfmeasurements(frames)
 ```
\ No newline at end of file
diff --git a/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl b/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl
index db9ea6d04..d65eb6e47 100644
--- a/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl
+++ b/ext/QuantumCliffordQuantikzExt/QuantumCliffordQuantikzExt.jl
@@ -19,8 +19,8 @@ function Quantikz.QuantikzOp(op::SparseGate)
         return Quantikz.MultiControlU([],[],op.indices) # TODO Permit skipping the string
     end
 end
-Quantikz.QuantikzOp(op::AbstractOperation) = Quantikz.MultiControlU(affectedqubits(op))
-Quantikz.QuantikzOp(op::PauliOperator) = Quantikz.MultiControlU("\\begin{array}{c}$(lstring(op))\\end{array}", affectedqubits(op))
+Quantikz.QuantikzOp(op::AbstractOperation) = Quantikz.MultiControlU(collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
+Quantikz.QuantikzOp(op::PauliOperator) = Quantikz.MultiControlU("\\begin{array}{c}$(lstring(op))\\end{array}", collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
 Quantikz.QuantikzOp(op::sId1) = Quantikz.Id(affectedqubits(op)...)
 function Quantikz.QuantikzOp(op::AbstractSingleQubitOperator)
     T = typeof(op)
@@ -40,37 +40,37 @@ Quantikz.QuantikzOp(op::sXCZ) = Quantikz.CNOT(reverse(affectedqubits(op))...)
 Quantikz.QuantikzOp(op::sCPHASE) = Quantikz.CPHASE(affectedqubits(op)...)
 Quantikz.QuantikzOp(op::sSWAP) = Quantikz.SWAP(affectedqubits(op)...)
 Quantikz.QuantikzOp(op::sZCZ) = Quantikz.CPHASE(affectedqubits(op)...)
-Quantikz.QuantikzOp(op::sXCX) = Quantikz.MultiControl([],[],affectedqubits(op),[])
+Quantikz.QuantikzOp(op::sXCX) = Quantikz.MultiControl([],[],collect(affectedqubits(op)),[]) # TODO make Quantikz work with tuples and remove the collect
 Quantikz.QuantikzOp(op::sZCY) = Quantikz.MultiControlU("Y",affectedqubits(op)...)
 Quantikz.QuantikzOp(op::sYCZ) = Quantikz.MultiControlU("Y",reverse(affectedqubits(op))...)
 function Quantikz.QuantikzOp(op::AbstractTwoQubitOperator)
     T = typeof(op)
     if T in (sXCY,sYCX,sYCY,sYCZ,sZCY)
-        return Quantikz.U(string(T),affectedqubits(op))
+        return Quantikz.U(string(T),collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
     else
-        return Quantikz.U("{U}",affectedqubits(op)...)
+        return Quantikz.U("{U}",collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
     end
 end
-Quantikz.QuantikzOp(op::BellMeasurement) = Quantikz.ParityMeasurement(["\\mathtt{$(string(typeof(o))[3])}" for o in op.measurements], affectedqubits(op))
+Quantikz.QuantikzOp(op::BellMeasurement) = Quantikz.ParityMeasurement(["\\mathtt{$(string(typeof(o))[3])}" for o in op.measurements], collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
 Quantikz.QuantikzOp(op::NoisyBellMeasurement) = Quantikz.QuantikzOp(op.meas)
-Quantikz.QuantikzOp(op::ConditionalGate) = Quantikz.ClassicalDecision(affectedqubits(op),op.controlbit)
-Quantikz.QuantikzOp(op::DecisionGate) = Quantikz.ClassicalDecision(affectedqubits(op),Quantikz.ibegin:Quantikz.iend)
+Quantikz.QuantikzOp(op::ConditionalGate) = Quantikz.ClassicalDecision(affectedqubits(op),op.controlbit) # TODO make Quantikz work with tuples and remove the collect
+Quantikz.QuantikzOp(op::DecisionGate) = Quantikz.ClassicalDecision(affectedqubits(op),Quantikz.ibegin:Quantikz.iend) # TODO make Quantikz work with tuples and remove the collect
 #Quantikz.QuantikzOp(op::DenseGate) = Quantikz.MultiControlU(affectedqubits(op))
-Quantikz.QuantikzOp(op::PauliMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",affectedqubits(op),op.bit)
-Quantikz.QuantikzOp(op::sMX) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{X}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit)
-Quantikz.QuantikzOp(op::sMY) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Y}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit)
-Quantikz.QuantikzOp(op::sMZ) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Z}\\end{array}",affectedqubits(op), iszero(op.bit) ? nothing : op.bit)
+Quantikz.QuantikzOp(op::PauliMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",collect(affectedqubits(op)),op.bit) # TODO make Quantikz work with tuples and remove the collect
+Quantikz.QuantikzOp(op::Union{sMX,sMRX}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{X}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect
+Quantikz.QuantikzOp(op::Union{sMY,sMRY}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Y}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect
+Quantikz.QuantikzOp(op::Union{sMZ,sMRZ}) = Quantikz.Measurement("\\begin{array}{c}\\mathtt{Z}\\end{array}",collect(affectedqubits(op)), iszero(op.bit) ? nothing : op.bit) # TODO make Quantikz work with tuples and remove the collect
 #Quantikz.QuantikzOp(op::SparseMeasurement) = Quantikz.Measurement("\\begin{array}{c}$(lstring(op.pauli))\\end{array}",affectedqubits(op),op.bit)
 Quantikz.QuantikzOp(op::NoisyGate) = Quantikz.QuantikzOp(op.gate)
-Quantikz.QuantikzOp(op::VerifyOp) = Quantikz.MultiControlU("\\begin{array}{c}\\mathrm{Verify:}\\\\$(lstring(op.good_state))\\end{array}",affectedqubits(op))
+Quantikz.QuantikzOp(op::VerifyOp) = Quantikz.MultiControlU("\\begin{array}{c}\\mathrm{Verify:}\\\\$(lstring(op.good_state))\\end{array}",collect(affectedqubits(op))) # TODO make Quantikz work with tuples and remove the collect
 function Quantikz.QuantikzOp(op::Reset) # TODO This is complicated because quantikz requires $$ for some operators but not all of them... Fix in Quantikz.jl
     m,M = extrema(op.indices)
     indices = sort(op.indices)
     str = "\\begin{array}{c}\\\\$(lstring(op.resetto))\\end{array}"
     if collect(m:M)==indices
-        Quantikz.Initialize("\$$str\$",affectedqubits(op))
+        Quantikz.Initialize("\$$str\$",affectedqubits(op)) # TODO make Quantikz work with tuples and remove the collect
     else
-        Quantikz.Initialize("$str",affectedqubits(op))
+        Quantikz.Initialize("$str",affectedqubits(op)) # TODO make Quantikz work with tuples and remove the collect
     end
 end
 Quantikz.QuantikzOp(op::NoiseOp) = Quantikz.Noise(op.indices)
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 32831bb24..efbc87a14 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -27,7 +27,7 @@ export
     stabilizerview, destabilizerview, logicalxview, logicalzview, phases,
     bitview, quantumstate, tab,
     BadDataStructure,
-    affectedqubits,
+    affectedqubits, #TODO move to QuantumInterface?
     # GF2
     stab_to_gf2, gf2_gausselim!, gf2_isinvertible, gf2_invert, gf2_H_to_G,
     # Canonicalization
@@ -91,6 +91,12 @@ function __init__()
     BIG_INT_FOUR[] = BigInt(4)
 end
 
+const MINBATCH1Q = 100
+const MINBATCH2Q = 100
+const MINBATCHDENSE = 25
+
+const NoZeroQubit = ArgumentError("Qubit indices have to be larger than zero, but you attempting are creating a gate acting on a qubit with a non-positive index. Ensure indexing always starts from 1.")
+
 # Predefined constants representing the permitted phases encoded
 # in the low bits of UInt8.
 const _p  = 0x00
diff --git a/src/affectedqubits.jl b/src/affectedqubits.jl
index 712a7f712..f27bc2469 100644
--- a/src/affectedqubits.jl
+++ b/src/affectedqubits.jl
@@ -1,12 +1,16 @@
 """A method giving the qubits acted upon by a given operation. Part of the Noise interface."""
 function affectedqubits end
-affectedqubits(g::AbstractSingleQubitOperator) = [g.q,]
-affectedqubits(g::AbstractTwoQubitOperator) = [g.q1, g.q2]
+affectedqubits(g::AbstractSingleQubitOperator) = (g.q)
+affectedqubits(g::AbstractTwoQubitOperator) = (g.q1, g.q2)
 affectedqubits(g::NoisyGate) = affectedqubits(g.gate)
 affectedqubits(g::SparseGate) = g.indices
-affectedqubits(b::BellMeasurement) = [m.qubit for m in b.measurements]
+affectedqubits(b::BellMeasurement) = map(m->m.qubit, b.measurements)
 affectedqubits(r::Reset) = r.indices
 affectedqubits(n::NoiseOp) = n.indices
 affectedqubits(g::PauliMeasurement) = 1:length(g.pauli)
 affectedqubits(p::PauliOperator) = 1:length(p)
-affectedqubits(m::AbstractMeasurement) = [m.qubit]
+affectedqubits(m::Union{AbstractMeasurement,sMRX,sMRY,sMRZ}) = (m.qubit,)
+
+affectedbits(o) = ()
+affectedbits(m::sMRZ) = (m.bit,)
+affectedbits(m::sMZ) = (m.bit,)
diff --git a/src/dense_cliffords.jl b/src/dense_cliffords.jl
index 7c810fcb5..a67bd2013 100644
--- a/src/dense_cliffords.jl
+++ b/src/dense_cliffords.jl
@@ -127,7 +127,7 @@ function _apply!(stab::AbstractStabilizer, c::CliffordOperator; phases::Val{B}=V
     nqubits(stab)==nqubits(c) || throw(DimensionMismatch("The tableau and the Clifford operator need to act on the same number of qubits. Consider specifying an array of indices as a third argument to the `apply!` function to avoid this error."))
     s_tab = tab(stab)
     c_tab = tab(c)
-    @batch minbatch=25 threadlocal=zero(c_tab[1]) for row_stab in eachindex(s_tab)
+    @batch minbatch=MINBATCHDENSE threadlocal=zero(c_tab[1]) for row_stab in eachindex(s_tab)
         zero!(threadlocal) # a new stabrow for temporary storage
         apply_row_kernel!(threadlocal, row_stab, s_tab, c_tab, phases=phases)
     end
diff --git a/src/experimental/NoisyCircuits.jl b/src/experimental/NoisyCircuits.jl
index 3cee5d18d..8dc61c64f 100644
--- a/src/experimental/NoisyCircuits.jl
+++ b/src/experimental/NoisyCircuits.jl
@@ -64,7 +64,7 @@ end
 """Compute all possible new states after the application of the given noise model. Reports the probability of each one of them. Deterministic, part of the Noise interface."""
 function applynoise_branches end
 
-function applynoise_branches(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,indices::AbstractVector{Int}; max_order=1)
+function applynoise_branches(s::AbstractStabilizer,noise::UnbiasedUncorrelatedNoise,indices; max_order=1)
     n = nqubits(s)
     l = length(indices)
     infid = noise.errprobthird
diff --git a/src/mctrajectory.jl b/src/mctrajectory.jl
index bc3ff6157..89b474669 100644
--- a/src/mctrajectory.jl
+++ b/src/mctrajectory.jl
@@ -56,7 +56,9 @@ function countmap(samples::Vector{CircuitStatus}) # A simpler faster version of
     Dict(CircuitStatus(i-1)=>counts[i] for i in eachindex(counts))
 end
 
-"""Run multiple Monte Carlo trajectories and report the aggregate final statuses of each."""
+"""Run multiple Monte Carlo trajectories and report the aggregate final statuses of each.
+
+See also: [`pftrajectories`](@ref), [`petrajectories`](@ref)`]"""
 function mctrajectories(initialstate,circuit;trajectories=500)
     counts = countmap([mctrajectory!(copy(initialstate),circuit)[2] for i in 1:trajectories]) # TODO use threads or at least a generator, but without breaking Polyester
     return counts
diff --git a/src/pauli_frames.jl b/src/pauli_frames.jl
index b41696a66..fa6c96f7d 100644
--- a/src/pauli_frames.jl
+++ b/src/pauli_frames.jl
@@ -6,15 +6,16 @@ although it does conjugate the same under Clifford operations.
 Each row in the tableau refers to a single frame.
 The row represents the Pauli operation by which the frame and the reference differ.
 """
-struct PauliFrame{T} <: AbstractQCState
-    frame::T # TODO this should really be a Tableau
-    measurements::Matrix{Bool}
+struct PauliFrame{T,S} <: AbstractQCState
+    frame::T # TODO this should really be a Tableau, but now most of the code seems to be assuming a Stabilizer
+    measurements::S # TODO check if when looping over this we are actually looping over the fast axis
 end
 
 nqubits(f::PauliFrame) = nqubits(f.frame)
 Base.length(f::PauliFrame) = size(f.measurements, 1)
 Base.eachindex(f::PauliFrame) = 1:length(f)
 Base.copy(f::PauliFrame) = PauliFrame(copy(f.frame), copy(f.measurements))
+Base.view(frame::PauliFrame, r) = PauliFrame(view(frame.frame, r), view(frame.measurements, r, :))
 
 """
 $(TYPEDSIGNATURES)
@@ -124,6 +125,50 @@ function pftrajectories end
 """
 $(TYPEDSIGNATURES)
 
+The main method for running Pauli frame simulations of circuits.
+See the other methods for lower level access.
+
+Multithreading is enabled by default, but can be disabled by setting `threads=false`.
+Do not forget to launch Julia with multiple threads enabled, e.g. `julia -t4`, if you want
+to use multithreading.
+
+Note for advanced users: Much of the underlying QuantumClifford.jl functionaly is capable of
+using Polyester.jl threads, but they are fully dissabled here as this is an embarassingly
+parallel problem. If you want to use Polyester.jl threads, use the lower level methods.
+The `threads` keyword argument controls whether standard Julia threads are used.
+"""
+function pftrajectories(circuit;trajectories=5000,threads=true)
+    Polyester.disable_polyester_threads() do
+        _pftrajectories(circuit;trajectories,threads)
+    end
+end
+
+function _pftrajectories(circuit;trajectories=5000,threads=true)
+    ccircuit = if eltype(circuit) <: CompactifiedGate
+        circuit
+    else
+        compactify_circuit(circuit)
+    end
+    qmax=maximum((maximum(affectedqubits(g)) for g in ccircuit))
+    bmax=maximum((maximum(affectedbits(g),init=1) for g in ccircuit))
+    frames = PauliFrame(trajectories, qmax, bmax)
+    nthr = min(Threads.nthreads(),trajectories÷(MINBATCH1Q))
+    if threads && nthr>1
+        batchsize = trajectories÷nthr
+        Threads.@threads for i in 1:nthr
+            b = (i-1)*batchsize+1
+            e = i==nthr ? trajectories : i*batchsize
+            pftrajectories((@view frames[b:e]), ccircuit)
+        end
+    else
+        pftrajectories(frames, ccircuit)
+    end
+    return frames
+end
+
+"""
+$(TYPEDSIGNATURES)
+
 Evolve each frame stored in [`PauliFrame`](@ref) by the given circuit.
 """
 function pftrajectories(state::PauliFrame, circuit)
diff --git a/src/sumtypes.jl b/src/sumtypes.jl
index 2f6167196..a3824c9d6 100644
--- a/src/sumtypes.jl
+++ b/src/sumtypes.jl
@@ -114,6 +114,8 @@ function make_all_sumtype_infrastructure()
             (:applywstatus!, (:(s::Register),), ()),
             (:apply!, (:(s::PauliFrame),), ()),
             (:applywstatus!, (:(s::PauliFrame),), ()),
+            (:affectedqubits, (), ()),
+            (:affectedbits, (), ()),
         ]
     ) |> eval
 end
diff --git a/src/symbolic_cliffords.jl b/src/symbolic_cliffords.jl
index 8972d7e0a..8f84fd5a5 100644
--- a/src/symbolic_cliffords.jl
+++ b/src/symbolic_cliffords.jl
@@ -9,9 +9,6 @@ abstract type AbstractTwoQubitOperator <: AbstractSymbolicOperator end
 """Supertype of all symbolic single-qubit measurements."""
 abstract type AbstractMeasurement <: AbstractOperation end
 
-const MINBATCH1Q = 100
-const MINBATCH2Q = 100
-
 # Stim has a good list of specialized single and two qubit operations at https://github.com/quantumlib/Stim/blob/e51ea66d213b25920e72c08e53266ec56fd14db4/src/stim/stabilizers/tableau_specialized_prepend.cc
 # Note that their specialized operations are for prepends (right multiplications), while we implement append (left multiplication) operations.
 
@@ -67,6 +64,7 @@ macro qubitop1(name, kernel)
     quote
         struct $(esc(prefixname)) <: AbstractSingleQubitOperator
             q::Int
+            $(esc(prefixname))(q) = if q<=0 throw(NoZeroQubit) else new(q) end
         end
         @doc $docstring $prefixname
         @inline $(esc(:qubit_kernel))(::$prefixname, x, z) = $kernel
@@ -132,6 +130,7 @@ struct SingleQubitOperator <: AbstractSingleQubitOperator
     zz::Bool
     px::Bool
     pz::Bool
+    SingleQubitOperator(q,args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 function _apply!(stab::AbstractStabilizer, op::SingleQubitOperator; phases::Val{B}=Val(true)) where B # TODO Generated functions that simplify the whole `if phases` branch might be a good optimization, but a quick benchmakr comparing sHadamard to SingleQubitOperator(sHadamard) did not show a worthwhile difference.
     s = tab(stab)
@@ -249,6 +248,7 @@ macro qubitop2(name, kernel)
         struct $(esc(prefixname)) <: AbstractTwoQubitOperator
             q1::Int
             q2::Int
+            $(esc(prefixname))(q1,q2) = if q1<=0 || q2<=0 throw(NoZeroQubit) else new(q1,q2) end
         end
         @doc $docstring $prefixname
         @inline $(esc(:qubit_kernel))(::$prefixname, x1, z1, x2, z2) = $kernel
@@ -363,18 +363,21 @@ end
 struct sMX <: AbstractMeasurement
     qubit::Int
     bit::Int
+    sMX(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 """Symbolic single qubit Y measurement. See also [`Register`](@ref), [`projectYrand!`](@ref), [`sMX`](@ref), [`sMZ`](@ref)"""
 struct sMY <: AbstractMeasurement
     qubit::Int
     bit::Int
+    sMY(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 """Symbolic single qubit Z measurement. See also [`Register`](@ref), [`projectZrand!`](@ref), [`sMX`](@ref), [`sMY`](@ref)"""
 struct sMZ <: AbstractMeasurement
     qubit::Int
     bit::Int
+    sMZ(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 sMX(i) = sMX(i,0)
@@ -471,6 +474,7 @@ See also: [`Reset`](@ref), [`sMZ`](@ref)"""
 struct sMRZ <: AbstractOperation
     qubit::Int
     bit::Int
+    sMRZ(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 """Measure a qubit in the X basis and reset to the |+⟩ state.
@@ -479,6 +483,7 @@ See also: [`sMRZ`](@ref), [`Reset`](@ref), [`sMZ`](@ref)"""
 struct sMRX <: AbstractOperation
     qubit::Int
     bit::Int
+    sMRX(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 """Measure a qubit in the Y basis and reset to the |i₊⟩ state.
@@ -487,6 +492,7 @@ See also: [`sMRZ`](@ref), [`Reset`](@ref), [`sMZ`](@ref)"""
 struct sMRY <: AbstractOperation
     qubit::Int
     bit::Int
+    sMRY(q, args...) = if q<=0 throw(NoZeroQubit) else new(q,args...) end
 end
 
 sMRX(i) = sMRX(i,0)
diff --git a/test/test_pauliframe.jl b/test/test_pauliframe.jl
index 14d0ccd8d..0669f2ac2 100644
--- a/test/test_pauliframe.jl
+++ b/test/test_pauliframe.jl
@@ -10,11 +10,15 @@ using Test
     ]
     frame = PauliFrame(100, 5, 2)
     frame = pftrajectories(frame, circuit)
-    m = pfmeasurements(frame)
+    frame1 = pftrajectories(circuit; trajectories=100, threads=false)
+    frame2 = pftrajectories(circuit; trajectories=100, threads=true)
     # If the x component is set on the second qubit, then the fourth and fifth qubits should also have it set
-    for (i, row) in enumerate(frame.frame)
-        if row[2] == (true, false)
-            @test row[4][1] && row[5][1] && m[i,1] && m[i,2]
+    for f in [frame, frame1, frame2]
+        m = pfmeasurements(f)
+        for (i, row) in enumerate(f.frame)
+            if row[2] == (true, false)
+                @test row[4][1] && row[5][1] && m[i,1] && m[i,2]
+            end
         end
     end
 end
@@ -24,44 +28,65 @@ end
         sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state
         sMZ(1,1), sMZ(2,2), sMZ(3,3) # measure each qubit
     ]
-    frame = PauliFrame(10^6, 3, 3)
+    n = 10^6
+    frame = PauliFrame(n, 3, 3)
     f = pftrajectories(frame, ghz_circuit)
     m = pfmeasurements(f)
-    rowtotal_1s = sum(m, dims=2)[:,1]
-    @test all(rowtotal_1s .% 3 .== 0)
-    fractotal_1s = sum(rowtotal_1s)/3 / 10^6
-    @test (fractotal_1s > 0.49) && (fractotal_1s < 0.51)
+    frame1 = pftrajectories(ghz_circuit; trajectories=n, threads=false)
+    m1 = pfmeasurements(frame1)
+    frame2 = pftrajectories(ghz_circuit; trajectories=n, threads=true)
+    m2 = pfmeasurements(frame2)
+    for _m in [m, m1, m2]
+        rowtotal_1s = sum(m, dims=2)[:,1]
+        @test all(rowtotal_1s .% 3 .== 0)
+        fractotal_1s = sum(rowtotal_1s)/3 / n
+        @test (fractotal_1s > 0.49) && (fractotal_1s < 0.51)
+    end
 end
 
-@testset "sMZ vs sMZR - mctrajectories vs pftrajectories" begin
+@testset "sMZ vs sMRZ - mctrajectories vs pftrajectories" begin
+    n = 2000
     state = Register(one(MixedDestabilizer, 3), 6)
-    frame = PauliFrame(100, 3, 6)
+    frame = PauliFrame(n, 3, 6)
 
     ghz_circuit1 = [
         sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state
         sMZ(1,1), sMZ(2,2), sMZ(3,3), # measure each qubit
         sMZ(1,4), sMZ(2,5), sMZ(3,6)  # measure each qubit again
     ]
-    ms1 = stack([bitview(mctrajectory!(copy(state), ghz_circuit1)[1]) for i in 1:100], dims=1)
-    mf1 = pfmeasurements(pftrajectories(copy(frame), ghz_circuit1))
-    @test all(0.25 .< sum(ms1, dims=1)./100 .< 0.75)
-    @test all(0.25 .< sum(mf1, dims=1)./100 .< 0.75)
+    for m in [
+        stack([bitview(mctrajectory!(copy(state), ghz_circuit1)[1]) for i in 1:n], dims=1),
+        pfmeasurements(pftrajectories(copy(frame), ghz_circuit1)),
+        pfmeasurements(pftrajectories(ghz_circuit1;trajectories=n,threads=false)),
+        pfmeasurements(pftrajectories(ghz_circuit1;trajectories=n,threads=true)),
+    ]
+        @test all(0.25 .< sum(m, dims=1)./n .< 0.75)
+    end
 
     ghz_circuit2 = [
         sHadamard(1), sCNOT(1,2), sCNOT(1,3), # prepare a GHZ state
         sMRZ(1,1), sMZ(2,2), sMZ(3,3), # measure and reset each qubit
         sMZ(1,4), sMZ(2,5), sMZ(3,6)  # measure each qubit again
     ]
-    ms2 = stack([bitview(mctrajectory!(copy(state), ghz_circuit2)[1]) for i in 1:100], dims=1)
-    mf2 = pfmeasurements(pftrajectories(copy(frame), ghz_circuit2))
-    @test all(0.25.*[1 1 1 0 1 1] .<= sum(ms2, dims=1)./100 .<= 0.75.*[1 1 1 0 1 1])
-    @test all(0.25.*[1 1 1 0 1 1] .<= sum(mf2, dims=1)./100 .<= 0.75.*[1 1 1 0 1 1])
+    for m in [
+        stack([bitview(mctrajectory!(copy(state), ghz_circuit2)[1]) for i in 1:n], dims=1),
+        pfmeasurements(pftrajectories(copy(frame), ghz_circuit2)),
+        pfmeasurements(pftrajectories(ghz_circuit2;trajectories=n,threads=false)),
+        pfmeasurements(pftrajectories(ghz_circuit2;trajectories=n,threads=true)),
+    ]
+        @test all(0.25.*[1 1 1 0 1 1] .<= sum(m, dims=1)./n .<= 0.75.*[1 1 1 0 1 1])
+    end
 
     noncom_circuit = [
         sHadamard(1), sMRZ(1,1), sX(1), sMZ(1,2), sMRZ(1,3), sMRZ(1,4), sHadamard(1), sMZ(1,5)
     ]
-    ms3 = stack([bitview(mctrajectory!(copy(state), noncom_circuit)[1]) for i in 1:100], dims=1)
-    mf3 = pfmeasurements(pftrajectories(copy(frame), noncom_circuit))
-    @test all(0.25.*[1 4 4 0 1 0] .<= sum(ms3, dims=1)./100 .<= 0.75.*[1 2 2 0 1 0])
-    @test all(0.25.*[1 0 0 0 1 0] .<= sum(mf3, dims=1)./100 .<= 0.75.*[1 0 0 0 1 0])
+    ms3 = stack([bitview(mctrajectory!(copy(state), noncom_circuit)[1]) for i in 1:n], dims=1)
+    @test all(0.25.*[1 4 4 0 1 0] .<= sum(ms3, dims=1)./n .<= 0.75.*[1 2 2 0 1 0])
+    for m in [
+        pfmeasurements(pftrajectories(copy(frame), noncom_circuit)),
+        pfmeasurements(pftrajectories(noncom_circuit;trajectories=n,threads=false)),
+        pfmeasurements(pftrajectories(noncom_circuit;trajectories=n,threads=true)),
+    ]
+        @test all(0.25.*[1 0 0 0 1] .<= (sum(m, dims=1)[:,1:5])./n .<= 0.75.*[1 0 0 0 1])
+    end
 end