qsim — a tiny, optimised quantum simulator

A complete state-vector quantum computer simulator in two forms:

Part	File	Size	Dependencies
Python engine	qsim.py	~400 lines	numpy only
Browser engine	web/engine.js	~130 lines	none
Interactive UI	web/index.html, web/app.js	~500 lines	none
Tests	test_qsim.py	19 tests	none beyond numpy
Benchmark	bench.py	—	—

Both engines implement the same algorithm and are cross-verified to agree to 2.2 × 10⁻¹⁶ (machine epsilon) on 10 reference circuits.

---

Quick start

pip install numpy
cd qsim
python test_qsim.py        # 19 tests, ~1 s
python bench.py            # performance table

from qsim import Circuit, ghz, qft, grover

# Bell pair
c = Circuit(2).h(0).cx(0, 1)
print(c.probs_dict())      # {'00': 0.5, '11': 0.5}
print(c.sample(1000))      # {'11': 508, '00': 492}

# 20-qubit GHZ state — 16 MiB, well under a second
g = ghz(20)
print(g)                   # <Circuit n=20 dim=1048576 mem=16384.0KiB ops=20>

# Grover search over 2^10 items
r = grover(10, marked=417)
print(r.probabilities().argmax(), r.probabilities().max())   # 417 0.9995

Open web/index.html in any browser for the interactive circuit builder — no build step, no server, no install.

---

1. What a state-vector simulator is

An n-qubit pure state is a unit vector in (\mathbb{C}^{2^n}):

[ |\psi\rangle = \sum_{i=0}^{2^n-1} \alpha_i,|i\rangle,\qquad \sum_i |\alpha_i|^2 = 1 ]

A gate is a unitary matrix. Applying it naively means building a (2^n \times 2^n) matrix and doing a matrix–vector product:

• memory: (2^{2n}) complex numbers — 256 MiB at just 12 qubits
• time: (O(4^n))

Nobody does this. The whole art of a fast simulator is applying gates without ever forming that matrix. That is what this code does.

Qubit ordering. Qubit 0 is the least-significant bit, so basis index i has qubit k equal to (i >> k) & 1, and the printed bitstring q_{n-1}…q_1 q_0 reads left-to-right from the highest qubit. This is the same little-endian convention Qiskit uses, so exported OpenQASM behaves identically there.

---

2. The core optimisation: reshape, don't multiply

A single-qubit gate on qubit q only ever mixes pairs of amplitudes whose indices differ in bit q. If we view the flat state array as a 3-D array

state.reshape(2**(n-q-1),  2,  2**q)
        ^ higher bits    ^ q  ^ lower bits

then the middle axis is qubit q, and the gate is just

v0, v1 = v[:, 0, :], v[:, 1, :]        # two strided views, no copy
v0, v1 = m00*v0 + m01*v1,  m10*v0 + m11*v1

• Time: (O(2^n)) — one pass over the state, fully vectorised by numpy.
• Memory: one (2^{n-1}) scratch buffer. No (2^n \times 2^n) operator.
• Any qubit: the same code works for every q; only the reshape changes.

Measured effect (bench.py, column "speedup" = naive ÷ optimised):

qubits	state vector	naive operator	optimised	naive	speedup
8	4 KiB	1 MiB	0.010 ms	1.70 ms	165×
10	16 KiB	16 MiB	0.011 ms	23.5 ms	2 079×
12	64 KiB	256 MiB	0.021 ms	228 ms	10 987×

The gap grows like (2^n); at 20 qubits the naive operator would need 16 TiB and is simply not representable.

---

3. Second optimisation: structure-aware gate kernels

Most useful gates are not dense 2×2 matrices. Everything funnels through one kernel, _mix(v0, v1, m), which branches on the matrix's structure:

Structure	Gates	Work done
Diagonal (m01 = m10 = 0)	Z, S, S†, T, T†, P, RZ, CZ, CP, CRZ	in-place scalar scaling — zero copies, zero adds
Anti-diagonal (m00 = m11 = 0)	X, Y, CNOT, CY, Toffoli	slab exchange — zero adds
Dense	H, SX, RX, RY, U	1 scratch buffer + 2 fused writes

This is not a micro-optimisation. At 20 qubits:

gate	before	after	gain
CNOT	19.6 ms	1.8 ms	11×
CZ	19.6 ms	0.6 ms	33×
T	17.9 ms	1.5 ms	12×
GHZ-20 (20 gates)	361 ms	53 ms	6.8×

---

4. Third optimisation: controlled gates without index arrays

The obvious way to do a controlled gate is to build a boolean mask over all (2^{n-1}) index pairs and gather/scatter. That allocates two int64 arrays and a bool array per gate — more traffic than the gate itself.

For the overwhelmingly common case of exactly one control, this code instead does a five-axis reshape that isolates both bits at once:

state.reshape(high, 2, mid, 2, low)     # axes 1 and 3 are the two qubits

