Build a single-qubit quantum simulator in Python using NumPy. You will represent quantum states as vectors, apply gates as matrix multiplications, and simulate measurement with probabilistic collapse. By the end, you will see the Hadamard gate create a 50/50 superposition and verify it with a histogram of 1000 measurements.
You need Python 3.8+ and NumPy installed.
python3 --version
pip3 install numpy
Create a single file called qubit.py. All code goes in this file.
A single qubit lives in a 2D complex vector space. The two basis states are:
Any qubit state is a linear combination: a|0> + b|1> where |a|^2 + |b|^2 = 1.
Start qubit.py with:
import numpy as np
# Basis states
ket_0 = np.array([1, 0], dtype=complex)
ket_1 = np.array([0, 1], dtype=complex)
def display_state(state, label="State"):
"""Print the state vector with amplitudes and probabilities."""
print(f"\n{label}:")
print(f" |0> amplitude: {state[0]:.4f} (prob: {abs(state[0])**2:.4f})")
print(f" |1> amplitude: {state[1]:.4f} (prob: {abs(state[1])**2:.4f})")
display_state(ket_0, "|0>")
display_state(ket_1, "|1>")
Run it:
python3 qubit.py
|0>:
|0> amplitude: 1.0000+0.0000j (prob: 1.0000)
|1> amplitude: 0.0000+0.0000j (prob: 0.0000)
|1>:
|0> amplitude: 0.0000+0.0000j (prob: 0.0000)
|1> amplitude: 1.0000+0.0000j (prob: 1.0000)
The state |0> has 100% probability of measuring 0 and 0% probability of measuring 1. This matches classical behavior – no superposition yet.
Quantum gates are unitary matrices. Add these three fundamental gates:
# Pauli-X gate (quantum NOT): flips |0> to |1> and vice versa
X = np.array([[0, 1],
[1, 0]], dtype=complex)
# Hadamard gate: creates equal superposition
H = (1 / np.sqrt(2)) * np.array([[1, 1],
[1, -1]], dtype=complex)
# Phase gate (S gate): adds a phase of i to |1>
S = np.array([[1, 0],
[0, 1j]], dtype=complex)
def apply_gate(gate, state):
"""Apply a gate (matrix) to a state (vector)."""
return gate @ state
Test the X gate – it should flip |0> to |1>:
flipped = apply_gate(X, ket_0)
display_state(flipped, "X|0> (should be |1>)")
X|0> (should be |1>):
|0> amplitude: 0.0000+0.0000j (prob: 0.0000)
|1> amplitude: 1.0000+0.0000j (prob: 1.0000)
The Hadamard gate is the key to quantum behavior. Apply it to |0>:
superposed = apply_gate(H, ket_0)
display_state(superposed, "H|0> (superposition)")
H|0> (superposition):
|0> amplitude: 0.7071+0.0000j (prob: 0.5000)
|1> amplitude: 0.7071+0.0000j (prob: 0.5000)
The amplitudes are both 1/sqrt(2) ~ 0.7071. Squaring gives 0.5 – a 50/50 chance of measuring either outcome. The qubit is in a genuine superposition, not “we don’t know which state it’s in.” Both states coexist until measurement.
Measurement collapses the superposition. The probability of each outcome is the squared amplitude. After measurement, the state is the basis state corresponding to the result.
def measure(state, num_shots=1):
"""Simulate quantum measurement.
Returns a list of measurement results (0 or 1).
Each measurement collapses the state probabilistically.
"""
prob_0 = abs(state[0]) ** 2
results = []
for _ in range(num_shots):
outcome = 0 if np.random.random() < prob_0 else 1
results.append(outcome)
return results
Test with a single measurement of the superposed state:
result = measure(superposed, num_shots=1)
print(f"\nSingle measurement of H|0>: {result[0]}")
print("(Run again -- you may get a different result)")
Single measurement of H|0>: 0
(Run again -- you may get a different result)
The output will be either 0 or 1 with roughly equal probability across multiple runs.
One measurement tells you little. Run 1000 measurements to see the probability distribution:
def histogram(results):
"""Print a text-based histogram of measurement results."""
counts = {0: 0, 1: 0}
for r in results:
counts[r] += 1
total = len(results)
print(f"\nMeasurement results ({total} shots):")
for outcome in [0, 1]:
count = counts[outcome]
frac = count / total
bar = "#" * int(frac * 50)
print(f" |{outcome}>: {bar} {count} ({frac:.1%})")
# Measure the superposition 1000 times
results = measure(superposed, num_shots=1000)
histogram(results)
Measurement results (1000 shots):
|0>: ######################### 503 (50.3%)
|1>: ######################## 497 (49.7%)
The exact counts will vary, but they will cluster around 500/500. This is the Born rule in action: the probability of outcome k is |amplitude_k|^2.
Compare with a state that is not in superposition. Measuring |0> directly should always give 0:
print("\n--- Control experiment: measuring |0> ---")
results_zero = measure(ket_0, num_shots=1000)
histogram(results_zero)
print("\n--- X gate then measure: measuring |1> ---")
results_one = measure(apply_gate(X, ket_0), num_shots=1000)
histogram(results_one)
--- Control experiment: measuring |0> ---
Measurement results (1000 shots):
|0>: ################################################## 1000 (100.0%)
|1>: 0 (0.0%)
--- X gate then measure: measuring |1> ---
Measurement results (1000 shots):
|0>: 0 (0.0%)
|1>: ################################################## 1000 (100.0%)
A single qubit cannot demonstrate entanglement. Extend the simulator to two qubits. A two-qubit state is a 4-element vector in the tensor product space.
print("\n=== Two-Qubit Extension ===")
# Two-qubit basis states
ket_00 = np.kron(ket_0, ket_0) # [1, 0, 0, 0]
ket_01 = np.kron(ket_0, ket_1) # [0, 1, 0, 0]
ket_10 = np.kron(ket_1, ket_0) # [0, 0, 1, 0]
ket_11 = np.kron(ket_1, ket_1) # [0, 0, 0, 1]
# CNOT gate: flips the second qubit if the first qubit is |1>
CNOT = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 0, 1],
[0, 0, 1, 0]], dtype=complex)
# Create a Bell state: H on first qubit, then CNOT
# This produces (|00> + |11>) / sqrt(2) -- maximally entangled
H_on_first = np.kron(H, np.eye(2)) # H tensor I
bell_state = CNOT @ (H_on_first @ ket_00)
def display_two_qubit(state, label="State"):
"""Print a two-qubit state vector."""
labels = ["00", "01", "10", "11"]
print(f"\n{label}:")
for i, lbl in enumerate(labels):
amp = state[i]
prob = abs(amp) ** 2
if prob > 1e-10:
print(f" |{lbl}> amplitude: {amp:.4f} (prob: {prob:.4f})")
display_two_qubit(bell_state, "Bell state (|00> + |11>)/sqrt(2)")
def measure_two_qubit(state, num_shots=1000):
"""Measure a two-qubit state, returning pairs of results."""
probs = np.abs(state) ** 2
outcomes = np.random.choice(4, size=num_shots, p=probs)
labels = ["00", "01", "10", "11"]
counts = {lbl: 0 for lbl in labels}
for o in outcomes:
counts[labels[o]] += 1
print(f"\nTwo-qubit measurement ({num_shots} shots):")
for lbl in labels:
count = counts[lbl]
frac = count / num_shots
bar = "#" * int(frac * 50)
if count > 0:
print(f" |{lbl}>: {bar} {count} ({frac:.1%})")
measure_two_qubit(bell_state)
Bell state (|00> + |11>)/sqrt(2):
|00> amplitude: 0.7071+0.0000j (prob: 0.5000)
|11> amplitude: 0.7071+0.0000j (prob: 0.5000)
Two-qubit measurement (1000 shots):
|00>: ######################### 498 (49.8%)
|11>: ######################### 502 (50.2%)
You will only ever see |00> or |11>, never |01> or |10>. The qubits are entangled: measuring one instantly determines the other. This is not a classical correlation – it was created by quantum gates, and no local hidden variable can reproduce the statistics for all possible measurement bases.