Building a Small Simulator From Scratch

A state vector with a declared index map#

For n qubits, allocate $N=2^n$ complex amplitudes. This reference uses big-endian display: q0 is the most-significant bit, so basis string 101 is integer index five and q0 contributes weight four. The zero constructor sets amplitude zero to one and every other amplitude to zero. Pure states inhabit this exponentially growing complex vector space [Nielsen and Chuang].

The constructor should reject a vector unless

within absolute tolerance 10⁻¹². Silent normalization is tempting and wrong: it can turn a caller’s bug into a plausible state. Validate dimensions, finite values, and norm; report the faulty quantity.

Applying a one-qubit matrix without building 2^n by 2^n#

A gate on target qubit q mixes amplitude pairs whose indices differ only at bit q. With big-endian display, the target mask is 1 << (n - 1 - q). Iterate indices whose target bit is zero, copy the pair (a0,a1), and write (u00*a0+u01*a1, u10*a0+u11*a1). Copying both old values before either write prevents an in-place corruption that often passes basis-state tests and fails superpositions.

If $U^\dagger U=I$, then the pair norm is preserved; summing over disjoint pairs preserves the whole norm. Test the implementation anyway. A transposed gate, wrong mask, or aliased write can be algebraically unitary and operationally wrong.

Why the pair update is the tensor product#

Fix every bit except the target. The two corresponding basis states span a two-dimensional slice of the full vector, with amplitudes $a_0$ and $a_1$. Applying $U$ to that qubit applies the same $2\times2$ map to every such slice. This is exactly $I\otimes\cdots\otimes U\otimes\cdots\otimes I$, evaluated without materializing its mostly zero $2^n\times2^n$ matrix. There are $2^{n-1}$ disjoint pairs, so a one-qubit gate costs $O(2^n)$ arithmetic and a copied output vector costs $O(2^n)$ memory.

A diagnostic trace should print the target mask and the first few index pairs. For three qubits, target q1 must produce (0,2), (1,3), (4,6), and (5,7). If the trace instead begins (0,4), the implementation is applying the gate to q0. That one line often localizes an ordering defect faster than comparing eight complex amplitudes after a long circuit.

CNOT exposes every ordering decision#

CNOT flips the target bit exactly when the control bit is one. Do not start with Bell states. First assert the complete two-qubit truth table: $00\mapsto00$, $01\mapsto01$, $10\mapsto11$, $11\mapsto10$ for control $q_0$ and target $q_1$ under this book’s display. Then repeat with control and target exchanged. Those eight cases identify which integer mask belongs to which logical label.

Qiskit and Cirq expose their own ordering and output conventions [Qiskit documentation] [Cirq documentation]. An adapter may reverse a displayed string without changing the state. That is a translation, not a reason to make the core simulator ambiguous. OpenQASM defines serialized circuit semantics and QIR defines an intermediate representation, but neither settles a Python library’s in-memory endianness [OpenQASM 3 specification] [QIR Alliance specification].

The implementation can be read as a permutation. For each source index $i$, compute the control mask; if it is set, move the amplitude to $i\mathbin{\oplus}m_t$, otherwise leave it at $i$. Every destination must be written exactly once. Summing into a zeroed vector makes that bijection visible and avoids a second in-place hazard: swapping an amplitude twice while iterating both members of a pair.

Control equal to target is not a degenerate CNOT; it is an invalid request and should raise before mutation. Negative or out-of-range labels fail for the same reason. Public methods need these rejection paths because Python’s shifts and negative indices can otherwise turn a bad circuit into a different, executable circuit.

The Bell fixture has exact and stochastic halves#

Apply H to q0 and then CNOT(q0,q1). The exact vector is (1/√2, 0, 0, 1/√2). The Born rule maps an amplitude to its squared magnitude [Nielsen and Chuang], producing probabilities 0.5 on 00 and 11. These are deterministic numerical assertions at tolerance 10⁻¹².

Sampling is a different API. It consumes probabilities, a shot count, and an optional seed; it returns integer counts. Assert that counts total the shots and impossible outcomes remain absent for the ideal Bell state. Use the seed only to stabilize a regression transcript. An unseeded scientific check needs a declared binomial interval, not an exact count.

Construct the cumulative distribution in basis-index order and force its final boundary to one after checking the probability sum within tolerance. The correction is only for floating-point roundoff; it must not rescue a vector whose norm is materially wrong. Given a uniform variate $r\in[0,1)$, return the first index whose cumulative mass exceeds $r$. Tests should cover zero shots, a negative shot request, a basis state with probability one, and a distribution whose tiny tail would expose an off-by-one comparison.