Selecting control = 1 on its axis yields two writable views of the target bit, and _mix runs on them directly. No index arrays, no masks, no gather/scatter — pure strided views. Multi-control gates (Toffoli, MCX) fall back to the mask path, which is still a single pass.

SWAP skips gates entirely: it is one permutation of the state array computed with XOR arithmetic, not the usual three CNOTs.

---

5. Sampling and measurement

• sample(shots) — draws all shots in one multinomial call: (O(2^n + \text{shots})), not (O(\text{shots} \times n)). 65 536 shots on a 20-qubit state takes ~18 ms. The state is not collapsed, so you can keep simulating.
• measure(q) — a genuine projective measurement: computes (P(1)) from the same reshape, zeroes the rejected half, renormalises the survivor, in place.
• expectation_z(q), bloch(q) — read observables straight off the reshaped views, no copies. bloch returns the reduced-density-matrix vector, so its length shrinks below 1 exactly when the qubit is entangled with the rest.

---

6. Measured performance

Ordinary 2-vCPU container, numpy 2.5, complex128 (16 bytes/amplitude).

qubits	memory	1q gate	CNOT	GHZ build	100 k shots
12	0.06 MiB	0.019 ms	0.009 ms	0.17 ms	0.08 ms
16	1 MiB	0.42 ms	0.07 ms	1.9 ms	1.0 ms
20	16 MiB	17.9 ms	1.8 ms	53 ms	17.5 ms
22	64 MiB	115 ms	15.5 ms	354 ms	101 ms

Algorithms:

algorithm	qubits	gates	runtime	result
QFT	8	40	0.6 ms	matches numpy.fft to 1e-10
QFT	16	144	16.5 ms	—
Grover	8	658	5.4 ms	P(marked) = 0.986
Grover	10	1 560	16.0 ms	P(marked) = 0.9995
Grover	12	3 712	65.8 ms	P(marked) = 0.9999

Practical ceiling: memory is (16 \times 2^n) bytes, so 24 qubits = 256 MiB, 26 qubits = 1 GiB, 30 qubits = 16 GiB. The class caps n at 30. On typical hardware, 22–24 qubits is the comfortable working range.

---

7. Cost model — why this is cheap

Resource	This simulator
Install	numpy (browser version: nothing at all)
Runtime	your own CPU / your own browser tab
Server	none — the web build is three static files
Per-run cost	zero
Code to audit	~530 lines total across both engines

There is no cloud backend, no queue, no account, no API key. The web app runs entirely in the tab; nothing leaves your machine.

---

8. API reference (qsim.py)

Construction

Circuit(n_qubits, seed=None)   # starts in |00…0>, n ≤ 30
c.reset()                      # back to |00…0>, clears the op log

Single-qubit gates

i x y z h s sdg t tdg sx — all take (q).

Parametrised	Signature	Matrix
rx	rx(theta, q)	(e^{-i\theta X/2})
ry	ry(theta, q)	(e^{-i\theta Y/2})
rz	rz(theta, q)	(e^{-i\theta Z/2})
p	p(lam, q)	diag(1, e^{iλ})
u	u(theta, phi, lam, q)	general U3

Multi-qubit gates

Method	Meaning
cx(c, t) / cy / cz / ch	controlled Pauli / Hadamard
crz(theta, c, t), cp(lam, c, t)	controlled rotations
ccx(c1, c2, t)	Toffoli
mcx(ctrls, t)	n-controlled X
swap(a, b)	SWAP (single permutation)
barrier()	log-only marker

All gate methods return self, so they chain: Circuit(3).h(0).cx(0,1).cx(1,2).

Readout

Method	Returns
probabilities()	ndarray of length (2^n)
probs_dict(cutoff=1e-12)	{'011': 0.5, …} for non-negligible states
sample(shots=1024)	{'011': 512, …}, does not collapse
measure(q)	0 or 1, collapses the state
measure_all()	bitstring, collapses everything
expectation_z(q)	(\langle Z_q \rangle \in [-1, 1])
bloch(q)	(x, y, z) of the reduced state
statevector()	a copy of the amplitudes
memory_bytes()	exact footprint
draw()	text dump of non-zero amplitudes
c.ops	human-readable list of everything applied

Prebuilt circuits

bell()                              # (|00> + |11>)/√2
ghz(n=3)                            # (|0…0> + |1…1>)/√2
qft(circuit, qubits=None)           # in-place QFT, O(n²) gates
grover(n, marked, iterations=None)  # defaults to ⌊π/4·√(2ⁿ)⌉ iterations
deutsch_jozsa(n, oracle="balanced") # n inputs + 1 ancilla

---

9. The web app

Open web/index.html. Layout:

