Quantum Circuit Simulator — Bloch Gallery

Explanations

1. What a qubit register is, mathematically

A register of n qubits is a unit vector |ψ⟩ in ℂ²ⁿ. Nothing more. Index the basis vectors by integers 0 … 2ⁿ−1 and write the integer in binary; bit q of the index is the value of qubit q. This page uses qubit 0 as the least significant bit, so on 3 qubits the state "qubit 2 = 1, others 0" is basis index 4 = |100⟩. The complex coefficients a_i are called amplitudes. If you measure all qubits you get outcome i with probability |a_i|² (the Born rule); that is why the vector must have unit norm. A single qubit is the case n = 1: |ψ⟩ = α|0⟩ + β|1⟩ with |α|² + |β|² = 1.

Two registers side by side form the tensor (Kronecker) product of their spaces. That is where the exponential comes from: 24 qubits need 2²⁴ = 16,777,216 amplitudes. As two 32-bit floats each, that is 128 MiB, which is why this page stores the state on the GPU as two buffers, one for real parts and one for imaginary parts.

2. Gates are unitary matrices, and how the simulator applies them

Physics allows exactly the linear maps that keep the norm, the unitary matrices U^†U = I. A gate on one qubit is a 2×2 unitary. The ones used here:

H = (1/√2)[[1, 1],[1, −1]]      X = [[0,1],[1,0]]     Y = [[0,−i],[i,0]]     Z = [[1,0],[0,−1]]
S = diag(1, i)   T = diag(1, e^{iπ/4})   P(θ) = diag(1, e^{iθ})
RX(θ) = [[cos θ/2, −i sin θ/2],[−i sin θ/2, cos θ/2]]   RY(θ) = [[cos θ/2, −sin θ/2],[sin θ/2, cos θ/2]]   RZ(θ) = diag(e^{−iθ/2}, e^{iθ/2})

Applying a one-qubit gate M to qubit q of an n-qubit register means applying the 2ⁿ×2ⁿ matrix I ⊗ M ⊗ I. Never build that matrix. It is block-sparse: it couples index i only with index i XOR 2^q, the index that differs in bit q. So the whole gate is 2ⁿ⁻¹ independent 2×2 multiplications on pairs (a_i0, a_i1). Each pair is one GPU thread. That is the entire compute shader on this page: pick a pair, multiply by a 2×2 complex matrix, write back in place. Cost per gate is O(2ⁿ), and a circuit of g gates costs O(g·2ⁿ).

Controlled gates (CNOT = cx, CZ, controlled phase cp, Toffoli ccx) apply M to the target only on the basis states where every control bit is 1. In the pair picture: skip the pair unless the control bits of its indices are all set. A control mask handles any number of controls with no extra code, which is how the multi-controlled Z inside Grover is done. SWAP exchanges two qubits: it is the same kernel with the pair chosen as (index with bit a = 1, b = 0) and (bit a = 0, b = 1) and M = X. A controlled SWAP (Fredkin, cswap) adds a control mask. Every gate on this page reduces to that one kernel.

3. The Bloch sphere, derived

A single pure qubit α|0⟩ + β|1⟩ has 4 real parameters, minus 1 for the norm, minus 1 for the global phase (which no measurement can see): 2 real degrees of freedom, a point on a sphere. To make that exact, and to cover qubits that are part of a larger entangled register, use the density matrix ρ. For the whole register ρ = |ψ⟩⟨ψ|. For one qubit inside it, take the partial trace over the other qubits, which in components is just a sum over the other bits:

ρ₀₀ = Σ_rest |a_0,rest|², ρ₁₁ = Σ_rest |a_1,rest|², ρ₀₁ = Σ_rest a_0,rest · conj(a_1,rest)

