PRISM-Q runs quantum circuits fast by matching each one to the right simulation strategy. It dispatches across nine CPU backends plus optional CUDA and MPI paths, optimizes circuits through a multi-pass fusion pipeline, and uses AVX2, FMA, and BMI2 SIMD in the inner loop. Input is OpenQASM 3.0, with backward-compatible 2.0 syntax. The same library handles a two-qubit Bell pair and a thousand-qubit Clifford circuit.

#![allow(unused)]
fn main() {
use prism_q::CircuitBuilder;

let result = CircuitBuilder::new(2).h(0).cx(0, 1).run(42).unwrap();
let probs = result.probabilities.unwrap();
// |00> = 0.5, |11> = 0.5
}

What it does

• Eight CPU backends selected automatically per circuit: statevector, stabilizer (with factored and filtered variants), sparse, MPS, product state, tensor network, and dynamic factored split-state.
• Compiled shot samplers that sample without rebuilding the full statevector each shot, including noisy and detector/QEC paths.
• Clifford+T strategies (stabilizer rank, stochastic and deterministic Pauli propagation) for circuits beyond the reach of a dense statevector.
• Optional CUDA backend for statevector and stabilizer execution.

Where to go next

• Installation and Your First Circuit: a hands-on start.
• Architecture: backends, dispatch tree, fusion pipeline, and SIMD strategy.
• Glossary: terminology used across the documentation.
• API reference: generated Rust API documentation.
• Source and issues: the repository on GitHub.

Installation

PRISM-Q is a Rust library with optional Python bindings. For Python, see Python Bindings:

pip install prism-q

Add the Rust crate to a project with Cargo:

cargo add prism-q                          # Rayon parallelism + faer SVD (default)
cargo add prism-q --no-default-features    # single-threaded, minimal dependencies

cargo add writes the current release into Cargo.toml for you. To pin a version by hand instead, take it from crates.io.

Feature flags

Building from source

git clone https://github.com/AbeCoull/prism-q
cd prism-q
cargo build --release

Running the test suite

cargo nextest run --all-features                          # unit + integration tests
cargo test --doc --all-features                           # doctests
cargo clippy --all-targets --all-features -- -D warnings  # lint

Use cargo test --all-features if cargo-nextest is unavailable.

--all-features includes distributed-mpi, which needs a system MPI installation and libclang for its bindgen step. Without those, substitute the features you actually want, for example --features "parallel gpu".

Next: build Your First Circuit.

Your First Circuit

There are two ways to build a circuit: the fluent CircuitBuilder API, or by parsing OpenQASM text.

With the builder

CircuitBuilder chains gate calls and runs the result. This builds a Bell pair, the two-qubit entangled state (|00⟩ + |11⟩) / √2:

#![allow(unused)]
fn main() {
use prism_q::CircuitBuilder;

let result = CircuitBuilder::new(2)
    .h(0)
    .cx(0, 1)
    .run(42)                       // seed = 42
    .expect("simulation failed");

let probs = result.probabilities.expect("no probabilities");
for i in 0..probs.len() {
    let p = probs.get(i);
    if p > 1e-10 {
        println!("|{i:02b}> = {p:.4}");
    }
}
// |00> = 0.5000
// |11> = 0.5000
}

Measurement is Z-basis by default. measure_in_basis(qubit, axis, bit) measures along X, Y, or Z, and measure_pauli_product(&terms, bit) records the parity of a Pauli string such as X0 Z1. Both lower onto the gates the backends already run, so every backend accepts them. A basis measurement leaves the qubit in the Z eigenstate of the recorded bit rather than rotating back, and a Pauli product accumulates its parity on one extra qubit appended past the declared register:

#![allow(unused)]
fn main() {
use prism_q::{CircuitBuilder, PauliAxis, PauliTerm};

let circuit = CircuitBuilder::new_with_classical(2, 2)
    .h(0)
    .cx(0, 1)
    .measure_in_basis(0, PauliAxis::X, 0)
    .measure_pauli_product(&[PauliTerm::z(0), PauliTerm::z(1)], 1)
    .build();
assert_eq!(circuit.num_qubits, 3);
}