• Circuit grid — rows are qubits (q0 at the bottom), columns are time steps. Pick a gate from the palette, click a cell to place it, click it again to erase. Place • on another wire in the same column to make the gate controlled; two controls give a Toffoli. Two × marks in one column form a SWAP. RX/RY/RZ/P prompt for an angle and accept expressions like pi/2.
• State probabilities — every non-zero basis state, sorted.
• Stats — qubit count, amplitude count, exact memory, gate count, the norm (\lVert\psi\rVert) (a live correctness check — it must stay at 1.000000000), and the Shannon entropy of the outcome distribution.
• Bloch vectors — one sphere per qubit, drawn from the reduced density matrix, with the von Neumann entropy S underneath. A shrunken sphere and S > 0 mean that qubit is entangled with the rest.
• Amplitudes — real part, imaginary part, probability.
• Measurement outcomes — sampled shot counts, re-sampleable.
• Export — copies the circuit as OpenQASM 2.0 or as Python for qsim.py.

Presets: Bell pair, GHZ, QFT (4q), Grover (3q searching for |101⟩), Deutsch–Jozsa, and a teleportation setup.

Browser engine notes

engine.js stores amplitudes in two Float64Arrays (re, im) rather than objects, so a gate is a single flat loop with zero allocation per gate. The loop visits each index pair once by skipping indices where the target bit is already set, and a control mask short-circuits the rest — the same (O(2^n))-per-gate cost as the Python version.

---

10. Correctness — how it is verified

python test_qsim.py runs 19 tests:

Test	What it proves
test_single_qubit_against_dense	every 1q gate on every wire of a 4-qubit random state matches the brute-force Kronecker-product result to 1e-12
test_rotations_are_unitary_and_correct	RX/RY/RZ/U3 are unitary; RX(π) = −iX
test_norm_preserved_random_circuit	300 random gates on 8 qubits keep ‖ψ‖ = 1
test_cx_truth_table, test_toffoli_truth_table	all 4 / all 8 classical inputs
test_swap_permutation, test_swap_equals_three_cnots	the permutation shortcut equals CX·CX·CX
test_qft_matches_dft	QFT of a random 4-qubit state equals numpy.fft.ifft(norm="ortho")
test_grover_finds_marked_state	marked state has P > 0.9 and is the argmax
test_deutsch_jozsa	constant oracle → all-zero inputs; balanced → all-one
test_measure_collapses_and_is_consistent	Bell measurements are perfectly correlated and the state stays normalised
test_sampling_statistics	20 000 shots land within 2 % of 50/50
test_expectation_and_bloch	⟨Z⟩ signs, H → +x̂, HS → +ŷ, entangled qubit → zero-length Bloch vector
test_phase_gates_relations	T² = S, S·S† = I
test_cz_symmetry	CZ(a,b) = CZ(b,a)
test_input_validation	bad qubit counts and repeated qubit arguments raise

node tests/test.mjs runs 32 more checks on the browser engine:

• Presets are physically verified, not eyeballed. Bell gives exactly 50/50 on |00⟩,|11⟩ with S = 1; GHZ-4 gives |0000⟩,|1111⟩ with S = 1 on every qubit; QFT of |0000⟩ is perfectly uniform at 6.25 % across all 16 states; Grover finds |101⟩ with P = 94.53 %; Deutsch–Jozsa's balanced oracle leaves the input register at all-ones.
• Engine invariants: norm preserved over 20 random 60-gate circuits, HH = I, T² = S, SWAP is an involution, shot counts sum exactly to shots, and zero-probability states are never drawn.
• Exporters produce the right OpenQASM and Python text.
• Cross-engine: tests/make_ref.py dumps reference state vectors from the Python engine, and test.mjs replays the same 10 circuits in JavaScript. All agree to ≤ 2.2 × 10⁻¹⁶.

This test suite is what caught two real bugs during development: a control/target reversal in the Bell preset, and a Grover circuit that was silently truncated because it needed 26 columns and had been given 14.

---

11. Known limits

• Pure states only — no noise channels, no density-matrix mode.
• No mid-circuit classical feedback in the web UI (the Python measure() does collapse the state, so you can branch on it in Python).
• Memory grows as (2^n); this is fundamental to state-vector simulation, not a flaw in the implementation. Beyond ~26 qubits you need tensor-network or stabiliser methods instead.
• Multi-control gates (3+ controls) use the mask path and are ~3× slower than single-control gates.

12. File map

qsim/
├── qsim.py            Python engine + prebuilt algorithms
├── test_qsim.py       19 correctness tests
├── bench.py           performance + naive-vs-optimised comparison
├── README.md          this document
├── web/
│   ├── index.html     UI shell and styles
│   ├── app.js         DOM layer only (rendering, clicks)
│   ├── circuit.js     pure circuit model: grid -> ops -> simulate, presets, exporters
│   └── engine.js      the simulator, in plain JS
└── tests/
    ├── test.mjs       32 headless checks of circuit.js + engine.js
    └── make_ref.py    dumps ref.json from qsim.py for the cross-engine check

Run everything:

python test_qsim.py                       # Python engine  — 19 tests
cd tests && python make_ref.py && node test.mjs   # JS engine — 32 checks incl. cross-check
python bench.py                           # performance table

app.js holds no physics at all — every circuit decision lives in circuit.js, which is exactly what tests/test.mjs exercises. The UI and the test therefore cannot drift apart.