Any 2×2 Hermitian matrix with trace 1 can be written as ρ = ½(I + xX + yY + zZ) with real x, y, z. Reading off components: z = ρ₀₀ − ρ₁₁, x = 2 Re ρ₀₁, y = −2 Im ρ₀₁. The vector r = (x, y, z) is the Bloch vector, and tr ρ² = ½(1 + |r|²) ≤ 1 gives |r| ≤ 1. A pure qubit sits on the surface. |0⟩ is the north pole, |1⟩ the south pole, (|0⟩+|1⟩)/√2 is +x, (|0⟩+i|1⟩)/√2 is +y. A qubit entangled with others has |r| < 1; in a Bell pair each qubit is at the centre, r = 0, which is the picture of "the individual qubit has no state of its own". This is what the gallery shows after each gate. The reduction is done on the CPU after reading the state back from the GPU.

Rotation gates earn their names here: RX(θ) rotates r by θ around the x axis, likewise RY and RZ; H is a 180° rotation about the axis (x+z)/√2, which swaps the z and x axes.

4. Entanglement: Bell and GHZ

Preset "Bell": H on qubit 0 makes (|00⟩ + |01⟩)/√2 (qubit 0 is the right-hand bit). Then CNOT with control 0, target 1 flips qubit 1 on the second term: (|00⟩ + |11⟩)/√2. This vector cannot be written as a product (a|0⟩+b|1⟩) ⊗ (c|0⟩+d|1⟩), since that would need ad = bc = 0 and ac = bd ≠ 0. Watch the gallery: after H, qubit 0 points along +x; after CNOT both arrows collapse to the centre. GHZ does the same with a chain of CNOTs to give (|0…0⟩ + |1…1⟩)/√2 on any number of qubits.

5. The quantum Fourier transform (QFT)

On m qubits, Q = 2^m, the QFT is the ordinary discrete Fourier transform matrix applied to the amplitude vector:

QFT |x⟩ = Q^−1/2 Σ_y e^{2πi xy/Q} |y⟩

The DFT of a 2^m-vector costs O(m2^m) classically (FFT). The quantum circuit does it with m(m+1)/2 gates because the output factorises: writing y in binary, e^{2πi xy/Q} = Π_k e^{2πi x y_k 2k/Q}, so each output qubit is a product state whose phase depends on the low bits of x only. Circuit: for the most significant qubit down to the least: H, then controlled-phase gates cp(π/2^d) from each lower qubit at distance d. That leaves the output bit-reversed, so the circuit ends with SWAPs. The preset applies the QFT to a basis state |j⟩ (parameter j), which produces a flat probability distribution with phase winding at frequency j; the bar chart shows uniform bars, and the verification block checks every amplitude against the formula above. Mind the cost: the simulator does one pass over the state per gate, so the QFT costs the same O(m²2^m) as a naive DFT here. The speed-up is a statement about the quantum computer, not about simulating it.

6. Grover search, and why it is a rotation

Problem: among N = 2ⁿ items, one item w (the preset parameter) satisfies a predicate you can evaluate but not invert. Classically you expect N/2 evaluations; Grover needs about (π/4)√N. The circuit: start with H on every qubit, the uniform superposition |s⟩ = N^−1/2 Σ|i⟩. Then repeat two steps. The oracle flips the sign of the amplitude of |w⟩ only: X on the qubits where w has a 0 bit, a multi-controlled Z, and the X gates again. The diffusion step reflects the amplitudes about their mean: H on all, X on all, multi-controlled Z, X on all, H on all, which is −(2|s⟩⟨s| − I).

The state always stays in the real plane spanned by |w⟩ and its orthogonal complement inside |s⟩. Two reflections compose to a rotation by twice the angle between the mirrors. With sin θ = ⟨w|s⟩ = 1/√N, after k iterations the amplitude on |w⟩ is (−1)^k sin((2k+1)θ) (the sign is the global phase from the way diffusion is built). Success probability sin²((2k+1)θ) peaks at k ≈ π/(4θ) − ½. Overshooting past that makes it worse again, which the stepper shows if you keep playing. The status line prints the predicted value next to the simulated one at every iteration.

7. Shor's algorithm: factor 15 with a = 7, step by step