from math import sqrt
H = ((1/sqrt(2), 1/sqrt(2)), (1/sqrt(2), -1/sqrt(2)))
state = [1+0j, 0j, 0j, 0j]

def one_qubit(vector, gate, target, n=2):
    out = vector[:]
    mask = 1 << (n - 1 - target)
    for i in range(len(vector)):
        if i & mask: continue
        j = i | mask; a, b = vector[i], vector[j]
        out[i] = gate[0][0]*a + gate[0][1]*b
        out[j] = gate[1][0]*a + gate[1][1]*b
    return out

def cnot(vector, control, target, n=2):
    out = [0j] * len(vector)
    cm = 1 << (n - 1 - control); tm = 1 << (n - 1 - target)
    for i, amplitude in enumerate(vector): out[i ^ tm if i & cm else i] += amplitude
    return out

product_state = one_qubit(state, H, 0)
state = cnot(product_state, 0, 1)
missing_entangler = cnot(product_state, 1, 0)
expected = [1/sqrt(2), 0, 0, 1/sqrt(2)]
assert all(abs(a-b) <= 1e-12 for a,b in zip(state, expected))
assert abs(sum(abs(a)**2 for a in state)-1) <= 1e-12
assert missing_entangler == product_state and missing_entangler != state
print(f"PASS: 37 simulator Bell_support={[i for i,a in enumerate(state) if abs(a) > 1e-12]} missing_entangler_support={[i for i,a in enumerate(missing_entangler) if abs(a) > 1e-12]}")

Tests that catch plausible-looking wrong simulators#

Reject non-power-of-two dimensions and non-normalized input. Apply H twice to several states, not only $\lvert 0 \rangle$. Exercise CNOT on every basis state in both directions. Check a phase-sensitive superposition so a sign bug cannot hide behind computational-basis probabilities. Assert impossible Bell outcomes, exact probability sums, sample totals, and error messages for invalid qubit indices. Include a randomized unitary test with a fixed generator seed, but keep small analytic fixtures as the diagnostic core.

A controlled-Z extension is especially useful: its computational-basis probabilities do not change, yet the relative phase of $\lvert 11 \rangle$ does. The Bell state used above happens to acquire a relative minus sign, not a global phase; an X-basis measurement or later interference reveals it.

Test relative phase by composing operations, not by inspecting a comment. Prepare $(|00\rangle+|11\rangle)/\sqrt2$, apply CZ, then apply H to both qubits. The resulting probability distribution differs from the unmodified Bell path even though the distribution immediately after CZ does not. This fixture catches a simulator that stores only magnitudes, conjugates at the wrong point, or quietly rounds complex amplitudes to real values.

Property tests complement, rather than replace, exact cases. Generate a small normalized vector and a unitary matrix from a fixed seed; assert norm preservation and agreement with a slow reference matrix product. Then shrink a failure to the explicit vector, matrix, target, and mask. A seed without the generated inputs is not a reproducible counterexample.

Where the model ends#

Memory and gate work grow exponentially with qubit count. This implementation has no mid-circuit measurement, reset, classical control, density matrices, leakage, or device noise. It is an oracle for small pure-state circuits, not a performance competitor to optimized simulators and not a proxy for hardware. Its public API earns trust only through its declared domain and tests.

Those omissions are architectural, not missing convenience methods. Mid-circuit measurement changes a pure state conditionally and requires a classical outcome record; a density-matrix backend changes storage from $2^n$ amplitudes to $4^n$ complex entries; leakage requires states outside the qubit basis. Adding any one of them changes invariants and tests. A small simulator stays trustworthy by refusing unsupported operations rather than attaching device-sounding names to an ideal state vector.

Numerical comparison needs a policy too. Compare complex amplitudes with an absolute tolerance on small analytic fixtures and compare state vectors up to global phase only when the contract permits it. Probability agreement is weaker: states with different relative phase can share every computational-basis probability. The API should expose amplitudes for exact tests, probabilities for measurement oracles, and counts for sampling; converting all three into one dictionary would erase the evidence that locates a failure.

Every public method named here now has a corresponding rejection or invariant test: constructor dimension and norm, one-qubit target range and pair update, CNOT control/target validity and truth table, Born probability normalization, and sampler shot totals. That coverage is more meaningful than a line percentage because it follows the semantic contracts a caller can violate.