A larger structurally similar circuit, the 5-qubit GHZ state, renders like this (diagram generated by PRISM-Q's own SVG renderer):

From OpenQASM

The same Bell pair, written in OpenQASM 3.0 and parsed:

#![allow(unused)]
fn main() {
use prism_q::circuit::openqasm;
use prism_q::simulate;

let qasm = r#"
    OPENQASM 3.0;
    include "stdgates.inc";
    qubit[2] q;
    h q[0];
    cx q[0], q[1];
"#;

let circuit = openqasm::parse(qasm).expect("failed to parse QASM");
let result = simulate(&circuit).seed(42).run().expect("simulation failed");
}

run_qasm(qasm, seed) is a shortcut that parses and simulates in one call. See the OpenQASM Support guide for the supported subset.

Next: sample measurement outcomes in Shots and Sampling.

Shots and Sampling

Probabilities give you the exact distribution. Shots give you sampled measurement outcomes, the way real hardware reports results. PRISM-Q samples deterministically from a fixed seed.

Sampling shots

#![allow(unused)]
fn main() {
use prism_q::circuit::openqasm;
use prism_q::simulate;

let qasm = r#"
    OPENQASM 3.0;
    include "stdgates.inc";
    qubit[2] q;
    bit[2] c;
    h q[0];
    cx q[0], q[1];
    c[0] = measure q[0];
    c[1] = measure q[1];
"#;
let circuit = openqasm::parse(qasm).expect("failed to parse QASM");

let result = simulate(&circuit).seed(42).shots(1024).expect("shots failed");
print!("{result}");   // ShotsResult implements Display
}

The same seed always produces the same samples. Pass rand::random() as the seed for non-deterministic sampling.

Counts and marginals

For large shot counts, you usually want aggregates rather than raw shots:

#![allow(unused)]
fn main() {
// Frequency histogram: bitstring -> count
let counts = simulate(&circuit).seed(42).sample_counts(100_000).unwrap();

// Per-qubit P(measuring |1>), without the full joint distribution
let marginals = simulate(&circuit).seed(42).marginals().unwrap();
}

Noisy sampling

Attach a NoiseModel to sample under depolarizing or readout noise:

#![allow(unused)]
fn main() {
use prism_q::{simulate, BackendKind, NoiseModel};

let noise = NoiseModel::uniform_depolarizing(&circuit, 0.001);
let result = simulate(&circuit)
    .backend(BackendKind::Statevector)
    .noise(noise)
    .seed(42)
    .shots(1024)
    .unwrap();
}

The Noise and QEC guide covers noise models and detector sampling in depth.

Next: learn how PRISM-Q picks a representation in Choosing a Backend.

Choosing a Backend

By default PRISM-Q inspects the circuit and picks a backend for you. You only need to choose explicitly when you know something the auto-dispatcher cannot infer, or when you are benchmarking a specific representation.

Let it choose

#![allow(unused)]
fn main() {
use prism_q::simulate;

let result = simulate(&circuit).seed(42).run().unwrap();   // BackendKind::Auto
}

Auto-dispatch walks this decision tree:

flowchart TD
    A[Auto] --> E{Entangling gates?}
    E -- none --> PS[ProductState]
    E -- yes --> CL{All Clifford?}
    CL -- yes --> STB[Stabilizer]
    CL -- no --> MEM{Above memory limit?}
    MEM -- "yes, sparse-friendly" --> SPR[Sparse]
    MEM -- "yes, otherwise" --> MPS[MPS bond 256]
    MEM -- no --> IND{Partial independence?}
    IND -- yes --> FAC[Factored]
    IND -- no --> SV[Statevector]

Choose explicitly

#![allow(unused)]
fn main() {
use prism_q::{simulate, BackendKind};

let result = simulate(&circuit)
    .backend(BackendKind::Stabilizer)
    .seed(42)
    .run()
    .unwrap();
}

Symptom to backend

The Backends Deep Dive and the architecture reference cover each backend's internals.

Python Bindings

prism-q ships Python bindings built with PyO3. They are a thin wrapper over the Rust crate: the compiled extension is prism_q._prism_q and the pure-Python prism_q package re-exports it. Simulation runs in Rust with the GIL released, so the wrapper adds no per-gate overhead.

Wheels are abi3 for Python 3.11 and newer, so one wheel per platform covers every supported interpreter.

Install

pip install prism-q

NumPy is the only runtime dependency. Building from a source checkout needs maturin:

pip install maturin
maturin develop --manifest-path bindings/python/Cargo.toml

The bindings enable the parallel feature by default. The gpu feature is optional and off in the published wheels (see GPU backends); the distributed backend is not reachable from Python.

Quick start

from prism_q import CircuitBuilder, simulate

circuit = CircuitBuilder(2, 2).h(0).cx(0, 1).measure_all().build()
counts = simulate(circuit).seed(42).shots(1000).counts()
print(counts)          # {'00': 507, '11': 493}

Building circuits

CircuitBuilder is a fluent API. Every gate method returns the builder, and build() produces the Circuit that simulation consumes.

from prism_q import CircuitBuilder

circuit = (
    CircuitBuilder(3, 3)
    .h(0)
    .cx(0, 1)
    .rz(0.5, 2)
    .cphase(0.25, 1, 2)
    .measure_all()
    .build()
)

pauli_rotation(theta, factors) appends exp(-i * theta * P / 2) for the Pauli string P given as (qubit, axis) factors with axis one of "X", "Y", "Z"; identity factors are omitted. A weight-1 string lowers to rx, ry, or rz and a two-qubit ZZ string to rzz, so fusion and Clifford recognition keep firing on them. Circuit.add_pauli_rotation is the imperative spelling.

builder.pauli_rotation(0.4, [(0, "X"), (1, "Y"), (3, "Z")]).param(0)

cu and mcu take a 2x2 matrix as nested Python sequences of complex numbers. Out-of-range qubits raise PrismError at build time rather than at simulation time.

Three other routes produce a Circuit:

from prism_q import Circuit, circuits, parse_qasm

manual = Circuit(2, 2)                    # imperative, add_gate / add_measure / add_reset
ghz = circuits.ghz(10)                    # pre-built corpus
parsed = parse_qasm(qasm_source)          # OpenQASM 3.0, with 2.0 accepted

The circuits submodule mirrors the Rust builders documented in Circuit Builders: qft, ghz, w_state, random, hardware_efficient_ansatz, clifford_heavy, clifford_random_pairs, qaoa, single_qubit_rotation, clifford_t, quantum_volume, cz_chain, phase_estimation, independent_bell_pairs, independent_random_blocks, and local_clifford_blocks. Seeded builders default to seed 42.

Running a simulation

simulate(circuit) returns a Simulation you configure with .seed(), .backend(), and .noise(), then finish with a terminal method. The default seed is 42.

from prism_q import BackendKind, simulate

sim = simulate(circuit).seed(7).backend(BackendKind.statevector())
outcome = sim.run()

shots() and sample_counts() average trajectories on any backend holding a per-shot pure state. The three rows marked "density matrix only" read the exact mixed state instead, so they need .backend(BackendKind.density_matrix()); auto dispatch never selects it. There the mixture is evolved once and every terminal reads that one evolution, so the probabilities are seed independent and the observables carry no sampling error. Circuits with mid-circuit measurement or classical conditioning are rejected on that route, since the mixture holds every measurement branch at once.

Terminals that cannot honor a model raise PrismError naming the reason, rather than silently ignoring it. state_vector() always uses the statevector backend and density_matrix_expectation_values() always uses the density-matrix backend, both regardless of .backend(...).

ShotsResult and CountsResult both expose counts(), returning a dict keyed by bitstring.

Distributions too wide to write down

A circuit whose qubits fall into independent groups is answered per group, and RunOutcome keeps it that way: the dense vector is built only when probabilities is read. Fifteen independent Bell pairs span 30 qubits, whose dense form is 8 GB, and the blocks are 15 arrays of four entries.

outcome = simulate(circuits.independent_bell_pairs(15)).seed(42).run()

outcome.num_basis_states            # 2 ** 30, and nothing was materialized
for qubits, probs in outcome.probabilities_factored():
    print(qubits, probs)            # [0, 1] [0.5 0. 0. 0.5], ...

Each block is (qubits, probs) with qubits ascending, and probs indexed by those qubits packed in that order with qubits[0] in the least significant bit. The probability of a basis state is the product of one entry per block, which is what probabilities computes.

probabilities_factored() returns None when the run produced a dense distribution, which is the common case: the decomposed route needs the widest group several qubits narrower than the register, so two Bell pairs stay dense. num_basis_states is None when the backend exposed no distribution at all.

Result metadata

RunOutcome, ShotsResult, and CountsResult each carry a metadata object describing how the result was produced.

result = simulate(circuit).seed(42).run()
print(result.metadata.backend)               # 'Statevector'
print(result.metadata.engine)                # None unless samplers share the backend
print(result.metadata.is_exact)              # True
print(result.metadata.fidelity_lower_bound)  # None when exact
print(result.metadata.placement)             # 'host' or 'device'

is_exact is False when the engine that ran can discard state weight or estimate by sampling. It marks the route, not the run: an MPS whose bond dimension the circuit never fills reports is_exact == False with fidelity_lower_bound == 1.0, so the flag answers whether the answer could have been approximated and the bound answers whether it was.

Automatic dispatch sends a circuit past the statevector cap to a bounded-bond MPS, which is the only route those circuits have. That is taken by default and the result says so. .require_exact() rejects it instead, raising PrismError naming the engine it would have used.

simulate(big_circuit).seed(42).require_exact().marginals()  # raises

Starting from a state other than |0...0>

.initial_state(amplitudes) replaces the default all-zero start. It takes any sequence of complex numbers, a complex128 NumPy array included, indexed with qubit 0 in the least significant bit.

import math
import numpy as np
from prism_q import CircuitBuilder, simulate

theta = math.pi / 8
start = np.array([math.cos(theta), math.sin(theta)], dtype=np.complex128)
circuit = CircuitBuilder(1).h(0).build()
probs = simulate(circuit).initial_state(start).seed(42).run().probabilities

The vector must have 2 ** num_qubits entries and unit norm. A wrong length, a non-finite entry, or a norm off unity raises PrismError; an unnormalized vector is rejected rather than rescaled, so a mistake surfaces instead of becoming a silent factor on every amplitude.

A start state also narrows the route. Auto dispatch reads circuit structure, and its shortcuts (tableau, product state, subsystem decomposition, Pauli propagation) are only valid from |0...0>: a Clifford circuit produces a stabilizer state only when its input is one. So auto() resolves to the statevector, the GPU and distributed statevectors and density_matrix() are the only other backends that accept one, and every other choice raises PrismError naming itself. run(), shots(), sample_counts(), marginals(), expectation_values(), and state_vector() carry it; expectation_gradient() and density_matrix_expectation_values() reject it, as do shots() and sample_counts() with a noise model attached, since trajectory replay reinitializes a pure state per shot. To evolve a start state under noise, read the exact mixture with run(), marginals(), or expectation_values() on density_matrix().

Selecting a backend

BackendKind.auto() is the default and picks a backend from circuit structure. Pass an explicit one to override it.

The density-matrix backend stores 4^n amplitudes, so its qubit ceiling is about half the statevector cap; exceeding it raises PrismError naming the cap. The distributed statevector backend has no Python constructor: MPI_Init ownership between the interpreter, mpi4py, and the extension is unsettled.

GPU backends

The GPU constructors take a GpuContext, an opaque handle to one CUDA device and its compiled kernels. Build it once and reuse it: construction compiles the kernel module, and passing the same handle to several simulations shares that work.

from prism_q import BackendKind, GpuContext, circuits, simulate

context = GpuContext(0)
outcome = simulate(circuits.qft(16)).backend(BackendKind.auto_gpu(context)).seed(42).run()
print(outcome.probabilities)

GpuContext(device_id) is where a missing or unusable device is reported, and it raises PrismError rather than falling back. Past construction, routing is soft by design and matches the Rust API: statevector_gpu runs circuits below the crossover (PRISM_GPU_MIN_QUBITS, default 14) on the host, auto_gpu routes each block independently, and a block whose device allocation fails degrades to the host rather than erroring. A run that produces host results is therefore normal, not a failure signal.

stabilizer_gpu sets its crossover at 100000 qubits (PRISM_STABILIZER_GPU_MIN_QUBITS), so it runs on the host tableau unless that override is lowered. The device tableau is correct; the default stays high until benchmarks justify lowering it.

density_matrix_gpu(context) holds the exact mixed state on the device and is the one device kind with no crossover and no host fallback: auto_gpu never selects it, every noisy terminal that density_matrix() serves answers from the device buffer, and a width whose 4^n buffer does not fit in free device memory raises PrismError before anything is allocated (13 qubits on an 11 GiB card).

The published wheels are built without CUDA, because two of the three wheel targets have no CUDA toolkit and macOS has no CUDA at all. In those wheels the constructors still exist and GpuContext(...) raises PrismError naming the missing build feature, so code written against the GPU API fails with a message rather than an AttributeError. Two predicates separate the cases:

GpuContext.is_supported()   # was this build compiled with CUDA support
GpuContext.is_available()   # ... and is a usable device present

To get a build with CUDA support, install the CUDA toolkit (12.x or newer) and build from a checkout:

maturin develop --manifest-path bindings/python/Cargo.toml --features gpu

On Windows that build links the toolkit's NVRTC library (nvrtc64_120_0.dll for CUDA 12.x) from the toolkit bin directory, which Python does not search. The package adds it on import when CUDA_PATH is set, which the toolkit installer does; without it the import fails with DLL load failed while importing _prism_q.

Distributed backend

The distributed statevector shards the dense state across MPI ranks. It is reachable from Python through a DistributedContext, which attaches to an MPI that is already running. This extension never calls MPI_Init or MPI_Finalize: mpi4py does both, at its own import and at interpreter exit, and a handle whose refcount drop finalized MPI would make every later MPI call in the process erroneous.

from mpi4py import MPI  # starts MPI; import before touching the context
from prism_q import BackendKind, DistributedContext, circuits, simulate

context = DistributedContext()
outcome = (
    simulate(circuits.qft(24))
    .backend(BackendKind.statevector_distributed(context))
    .seed(42)
    .run()
)
print(context.rank, context.size, outcome.probabilities[:4])

Run it with mpiexec -n 4 python script.py. The rank count must be a power of two, and every rank needs enough local qubits after the shard bits are taken (PRISM_DIST_MIN_LOCAL_QUBITS).

The contract is SPMD, and it is enforced rather than assumed. Four ranks are four interpreters running the same source, and every collective inside the backend is entered by all of them, so a script that branches before the call

if comm.rank == 0:
    result = simulate(circuit).backend(...).run()   # deadlocks

hangs the job: rank 0 blocks in a collective the others never enter. Every rank calls run with the same circuit and seed; only what you do with the returned value may branch on rank. Two cross-checks turn the common violations into errors instead of hangs, each costing one collective at run entry: ranks disagreeing about the register shape, the seed, or the tuning knobs are rejected, and so are ranks handed different circuits.

Two surfaces are deliberately loud rather than convenient.

A world of one rank raises. MPI-2 and later make a singleton MPI_Init succeed, so a script launched without mpiexec would otherwise get a correct answer from one rank at single-host speed with no signal that nothing was distributed. Pass allow_single_rank=True when that is what you meant.

Constructing the context without MPI running raises rather than starting MPI, so a forgotten from mpi4py import MPI is reported at the point it happened.

The thread level mpi4py negotiated must be at least MPI_THREAD_FUNNELED: Rayon workers run beside MPI in the same process, and every MPI call stays on the thread that constructed the context. mpi4py requests MPI_THREAD_MULTIPLE by default; a script that lowers it (mpi4py.rc.thread_level = 'single' or mpi4py.rc.threads = False) gets an error from the constructor rather than a run that mixes threads with a single-threaded MPI.

The published wheels have no MPI support: mpi-sys runs bindgen and needs a system MPI at build time, and the extension has to link the same MPI implementation and ABI as mpi4py and as the launcher. Mixing two implementations in one process corrupts rather than failing loudly, which is why this is a from-source path:

maturin develop --manifest-path bindings/python/Cargo.toml --features distributed-mpi

DistributedContext.is_supported() says whether a build has MPI support, the same way GpuContext.is_supported() does for CUDA.

Sub-communicators, per-rank device placement, and any distributed path other than the statevector backend are out of scope: the context is the world communicator and nothing else.

Noise

Build a NoiseModel from a circuit, then attach it. A model is sized to the circuit it was built from, and using it with a different circuit raises PrismError.

from prism_q import NoiseChannel, NoiseModel, simulate

model = NoiseModel.uniform_depolarizing(circuit, 0.01)
counts = simulate(circuit).seed(42).noise(model).sample_counts(4000).counts()

NoiseModel.uniform_depolarizing(circuit, p) and NoiseModel.with_amplitude_damping(circuit, gamma) cover the common cases. For per-instruction control, start from NoiseModel.empty(circuit) and attach events:

model = NoiseModel.empty(circuit)
model.add_event(0, NoiseChannel.amplitude_damping(0.05), [0])
model.add_event(1, NoiseChannel.two_qubit_depolarizing(0.02), [0, 1])
model.with_readout_error(0.01, 0.01)
model.validate()

Channels are pauli(px, py, pz), depolarizing(p), amplitude_damping(gamma), phase_damping(gamma), thermal_relaxation(t1, t2, gate_time), two_qubit_depolarizing(p), and custom(kraus) for an explicit list of 2x2 Kraus operators. validate() checks probabilities and Kraus completeness; is_pauli_only() reports whether the model holds only single-qubit Pauli channels and no readout error. A model carrying readout error or two_qubit_depolarizing answers False there and still runs on the stabilizer samplers, which apply readout to the measurement record and sample the pair channel as one joint draw over its 15 branches.

Expectation values

An observable is a list of (qubit, axis) factors with axis one of "X", "Y", "Z". Identity factors are omitted, so [(0, "Z"), (2, "X")] means Z0 ⊗ I1 ⊗ X2. Both expectation terminals take a list of observables and return one float each.

observables = [[(0, "Z")], [(0, "Z"), (1, "Z")]]
values = simulate(circuit).seed(42).expectation_values(observables)

expectation_values requires a unitary circuit and gives ⟨ψ|P|ψ⟩. density_matrix_expectation_values evolves the density matrix through the circuit and any attached noise model and gives exact Tr(rho P), with measurements read off the final mixed state without collapse. It is the zero-variance analogue of averaging over trajectories:

model = NoiseModel.empty(circuit)
model.add_event(0, NoiseChannel.amplitude_damping(0.3), [0])
exact = simulate(circuit).seed(42).noise(model).density_matrix_expectation_values(
    [[(0, "Z")]]
)

Parameter sweeps

A variational loop rebinds angles while the gate sequence stays fixed. Parameters names the slots those angles land in, and PreparedCircuit holds the circuit across bindings so fusion and backend selection are settled once rather than per point.

from prism_q import CircuitBuilder, PreparedCircuit

builder = CircuitBuilder(4)
for q in range(4):
    builder.ry(0.1, q).param(q)
for q in range(3):
    builder.cx(q, q + 1)

prepared = PreparedCircuit(builder.build(), builder.parameters())
for values in [[0.1, 0.2, 0.3, 0.4], [0.5, 0.6, 0.7, 0.8]]:
    outcome = prepared.run(values, seed=42)

run(values, seed=42) returns the same RunOutcome as simulate(...).run(). bind(values) returns the bound Circuit instead, for handing to simulate(...) with different options or to any other consumer of a circuit.

Automatic dispatch reads the template, so build it at angles representative of the sweep. A template whose rotations are all zero reads as Clifford and settles on a backend that then rejects the bound circuit. Pass an explicit backend as the third argument to decide it yourself:

prepared = PreparedCircuit(circuit, params, BackendKind.statevector())

reuses_fusion_plan reports whether the recorded fused structure is being replayed. It is a performance fact, not an error: results agree either way.

Build a Parameters three ways. builder.parameters() returns the set recorded by param(slot); Parameters.all_rotations(circuit) gives every bindable gate its own slot in circuit order; Parameters(n) plus link(instruction, slot) declares the slots up front. Several gates may share a slot, in which case binding writes one angle to each.

A wrong-length value vector, a non-finite angle, and a link pointing at a gate that carries no angle all raise PrismError.

Gradients

expectation_gradient computes ⟨H⟩ and its exact gradient with respect to the bound parameters by the adjoint method, at a cost independent of the parameter count. Mark parameters with param(slot) while building, then pass parameter_links() through.

import numpy as np
from prism_q import CircuitBuilder, simulate

builder = CircuitBuilder(2)
builder.ry(0.3, 0).param(0)
builder.cx(0, 1)
builder.rz(0.7, 1).param(1)
circuit, links = builder.build(), builder.parameter_links()

hamiltonian = [(1.0, [(0, "Z")]), (0.5, [(0, "Z"), (1, "Z")])]
value, gradient = simulate(circuit).seed(42).expectation_gradient(hamiltonian, links)

A Hamiltonian term is (coefficient, observable). Several gates may share a slot, in which case their gradients accumulate. param() rejects anything but a differentiable gate (rx, ry, rz, rzz, p, pauli_rotation), and the circuit must be unitary.

Quantum error correction

QecProgram exposes the native QEC IR: reset, measure, detector, observable_include, postselect, and noise, with QecBasis, QecNoise, and RecordRef as the supporting types. run() returns a QecResult carrying detector, observable, and measurement arrays as NumPy bool_ matrices, plus logical_error_rates() and survivor_rate(). Programs can also be parsed from text with QecProgram.from_text. See the Noise and QEC guide for the model itself.

detector_error_model() derives the program's DetectorErrorModel for decoding: probabilities() (float64), detector_matrix(), and observable_matrix() (bool, detectors or observables by mechanisms) feed check-matrix decoder constructors directly, detector_coords() carries the per-detector coordinates, decompose_graphlike() returns the form matching decoders need (at most two detectors per mechanism), and to_text() writes the common detector error model text format for file-based decoders:

dem = qp.detector_error_model()
H, p, L = dem.detector_matrix(), dem.probabilities(), dem.observable_matrix()
with open("memory_d3.dem", "w") as f:
    f.write(dem.to_text())

Decoder runs the built-in union-find decoder over a graphlike model, so a memory experiment's logical error rate needs no external decoder. decode takes the (shots, num_detectors) bool detector array and returns the (shots, num_observables) predicted observable flips:

decoder = prism_q.Decoder(dem.decompose_graphlike())
res = qp.run()
predicted = decoder.decode(res.detectors)
failures = (predicted[:, 0] != res.observables[:, 0]).sum()

Errors and typing

Every failure surfaces as prism_q.PrismError, carrying the message from the Rust error. Backend limits, unsupported operations, and invalid arguments all raise it rather than panicking:

import prism_q

try:
    simulate(huge).density_matrix_expectation_values([[(0, "Z")]])
except prism_q.PrismError as exc:
    print(exc)   # backend `density_matrix` is incompatible: circuit has ... qubits

The package ships type stubs (prism_q/_prism_q.pyi) and a py.typed marker, so mypy and Pyright see the full surface.

Backends Deep Dive

PRISM-Q does not have one simulation algorithm. It has nine, each optimal for a different class of circuit. This guide is the task-oriented companion to the architecture reference: it focuses on scaling and when to reach for each one. To select a backend in code, see Choosing a Backend. For which CPU and GPU architectures each backend supports, see the Capability and Support Matrix.

Scaling at a glance

The distributed statevector backend (behind the distributed feature) shards the dense state across MPI ranks; see the Capability and Support Matrix for its status.

Statevector

The default for dense circuits. Exact, fully general, and the fastest option whenever the state fits in RAM. The memory cap is derived from system RAM (overridable with PRISM_MAX_SV_QUBITS). Above it, auto-dispatch falls back to Sparse or MPS.

Stabilizer

If your circuit uses only Clifford gates (H, S, Sdg, SX, SXdg, X, Y, Z, Id, CX, CZ, SWAP, measurement), the stabilizer tableau simulates it in $O(n^2)$ and scales to thousands of qubits. Auto-dispatch selects it whenever the circuit is Clifford-only.

Sparse, MPS, Product, Tensor Network, Factored, Density Matrix

• Sparse wins when the state stays concentrated in a handful of computational-basis states (amplitude pruning keeps the map small).
• MPS trades exactness for polynomial memory in the bond dimension. Ideal for low-entanglement circuits over many qubits.
• Product is the degenerate, entanglement-free case: $O(n)$ memory, $O(1)$ per 1q gate.
• Tensor Network defers contraction until measurement, useful for shallow or structured circuits.
• Factored detects partial independence and simulates sub-registers separately, merging lazily via a Kronecker product computed on demand.
• Density Matrix evolves the full mixed state exactly, for noise studies below the $4^n$ memory ceiling. Explicit dispatch only; Auto never selects it.

For the internal kernels behind each of these, read the architecture reference. For raw speed mechanics, see Performance and SIMD.

Capability and Support Matrix

This page records which CPU and GPU architectures each PRISM-Q backend supports, and where distributed execution stands. CPU backends are written in portable Rust and run on every supported architecture; SIMD acceleration (AVX2/FMA/BMI2 on x86-64, NEON on ARM64) is selected at runtime where a kernel exists, otherwise a scalar path is used.

Legend

Backend support by architecture

The nine CPU backends implement the Backend trait; the distributed statevector backend is a tenth, feature-gated implementation covered by the Distributed column. Planned marks only work the roadmap carries: a ROCm port of the existing CUDA kernels. Backends without a CUDA kernel have nothing to port, so their ROCm cell is No, and the roadmap carries no distributed execution for any backend other than the statevector.

The Clifford+T engines below are not Backend implementations; they serve probability, shot, and observable queries through their own routes (see Clifford+T Simulation).

Notes:

• AVX2/FMA/BMI2 is the x86-64 SIMD tier. The active tier is chosen at runtime (AVX2+FMA, then FMA, then SSE2 baseline). See Threading, SIMD, and Memory Layout.
• NEON is the ARM64 SIMD tier. Backends marked SIMD carry a NEON kernel that mirrors the x86-64 path; the rest fall back to scalar code on ARM64.
• CUDA covers the optional gpu feature. Only the statevector and stabilizer paths have device kernels; every other backend runs on CPU. See the GPU Backend guide.
• Distributed covers the optional distributed and distributed-mpi features. The statevector backend splits the state across MPI ranks with exact results, including gates, measurement, reset, and multi-shot sampling without gathering the dense state. Use simulate(&circuit).distributed(context). A run can start from an injected amplitude vector through .initial_state(...): every rank receives the full 2^n vector and keeps only its own slice.

Shot and observable queries above the dense cap

simulate(...).shots(n), .sample_counts(n), and .expectation_values(...) answer from a dense 2^n vector unless the backend carries its own path. The dense route is capped by system memory (roughly 29 qubits on a 16 GiB host; see PRISM_MAX_SV_QUBITS). Backends marked Native below answer without it and are bounded only by their own representation.

Native sampling is deterministic from the seed alone: the same seed and shot count reproduce the same bitstrings. It is not shot-for-shot identical to the dense route, which consumes its randomness on a different schedule; the distributions agree.

Backends without an observable path return BackendUnsupported naming themselves, so a rejected request says which engine could not serve it rather than blaming the route that selected it.

simulate(...).marginals() reads per-qubit Z expectations rather than a distribution when the resolved backend has an observable path and the circuit routes straight to it. Sparse, MPS, product state, factored, tensor network, distributed, the density matrix, and both stabilizer backends have one, and the dense output cap does not apply on that route; a Clifford circuit with measurements reads its marginals off the tableau at any width. It falls back to the dense distribution on backends without an observable path, on a circuit that splits into independent blocks unless those blocks run as product states, and under a noise model, where the mixture is read densely and the density-matrix memory limit applies instead.

Under a noise model run() and marginals() answer from the exact mixture, so both reject a model carrying readout error instead of serving one: readout acts on the measurement record rather than the state, and it is indexed by classical bit where a marginal is indexed by qubit. shots and sample_counts are the terminals that apply it.

simulate(...).run() is the one terminal that needs the whole distribution, so on the distributed backend it rejects a register past the dense cap up front rather than running first and answering with no distribution.

What kind of answer a result is

Every result type carries a metadata field describing how it was produced: the engine automatic dispatch resolved to, whether that engine can discard state weight, where the state lived, and the shot count for a sampled result.

Approximate marks the route rather than the run. An MPS at bond 256 on a circuit that never fills a bond truncates nothing and still reports Approximate, with a bound of 1.0: the variant answers whether the answer could have been approximated, the bound answers whether it was.

The bound describes the normalized state, which is what every read returns: a truncating MPS does not renormalize its chain, but expectation values, shot sampling, the probability vector and the exported statevector all rescale on read, so the probabilities sum to 1 at any bond cap. The discarded weight the bound reports is error in that state, not weight missing from it.

Auto sends a circuit past the statevector cap to an MPS at a bounded bond dimension, which is the only route those circuits have. It is taken by default and the result says so. simulate(...).require_exact() rejects that route instead, with an error naming the engine it would have used.

Simulate::expectation_values_reported returns the values with a standard error per value on a route that estimates rather than evaluates. An evaluated route reports Exact and no interval, so a caller distinguishes "converged" from "not estimated" without comparing a float against zero.

Not yet supported

These targets are listed so the matrix reflects the roadmap rather than hiding the gaps.

Performance and SIMD

Performance is the primary product requirement. This guide explains the mechanisms that make PRISM-Q fast and the knobs you can turn. The internals live in the architecture reference under Fusion Pipeline and Threading, SIMD, and Memory Layout.

The three levers

1. Fusion collapses many small gate passes into fewer, larger ones before execution, reducing memory traffic over the statevector. It is qubit-count gated and zero-cost when it does not apply.
2. Cache-resident tiling keeps batched gates (MultiFused, Multi2q) operating on L2/L3-sized tiles so repeated passes reuse hot data.
3. SIMD vectorizes the inner complex-arithmetic loop with AVX2+FMA, FMA, and BMI2, with a scalar fallback on non-x86_64.

The levers are ordered, and lever 3 comes with a prior question: can the arithmetic be removed rather than issued faster? A kernel whose operations an algebraic identity or an operator structure deletes is bounded by its memory floor; vectorizing what remains is bounded by the complex-arithmetic issue ceiling, near 23% of FMA peak in the interleaved layout.

Threading

Rayon parallel kernels engage at ≥14 qubits (below that, thread-pool overhead dominates), with MIN_PAR_ELEMS = 4096 per task. The pool defaults to all logical cores.

Determinism

Deterministic partitioning makes unitary evolution and seeded terminal sampling on the dense backends bitwise reproducible at any thread count. Parallel reductions (norms, collapse probabilities, expectation values) are stable to about 1e-12 but not bitwise, and the batched compiled sampler seeds one RNG stream per worker, so its shots reproduce only at a fixed thread count. The per-path contract is in Threading, SIMD, and Memory Layout.

Tuning environment variables

Benchmarking

cargo bench --bench circuits --features parallel       # circuit macrobenchmarks
cargo bench --bench bench_driver --features parallel   # gate microbenchmarks

For current wall-clock numbers across the circuit suite, see the Benchmarks page.

OpenQASM Support

PRISM-Q parses a practical subset of OpenQASM 3.0, with backward compatibility for common 2.0 syntax. The parser converts text directly to the Circuit IR with no intermediate AST.

Parsing and running

#![allow(unused)]
fn main() {
use prism_q::circuit::openqasm;
use prism_q::simulate;

let circuit = openqasm::parse(qasm_str).expect("parse error");
let result = simulate(&circuit).seed(42).run().unwrap();
}

run_qasm(qasm, seed) parses and simulates in one call.

Exporting

#![allow(unused)]
fn main() {
use prism_q::circuit::qasm_export;

let qasm = qasm_export::to_qasm3(&circuit).expect("export error");
}

Export inverts the parser: re-parsing the result gives back the same instruction stream, with gate matrices agreeing to floating-point round-off and inline angles (rx, rz, rzz, p) surviving exactly. Qubits come out as one qubit[n] q register, classical bits as one bit[m] c register, split only where a condition compares against a register narrower than the whole.

A circuit that has been through fusion is not exportable: fused blocks, tiled multi-gate passes, and diagonal batches carry matrices with no OpenQASM spelling, and to_qasm3 returns ExportUnsupported naming the instruction index. Export the circuit before fusing it, or the template a PreparedCircuit binds. QftBlock and PauliRot are the exceptions: export expands the first to its textbook Hadamard, controlled-phase, and swap sequence and the second to its CNOT-ladder lowering on the way out.

Declarations and measurement

OPENQASM 3.0;
include "stdgates.inc";
qubit[3] q;          // OpenQASM 3.0 register
bit[3] c;
h q[0];
cx q[0], q[1];
c[0] = measure q[0]; // OQ3 measurement

OpenQASM 2.0 syntax also works: qreg q[3]; / creg c[3]; declarations and measure q[0] -> c[0]; measurement.

output bit[3] c; declares the register and marks it as the program's result. Every classical bit is reported either way, so the marking costs nothing and changes nothing.

Input parameters

An input declaration names a parameter slot. openqasm::parse_parametric returns the template circuit alongside the Parameters that binds it, in declaration order and under the declared names:

OPENQASM 3.0;
input float[64] theta;
input float[64] phi;
qubit[2] q;
h q[0];
rx(theta) q[0];
cx q[0], q[1];
rz(phi) q[1];

#![allow(unused)]
fn main() {
let (template, params) = openqasm::parse_parametric(qasm)?;
let bound = params.bind(&template, &[0.41, 1.27])?;
let text = to_qasm3(&bound)?;   // angles written out, no `input` line
}

Several gates may read one input, which is the weight sharing Parameters already models: rx(theta) q; over a register links every gate it produces to the same slot, and binding writes one angle to each.

parse itself rejects a program that declares an input, because it has nowhere to take the value and a zero would be a quiet wrong answer. Feed those through parse_parametric, or through PreparedCircuit for a sweep.

An input binds an angle whole, so it may only be the entire angle argument of a directly named parametric gate at the top level. rx(2 * theta), an input on a gate carrying no rotation angle, one reaching a gate, def, for, or guarded body, and one on a modified gate all return UnsupportedConstruct naming the reason rather than binding something the source did not mean.

Supported gates

• Standard / aliases: x, y, z, h, s, sdg, t, tdg, sx, rx, ry, rz, p/phase, cx/CX/cnot, cy, cz, cp/cphase, crx, cry, crz, ch, swap, ccx/toffoli, cswap/fredkin, cu, u1, u2, u3/u/U.

• Qiskit / exporter: sxdg, cs, csdg, csx, ccz, r, xx_plus_yy, xx_minus_yy, ecr, iswap, dcx, c3x, c4x, mcx, rccx, rc3x/rcccx.

• Hardware-native: gpi, gpi2, ms, syc, sqrt_iswap, sqrt_iswap_inv.

• Pauli rotations: r followed by one Pauli letter per qubit argument. rxx, ryy, and rzz are the two-letter cases; rxyz(0.7) q[0], q[1], q[2]; is exp(-i * 0.7 * (X⊗Y⊗Z) / 2) with x on q[0]. rzz resolves to the native two-qubit rotation and a one-letter name to rx/ry/rz; wider strings build the native multi-qubit gate, which the statevector applies in one pass and every other backend receives as its CNOT-ladder lowering.

  The wider spelling is a PRISM-Q extension rather than standard OpenQASM, and to_qasm3 emits it so a round trip preserves the gate instead of a lowering of it. For output another toolchain reads, run circuit::expand_pauli_rotations before exporting.

Other supported constructs

• Gate modifiers: ctrl @, inv @, pow(k) @. inv @ and pow(k) @ also apply to a user gate, a def call, and a gate that lowers to a sequence (u3, ecr, iswap, and peers), reversing or repeating the expanded body. ctrl @ on those is an error, since an expanded body has no controlled form.
• User-defined gate blocks.
• Classical if conditionals, guarding either a single statement or a braced body. A braced body admits any supported statement, measure and reset included, and may nest.
• else and else if arms, and switch with case and default arms. Both lower to guards on the existing condition language rather than new syntax in the IR.
• Parity conditions, if (c[0] ^ c[2]) or if ((c[0] ^ c[2]) == 0).
• Multi-register broadcast, barrier, and an expression evaluator with math functions.

bit[2] c;
qubit[3] q;
c[0] = measure q[0];
if (c[0]) {
  x q[1];
  c[1] = measure q[1];
  if (c[1]) { reset q[2]; }
}

Clifford+T Simulation

Circuits that mix Clifford gates with a modest number of T gates sit between the efficient stabilizer regime and the exponential statevector regime. PRISM-Q offers three strategies. The right one depends on your T-count, qubit count, and whether you need exact answers or can tolerate Monte Carlo error.

These route through the Clifford+T strategies before the standard dispatch tree when the T-count permits.

Accepted gate forms

All three strategies take Clifford gates, T, Tdg, Rz, P, and Rzz as they are, and lower Rx, Ry, and PauliRot to Clifford conjugation around one Rz. A Fused matrix (the parser's u, u3, u2, r, gpi, and gpi2) lowers to the named gate or Rz it equals up to a global phase, or else to its Rz-Ry-Rz Euler triple; a cu lowers when its target is diagonal or a Pauli up to phase (cy, cp, crz, cs), and is rejected otherwise, as are ccx and the other multi-controlled forms. The Pauli engines take any rotation angle. Stabilizer rank needs every lowered Z rotation on the pi/4 grid, and Rz(pi/4), P(pi/4), and their odd multiples count as one T each, so is_clifford_plus_t and the automatic route treat them like T.

Stabilizer rank (src/sim/stabilizer_rank.rs)

Exact probability output remains capped because it returns a dense vector with 2^n entries. Shot sampling uses coherent weighted MPS branches instead of a dense statevector fallback. Clifford gates mutate each branch state, T and Tdg split branches, and measurement computes outcome probabilities from the weighted branch ensemble before projecting every branch to the sampled outcome. This removes the hard qubit-count cap from run_stabilizer_rank_shots; practical scaling is governed by branch count, MPS bond growth, and measurement count.

The dense probability path maintains a weighted sum of stabilizer states. Each T gate doubles the term count via the T = alpha*I + beta*Z decomposition. Clifford gates are O(n²) per term and weighted amplitudes are accumulated for exact probabilities.

run_stabilizer_rank_approx truncates deterministically. After each T it drops the smallest-magnitude terms in excess of the budget, renormalizes the probabilities it returns, and reports the summed 1-norm of what it dropped as discarded_weight. Writing that as d, the returned distribution differs from the exact one by at most 4 * d, and StabRankResult::fidelity_bound converts d to a fidelity floor. The bound assumes worst-case interference between the non-orthogonal branch states, so an aggressive budget drives it to zero. Auto never selects this function: the automatic probability route runs the exact expansion when the T count fits the size-derived budget and falls through to the dispatch tree otherwise, so the pruned expansion is reached only by calling it directly.

Stochastic Pauli Propagation (src/sim/unified_pauli.rs)

Backward-propagates measurement observables as Pauli strings. Clifford gates conjugate in O(1). T gates branch stochastically into two Pauli paths with appropriate weights. Per-path cost O(d×n/64), independent of T-gate count. Returns marginal probabilities via Monte Carlo estimation.

#![allow(unused)]
fn main() {
run_spp(circuit, num_samples, seed) // -> SppResult
}

Deterministic Sparse Pauli Dynamics (src/sim/unified_pauli.rs)

Backward-propagates as a weighted sum of Pauli strings stored in a HashMap. T gates deterministically branch X/Y terms. Identical strings auto-merge. Optional ε-truncation for approximate mode. Exact for small T-counts, approximate with bounded error for larger ones.

#![allow(unused)]
fn main() {
run_spd(circuit, epsilon, max_terms) // -> SpdResult
}

Pauli path propagation under noise (src/sim/unified_pauli.rs)

The same weighted Pauli sum, carried through a noise model. Select BackendKind::PauliPath { epsilon, max_terms } and attach a noise model; the engine answers expectation_values and observable_expectation, and nothing else.

#![allow(unused)]
fn main() {
let noise = NoiseModel::uniform_depolarizing(&circuit, 0.01);
let values = simulate(&circuit)
    .backend(BackendKind::PauliPath { epsilon: 1e-8, max_terms: 1 << 16 })
    .noise(&noise)
    .expectation_values(&observables)?;
}

The error model is worth stating plainly, because it has two independent parts and only one of them is approximate.

The channels are exact. Each one is applied as its adjoint on the Pauli basis rather than as a twirl, so depolarizing, dephasing, thermal relaxation, and amplitude damping all reproduce the density matrix to machine precision at max_terms = 0. Amplitude damping is the case worth naming: it is not unital, and the identity term its adjoint produces from Z is carried rather than dropped. A channel with no Pauli-basis form (custom Kraus, two-qubit Kraus, readout error) is rejected rather than approximated.

The truncation is the approximate part, and it reports its own bound. With max_terms = 0 nothing is dropped and the value is exact. With a budget set, terms whose coefficient magnitude falls below epsilon are dropped once the sum exceeds the budget, and the total dropped magnitude bounds the error in the returned value. That bound is a worst case rather than an estimate: it holds because every remaining operation is a contraction in the Pauli 1-norm.

What decides whether the engine is usable is the term count, not the qubit count. Every non-Clifford rotation in the observable's backward light cone can double the sum; every channel shrinks it. Circuits where noise wins stay cheap at widths no dense representation reaches, and circuits where it does not will hit the budget and report a large discarded mass, which is the signal to use the density matrix or trajectory averaging instead.

In practice the observable's weight is what moves that count, ahead of width and depth, because the backward light cone opens from every letter it starts with. On a two-layer hardware-efficient ansatz under 1% depolarizing, Z on one qubit stays at 11 terms whether the register is 20 qubits or 100; Z on two adjacent qubits fills a 16384-term budget by 30 qubits; a Z on every qubit fills it at 20 and returns a discarded mass larger than the observable's own norm, which is the engine saying the answer is not usable rather than returning a wrong one quietly. Check the reported discarded mass against the precision the caller needs before trusting a truncated run.

You can also build Clifford+T test circuits directly with clifford_t_circuit.

Noise and QEC

PRISM-Q models noise and quantum error correction without falling back to a dense statevector per shot. The machinery is the compiled samplers and the native QEC program IR; this guide shows how they fit together.

Noisy shot sampling

Attach a NoiseModel and sample:

#![allow(unused)]
fn main() {
use prism_q::{simulate, BackendKind, NoiseModel};

let noise = NoiseModel::uniform_depolarizing(&circuit, 0.001);
let result = simulate(&circuit)
    .backend(BackendKind::Statevector)
    .noise(noise)
    .seed(42)
    .shots(1024)
    .unwrap();
}

NoiseModel carries per-instruction depolarizing channels (NoiseOp { qubit, px, py, pz }) and supports readout error and amplitude damping. For Clifford circuits, the noisy compiled sampler propagates noise sensitivity rows and XORs fired channels into each sample, avoiding per-shot state evolution entirely.

Detector sampling

For repeated syndrome extraction, compile_detector_sampler compiles a Clifford circuit with measurement and reset reuse into a packed sampler, then derives detector and observable records as parity rows over the measurement record. Reset reuse becomes fresh qubit aliases, so there is no per-shot tableau replay.

Native QEC programs

When you need detectors, logical observables, postselection, and Pauli-noise annotations as first-class constructs, use the native QEC program IR rather than a Circuit:

#![allow(unused)]
fn main() {
use prism_q::{parse_qec_program, run_qec_program};

let program = parse_qec_program(qec_text).unwrap();
let result = run_qec_program(&program).unwrap();
}

run_qec_program lowers Clifford-compatible programs into the packed compiled sampler. run_qec_program_reference is the per-shot statevector oracle for validating small programs.

Detector error model export

Matching and belief-propagation decoders consume an error model, not raw detector samples. QecProgram::detector_error_model derives one from the program's noise annotations, detectors, and observables, and to_text renders it in the common detector error model text format that external decoders read:

#![allow(unused)]
fn main() {
let model = program.detector_error_model().unwrap();
std::fs::write("memory_d3.dem", model.to_text()).unwrap();
}

Each mechanism carries a probability and the detector and observable indices it flips; detector coordinates pass through from the program. In Python the model also exposes probabilities(), detector_matrix(), and observable_matrix(), the check-matrix triple that in-process decoder libraries accept directly. Matching decoders need at most two detectors per mechanism: decompose_graphlike returns that form, splitting each hypergraph mechanism across existing graphlike ones and erroring loudly when no split exists. See QEC program execution for the derivation semantics and the emitted grammar.

Decoding

UnionFindDecoder decodes sampled detectors against a graphlike model in-process, so the logical error rate of a memory experiment never leaves the tool:

#![allow(unused)]
fn main() {
use prism_q::{UnionFindDecoder, run_qec_program};

let model = program.detector_error_model()?.decompose_graphlike()?;
let decoder = UnionFindDecoder::from_model(&model)?;
let result = run_qec_program(&program)?;
let predicted = decoder.decode_packed(&result.detectors)?;
let failures = (0..result.total_shots)
    .filter(|&shot| predicted.get_bit(shot, 0) != result.observables.get_bit(shot, 0))
    .count();
}

The decoder is weighted union-find with peeling: edges weigh ln((1-p)/p), one-detector mechanisms are boundary edges, and mechanisms flipping no detector bound the achievable logical error rate from below. Construction rejects hypergraph models with a pointer to decompose_graphlike. Decoding is deterministic and allocation-free per shot; large batches decode in parallel. See the decoding section of QEC program execution for the growth and peeling semantics and the validation against the exact ML rate.

Homological sampling

run_shots_homological and ErrorChainComplex model the GF(2) chain complex over noise locations, identifying undetectable error cycles. noisy_marginals_analytical computes marginals in closed form from the parity matrix and noise rates, with no Monte Carlo.

GPU Backend

cargo build --release --features "parallel gpu"
cargo nextest run --features "parallel gpu" --test golden_gpu --test golden_gpu_density_matrix

Constructing a GpuContext compiles the CUDA kernels through NVRTC, one to three seconds on a GTX 1080 Ti. The PTX is shared by every context in the process and cached on disk in prism-q-ptx under the user cache directory (XDG_CACHE_HOME, else LOCALAPPDATA, else HOME/.cache, else the OS temp directory), keyed by device arch, crate version, and a hash of the kernel source, so later processes skip the compile. A missing, unreadable, or corrupt cache file only costs a recompile; delete the directory to force one.

CUDA acceleration covers statevector execution, stabilizer execution, density-matrix execution, and compiled BTS sampling. Seven entry points are available:

• BackendKind::AutoGpu { context } (simulate(circuit).gpu_auto(ctx)). Automatic backend selection with the device opted in. The shape-based decision tree runs unchanged; a selected statevector or stabilizer workload that clears the family's qubit crossover and fits in VRAM runs on the device. Everything else, including a device allocation that fails at init, takes the identical CPU path (the soft VRAM fallback).
• BackendKind::StatevectorGpu { context }. Public dispatch path for statevector GPU execution. It routes through simulate(circuit).backend(kind).seed(seed).run(), keeps fusion and subsystem decomposition, and uses crate::gpu::min_qubits() (default 14, PRISM_GPU_MIN_QUBITS override) to keep small sub-circuits on CPU.
• BackendKind::StabilizerGpu { context }. Public dispatch path for stabilizer GPU execution. Gate application uses a device tableau and one word-grouped batched Clifford kernel (stab_apply_word_grouped). Measurement and reset keep pivot search, row cascade, phase fixup, and deterministic outcomes on the device. The default crossover stays conservative (STABILIZER_MIN_QUBITS_DEFAULT = 100_000, PRISM_STABILIZER_GPU_MIN_QUBITS override) until benchmarks justify lowering it. Direct backend benchmarks should use StabilizerBackend::with_gpu(ctx) to exclude diagnostic readbacks from probabilities(), export_tableau(), and export_statevector(). Golden tests cover every kernel path, including 500q GHZ measure-all.
• BackendKind::DensityMatrixGpu { context }. The exact mixture held in device memory. Explicit only: neither Auto nor AutoGpu selects it, there is no crossover, and there is no host fallback. init budgets the 4^n buffer against the free VRAM and errors before allocating when it does not fit, so an 11 GiB card holds 13 qubits (1 GiB at 13, 4 GiB at 14 plus scratch). The unitary half reuses the dense statevector kernels on the embedded 2n-qubit buffer, and every channel, measurement, and readout sweep runs as a kernel of its own. The noisy Simulate terminals answer from the device mixture exactly as DensityMatrix does. DensityMatrixBackend::new(seed).with_gpu(ctx) is the direct form.
• StatevectorBackend::new(seed).with_gpu(ctx). Direct statevector GPU opt-in. Every instruction routes to CUDA after the context is attached. No crossover or subsystem decomposition applies.
• StabilizerBackend::new(seed).with_gpu(ctx). Direct stabilizer GPU opt-in for kernel benchmarks and targeted correctness tests.
• run_shots_compiled_with_gpu (or CompiledSampler::with_gpu(ctx)). GPU BTS sampling for flat sparse parity. The path launches one kernel per 65_536-shot chunk, uses random bits generated on the host, and preserves the CPU sample_bts_meas_major layout. The sampler caches sparse parity CSR arrays, packed reference bits, and reusable scratch on the device. It is active only when num_shots >= BTS_MIN_SHOTS_DEFAULT (131_072 by default, PRISM_GPU_BTS_MIN_SHOTS override). sample_bulk_packed_device returns a DevicePackedShots handle. Marginals reduce to one counter per measurement row on the device. Exact counts use a bounded device hash reduction for up to 8 packed measurement words when the compact result is cheaper to transfer than the full shot matrix. Otherwise the API uses a host copy for correctness.

When a GPU context is attached, Backend::init allocates state on the device instead of a host Vec<Complex64> and every instruction routes to a CUDA kernel. On the hard statevector path (StatevectorGpu, with_gpu), a state that does not fit the currently free VRAM is rejected at init with an error naming the requested and free device memory; GpuContext::max_qubits_for_statevector reports the advisory cap from free memory.

The four BackendKind entry points are also reachable from Python, from a build carrying the gpu feature. See Python Bindings.

Module layout (src/gpu/)

Kernel coverage

Every variant in the Gate enum has a dedicated kernel. Batched variants (BatchPhase, BatchRzz, DiagonalBatch, MultiFused { all_diagonal: true }) use LUT kernels that consume the same host table builders as the CPU path. Non-diagonal MultiFused uses a shared memory tiled kernel (apply_multi_fused_tiled, TILE_Q = 10, TILE_SIZE = 1024) over a chosen set of ten qubits per pass: the five lowest qubits, which keep a warp's loads contiguous, plus up to five of the sub-gates' higher targets. A MultiFused over n qubits therefore takes about (n - 5) / 5 passes, each applying its sub-gates in shared memory. A pass with fewer than three sub-gates falls back to per gate launches. Multi2q still launches once per sub-gate; rare in practice.

PTX template substitution: the CUDA C source lives in .cu files beside the Rust launchers and reaches KERNEL_SOURCE_TEMPLATE through include_str!. kernels/dense.cu carries placeholders such as {{BP_TABLE_SIZE}} and {{TILE_Q}}. The kernel_source() function substitutes them at device construction from the Rust constants in src/backend/statevector/kernels.rs, keeping CPU and GPU in sync.

Correctness

tests/golden_gpu.rs compares GPU amplitudes against the CPU statevector within 1e-12 for every gate variant, the fusion paths, and the BackendKind::StatevectorGpu public dispatch path at the crossover boundary; tests/golden_gpu_density_matrix.rs does the same for the full density matrix buffer. Both suites skip when no device opens, and no CI job opens one. scripts/test-gpu.ps1 runs them with PRISM_REQUIRE_GPU=1 so a missing or unusable device fails the run.

Shot reproducibility

Two limits bound what a seed guarantees, and neither is visible from the golden equality tests.

• CPU against GPU: agreement in distribution, not bit for bit. Both paths draw the same RNG stream for a given shot seed, but the device reduction that produces a measurement probability sums in tree order with FMA contraction, so it can differ from the host sum in the last ulp, and a uniform draw landing between the two flips that outcome and every outcome after it. Where the probability is a dyadic rational (0.5, 1.0, and the amplitudes reachable from Clifford gates) both sums are exact and the shots do match, which is what statevector_gpu_mid_measure_shots_match_cpu pins. statevector_gpu_shot_frequencies_match_cpu_off_dyadic pins the general case: equal frequencies within 5 sigma after Rx(0.3).
• GPU BTS sampling: reproducible at a fixed Rayon thread count. Above MIN_PAR_DRAWS random bits per chunk, fill_random_bits seeds one stream per worker and partitions the draws by rayon::current_num_threads(), so the same seed on a host with a different worker count produces different shots. Below that threshold the serial single-stream path runs and the seed reproduces outright. Pin RAYON_NUM_THREADS when byte-identical shot payloads matter across machines.

Current limits

• Device placement is silent. Circuits below the crossover run on the host, and the AutoGpu soft VRAM fallback degrades to host execution without a report; nothing user-visible says whether a run executed on the device.
• Stabilizer probabilities(), export_tableau(), and export_statevector() read back to the CPU.
• Every trajectory shot rebuilds the backend, reallocating the device buffer (measured at 0.1 ms per shot at 20 qubits, so not a practical cost). Custom Kraus branch probabilities come from an on-device reduced-density-matrix reduction (rdm_qubit), not a full-state readback.
• Kernel design and crossover analysis live in the module docstrings on src/gpu/kernels/dense.rs.

Architecture: Overview and Layered Design

This section is the technical reference for how PRISM-Q is built. See the Glossary for definitions of terms used throughout.

Goals

• Primary: Fastest practical quantum circuit simulation in Rust.
• Correct simulation of supported gate sets across multiple backend strategies.
• Clean backend plugin model. New simulation strategies can be added without touching the core.

Non-goals

• Full OpenQASM 3.0 compliance (supports a practical subset).
• GUI or notebook integration (library-first).
• Hardware backend / QPU connectivity.

Layered design

A circuit flows top to bottom: text is parsed into a backend-agnostic IR, optimized by the fusion pipeline, then dispatched by the simulation engine to one of the backends or a compiled sampler.

flowchart TD
    U[User / Application]
    API["Public API — run_qasm, simulate (src/lib.rs)"]
    P["OpenQASM 3.0 Parser — &str to Circuit IR (src/circuit/openqasm.rs)"]
    IR["Circuit IR — gates, measures, barriers, conditionals (src/circuit/mod.rs)"]
    F["Fusion Pipeline — cancel, fuse, reorder, batch (src/circuit/fusion.rs)"]
    E["Simulation Engine — dispatch, decompose, execute (src/sim/mod.rs)"]
    U --> API --> P --> IR --> F --> E
    E --> B[Backends]
    E --> C["Compiled Samplers — shot-based (src/sim/compiled, noise.rs, homological.rs)"]
    B --> SV[Statevector]
    B --> TN[Tensor Network]
    B --> MPS[MPS]
    B --> SP[Sparse]
    B --> PR[Product]
    B --> ST[Stabilizer]
    B --> FA[Factored]

The remaining pages in this section follow that flow: the parser and circuit IR, the fusion pipeline, the simulation engine and dispatch, the individual backends, the compiled samplers, the native QEC program IR and its execution path, the threading, SIMD, and memory layout, and the error model and public API surface.

Parser and Circuit IR

Parser

Handwritten parser targeting a practical OpenQASM 3.0 subset. It processes input line by line and converts &str directly to Circuit IR with no intermediate AST.

Supported: qubit/bit declarations, OpenQASM standard gates and aliases (x, y, z, h, s, sdg, t, tdg, sx, rx, ry, rz, p/phase, cx/CX/cnot, cy, cz, cp/cphase, crx, cry, crz, ch, swap, ccx/toffoli, cswap/fredkin, cu, u1, u2, u3/u/U), Qiskit and exporter gates (sxdg, cs, csdg, csx, ccz, r, rzz, rxx, ryy, xx_plus_yy, xx_minus_yy, ecr, iswap, dcx, c3x, c4x, mcx, rccx, rc3x/rcccx), hardware-native gates (gpi, gpi2, ms, syc, sqrt_iswap, sqrt_iswap_inv), gate modifiers (ctrl @, inv @, pow(k) @), user-defined gate blocks, classical if conditionals with a single statement or a braced body, multi-register broadcast, measure, barrier, expression evaluator with math functions. OpenQASM 2.0 backward compatibility (qreg/creg, measure q -> c syntax).

Unsupported: for/while loops, subroutines, classical expressions beyond if.

See the OpenQASM Support guide for a user-facing walkthrough.

Circuit IR

Circuit holds num_qubits, num_classical_bits, and Vec<Instruction>. Instructions are an enum:

Targets use SmallVec<[usize; 4]>, inline storage for up to 4 qubits, no heap allocation for typical gates.

Guarded regions

Region carries a condition and a body that runs, in order, only when the condition holds against the classical bits as they stand when control reaches it. The body admits any instruction, measurement and reset included, and may nest to MAX_REGION_DEPTH. The body is boxed, so Instruction stays at 96 bytes.

GuardedRegion caches the sorted union of the qubits its body touches, nested bodies included. Passes read that set rather than re-walking the body: fusion flushes exactly those qubits at the region boundary and is transparent elsewhere, and no pass fuses across the boundary because a region is not an Instruction::Gate.

The body itself is fused, by fuse_region_bodies running the same pipeline over it as an ordinary instruction list on the same register, nested bodies included. That is body-local only and does not weaken the boundary above: a body is fused as a unit, which is sound because it executes as a unit. Leaving it unfused cost 73.5% of dynamic/guarded_region/20.

Conditional is the single-gate lowering of the same construct. if (c) x q[0]; keeps that form, so the common guarded gate costs no allocation; anything else becomes a Region. Build either through circuit::guarded, which picks the form and returns None for an empty body.

A circuit holding a region runs once per shot: measurement-conditioned execution has no single evolved distribution to sample, so the compiled samplers reject it and the run falls back to replay. Routes requiring a unitary circuit (adjoint gradients, exact expectation values, Pauli propagation, stabilizer-rank probabilities) reject a region for the same reason they reject a bare conditional. Noise models index one event slot per instruction, so they reject a region rather than leave its body noiseless.

Circuit::fold_static_guards runs once per shots or counts call and resolves the guards that cannot depend on a measurement. A condition reading only bits no preceding measurement writes is a function of the initial classical state, which every backend zeroes, so the guard is statically dead or statically taken and is dropped or inlined. That returns a circuit whose only guard can never fire to the terminal-measurement sampling path, worth 252x on dynamic/dead_region/16 at 1000 shots. A circuit with no guard borrows through unchanged.

Condition language

ClassicalCondition is a pure function of a bit slice: BitIsOne, BitIsZero, RegisterEquals and RegisterNotEquals over a contiguous range read as u64, and Parity over an arbitrary bit set compared against an expected value. The bit set is boxed, so the enum stays at the size the register variants set. The language is closed under negation (ClassicalCondition::negate), which is what lets else lower rather than grow the instruction.

else emits the then guard followed by a second guard carrying the negated condition, and switch emits one RegisterEquals guard per case label with the default arm nested once per label to spell the conjunction of their negations. Both lowerings re-read the classical bits after an earlier body has run, so the parser rejects a source whose body measures into a bit its own guard reads.

Gate enum

Gate is a Clone enum kept at 16 bytes. Simple variants carry parameters inline. Composite variants use Box to stay within the 16-byte budget for cache-friendly dispatch.

Fusion Pipeline

Gate optimizations before execution, gated by qubit count thresholds. Every pass returns Cow<Circuit>. Borrowed when no optimization applies, so circuits that do not benefit pay zero overhead.

flowchart TD
    IN[Input Circuit] --> PR["fuse_region_bodies (always): each guarded body through this same pipeline"]
    PR --> P0["cancel_self_inverse_pairs (always)"]
    P0 --> P0r["fuse_rzz (always): CX·Rz·CX to Rzz"]
    P0r --> P0b["fuse_batch_rzz (>=16q): N×Rzz to BatchRzz"]
    P0b --> G{"qubits >= MIN_QUBITS_FOR_FUSION (10)?"}
    G -- no --> OUT[Output Circuit]
    G -- yes --> P1["fuse_single_qubit_gates (>=10q)"]
    P1 --> P1r["reorder_1q_gates (>=10q)"]
    P1r --> P1c["cancel_self_inverse_pairs (>=10q)"]
    P1c --> P1f["fuse_single_qubit_gates re-fuse (>=10q)"]
    P1f --> P2q["fuse_2q_gates (>=12q): CX/CZ + adjacent 1q to Fused2q"]
    P2q --> P2qb["fuse_same_pair_2q_blocks (>=12q)"]
    P2qb --> P2["fuse_multi_1q_gates (>=14q) to MultiFused"]
    P2 --> P2qr["reorder_disjoint_fused2q (>=12q)"]
    P2qr --> Pm2q["fuse_multi_2q_gates (>=12q) to Multi2q"]
    Pm2q --> Pcp["fuse_controlled_phases (>=16q) to BatchPhase"]
    Pcp --> Pdb["fuse_diagonal_batch (>=16q) to DiagonalBatch"]
    Pdb --> Ppp["batch_post_phase_1q (>=18q)"]
    Ppp --> OUT

Threshold constants

Payload capacities

The batched gates carry a lookup table sized at compile time, so the pass that emits them is what keeps the payload inside it. Both caps are declared on the gate payload (BatchRzzData::MAX_EDGES, BatchPhaseData::MAX_PHASES) and pinned to the kernel table shape by a compile-time assertion; the kernels assert on entry in release builds as well, so a producer that outgrows a table fails loudly instead of dropping work.

Splitting is sound because both payloads hold mutually commuting diagonal terms. A repeated (control, target) pair folds into the entry already present rather than adding a second one, which both keeps the two paths in agreement (the BMI2 kernel indexes one bit per distinct qubit, so a repeated target has no bit of its own) and bounds a chain by the qubit count.

DiagonalBatch instead declines at the kernel: build_diagonal_batch_tables returns None when the grouping does not fit and the backend runs the per-element path.

Plan capture and replay

A variational sweep holds one gate sequence and varies only the angles. Fusion decides the same block structure at every point, so PreparedCircuit settles it once and rebinds against it.

What is reusable is the plan, not the matrices: a changed angle changes every fused matrix it feeds. FusionPlan therefore records a recipe per angle-derived payload, a list of template instructions and how each one's matrix enters the product. Replay recomputes the products; it never caches them. Nested recipes splice by rewriting a placement flag rather than by materializing the inner product, which is sound because both widening to a pair and SWAP conjugation are multiplicative.

The passes record this under a Tracer that is inactive on the ordinary path, so a fusion outside the prepared form allocates what it always did.

A few decisions read a matrix rather than the gate sequence, and those the plan cannot assume:

Each guard is re-tested per binding, and a binding that flips one falls back to running the pipeline, so the stream matches an independently fused circuit either way. The 1q diagonality guard has to live on the run rather than on a payload site: the block is often absorbed into a Fused2q and keeps no 2x2 of its own, while the reorder it steered has already happened.

DiagonalBatch and BatchPhase have no recipe, so a template reaching either declines capture and re-fuses per binding.

A gate the passes leave untouched still needs a site when it carries an angle. PauliRot holds its angle in a boxed payload rather than in the enum slot, and replay patches it exactly as it patches an Rz; without that site a bound rotation would keep the template's angle while every other payload moved.

Fusion cost against apply cost

Fusion cost tracks instruction count and is close to flat in qubit count, while gate application is 2^n. The ratio therefore moves by an order of magnitude across the useful range: for hardware_efficient_ansatz(n, 5) on this project's reference host, fusion is about 22% of a run at 12 qubits and about 0.3% at 20.

Simulation Engine and Dispatch

Backend trait

#![allow(unused)]
fn main() {
pub trait Backend {
    fn name(&self) -> &'static str;
    fn init(&mut self, num_qubits: usize, num_classical_bits: usize) -> Result<()>;
    fn apply(&mut self, instruction: &Instruction) -> Result<()>;
    fn classical_results(&self) -> &[bool];
    fn probabilities(&self) -> Result<Vec<f64>>;
    fn num_qubits(&self) -> usize;

    // Optional overrides:
    fn apply_instructions(&mut self, instructions: &[Instruction]) -> Result<()>;  // batch apply
    fn supports_fused_gates(&self) -> bool;   // false for symbolic backends (stabilizer)
    fn export_statevector(&self) -> Result<Vec<Complex64>>;  // for backend transitions
}
}

Contract: init before apply. Instructions arrive in circuit order. Measurement is destructive. Deterministic given same RNG seed.

Entry points

Orchestration layer in src/sim/mod.rs.

RunOutcome::probabilities is None only when the selected backend has no dense probability terminal for the requested circuit: the backends built to run past the dense cap (sparse, MPS, stabilizer, factored, product) above PRISM_MAX_PROB_QUBITS, or a register too wide to index. A statevector or tensor-network run above that cap is an error naming it, since the register already fits the state, and other probability extraction failures propagate as errors too. marginals() requires either a direct Pauli marginal route or backend probability output; it returns BackendUnsupported instead of fabricating uniform marginals when neither path is available. Stochastic and deterministic Pauli marginal backends accept only unitary circuits of Clifford gates and Pauli rotations without measurement, reset, or conditional instructions: T, Tdg, Rz, P, and the two-qubit Rzz branch natively, while Rx, Ry, and multi-qubit PauliRot strings lower to Clifford conjugation around one Rz before the run, Fused matrices lower to the named gate they equal up to phase or to an Euler triple, and a Cu with a diagonal or Pauli target lowers to Z rotations and Cliffords. Automatic dispatch still routes only Clifford+T circuits to SPD, with rotations at multiples of pi/4 counted as Clifford or T; a circuit carrying arbitrary rotation angles reaches the Pauli engines by explicit backend selection.

Noise across the terminals

A noise model reaches a terminal by one of two routes, and which one applies is fixed by the selected backend rather than by the terminal.

Backends holding a per-shot pure state average trajectories: each shot re-evolves the circuit with the channels sampled, so a distribution converges as 1/sqrt(shots). Only shots and sample_counts take that route, since a single trajectory is not an answer to run, marginals, or expectation_values.

The density matrix holds the mixture instead of a trajectory, so shots cannot mean "replay the circuit per shot". It means one exact evolution followed by a draw per shot from the resulting distribution. Every terminal reads that one evolution: run and marginals return the exact noisy distribution, expectation_values returns the exact Tr(rho P), and shots and sample_counts carry sampling noise but no trajectory variance. Readout error is applied to the drawn outcomes rather than to the state, on an RNG stream of its own, which is why run and marginals reject a model carrying it rather than returning a state distribution that the sampled terminals would contradict.

The mixture holds every measurement branch at once, which is what makes it exact and also what it cannot undo. A circuit with mid-circuit measurement or classical conditioning is rejected on this route, with or without a noise model attached: the outcome that a later gate would have been conditioned on was never fixed. The rejection sits on the evolution itself, so density_matrix_expectation_values refuses the same circuits the Simulate terminals do. Those circuits stay on trajectory averaging. This is the same property that makes the density matrix the mixture oracle rather than a comparable participant in the branching families of tests/conformance_matrix.rs.

expectation_gradient rejects a noise model on every backend, because the adjoint method backpropagates through a pure state. expectation_gradient_shift accepts one on the density-matrix kinds: the channels do not depend on the shifted angle, so the shift rule holds and each of the 1 + 2 * links evaluations reads the exact mixture. Every other backend rejects the pair, naming the density matrix.

Result provenance

Every terminal returns a result carrying a RunMetadata: the resolved engine, whether that engine is exact, where the state lived, and the shot count for a sampled result. Auto selecting an approximate backend is disclosed by the result, which is what makes require_exact() an opt-out: rejecting by default would remove the only route an oversize non-sparse circuit has.

require_exact() resolves the route from the circuit and errors before allocating, so it does not pay for state it would discard. Sparse Pauli dynamics truncates on coefficient magnitudes it only learns while propagating, so that route cannot be decided in advance and is caught by a second check on the finished result. Both checks run in every terminal.

Auto-dispatch decision tree

flowchart TD
    A[BackendKind::Auto] --> E{Entangling gates?}
    E -- none --> PS["ProductState (O(n))"]
    E -- yes --> CL{All Clifford?}
    CL -- yes --> STB["Stabilizer (O(n^2))"]
    CL -- no --> MEM{Above memory limit?}
    MEM -- "yes, sparse-friendly" --> SPR["Sparse (O(k))"]
    MEM -- "yes, otherwise" --> MPS["MPS (bond dim 256)"]
    MEM -- no --> IND{Partial independence?}
    IND -- yes --> FAC["Factored (split-state)"]
    IND -- no --> SV["Statevector (exact)"]

Memory limit is dynamically computed from available system RAM (50% budget, capped at 33 qubits). Overridable via PRISM_MAX_SV_QUBITS environment variable. Falls back to 28 qubits (4 GB) when detection unavailable.

For a user-facing version of this decision, see Choosing a Backend.

Start states other than |0...0>

Simulate::initial_state bypasses the tree above entirely. Every branch of it reads circuit structure alone and is sound only from the all-zero start: a Clifford circuit yields a stabilizer state when its input is one, the product state and the subsystem split assume an unentangled input, and the Pauli engines propagate observables back to |0...0>. Picking one of them for an arbitrary start state returns a wrong answer rather than an error, so initial_state_plan in src/sim/dispatch.rs constrains the route instead of consulting it: Auto resolves to the statevector, StatevectorDistributed starts the sharded statevector from it (every rank receives the full vector and keeps the slice its rank bits select), DensityMatrix accepts one as the pure mixture |psi><psi|, and every other kind returns IncompatibleBackend. Auto needs no memory check on that path, since a caller holding 2^n amplitudes can already afford the dense state.

The amplitude vector is validated before the run: 2^n entries for the circuit's n qubits, every component finite, and unit norm to 1e-9. An unnormalized vector is rejected rather than rescaled, because the statevector's deferred-normalization factor is reset to 1 by the load and a silent rescale would hide the error inside it.

Subsystem decomposition

Union-find detects independent qubit groups in O(n·α(n)). Each block runs separately with per-block Auto dispatch. Results merge lazily via Probabilities::Factored, a Kronecker product computed on demand per element in O(K), avoiding the O(2^N) dense materialization unless explicitly requested.

Block-level Rayon parallelism when all blocks are <14 qubits (avoids oversubscription with block-internal parallelism).

Temporal Clifford decomposition

For Clifford+T circuits: Clifford prefix runs on the Stabilizer backend, state is exported to Statevector for the non-Clifford tail. Saves exponential memory for circuits with a long Clifford preamble.

Expectation-value gradients

Two methods, chosen by the caller rather than by the engine: the adjoint method is the default and the reference, parameter shift is the fallback for what the adjoint declines. Selection is explicit because the two differ by a factor of 2 * links in evaluation count, which is not a difference a caller should discover from a wall clock.

Adjoint method

run_expectation_gradient(circuit, hamiltonian, params, seed) and simulate(circuit).seed(seed).expectation_gradient(hamiltonian, params) compute ⟨H⟩ and the exact gradient d⟨H⟩/dθ for a weighted Pauli-sum Hamiltonian H = Σ c_k P_k, at a cost independent of the parameter count. Implementation in src/sim/gradient.rs.

The method back-propagates two statevectors. With U = U_L…U_1 and |φ⟩ = U|0⟩:

1. Forward pass (unfused) keeps |φ⟩. Build |λ⟩ = H|φ⟩; the value is Re⟨φ|λ⟩.
2. Sweep i = L…1. For a trainable gate with generator G_i, accumulate Im⟨λ|G_i|φ⟩ (projector form for P), then step both states back through U_i† (Gate::inverse()).

The ⟨λ|G|φ⟩ sandwich generalizes the forward pauli_expectation_from_masks kernel to two vectors (pauli_sandwich, Rayon-parallel at 16+ qubits).

Differentiable gates are Rx, Ry, Rz, Rzz, P, and PauliRot (identified by Gate::pauli_generator, a method, so Gate stays 16 bytes). A multi-qubit Pauli rotation is exp(-iθP/2) like the named rotations, so its generator is the string itself: GeneratorKind::RotPauli borrows the letters off the gate, and the sandwich takes the masks they imply rather than a per-variant special case. Trainable links on other gates, non-unitary instructions, and QftBlock are rejected. Parameter identity is an index-based side table (Parameters, instruction→slot links recorded by CircuitBuilder::param); many gates may share a slot.

Differentiation runs on the unfused instruction stream so each gate keeps a 1:1 correspondence with its generator (fusion would erase both the stored angle and that correspondence). Two prunings cut work without changing results: the sweep stops at the earliest in-cone trainable gate (a non-trainable prefix costs no inverse applications), and a trainable gate outside the Hamiltonian's inverse light cone has a provably zero gradient, so its sandwich is skipped.

Memory is two statevectors, so the qubit ceiling is about one below a single run. Only the statevector backend is supported.

Parameter shift (fallback)

run_expectation_gradient_shift(circuit, hamiltonian, params, seed) and simulate(circuit).backend(kind).seed(seed).expectation_gradient_shift(hamiltonian, params) evaluate d⟨H⟩/dθ = (f(θ+π/2) - f(θ-π/2)) / 2 per trainable gate, exactly, from forward expectation_values calls alone. That is what makes it the fallback: it inherits whatever the selected backend can represent, so it serves Sparse, MPS, Factored, ProductState, DensityMatrix, and Distributed, widths past the statevector cap, and circuits containing QftBlock. A backend with no native observable path reports BackendUnsupported naming itself.

The rule is exact because each differentiable gate is exp(-iθG/2) with G of eigenvalues ±1, making ⟨H⟩ a degree-1 trigonometric polynomial in that angle. P(θ) is the apparent exception, its generator being the projector |1⟩⟨1| with eigenvalues {0, 1}, but P(θ) = e^{iθ/2} Rz(θ) and that scalar cancels against its conjugate in ⟨ψ|H|ψ⟩ wherever the gate sits, so the same shift applies.

The differentiable gate set is the same one the adjoint takes: Rx, Ry, Rz, Rzz, P, and PauliRot are the Gate variants carrying a rotation angle, so parameter shift reaches no gate the adjoint rejects. Its reach is backends and circuit shapes, not gates. A PauliRot on a backend without the native kernel is shifted in the same place, because the ladder expansion happens below the forward evaluation the rule calls.

Gates sharing a parameter slot are shifted one at a time and their contributions summed. Shifting them together is a different quantity: two Rx(θ) on one qubit under ⟨Z⟩ give cos 2θ, whose joint ±π/2 shift is zero rather than -2 sin 2θ.

Cost is 1 + 2 * links circuit evaluations against the adjoint's one, and there is no light-cone pruning to recover any of it: a trainable gate with a provably zero gradient is still evaluated twice. Evaluations run in sequence, because each already drives a Rayon-parallel backend run and holding several in flight would multiply peak state memory by the parameter count, which is the resource this path exists to stay under.

Backend dispatch variants

All BackendKind variants:

Backends

PRISM-Q ships nine CPU backends, an optional CUDA path attached to the statevector and stabilizer backends, and a feature-gated distributed statevector backend that shards the dense state across MPI ranks. The simulation engine picks a backend automatically (the density matrix, tensor network, and distributed backends are explicit-dispatch only), or you can select explicitly. For a task-oriented version of this material, see the Backends Deep Dive guide.

The diagrams below are rendered directly from PRISM-Q's own SVG circuit renderer.

Reset semantics

reset is the channel rho -> |0⟩⟨0| ⊗ tr_q rho on every backend: the qubit is traced out and replaced by |0⟩, leaving the rest of the register in the mixture the trace produces. Projecting onto |0⟩ and renormalizing is not equivalent. The two agree only when the reset qubit is unentangled; when it is entangled, projection also collapses its partners into the branch correlated with the |0⟩ outcome. Resetting qubit 1 of a Bell pair leaves ⟨Z0⟩ = 0 under the channel and ⟨Z0⟩ = 1 under projection.

A backend holding a single pure state cannot represent the resulting mixture, so it runs one trajectory of the channel: sample the measurement outcome, collapse onto it, and apply X when the outcome is 1. Averaged over shots that reproduces the channel, and a reset consumes one draw from the backend's RNG stream. The density-matrix backend holds the mixture and applies the channel directly, with no draw. tests/reset_channel.rs pins the contract across backends against the density-matrix oracle.

Memory budget

A circuit that does not fit in memory is an error, not a fallback. No backend silently hands the work to a different one when its state would not fit: it returns PrismError::IncompatibleBackend naming itself, the qubit count, the cap, and the environment variable that overrides it. Choosing a different backend is the caller's decision, and BackendKind::Auto makes it from circuit structure before any backend is constructed.

The check lives in Backend::init, which is the one point every execution path passes through before reserving its state. Putting it there means a caller that drives a backend directly, through run_on rather than simulate, gets the same guard as one that goes through dispatch.

The five growth caps are deliberately independent of PRISM_MAX_SV_QUBITS: the sparse, factored, MPS, tensor network, and factored stabilizer backends exist to run above the statevector cap, so lowering that cap to steer routing must not shrink what they may hold. Their defaults come from the same detected-memory budget.

The factored stabilizer cap is the one that is not a 2^n amplitude count. A stabilizer cluster costs O(n^2 / 64) words, so a dense cap is the wrong scale here: it would hold a cluster to the dense backends' qubit ceiling, far below the widths reached by the Clifford circuits at 128 qubits and above that dispatch selects that backend for.

The density-matrix cap is the tighter of its own override and half the statevector cap, computed in one place so dispatch-time validation and the backend's init guard cannot disagree about where the ceiling is. Raising PRISM_MAX_DM_QUBITS past that bound needs PRISM_MAX_SV_QUBITS raised with it, which is what the rejection says: the backend reports it itself rather than surfacing an error naming the statevector it allocates internally. When physical memory cannot be detected the caps are disabled and a warning is printed, because guessing a budget is worse than saying the budget is unknown.

Parallel noisy trajectories are the one path holding more than one state at a time: each Rayon thread runs its own backend, so peak memory is threads * state(n). That path is restricted to circuits below 14 qubits, where a statevector replica is 256 KiB and a full thread pool stays in the tens of megabytes. Above it trajectories run serially with one live backend, bounded by the ordinary state cap.

The three growth paths that once sat outside this contract are bounded at their growth events, each rejecting with an error naming its own backend before the allocation: a factored sub-state merge checks the merged block width against the statevector cap, the sparse map checks a branching gate's worst-case fan-out against the entry cap (so a rejection can fire one gate early on a state that would have deduplicated below it), and MPS gate application checks its live contraction-buffer total against the statevector budget. The MPS check bounds the workspace, not max_bond_dim itself: a large cap on a circuit whose bonds stay small is fine, and the rejection fires only when the bonds actually grow past what memory holds.

Statevector

Full-state simulation in a flat Vec<Complex64> of 2^n amplitudes. The primary backend for circuits up to ~28 qubits.

Gate kernels use enum dispatch with specialized routines for CX, CZ, SWAP, Cu, MCU, Rzz, BatchRzz, BatchPhase, DiagonalBatch, MultiFused, and PauliRot (one pass over (j, j ^ xmask) pairs for exp(-i θ P / 2), a parity-phase sweep when the string is Z-only; backends without the kernel receive the CNOT-ladder lowering from expand_pauli_rotations). Single-qubit gates go through PreparedGate1q with FMA-vectorized SIMD. MultiFused gates use a three-tier tiled kernel (L2 16K / L3 131K / individual passes) for cache locality. MultiFused batches where all gates are diagonal dispatch to a dedicated fast path (1 complex multiply/element vs 4+2 for full 2×2).

Rayon parallelism at ≥14 qubits with par_chunks_mut and MIN_PAR_ELEMS = 4096 per task. BMI2 _pext_u64 accelerates BatchPhase, BatchRzz, and DiagonalBatch LUT indexing.

Deferred measurement normalization: pending_norm accumulates normalization factors without full-state scaling passes. Zero-cost for circuits without measurements.

The Quantum Fourier Transform is a representative statevector workload, dense with controlled-phase gates that the fusion pipeline batches:

Stabilizer

Aaronson-Gottesman bit-packed tableau for Clifford circuits. O(n²) time and space. Scales to thousands of qubits. Gate kernels use wordwise bitwise ops and popcount for phase computation. Supports H, S, Sdg, SX, SXdg, X, Y, Z, Id, CX, CZ, SWAP, plus measurement, reset, and classical conditionals.

Word-group batching fuses multiple 1q gate flushes into single tableau passes. Type-grouped masks apply all gates of the same Pauli type with one wordwise op instead of per-gate dispatch. Sparse Generator Indexing (SGI) tracks per-qubit active generator lists, enabling targeted row operations instead of full-tableau scans. Lazy destabilizer materialization defers destabilizer rows until probabilities are requested.

Probability extraction uses coset-based enumeration with GF(2) Gaussian elimination. O(2^k) where k is the number of non-diagonal generators, rather than O(2^n).

Factored Stabilizer (FactoredStabilizerBackend): Per-cluster tableaux with dynamic merging. Starts with one qubit per cluster. Cross-cluster 2q gates merge tableaux. Measurement and reset can split independent sub-tableaux again. Independent subsystems avoid full-tableau work when product structure is preserved.

Sparse

HashMap<usize, Complex64> for states with few non-zero amplitudes. O(k) memory. Entries at or below a pruning threshold on |a|² (1e-16) are removed after gates that can shrink or cancel amplitudes, and the kept entries are rescaled so the state keeps its norm; a run whose threshold was raised above the default reports Approximate with a fidelity bound derived from the dropped weight. Best for circuits whose support stays concentrated in computational-basis states at large qubit counts.

The map's per-entry gate cost is about 16x the statevector's per-amplitude cost on a mixed diagonal and permutation workload (the sparse/densify bench rows), so a state that densifies past roughly 1/16 load factor runs slower than a dense vector at the same width would. There is deliberately no mid-run handoff to the statevector: automatic dispatch selects this backend only above the statevector memory cap, where the dense state exceeds the memory budget, and a run that branches past the entry cap rejects the gate rather than degrading silently. The map is keyed by a usize basis index, so a circuit wider than usize::BITS qubits is rejected at init and automatic dispatch sends it to MPS instead. An explicitly selected sparse run on a densifying circuit degrades in place, measured at up to 24x the dense cost when fully dense at 20 qubits.

MPS (Matrix Product State)

Chain of rank-3 tensors with adaptive bond dimension (default max 256). O(n·χ²) memory. Single-qubit gates absorb via FMA-vectorized SIMD over bond-dimension slices. Two-qubit gates contract adjacent sites, apply the gate, then SVD-truncate back. Non-adjacent gates route through SWAP chains.

Hybrid SVD dispatch: faer (bidiag+D&C) for matrices with m×n ≥ 256, hand-rolled Jacobi for small matrices.

Product State

Per-qubit [Complex64; 2] storage. O(n) memory, O(1) per single-qubit gate. Rejects entangling gates. Selected automatically for circuits with no 2q gates.

Shots and Pauli expectations answer from the per-qubit states rather than the 2^n probability vector, so both stay O(n) and the backend runs queries at widths no dense route reaches. See Sampling Architecture.

Tensor Network

Deferred contraction planned on metadata: a greedy min-size pass picks the pair order, seeded noisy restarts rerun it when the greedy tree's peak intermediate grows large, and the kernel replays the winner. Gates append tensors; contraction happens lazily at probability extraction, where the PRISM_MAX_PROB_QUBITS cap guards the dense readout and an explicit run past it errors naming the cap rather than reporting probabilities: None. Every contraction, dense or doubled, checks its planned peak intermediate against PRISM_MAX_TN_PEAK_QUBITS before allocating.

Measurement and reset do not contract to the dense state: the outcome draws from the single-qubit reduced density matrix and the renormalizing projector is absorbed into the tensor holding the measured qubit's output leg, so the network keeps its deferred form, mid-circuit measurement carries no width ceiling, and the tensor count does not grow across measurements.

Two further queries stay off the dense route by contracting the network against its conjugate. The bra copy's legs are shifted clear of the ket index space, and each qubit's boundary is either closed against its twin, which is a trace, or joined through an operator. A one-qubit reduced density matrix leaves that qubit's ket and bra indices open and returns a 2x2; a Pauli expectation joins every non-identity factor through its operator and contracts to a scalar. Both follow the doubled network's cost rather than the qubit count. An identity factor is a closed leg rather than an appended tensor, so a weight-k observable adds k tensors and not n.

Nothing about the planner changed: an index is open when exactly one tensor holds it, which the greedy ordering already carries through to its result.

The reduced density matrix is the half of general-noise support the backend was missing, so the trajectory engine now runs amplitude damping, phase damping, thermal relaxation, and custom Kraus channels here under explicit dispatch.

Factored

Dynamic split-state simulation. Starts with n independent 1-qubit states, merges via tensor product only when 2q gates bridge groups. Parallel kernels match statevector patterns for sub-states ≥14 qubits. Selected when subsystem decomposition detects partial independence.

Density Matrix

Exact mixed-state evolution. Stores the full density operator rho for n qubits as a 4^n Complex64 buffer laid out row-major: index (r << n) | c holds ⟨r|rho|c⟩. That layout is isomorphic to a 2n-qubit statevector whose high n qubits index the ket (row) and low n qubits index the bra (column), so gate application reuses the statevector kernels. A unitary U on the ket register gives the left product U rho; the right product rho U^dagger takes the gate's conjugate form on the bra register where one exists, and otherwise conjugates the buffer around the pass. So U rho U^dagger costs two statevector passes, plus two conjugations only for the variants with no conjugate form. Rzz is the exception that carries gate math of its own: both factors are diagonal, so the ket and bra phases cancel wherever the two registers agree on the target pair's parity, and the sandwich collapses to a single pass over a combined table.

Memory is 16 * 4^n bytes, so the ceiling is about 14 qubits on a 16 GiB host and 15 on 32 GiB (PRISM_MAX_DM_QUBITS moves it within the statevector budget). With a device attached (BackendKind::DensityMatrixGpu, or with_gpu on the backend) the buffer lives in VRAM, budgeted against free memory at init: an 11 GiB card holds 13 qubits, so the device lifts the ceiling by about one qubit and is a throughput arm rather than a width arm. On the device the unitary half still runs the dense statevector kernels over the embedded buffer, while the channels, projection, reset, and the diagonal and Pauli readouts have kernels of their own. Both kinds are explicit-dispatch only; Auto and AutoGpu never select them, and the device kind has no host fallback.

Selecting it with a noise model attached is the exact route for every Simulate terminal except the adjoint gradient: the mixture is evolved once and observables, marginals, probabilities, and shots all read that one evolution, and the parameter-shift gradient evaluates the mixture once per shifted angle. Readout error is the one part of a model that no evolution holds, so run and marginals reject a model carrying it and point at sample_counts. The adjoint stays excluded because it backpropagates against a pure state and a channel has no reverse evolution to walk. See Noise across the terminals for what that route accepts and what stays on trajectory averaging.

Pauli Path

Heisenberg propagation of an observable through a noisy circuit, as a weighted sum of Pauli strings. It holds no state of any kind, so it is not a Backend: like SPP and SPD it is an engine the dispatcher reaches directly, and it serves expectation_values and observable_expectation (and the parameter-shift gradient built on them) and nothing else. run, shots, marginals, and probabilities are rejected naming the two terminals it does serve.

The observable starts as one term and propagates backward. Clifford gates conjugate it one term at a time. Rz and Rzz split every anticommuting term into a cos branch and a -i sin branch, which is what grows the sum. A noise channel scales each term by the channel's action on the Pauli letters it touches, which is what shrinks it. The engine is polynomial exactly where the shrinking wins, so the useful regime is a circuit whose noise rate outpaces its density of non-Clifford rotations, and the term count at a given width is the thing to watch rather than the width itself.

Observable weight dominates that term count, ahead of both width and depth. On a two-layer hardware-efficient ansatz under 1% depolarizing, Z on one qubit holds the sum at 11 terms from 20 qubits to 100, while Z on two adjacent qubits reaches a 16384-term budget by 30 qubits and a full-width Z chain reaches it at every width. term_count_is_width_independent_at_unit_weight pins both halves.

Channels enter through their adjoint on the Pauli basis, not through a twirl, so nothing is approximated at the channel. A unital Pauli channel (Pauli, Depolarizing, PhaseDamping, TwoQubitDepolarizing) scales each letter by one eigenvalue. AmplitudeDamping and ThermalRelaxation are not unital: their adjoint sends Z to (1 - gamma) Z + gamma I, and the sum carries that identity branch as a second term. Custom Kraus, Kraus2q, and readout error have no Pauli-basis form and are rejected naming the density matrix.

With max_terms = 0 the run is exact and errors at the shared term ceiling rather than truncating silently. With a budget set, terms below epsilon are dropped once the sum exceeds it, and the discarded coefficient mass bounds the error: every channel and every Clifford conjugation is a contraction in the Pauli 1-norm, so a dropped term contributes at most its own magnitude to the terminal value. A run that truncated nothing reports itself exact whatever budget it was given.

State diagnostics

Simulate::reduced_density_matrix, Simulate::entanglement_entropy and Simulate::overlap read the output state once the circuit has been applied, so all three require a unitary circuit: a measurement, reset or conditional leaves one seeded branch of several, not the state the diagnostic is defined on. Each resolves to a single backend, as the native expectation path does, and asks that backend for the answer in its own representation; an explicitly selected backend with no kernel for one of them reports BackendUnsupported naming itself and the diagnostic rather than falling back to a dense export.

Under BackendKind::Auto the route is the dispatcher's choice, not the caller's, so a resolved backend that cannot answer is replaced by the statevector while the circuit fits its cap. A partially independent circuit routed to the factored backend and a sparse-friendly one both read their entropy that way, and the same circuits decline when the backend is named explicitly. A Clifford circuit keeps its tableau, which answers the entropy and the marginal without expanding anything.

The entropy is the von Neumann entropy in nats, so a Bell pair reads ln 2, and the Schmidt values come back descending with their squares summing to one whatever norm the representation carried. A stabilizer cut of rank r has 2^r equal weights, so its tableau reads the rank off one elimination and builds the list from it; past the dense export cap those values no longer fit while the rank still does, and the entropy comes back alone with EntropyResult::schmidt_values at None. That is the width where the old fallback to the statevector could not answer at all. The reduced density matrix is row major with side 2^k and trace one, and its 4^k entries are priced as a 2k-qubit statevector against the dense export cap. A noise model sends the marginal and the entropy to the density matrix, which answers the marginal of the exact mixture and declines the entropy.

Simulate::overlap takes a second seeded builder, so each side carries its own backend, seed and start state, and the two circuits must declare the same width. The result is the modulus squared of the inner product over the two normalized states; the amplitude is not reported, since a tableau keeps no global phase and every MPS truncation moves one. Normalization divides by both norms on every route, so an unnormalized chain answers the same as a unit-norm state. Two states in the same representation take the native route at any width, and every other pair is served by a dense export of both, which reaches exactly as far as the export cap does. A noise model on either side is rejected, since the fidelity of two mixtures is a different computation.

What a backend reports about its own result

Three Backend methods carry provenance onto every result: resolved names the engine, exactness says whether its representation can discard state weight and how much this run discarded, and placement says whether the state lived on the device. All three have defaults, so an out-of-tree backend compiles unchanged and is named by Backend::name.

These are reports, not predictions. exactness is read after the circuit has been applied, so the MPS bound reflects the singular values this run actually discarded, and placement reflects where the amplitudes ended up after any device fallback. The MPS accumulates discarded weight per SVD and returns 1 - total as a fidelity lower bound; the sum is over relative discarded weights, so the bound is conservative.

The decomposed route runs one backend per independent block and merges: its exactness is the weakest of the parts, its fidelity bound is the product, and its placement is Device only when every block was. Per-shot routes evolve one state per shot and keep the weakest claim across them, with the bound a minimum rather than a product.

The GPU backend is documented as a user guide. The distributed statevector backend is covered in the Capability and Support Matrix.

Compiled Samplers

For multi-shot sampling without materializing the full statevector on every shot.

Noiseless compiled sampler (src/sim/compiled/)

Backward path (compile_measurements): Propagates Pauli Z observables backward through the circuit. Each measurement qubit becomes a row in a GF(2) parity matrix M. Clifford gates conjugate Pauli strings in O(1). The resulting M encodes which input qubits each measurement depends on.

Forward path (compile_forward): Tracks stabilizer generator dependencies forward through the circuit. Produces the same parity matrix via dependency tracking.

Sampling: Random bits for independent generators, then XOR-cascade through the parity matrix. Multiple dispatch tiers:

ShotAccumulator trait: Pluggable result collection.

PackedShots raw format: PackedShots::RAW_FORMAT_VERSION is the replay contract for raw_data() and into_data(). Version 1 stores little-endian bit order within each u64. ShotMajor stores one row per shot with m_words() = ceil(num_measurements / 64). MeasMajor stores one row per measurement with s_words() = ceil(num_shots / 64). The checked try_from_shot_major and try_from_meas_major constructors reject shape mismatches and non-zero semantic padding. Histograms, marginals, parity rows, and accumulators mask only semantic padding: measurement-tail bits in shot-major data and shot-tail bits in measurement-major data.

Detector sampler (compile_detector_sampler): Compiles Clifford circuits with measurement and reset reuse into the same packed measurement sampler, then derives detector and observable records as packed parity rows over measurement record indices. Reset reuse is represented by fresh qubit aliases, so repeated syndrome extraction avoids per-shot tableau replay. The sampler can return packed measurements, packed detectors, packed observables, detector counts, or feed packed detector chunks into any ShotAccumulator.

Native backend sampling (Backend::sample_basis_states)

The compiled samplers above cover Clifford circuits. Everything else used to funnel through Backend::probabilities(), a dense 2^n allocation, and sample from that, which put a hard qubit ceiling on shots for backends whose state is polynomial.

Two Backend hooks lift it. supports_native_sampling declares that a backend draws outcomes from its own representation, and sample_basis_states(num_shots, seed) returns packed per-qubit outcomes as BasisSamples (ceil(n / 64) words per shot). Seeding is from the argument, not the backend's RNG, so a shot request replays exactly, and the call does not collapse the state. It takes &mut self because a backend may first have to reorganize its own storage: the distributed backend restores its qubit map, a collective, so that each rank owns a contiguous slice in circuit order. supports_pauli_expectation and pauli_expectations(observables) are the observable-side pair, normalization independent so a truncated MPS is divided by ⟨ψ|ψ⟩ rather than assumed unit.

run_shots_with picks the native path through try_native_terminal_backend, which requires the route to land on a single backend and probes the capability before init, so a backend without one costs an allocation and nothing else. run_counts_with needs no separate path: its tail is run_shots_with(..).counts().

The product state is the one backend taken past subsystem decomposition. It already stores one factor per qubit, so splitting a non-entangling circuit into independent blocks pays a backend, a partition, and a merge per block to rebuild what one native draw reads off the state, and past 64 qubits the merged block distribution has no representation at all. Every other backend keeps the block split it had before, which only_the_product_state_takes_the_native_sampler_past_decomposition (src/sim/mod.rs) pins from both sides.

MPS records each bit against the logical qubit currently hosted at a site rather than the site index, so a layout permuted by SWAP routing needs no canonicalization pass. tests/native_sampling.rs pins that case; the exact check that the conditional decomposition reproduces the dense vector to 1e-12 lives in mps_conditional_path_probabilities_match_the_dense_vector (src/backend/mps.rs), and the corpus-wide comparison is the query matrix in tests/conformance_matrix.rs.

The tensor network answers below its dense ceiling from one contraction of the full distribution, and past it samples qubit by qubit: each bit draws from the conditioned single-qubit marginal, contracted on the doubled network, and the outcome projector is absorbed before the next qubit's marginal. A sweep shot costs one doubled contraction per qubit with peaks set by treewidth rather than 2^n, which is what carries sampling past the ceiling; the measured crossover that put the dense arm below it is recorded in the backend's module docstring.

Weighted observables and commuting-set grouping (src/sim/observable.rs)

PauliObservable carries a weighted Pauli sum H = Σ c_k P_k in the same (f64, Vec<PauliTerm>) term shape the gradient surface takes. Terms are kept canonical (factors sorted by qubit, identical strings merged) and the qubit-wise-commuting grouping is computed by greedy first-fit-decreasing coloring, cached inside the observable, and invalidated on mutation. Two strings qubit-wise commute when every shared qubit carries the same axis, so a group has one well-defined axis per qubit it touches.

observable_expectation computes Var(H_g) = ⟨H_g²⟩ - ⟨H_g⟩² per group on the statevector family, with the mean and most group variances served by one shared batched traversal (pauli_expectations_from_masks). Two strings in a QWC group multiply phase-free: shared qubits carry equal axes and cancel to identity, so each P_i P_j is just another Pauli string, and a small group's ⟨H_g²⟩ expands into pairwise product masks appended to the same traversal that serves the term means. A group past the pair budget (MAX_PAIR_MASKS_PER_GROUP, set where the quadratic expansion would cost more than a state sweep) takes a dedicated single-pass moment accumulation instead: on the state as run when the group is Z-only, otherwise on a copy rotated by H on X-assigned qubits and Sdg then H on Y-assigned qubits, after which members are plus-sign Z strings, and h(j) accumulates per element before squaring so both moments come from the one pass.

A per-group state pass is the fallback rather than the default because it loses at molecular shapes: on the 2000-string Jordan-Wigner bench fixture the grouping yields about 850 groups of mean size 2.3, and an engine paying copy, rotate, and sweep per group measured 8.3-8.5x slower than the ungrouped batched traversal. That traversal already made the mean one pass regardless of grouping, so grouping buys no mean throughput; what it buys is the variance, priced at the pair-mask expansion.

The reported variance is Σ_g Var(H_g): the variance of a grouped measurement estimate drawing one shot per group, and the input to shot allocation via the per-group vector. It excludes cross-group covariances, so it equals Var(H) of the full operator only when a single group covers every term; including them would cost O(M²) Pauli products at molecular term counts. Routes without the grouped evaluator (Clifford/SPD, the per-backend native paths, a run with a noise model or start state) report the weighted mean with no variance.

Grouping cost is O(M · G) word operations for M terms and G groups, sub-millisecond at thousands of terms; stronger colorings were declined because fewer groups would shrink neither the mean's single traversal nor the pair expansion, which scales with group size rather than group count.

Noisy compiled sampler (src/sim/noise.rs)

Backward Pauli propagation through circuit + noise sensitivity analysis. Each noise location gets an X-flip and Z-flip sensitivity row. During sampling, Bernoulli coin flips determine which noise channels fire, then XOR the sensitivity rows into the sample.

NoiseModel: per-instruction noise events. Pauli and depolarizing channels are supported by every noisy engine, and two-qubit depolarizing by every one but the homological sampler. Amplitude damping, phase damping, thermal relaxation, and one- and two-qubit custom Kraus operators require the trajectory engine. Readout error is separate: it acts on the measurement record rather than the state, so the route is chosen on the channels alone and each engine applies readout itself.

NoiseBuilder (src/sim/noise_builder.rs) compiles declarative rules into that same per-instruction vector: per-gate-type and per-qubit rates, idle decoherence against circuit layers, crosstalk through a coupling map, coherent over-rotation proportional to a rotation gate's own angle, reset error, pre-measurement error, and per-bit readout. Every rule is evaluated once at build time, so nothing it expresses reaches a per-shot or per-instruction loop.

Two-qubit Kraus operators are indexed K[t][t'] with t = 2*bit(q0) + bit(q1), the packing Gate::matrix_4x4 uses. The density matrix compiles the set into a 16x16 block superoperator (apply_2q_kraus, which apply_2q_depolarizing lowers onto). The trajectory engine draws a branch from Tr(Kdagger K rho) over Backend::reduced_density_matrix_2q and applies the normalized operator as a Fused2q. Only the host statevector implements that reduction, and run_shots_with_noise checks Backend::supports_two_qubit_kraus before the first shot, so an Auto route that picked another backend is named at dispatch rather than part way through a trajectory.

Noisy engine routing and the observable-result contract

Noisy sampling is reachable through simulate(...).noise(...).shots(n) / .sample_counts(n); the other builder terminals reject an inline noise model. run_shots_with_noise (src/sim/mod.rs) routes Pauli-only noise on Clifford circuits with terminal measurements (no resets, no classical conditionals) to the compiled family when the backend is Auto or stabilizer-family; every other accepted combination runs the trajectory engine over the resolved backend. Within the compiled family, run_shots_noisy (src/sim/noise.rs) picks one engine per call:

The last three, and the noiseless compiled sampler, all resolve to ResolvedBackend::CompiledStabilizer. RunMetadata::engine names which one ran, so a test pins a route off the result rather than off the predicates that picked it.

The trajectory engine (src/sim/trajectory.rs) covers everything the compiled family rejects: non-Pauli channels, mid-circuit measurement, reset, classical conditionals, and non-Clifford gates, at per-shot state evolution cost. Distributed backends reject noisy sampling entirely; per-shot trajectories cannot keep rank collectives in lockstep.

A two-qubit channel is sampled as one joint draw over its 15 non-identity Pauli products, never as two single-qubit draws: the true probability that both letters move is 0.6p, where independent draws give (0.8p)^2. The compiled sampler stores the propagated X and Z components of both targets, four rows to an event, in a table parallel to the single-qubit one, and runs them as a second pass so the single-qubit rows keep the positions the flip LUT and the device buffers were built against. Component rows over precomputed branches is a storage decision: 15 rows an event against 4. A pair whose targets fall in different subsystem blocks takes the monolithic compile, since the block-filtered one holds only one block's propagated masks at a time.

A non-empty pair table keeps noise application on the host: the device noise kernel derives each thread's bit from (seed, event, batch) for the one measurement row it owns and would draw the two qubits independently. The noiseless parity sample and the bit transpose still run on the device, so what the pair table costs is the fused device noise kernel and the counts and marginals reductions that sit behind it, not GPU sampling as a whole.

The frame and compiled samplers apply readout error to the packed measurement record, after the reference outcomes are folded in so that a set bit is a measured one. Records are walked one at a time and thinned at max(p01, p10), then a candidate is accepted at the rate its live bit selects, which is what asymmetric rates cost over a single flip mask. Brute-force replay instead flips the unpacked record of each shot on that shot's own stream, as the trajectory engine does.

Every entry point that draws shots calls NoiseModel::validate_for against the circuit before allocating state: one event slot per instruction, channel parameters in range, distinct targets on a two-qubit channel, and every target inside the register. Bounds cannot be checked from the model alone, and a target outside the register reaches kernels that index amplitudes without one. The analytic noisy_marginals_analytical is the exception, having no per-shot state to allocate.

A guarded region (Instruction::Region) is rejected whenever the model carries at least one quantum event: slots are indexed per top-level instruction, so a region body has none and would run noiselessly. A readout-only model has nothing to lose there and is accepted. Reaching noise inside a region body needs the event stream keyed by something other than a top-level index, which the compiled sampler, the homological builder, and the density-matrix evolution all walk today.

Custom Kraus sets must be trace preserving, sum_k Kdagger_k K_k = I to 1e-9. The exact route applies a declared set literally while the trajectory route normalizes its branch probabilities, so a set that is not trace preserving would mean two different things depending on which engine ran it.

ThermalRelaxation { t1, t2, gate_time } is amplitude damping composed with pure dephasing on both routes, at rates chosen so populations decay as exp(-gate_time/t1) and coherences as exp(-gate_time/t2). A mixture of reset and Z reproduces the population decay but reaches the coherence decay only for t2 <= t1, and needs a negative dephasing probability above it.

All engines sample from the same measurement-record distribution for the noise models and circuits they accept. The equivalence is statistical, not shot-for-shot: engines consume independent RNG streams, so the same seed produces different shots with matching observable statistics (marginals, correlators, histograms). Cross-engine tests pin every engine to the analytic marginals from noisy_marginals_analytical and to each other's correlator statistics: pauli_engines_share_observable_statistics (src/sim/noise.rs), trajectory_pauli_matches_brute_force (src/sim/trajectory.rs), and the channel-level analytic checks in tests/trajectory_correctness.rs.

GPU reductions (gpu feature): with a context attached via with_gpu, the noisy compiled sampler can sample, apply noise, and reduce counts or marginals on the device. Device noise masks come from a device-seeded RNG stream, so GPU output matches CPU output statistically, not bit for bit. On-device counts are limited to 512 measurements (8 packed words); larger circuits fall back to the CPU reduction. Golden test: noisy_compiled_gpu_reductions_match_cpu_statistics (tests/golden_gpu.rs).

Homological sampler (src/sim/homological.rs)

ErrorChainComplex: GF(2) chain complex over the circuit's noise locations. Computes the kernel (null space) of the boundary map to identify error cycles that are undetectable by syndrome measurements. HomologicalSampler uses this for sampling with topological error correction awareness.

noisy_marginals_analytical: Closed-form marginal computation using the parity matrix and noise rates. Avoids Monte Carlo sampling entirely.

See the Noise and QEC guide for how these fit together in practice.

Native QEC Program IR

QecProgram in src/qec/mod.rs is a measurement-record IR for QEC workloads that need detectors, logical observables, postselection, expectation metadata, and Pauli-noise annotations before sampler lowering. It is separate from Circuit so measurement-record programs do not need to fit final-measurement OpenQASM semantics.

QecOp stores gates, basis measurements, MPP-style Pauli-product measurements, resets, detector rows, observable includes, expectation-value metadata, postselection predicates, feed-forward corrections, noise annotations, and tick separators. Record references can be absolute indices or rec[-k] style lookbacks. Construction validates qubit bounds, gate arity, finite coordinates and coefficients, finite probabilities, and measurement-record scope. Detector, observable, and postselection rows can be resolved to absolute measurement indices for later compilation into packed samplers.

Feed-forward

QecOp::Feedforward conditions a body on the parity of a record list against an expected value, the shape a detector already has. It is a consumer of the guarded-region contract rather than a second conditional mechanism: the QEC record space and the classical bit vector are one address space, because the reference runner writes record i to classical bit i, so a resolved record index is directly the bit a ClassicalCondition::Parity reads. The op executes as the guarded instruction circuit::guarded picks for its body: a Region through Backend::apply_region, or a Conditional when the body is one gate.

The body admits gates and resets only. Detectors and observables index the record space absolutely, so a measurement whose execution depended on a record would make every later index depend on the shot.

The compiled QEC sampler evaluates a static affine map from random bits to outcomes, which a record-conditioned branch makes depend on the sample. It rejects by name and points at run_qec_program_reference, as do the deferred lowering and the density-matrix estimator. The detector-error-model derivation inherits the deferred lowering's rejection rather than carrying its own. The op is built through QecProgram::feedforward; the native text format does not spell it.

Parsing

parse_qec_program and QecProgram::from_text parse the native QEC text subset used by current benchmark planning: H, S, S_DAG, T, T_DAG, CX, CZ, R/RX/RY, M/MX/MY, MR variants, MPP, DETECTOR, OBSERVABLE_INCLUDE, POSTSELECT, EXP_VAL, Pauli-noise instructions, TICK, QUBIT_COORDS, SHIFT_COORDS, and flattened REPEAT blocks. The parser resolves rec[-k] references while building the program. Numeric arguments on basis measurements, such as M(0.001), lower to pre-measurement Pauli flips that affect the measurement record.

Lowering

compile_qec_program_rows lowers basis measurements and MPP records into the same packed X/Z Pauli row representation used by the compiled sampler internals. It also carries detector, observable, and postselection rows forward as absolute measurement-record indices. Detector, observable, and postselection projection uses PackedShots::parity_rows, so the QEC layer reuses the existing packed parity engine instead of maintaining a second one. This is a sampler-lowering artifact, not an execution engine. Gate, reset, and noise execution lives in run_qec_program; EXP_VAL has no packed-row representation, so the row compiler rejects it and run_qec_program routes such programs to the estimator paths described below.

Execution

run_qec_program lowers Clifford-compatible programs into the packed compiled sampler, compiles Pauli-noise annotations into sensitivity rows XORed into the records, and routes EXP_VAL programs to estimator paths. run_qec_program_reference is the per-shot state-vector correctness oracle. The runner routing, the compiled and noisy sampling paths, the circuit lowerings, and the result shape are covered in QEC program execution.

Expectation values

EXP_VAL(c) P1*...*Pk estimates c * <P> for the Pauli product P in the program's final state and returns one QecObservableEstimate per op, in op order, in QecSampleResult::expectation_values. The placement rules, the estimator paths, and the analytical strategy ladder are defined in QEC program execution.

QEC Program Execution

This page covers the execution architecture behind native QEC programs: the runner routing, the compiled row machinery, the circuit lowerings, the noisy data flow, the expectation-value estimator paths, and the result shape. The data model and text format are defined in the native QEC program IR page; workflow examples live in the noise and QEC guide. The public entry points are run_qec_program, run_qec_program_reference, run_qec_program_with_strategy, run_qec_program_spd_rerouted, and compile_qec_program_rows.

Program, ops, and options

QecProgram holds a qubit count, an ordered op list, and a QecOptions value. The typed builders (push_gate, measure, measure_pauli_product, reset, detector, observable_include, expectation_value, postselect, noise) validate each op as it is appended; measure and measure_pauli_product return the new measurement-record index, and detector returns the detector index. Validation covers gate arity, qubit bounds, finite coordinates, coefficients, and probabilities, record-reference scope, duplicate qubits in a Pauli product, and DEPOLARIZE2 target pairing.

Record references are QecRecordRef::Absolute indices or QecRecordRef::Lookback distances (distance 1 is the most recent record, distance 0 is rejected), resolved against the records that exist when the referencing op is appended. The text parser resolves rec[-k] references to absolute indices while parsing.

Runner routing

run_qec_program picks one of six paths:

flowchart TD
    A[run_qec_program] --> B{EXP_VAL ops present?}
    B -- yes --> C[validate placement rules]
    C --> D{active noise?}
    D -- yes --> M{density-matrix eligible?}
    M -- yes --> X[density-matrix estimator]
    M -- no --> R[reference runner]
    X -- lowering or oracle rejects --> R
    D -- no --> E{detectors present?}
    E -- yes --> S[two-run split]
    E -- no --> L[analytical Auto ladder]
    S -- either half fails --> R
    B -- no --> F{active noise?}
    F -- yes --> N[noisy compiled sampler]
    F -- no --> P[clean compiled sampler]

Programs containing EXP_VAL ops are routed instead of packed-sampled. The placement rules are validated first, then active noise sends the program to the density-matrix estimator or run_qec_program_reference, detectors send it to the two-run split, and the remaining noiseless case runs the analytical Auto ladder described under Expectation values.

Programs without EXP_VAL ops take the packed compiled path. Validation rejects non-Clifford gates and reports whether active noise is present. Programs with no measurements return an empty record buffer without compiling a sampler. Noisy Clifford programs compile a noise-aware sampler; clean Clifford programs lower to a circuit and compile a detector sampler. Both honor QecOptions::chunk_size: sampling proceeds in batches of at most chunk_size shots, and detector, observable, postselection, and logical-error accounting runs per chunk, so peak memory stays at one chunk of measurement records when raw measurements are not kept.

The two-run split serves noiseless EXP_VAL programs with detectors. Detectors are record metadata with no effect on the state, so the halves compose: the packed sampler runs the program without its EXP_VAL ops (real sampled measurement, detector, and observable records), and the analytical ladder runs the program without its detectors (exact estimates attached to the sampled result). When either half cannot run, for example non-Clifford gates on the packed half, the whole program falls back to the reference runner.

The density-matrix estimator serves noisy EXP_VAL programs whose noisy ensemble the mixed state can carry in full, where Tr(rho P) is the exact value the reference runner approximates by averaging per-shot statevector expectations. Eligibility:

• No measurement records. M and MPP collapse the state per shot and feed the measurement, detector, and observable rows of the result; the mixed state holds no record stream, so those programs need real sampling.
• No postselection predicate, which would condition the estimate on an accepted subensemble that Tr(rho P) does not express.
• Width within the density-matrix cap (PRISM_MAX_DM_QUBITS, default half the statevector cap), since the backend stores 4^n amplitudes.

R is eligible. Both paths implement the reset channel rho -> |0><0| (x) tr_q rho per the reset contract on the backends page: the density matrix applies it directly, and the reference runner samples one trajectory of it per shot, so the shot mean still converges to Tr(rho P).

Gates and channels the density-matrix path rejects surface as an error from the lowering or the oracle, which also falls back to the reference runner. The lowering maps gates one to one, expands a basis reset into Reset plus its Z-to-basis rotation, and turns each Pauli-noise annotation into NoiseModel events on the instruction it follows, applied through the backend's exact one-qubit Kraus and two-qubit depolarizing channels.

run_qec_program_reference is the correctness oracle: one statevector simulation per shot, O(shots * 2^n). It executes ops in order, samples Pauli noise stochastically from a dedicated RNG stream derived from the seed, lowers MPP onto one scratch qubit at index num_qubits, and evaluates postselection parities per shot.

Compiled rows

compile_qec_program_rows is the public sampler-row primitive. It lowers basis measurements and MPP ops into QecCompiledRows: one packed X/Z Pauli row per measurement record (bitmask words over the qubits), plus detector, observable, and postselection rows carried as absolute record indices with their expected values. Parity projection delegates to PackedShots::parity_rows, so the QEC layer reuses the packed parity engine of the compiled sampler. The row compiler is a lowering artifact, not an execution path: it rejects programs containing gates, resets, active noise, or EXP_VAL.

The clean compiled path builds its sampler from the lowered Clifford circuit instead. The compiled sampler backward-propagates each measurement observable through the circuit into a Pauli sensitivity row, reduces the rows by Gaussian elimination into a set of independent flip rows plus a deterministic reference outcome, and samples shots by XORing a random subset of flip rows into the reference bits. Detector and observable rows are applied afterward as parities over the sampled records. The propagation and sampling machinery is described on the compiled samplers page.

Lowering

Two lowerings turn a QEC program into a circuit the samplers accept, sharing the same helpers for basis rotations (append_basis_to_z_rotation and its inverse), MPP parity accumulation (append_mpp_parity_rotations: rotate each term into the Z basis, accumulate parity on a scratch qubit via CX, undo the rotations, measure the scratch), and a final record-count check that the lowering emitted exactly num_measurements records.

The clean Clifford lowering emits measurements in place: rotate the measured qubit into the Z basis, measure into the next record, rotate back. MPP uses one scratch qubit at index num_qubits, reset between uses. Resets emit a reset followed by the basis rotation.

The deferred lowering backs the noisy path and the analytical estimator prelude. A reset assigns the program qubit a fresh circuit alias, so every measurement can be deferred to a terminal record of the lowered circuit; reusing a measured qubit without a reset is rejected. The final alias of each program qubit is recorded so terminal EXP_VAL terms translate onto the lowered circuit. Noise ops are not applied to the circuit; they are recorded as positioned events for sensitivity compilation.

Noisy data flow

Noise never executes during compiled sampling; it compiles into record flips.

At compile time, the deferred lowering produces the noiseless circuit plus the positioned noise events. Backward Pauli propagation then walks the circuit in reverse, carrying each measurement's observable to each event position, and converts every event into sensitivity rows: for each noise branch, the set of measurement records whose propagated Pauli anti-commutes with the injected error, packed as flip masks. The noiseless sampler is compiled on the deferred circuit.

At sample time, the noiseless records are sampled first, then each noise event stochastically XORs its branch flip masks into the shot-major record buffer from a noise RNG stream derived from the seed. Small-probability events skip between firing shots with geometric sampling; dense events (probability at or above 0.5, or fewer than 32 shots) iterate every shot. DEPOLARIZE2 precomputes the flip masks of all 15 non-identity two-qubit Pauli branches and picks one uniformly per firing.

Supported channels are X_ERROR, Z_ERROR, DEPOLARIZE1, and DEPOLARIZE2. Noise on an already-measured target is dropped (it can no longer affect any record), and a DEPOLARIZE2 pair with one measured target degrades to DEPOLARIZE1 at p * 0.8 on the survivor, preserving the marginal error rate. The reference runner instead applies the same channels stochastically to the per-shot state, and the density-matrix estimator applies them exactly: X_ERROR and Z_ERROR become one-axis Pauli channels, DEPOLARIZE1 a symmetric one-qubit depolarizing channel, and DEPOLARIZE2 a two-qubit depolarizing channel on each target pair. How the QEC noise path relates to the circuit-level noisy engines is covered by the noisy engine routing section of the compiled samplers page.

Detector error model export

QecProgram::detector_error_model derives a DetectorErrorModel: the set of independent error mechanisms implied by the program's noise annotations, detectors, and observables. Matching and belief-propagation decoders consume this model rather than raw detector samples.

The derivation reuses the noisy data flow above. The deferred lowering produces the noiseless circuit and the positioned noise events, and the same backward Pauli propagation supplies, at every event position, the set of measurement records each single-Pauli fault flips. Every annotation then expands into its fault branches (one per target for X_ERROR and Z_ERROR, three per target for DEPOLARIZE1, fifteen per pair for DEPOLARIZE2), and each branch's record mask is projected through the detector and observable rows to its symptom: the detectors and observables it flips.

Branches merge into mechanisms under two rules, at fault-site granularity (one target of a single-qubit annotation, or one target pair of DEPOLARIZE2):

• Branches at one fault site are mutually exclusive, so branches with the same symptom sum. A DEPOLARIZE1 site whose X and Y branches flip the same records yields one mechanism at exactly 2p/3.
• Distinct fault sites are independent, distinct targets of one annotation included, so mechanisms with the same symptom compose as p = p1(1-p2) + p2(1-p1).

Faults that flip no detector and no observable are omitted. Mechanisms keep program order (the position of the annotation that first produced each symptom). Mechanisms are independent in the model even where the underlying branches were exclusive, so model statistics agree with the sampler to second order in the branch probabilities; validation compares at tolerances, never exactly.

Consequences of the sampler semantics carry over unchanged: a measurement error argument (M(p)) is already a pre-measurement Pauli fault, so it appears as an ordinary mechanism on that record's detectors; noise on an already-measured qubit contributes nothing; and a DEPOLARIZE2 pair with one measured target enters as the exact DEPOLARIZE1(0.8p) marginal on the survivor. Non-Clifford gates and reuse of a measured qubit without reset are rejected, as on the compiled sampling path.

Hypergraph mechanisms (more than two detectors, as DEPOLARIZE2 produces) stay intact in the derived model; decoders that accept a check matrix consume them directly. A matching decoder needs the graphlike form, which is opt-in: decompose_graphlike returns a new model in which every hypergraph mechanism is replaced by two or more existing graphlike mechanisms whose non-empty detector sets partition its detectors and whose observable XOR matches it, with the mechanism's probability composed into every component. A DEPOLARIZE2 Y-Y branch, for example, splits into the same channel's single-qubit branches. Cross-component correlations are lost; single-detector marginals are unchanged, since each hyperedge composes into exactly one component containing any given detector. A mechanism with no such cover is an error naming its symptom; decomposition is never applied silently, and the writer serializes whichever model it is given.

Text format

DetectorErrorModel::to_text renders the model in the common detector error model text format that external matching and belief-propagation decoders read. The emitted subset:

error(<p>) D<i> ... L<j> ...
detector D<i>
detector(<c0>, <c1>, ...) D<i>
logical_observable L<j>

One error line per mechanism in mechanism order, carrying its probability and the D-prefixed detector indices and L-prefixed observable indices it flips, both ascending. One detector line per detector in index order, with the coordinates of the program's detector op when present. One logical_observable line per observable slot. Indices are zero-based and dense; probabilities print with enough digits to round-trip exactly. Flat models only: no repeat blocks, no coordinate shifts, no decomposition suggestions.

Example, one syndrome round of the three-qubit repetition memory under X_ERROR(0.05) on the data qubits:

error(0.05) D0
error(0.05) D0 D1
error(0.05) D1
detector D0
detector D1

In Python, QecProgram.detector_error_model() returns the model with probabilities() (float64), detector_matrix() and observable_matrix() (bool, detectors or observables by mechanisms), detector_coords(), and to_text(). The matrix triple feeds check-matrix decoder constructors directly, with no file in between.

Decoding

UnionFindDecoder closes the pipeline in-tool: sample, derive, decompose, decode, logical error rate, one API. The decoder family is union-find with peeling (Delfosse and Nickerson, arXiv:1709.06218), chosen for its almost-linear decode cost; a minimum-weight perfect matching decoder remains a possible second family behind the same model input if accuracy on hard workloads ever justifies its cost.

UnionFindDecoder::from_model compiles a graphlike model: detectors become vertices, a two-detector mechanism an internal edge, a one-detector mechanism a boundary edge, each weighted ln((1-p)/p) clamped at zero. A mechanism with more than two detectors is rejected with a pointer to decompose_graphlike. Mechanisms flipping no detector cannot enter the graph; their probability mass is a floor under the logical error rate of any decoder over the model. Mechanisms sharing one detector set collapse to the most probable of them, and the mass of single faults this misroutes is measured (not assumed away) by the enumeration tests in tests/qec_decoder.rs.

decode_packed maps a PackedShots of detector samples (either layout) to shot-major predicted observable flips, one bit per observable per shot. Per shot, clusters grow from the defects in the ln((1-p)/p) metric: each round adds the minimum slack over the active clusters' unsaturated incident edges, so at least one edge saturates per round; a cluster becomes inactive when its parity is even or it touches a boundary edge. Peeling then walks a spanning forest of each grown cluster, rooted at the boundary contact when one exists, and XORs the observable mask of every selected edge into the prediction. A component with odd parity and no boundary edge is impossible under the model and rejects the batch, naming the shot. Decoding is deterministic: no randomness anywhere, ties break by ascending edge index in mechanism order, and results are identical on the serial and Rayon paths. Shots are independent: under the parallel feature, batches of at least 1024 shots decode in parallel over 256-shot chunks with per-chunk scratch reuse and no per-shot allocation; a model with no observables decodes serially at any batch size.

Validation ties the decoder to the exact ML lookup rate: at distance 3 the full syndrome set is enumerable, so the tests compute the exact expected union-find failure rate alongside the exact ML rate and assert rate_ML <= rate_UF <= P(two or more faults) + measured single-fault misses, pin the analytic rate at 1e-12, and hold fixed-seed golden decode counts. On repetition memory at p=0.02 (20k shots, seed 42, 3 rounds) the decoded rate falls from distance 3 to distance 5 and both sit below the physical rate. Python exposes the same surface as Decoder(model) with decode(detectors) -> (shots, num_observables) over numpy bool arrays.

Expectation values

EXP_VAL(c) P1*...*Pk estimates c * <P> for the Pauli product P in the program's final state and returns one QecObservableEstimate per op, in op order, in QecSampleResult::expectation_values. Two placement rules make "final state" well defined on every path:

• Terminal placement: no gate, measurement, reset, or active noise may follow an EXP_VAL op. Detector, observable, postselection, and tick metadata may.
• Live qubits only: a term may not reference a qubit that was single-qubit-measured after its last reset. Under this rule the Pauli commutes with every prior measurement projector, so the sampled post-measurement average equals the measurement-stripped pure-state expectation the analytical strategies compute. MPP does not affect liveness: the deferred lowering measures a scratch alias, and the projected cross terms cancel exactly.

Estimator paths:

The reference runner precomputes Pauli masks per observable and evaluates c * <P> on the final statevector of every accepted shot, reporting the sample mean and unbiased sample variance with num_shots equal to the accepted count. Programs wider than 64 qubits (including the MPP scratch qubit) are rejected, since the mask reduction packs X and Z masks into 64-bit words. With noise annotations the per-shot trajectories average to the mixed-state expectation Tr(rho P).

The analytical strategies evaluate <0|U^dag P U|0> on the deferred lowering with the trailing measurements stripped, translating record observables into Z strings and EXP_VAL terms through the final qubit aliases. They reject detectors and active noise; run_qec_program composes those cases through the split and reference routes instead. QecTStrategy selects the path:

• Auto, the production ladder. Light-cone SPD runs first (truncation tolerance 1e-10, term cap 16384) and its result is accepted when every estimate reports no truncation; SPD encodes truncation as variance = total_discarded^2, and variance at or below 1e-12 counts as exact. Otherwise CAMPS runs (bond dimension cap 256), which hard-errors when SVD truncation discards weight above 1e-12 so the ladder falls through to the exact tensor-network scalar fallback. When all three fail, the combined error reports each stage's reason.
• Reference, the per-shot oracle.
• Spd and Camps run their stage directly.

CAMPS evaluates arbitrary X/Y/Z strings by conjugating each letter through the signed Clifford prefix (Y = i * X * Z composes the two inverse-tableau rows). Postselection composes on both paths: the reference runner averages over accepted shots, and the analytical combiner conditions via <O * Pi> / <Pi> evaluated over the projector subsets, capped at 12 postselection predicates; the projector Z strings live on measured aliases and so never overlap EXP_VAL terms.

run_qec_program_spd_rerouted accepts caller-supplied Z stabilizers per observable and evaluates each rerouted observable on the XOR-equivalent support with the smallest inverse light cone, verifying first (via SPD) that the substituted stabilizer product holds <S> = +1 in the lowered state; observables without a reroute evaluate on their original support. The path rejects EXP_VAL ops, postselection, and resets (resets relabel qubits, making stabilizer indices ambiguous).

Routing and placement are pinned by tests/qec_exp_val.rs; tests/qec_e2e_d3.rs checks distance-3 repetition-code fixtures against closed-form expectations on the compiled, analytical, and reference paths.

Result shape

Every runner returns a QecSampleResult. The type is #[non_exhaustive]: construct through new, new_with_total_shots, or empty, and match fields with a trailing ...

QecObservableEstimate carries mean, variance, and num_shots (the shots that contributed, excluding postselection rejections). Zero accepted shots yield {mean: 0.0, variance: 0.0, num_shots: 0}.

Analytical strategies have no per-shot stream, so they synthesize the packed observables records to match the logical_errors counts: the one-bits occupy positions [0, accepted_shots) and the remainder is inert padding. Derive rates with logical_error_rates (denominator accepted_shots); do not align these rows shot-for-shot with detector rows on the analytical path.

The compiled runner delivers shot-major PackedShots record buffers: one row of packed record bits per shot, record j at bit j % 64 of word j / 64 of that row. The Clifford+T sampler emits measurement-major records instead, so consumers branch on layout() or read through get_bit. Detector, observable, and postselection parities are XORs of the referenced measurement records. Summary statistics are available as survivor_rate, logical_error_rates, and their Wilson-interval variants.

Threading, SIMD, and Memory Layout

For which SIMD tiers and architectures each backend supports, see the Capability and Support Matrix.

Memory layout

Threading

Gate kernels have _par variants using par_chunks_mut for safe Rayon parallelism (behind the parallel feature flag):

• <14 qubits: Single-threaded. Thread-pool overhead exceeds computation.
• ≥14 qubits: Rayon parallel iterators with MIN_PAR_ELEMS = 4096 (64KB per task).

Thread pool defaults to all logical cores (HT helps at 24q+ by hiding memory latency). Overridable via RAYON_NUM_THREADS.

The default is the process-wide Rayon pool, sized on the first simulation call. An application that owns that pool can hand PRISM-Q a bounded one instead: ThreadPool::with_threads(n) builds it and install runs a closure on it, leaving the global pool unbuilt and unresized. Calls made outside install take the global path. Pool width is not free of consequences for results; see the determinism contract below.

SIMD

Complex64 maps to 128-bit SIMD naturally. Single-qubit gate kernels use PreparedGate1q with runtime CPU detection and tiered dispatch:

1. AVX2+FMA (256-bit): 2 complex pairs per iteration. Gated by MAX_AVX2_STATE for full-state passes (Skylake frequency throttling), but used freely within MultiFused L2 tiles where data is cache-resident.
2. FMA (128-bit): Default for larger states. 3-op complex multiply (permute + mul + fmaddsub).
3. BMI2: _pext_u64 for BatchPhase, BatchRzz, and DiagonalBatch LUT indexing. One BMI2 bit extraction replaces loops with repeated shifts and ORs.
4. Scalar fallback: No intrinsics. All SIMD functions have a #[cfg(not(target_arch = "x86_64"))] fallback.

Two key SIMD structs hoist matrix broadcast at construction time, avoiding per-element dispatch:

• PreparedGate1q: Broadcasts 2×2 matrix into SIMD registers. Methods: apply_full_sequential (full state), apply_tiled (cache-resident tile, no AVX2 throttle guard), apply_slice_pairs (MPS bond-dimension slices), apply_pair_ptr (Cu/Mcu parallel).
• PreparedGate2q: Broadcasts 4×4 matrix. Methods: apply_full (mask-based iteration), apply_tiled (cache-resident Multi2q tiles, AVX2 paired-group kernel when available), apply_group_ptr (4 scattered indices).

The 2q tiled AVX2 path processes paired k and k + 1 groups when the lower target qubit is above 0, which makes each row load contiguous. It falls back to the 128-bit FMA kernel for lo == 0 and when AVX2+FMA is unavailable. Set PRISM_NO_AVX2_2Q to compare against the 128-bit FMA path, or PRISM_NO_REORDER to disable disjoint Fused2q tier grouping for A/B timing.

Determinism

Reproducibility is a per-path contract, not a blanket guarantee. Gate application never reduces across tasks, so it is exactly reproducible; everything that sums floating-point values in parallel is reproducible to the last ulp only; the batched compiled sampler is reproducible only at a fixed thread count. tests/determinism.rs pins the dense unitary, terminal sampling, reduction, and compiled-sampler claims by running the same seeded circuits in scoped 1-thread and 4-thread pools; the trajectory, SPD, and stabilizer bullets stand on the mechanisms they state.

What every clause below turns on is thread count, not which pool supplied the threads. A caller-supplied ThreadPool therefore moves exactly the results a different RAYON_NUM_THREADS would: the bitwise claims survive it, the reduction bound holds across it, and compiled shot payloads change with it whenever the pool is narrower or wider than the one the comparison run used.

• Unitary evolution on the dense kernels: bitwise, at any thread count. Gate kernels partition the state into index-derived disjoint ranges (fixed chunk boundaries, index bijections for the SendPtr kernels) and write elementwise, so the schedule cannot reach any value. Pinned bitwise for representative fused circuits and a QFT.
• Terminal sampling on the dense route: bitwise for a given seed, at any thread count, when the measurement map covers the register directly. Shot thresholds are drawn and sorted up front from the seeded generator, the probability vector is an elementwise transform, and the CDF walk is sequential. A measurement map that compacts into fewer outcomes builds its histogram through a parallel reduction and moves to the ulp-stable class below.
• Noisy trajectory shots: bitwise for a given seed, at any thread count. Each shot's generator is seeded from the shot index, not the worker, and results are collected in shot order.
• Parallel reductions: stable to about 1e-12, not bitwise. Norms, measurement collapse probabilities, reduced density matrices, and expectation values sum deterministic per-chunk partials in Rayon's combine order, which varies with pool width and work stealing. A mid-circuit measurement compares a seeded draw against such a sum, so an outcome flip is possible in principle when the draw lands inside the ulp gap.
• Compiled (BTS) sampling: reproducible at a fixed thread count only. The batched sampler derives one RNG stream per worker and splits shots by rayon::current_num_threads(), so a different pool width yields a different, equally distributed shot set. Pin the width when byte-identical shot payloads matter across machines, through RAYON_NUM_THREADS on the global path or through the width passed to ThreadPool::with_threads on the scoped one. The GPU analogue is documented in the GPU guide.
• SPD analytic estimates: stable to about 1e-12 between runs. Hash-order term accumulation moves the last ulp even at a fixed thread count; tests carry that tolerance.
• Stabilizer tableau: bitwise, at any thread count. Row operations are integer and exactly associative. MPS carries no bitwise claim: truncation runs through faer's threaded SVD.

Error Model and Public API

Error model

Fallible public APIs return Result<T, PrismError>. Error variants:

Public API surface

Top-level re-exports from src/lib.rs. The full generated documentation is on docs.rs.

Simulation: simulate, Simulate, Unseeded, Seeded, run_on, run_on_state, run_qasm, run_expectation_values, run_observable_expectation, PauliObservable, ObservableExpectation, bitstring

State diagnostics: Simulate::reduced_density_matrix returns a ReducedDensityMatrix (row major, side 2^k, qubits[0] the lowest bit of the row index) with a purity method for Tr(rho^2); Simulate::entanglement_entropy returns an EntropyResult carrying the von Neumann entropy in nats and the descending Schmidt spectrum, which is None where the backend holds the entropy without the spectrum behind it; Simulate::overlap takes a second seeded builder over a circuit of the same width and returns an OverlapResult carrying the squared inner product and the provenance of both runs. All three require a unitary circuit, and under BackendKind::Auto a route that cannot answer falls back to the statevector. Which backends answer each is tabulated in Backends.

Gradients: run_expectation_gradient, run_expectation_gradient_shift, ExpectationGradient

Parameters and binding: Parameters, ParamLink, PreparedCircuit. One parameter model serves both consumers: the gradient path reads the links, and binding writes through them.

Compiled sampling: compile_measurements, compile_forward, compile_detector_sampler, compile_noisy, run_shots_compiled, run_shots_noisy, run_shots_homological, noisy_marginals_analytical, density_matrix_expectation_values; with the gpu feature: run_shots_compiled_with_gpu, DevicePackedShots

Native QEC: parse_qec_program, compile_qec_program_rows, run_qec_program, run_qec_program_reference, run_qec_program_with_strategy, run_qec_program_spd_rerouted, QecProgram, QecOp, QecOptions, QecSampleResult, QecBasis, QecPauli, QecRecordRef, QecNoise, QecMeasurementRow, QecCompiledRows, QecObservableEstimate, QecObservableReroute, QecTStrategy, DetectorErrorModel, ErrorMechanism, UnionFindDecoder

Clifford+T: run_stabilizer_rank, run_stabilizer_rank_approx, stabilizer_overlap_sq, stabilizer_inner_product, StabRankResult, run_spp, run_spp_observable, run_spd, run_spd_observable, run_spd_observable_light_cone, inverse_light_cone, PauliAxis, PauliTerm, SppResult, SppObservableResult, SpdResult, SpdObservableResult

Types: Circuit, CircuitBuilder, Instruction, ClassicalCondition, SvgOptions, TextOptions, Gate, GeneratorKind, BackendKind, RunOutcome, CountsResult, MarginalsResult, ReducedDensityMatrix, EntropyResult, Probabilities, FactoredBlock, ShotsResult, PrismError, Result, MultiFusedData, BatchPhaseData, McuData, Multi2qData

Backends: StatevectorBackend, StabilizerBackend, SparseBackend, MpsBackend, ProductStateBackend, TensorNetworkBackend, FactoredBackend, FactoredStabilizerBackend; with the distributed feature: DistributedStatevectorBackend, DistributedContext, RankComm, SerialComm; with the distributed-mpi feature: MpiComm

Threading: with the parallel feature: ThreadPool, a caller-supplied Rayon pool that simulation runs inside instead of sizing the process-wide one. See Threading, SIMD, and Memory Layout.

Accumulators: ShotAccumulator, HistogramAccumulator, MarginalsAccumulator, PauliExpectationAccumulator, CorrelatorAccumulator, NullAccumulator, PackedShots, ShotLayout, ParityStats

Data types: CompiledSampler, CompiledDetectorSampler, DetectorSampleBatch, NoisyCompiledSampler, NoiseChannel, NoiseEvent, NoiseModel, NoiseBuilder, GateFilter, ReadoutError, HomologicalSampler, ErrorChainComplex

Not re-exported at the root but part of the documented surface: the Backend trait and BasisSamples at prism_q::backend, the density matrix backend at prism_q::backend::density_matrix, and the accumulator chunk-size helpers (default_chunk_size, optimal_chunk_size) at prism_q::sim::compiled.

Benchmarks

Wall-clock simulation time for PRISM-Q on a fixed circuit suite built from the prism_q::circuits generators, pushed toward this machine's limits. Every number is reproducible with the command at the bottom of this page.

Setup

• Date: 2026-05-29
• CPU: Intel64 Family 6 Model 94 Stepping 3, GenuineIntel
• Threads available: 8
• PRISM-Q version: 0.16.0

Methodology

• Metric: median wall-clock over repeated runs after warmup, lower is better.
• Timed region is simulation only. Circuit construction happens once per circuit outside the timer.
• Timings use the full simulate().run() path, including the fusion pass and probability extraction.
• auto lets backend dispatch choose a specialized backend per circuit. The dense families (QFT, HEA, QV) are bounded by the statevector memory cap; GHZ is Clifford and runs into the thousands of qubits on the stabilizer backend. QV uses square depth, so its gate count grows with the qubit count and it is compute-bound earlier than the fixed-depth families.

GHZ

QFT

HEA

QV

Reproducing

cargo run --release --features parallel --example bench_suite

Times PRISM-Q on every circuit in the suite and rewrites this page.

Circuit Builders

Pre-built circuits for benchmarking and testing, in src/circuits.rs. Each returns a Circuit you can pass straight to simulate(&circuit) or any backend.

Example

#![allow(unused)]
fn main() {
use prism_q::circuits::qft_circuit;
use prism_q::simulate;

let circuit = qft_circuit(10);
let result = simulate(&circuit).seed(42).run().unwrap();
}

For hand-built circuits, use the CircuitBuilder fluent API instead.

The complete generated API documentation lives on docs.rs.

Architecture Glossary

This glossary defines key terms used throughout the PRISM-Q architecture and documentation.

Terms

Backend A simulation strategy implementing the Backend trait (e.g., statevector, stabilizer, MPS). Backends are swappable and selected automatically or explicitly per circuit.

Clifford A class of quantum gates (H, S, CNOT, etc.) that can be simulated efficiently on stabilizer tableaus in O(n²) time.

Count A frequency histogram of measured bitstrings across multiple shots, returned by simulate(circuit).seed(seed).sample_counts(shots).

Fusion The pre-simulation optimization pipeline that merges, cancels, and reorders gates to reduce execution cost.

Marginal Per-qubit probabilities of measuring |1⟩, extracted without computing the full multi-qubit distribution.

MPS (Matrix Product State) A tensor-network backend that stores the quantum state as a chain of rank-3 tensors with bounded bond dimension, trading exactness for polynomial memory.

Noise model A collection of probabilistic error channels (e.g., depolarizing noise) applied during noisy multi-shot simulation.

OpenQASM An assembly-like language for describing quantum circuits. PRISM-Q parses a practical subset of OpenQASM 3.0.

Shot A single execution of a quantum circuit from initialization to measurement, producing one bitstring sample.

State vector A flat vector of 2^n complex amplitudes representing the full quantum state, used by the statevector backend.

Stabilizer A representation of a Clifford state as a bit-packed tableau of Pauli generators, enabling simulation of thousands of qubits.

Tensor network A backend that represents gates as tensors and defers contraction until measurement or probability extraction.