Number theory first. To factor N, pick a coprime to N and find the period r of f(x) = a^x mod N. If r is even and a^r/2 ≠ −1 (mod N), then (a^r/2 − 1)(a^r/2 + 1) = a^r − 1 ≡ 0 (mod N) and neither factor is 0 mod N, so gcd(a^r/2 ± 1, N) are non-trivial factors. For N = 15, a = 7: powers of 7 mod 15 are 1, 7, 4, 13, 1, … so r = 4, 7² = 49 ≡ 4, and gcd(3, 15) = 3, gcd(5, 15) = 5. Everything hard is in finding r; the quantum part does exactly that.

Registers. A counting register of t = 8 qubits (qubits 0–7, Q = 256) and a work register of 4 qubits (qubits 8–11) holding numbers mod 15. The work register starts at |1⟩ (an X on qubit 8).

Stage 1. H on every counting qubit: Q^−1/2 Σ_x |x⟩|1⟩.

Stage 2, modular exponentiation. Apply |x⟩|w⟩ → |x⟩|w · 7^x mod 15⟩. Write x in binary, 7^x = Π_j (7^2j)^x_j, so counting qubit j controls "multiply by 7^2j mod 15". Those constants are 7, 4, 1, 1, 1, 1, 1, 1: only qubits 0 and 1 do anything, because r = 4 divides 2². The page still shows the identity stages, labelled, so the structure is visible. Multiplication by a constant mod 15 on 4 bits is a permutation of the 16 basis states, and for these constants it needs no arithmetic: 15 = 2⁴ − 1, so multiplying by 2 mod 15 is a cyclic bit rotation. ×4 is rotate by 2 = swap(w0,w2), swap(w1,w3). ×8 is rotate by 3, and ×7 = −8 mod 15 is that rotation followed by NOT on all four bits (valid on every nonzero residue, and the work register never holds 0). Each swap becomes a Fredkin gate and each NOT a CNOT, all controlled by the counting qubit. The verification block checks this multiplier on all eight reachable residues against the integer arithmetic. After this stage the state is Q^−1/2 Σ_x |x⟩|7^x mod 15⟩, and the counting register, seen alone, is a uniform mixture of four periodic combs, one per value of the work register; its Bloch vectors sit near the centre.

Stage 3, inverse QFT on the counting register. Each comb has period r = 4 in x, so its Fourier transform is concentrated on multiples of Q/r = 64. Because r divides Q exactly, the result is exact: probability ¼ on each of k ∈ {0, 64, 128, 192} and 0 elsewhere. The bar chart shows the marginal distribution of the counting register; the verification block checks the four peaks.

Classical post-processing. A measured k satisfies k/Q ≈ s/r for some integer s. The continued-fraction expansion of k/256 finds the fraction: 64/256 = 1/4 and 192/256 = 3/4 give r = 4, which works. 128/256 = 1/2 gives a candidate r = 2, and 7² = 49 ≡ 4 ≠ 1, so that run fails; k = 0 also fails. The algorithm is probabilistic: half of the runs succeed here, and a failure just means running it again. The status line walks through this for each of the four outcomes after the last stage.

8. Implementation notes (WebGPU)

State: two storage buffers of 2ⁿ f32 (real, imaginary). Uniform per gate: mode (target or swap pair), target indices, control mask, and the 2×2 complex matrix as 8 floats, packed 512 gates per command buffer at 256-byte offsets. One dispatch per gate, workgroup size 256, 2ⁿ⁻¹ invocations (2ⁿ⁻² for SWAP), which at 24 qubits is 32,768 workgroups, under the 65,535 limit. Reading the state back for the charts and Bloch vectors copies both buffers into a mapped staging buffer; at 24 qubits that is 128 MiB per readback, so the page reads back only when it needs to draw. The CPU backend runs the identical pair loop in JavaScript with 64-bit floats and is capped at 22 qubits. Float32 rounding on the GPU is visible in the verification as max |Δ| of order 10⁻⁶ against the 64-bit references.

Gate angle syntax accepts pi, e.g. rz(pi/4) 3, cp(2*pi/8) 0 1.