@@ -0,0 +1,6 @@
+# Programs for chapter 1
+
+Every number, table and figure of chapter 1 is produced by the programs in this folder (see `../../CLAUDE.md` and the `qc-book-code` skill). Keep this list and the `Programs for this chapter` section of `chapter1.tex` in sync.
+
+| Program | What it computes | Runtime |
+|---|---|---|
@@ -0,0 +1,6 @@
+# Programs for chapter 10
+
+Every number, table and figure of chapter 10 is produced by the programs in this folder (see `../../CLAUDE.md` and the `qc-book-code` skill). Keep this list and the `Programs for this chapter` section of `chapter10.tex` in sync.
+
+| Program | What it computes | Runtime |
+|---|---|---|
@@ -0,0 +1,6 @@
+# Programs for chapter 11
+
+Every number, table and figure of chapter 11 is produced by the programs in this folder (see `../../CLAUDE.md` and the `qc-book-code` skill). Keep this list and the `Programs for this chapter` section of `chapter11.tex` in sync.
+
+| Program | What it computes | Runtime |
+|---|---|---|
@@ -0,0 +1,15 @@
+# Programs for chapter 2
+
+Every number, table and figure of chapter 2 is produced by the programs in this folder together with the shared modules in `../common/` that this chapter introduces (see `../../CLAUDE.md` and the `qc-book-code` skill). Keep this list and the `Programs for this chapter` section of `chapter2.tex` in sync. Figures are written to `../../BookFigures/chapter2/` as PDF (for LaTeX) and PNG (for the notebook).
+
+| Program | What it computes | Runtime |
+|---|---|---|
+| `dirac_algebra.py` | Worked examples of Sec. 2.2: inner/outer products, projectors, Pauli algebra, tensor products, spectral decomposition | instant |
+| `bloch_sphere.py` | Bloch angles and vectors, Table 2.1, Fig. 2.1 (Pauli eigenstates, Larmor precession, rotations) | seconds |
+| `measurement.py` | Born rule, post-measurement states, X-basis measurement, marginals of a 3-qubit register, Fig. 2.2 (shot noise) | seconds |
+| `decoherence.py` | Dephasing and amplitude-damping channels (Kraus form), Fig. 2.3 | instant |
+| `entanglement.py` | Bell states, measurement correlations, partial traces, Schmidt decomposition via SVD, Table 2.2, Fig. 2.4 | instant |
+| `protocols.py` | Superdense coding (all four messages) and teleportation (200 random states) on the simulator | instant |
+| `two_level_system.py` | Two-level Hamiltonian in matrix and Pauli form, Table 2.3, grid search over the R_y R_x ansatz, Fig. 2.8 | seconds |
+
+Shared modules introduced by this chapter (in `../common/`): `gates.py`, `statevector.py`, `circuit.py`, `density_matrix.py`, `bookplot.py`. All programs run on `numpy` and `matplotlib` only; each runs as a script with `assert`-based checks.
@@ -0,0 +1,116 @@
+"""Chapter 2, Sec. 2.3: the Bloch sphere.
+
+Converts between the amplitudes of a qubit and its Bloch angles
+(Eq. 2.blochstate), computes the Bloch vector r = (<X>, <Y>, <Z>)
+(Eq. 2.blochvector), lists the Bloch angles of the six Pauli eigenstates
+(Table 2.paulieigen) and draws Fig. 2.bloch: the sphere with the Pauli
+eigenstates, the trajectory of a spin precessing in a magnetic field along
+z (Eq. 2.larmor) and the two-rotation path R_y(phi) R_x(theta)|0> used as a
+variational ansatz in Sec. 2.11.  Runs in a few seconds.
+"""
+import sys, pathlib
+root = pathlib.Path(__file__).resolve().parents[1]
+sys.path[:0] = [str(p) for p in root.iterdir() if p.is_dir()]
+
+import numpy as np
+import gates as G
+from bookplot import savefig, plt
+
+
+def bloch_state(theta, phi):
+    """|psi> = cos(theta/2)|0> + e^{i phi} sin(theta/2)|1>."""
+    return np.array([np.cos(theta / 2), np.exp(1j * phi) * np.sin(theta / 2)])
+
+
+def bloch_angles(psi):
+    """Inverse of bloch_state, removing the global phase first."""
+    psi = np.asarray(psi, dtype=complex)
+    psi = psi / np.linalg.norm(psi)
+    psi = psi * np.exp(-1j * np.angle(psi[0]))       # make the |0> amplitude real
+    theta = 2 * np.arccos(np.clip(psi[0].real, -1, 1))
+    phi = np.angle(psi[1]) if abs(psi[1]) > 1e-12 else 0.0
+    return theta, phi % (2 * np.pi)
+
+
+def bloch_vector(psi):
+    """r = (<X>, <Y>, <Z>) for a pure single-qubit state."""
+    psi = np.asarray(psi, dtype=complex)
+    return np.array([np.vdot(psi, P @ psi).real for P in (G.X, G.Y, G.Z)])
+
+
+PAULI_EIGENSTATES = {
+    "|0>": ([1, 0], "Z", +1), "|1>": ([0, 1], "Z", -1),
+    "|+>": ([1, 1] / np.sqrt(2), "X", +1), "|->": ([1, -1] / np.sqrt(2), "X", -1),
+    "|+i>": ([1, 1j] / np.sqrt(2), "Y", +1), "|-i>": ([1, -1j] / np.sqrt(2), "Y", -1),
+}
+
+
+def larmor_trajectory(theta0, phi0, omega, times):
+    """Points on the sphere for |psi(t)> = exp(-i H t) |psi(0)>, H = -omega Z/2."""
+    psi0 = bloch_state(theta0, phi0)
+    pts = []
+    for t in times:
+        U = G.Rz(-omega * t)                     # exp(-i H t) = exp(+i omega t Z / 2)
+        pts.append(bloch_vector(U @ psi0))
+    return np.array(pts)
+
+
+def draw_bloch(ax, elev=20, azim=35):
+    u, v = np.mgrid[0:2 * np.pi:60j, 0:np.pi:30j]
+    ax.plot_surface(np.cos(u) * np.sin(v), np.sin(u) * np.sin(v), np.cos(v),
+                    color="lightsteelblue", alpha=0.15, linewidth=0)
+    t = np.linspace(0, 2 * np.pi, 200)
+    for circ in ((np.cos(t), np.sin(t), 0 * t), (np.cos(t), 0 * t, np.sin(t)), (0 * t, np.cos(t), np.sin(t))):
+        ax.plot(*circ, color="gray", lw=0.5)
+    for axis, label in (((1.3, 0, 0), "x"), ((0, 1.3, 0), "y"), ((0, 0, 1.3), "z")):
+        ax.quiver(0, 0, 0, *axis, color="k", arrow_length_ratio=0.08, lw=0.8)
+        ax.text(*(1.4 * np.array(axis)), label)
+    ax.set_box_aspect((1, 1, 1)); ax.set_axis_off(); ax.view_init(elev=elev, azim=azim)
+    ax.set_xlim(-1.2, 1.2); ax.set_ylim(-1.2, 1.2); ax.set_zlim(-1.2, 1.2)
+
+
+if __name__ == "__main__":
+    np.set_printoptions(precision=4, suppress=True)
+    # Table 2.paulieigen
+    print(f"{'state':6s} {'operator':9s} {'eigenvalue':>10s} {'theta':>8s} {'phi':>8s}   Bloch vector")
+    for name, (psi, op, ev) in PAULI_EIGENSTATES.items():
+        th, ph = bloch_angles(psi)
+        print(f"{name:6s} {op:9s} {ev:>10d} {th/np.pi:7.2f}π {ph/np.pi:7.2f}π   {bloch_vector(psi)}")
+    # round trip amplitudes -> angles -> amplitudes
+    theta, phi = 2.0, 0.7
+    psi = bloch_state(theta, phi)
+    assert np.allclose(bloch_angles(psi), (theta, phi))
+    assert np.allclose(bloch_vector(psi), [np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta)])
+    # a global phase leaves the Bloch vector unchanged
+    assert np.allclose(bloch_vector(np.exp(0.3j) * psi), bloch_vector(psi))
+    print("\nBloch angles of psi =", psi, "->", f"theta={theta}, phi={phi}")
+
+    # Figure 2.bloch
+    fig = plt.figure(figsize=(9, 4.4))
+    ax = fig.add_subplot(1, 2, 1, projection="3d")
+    draw_bloch(ax)
+    for name, (v, op, ev) in PAULI_EIGENSTATES.items():
+        r = bloch_vector(v)
+        ax.scatter(*r, color="C3", s=25)
+        ax.text(*(1.12 * r), name, fontsize=8)
+    r = bloch_vector(psi)
+    ax.quiver(0, 0, 0, *r, color="C0", arrow_length_ratio=0.1, lw=2)
+    ax.text(*(1.1 * r), r"$|\psi\rangle$", color="C0")
+    ax.set_title("(a) Pauli eigenstates and a general state", fontsize=9)
+
+    ax = fig.add_subplot(1, 2, 2, projection="3d")
+    draw_bloch(ax)
+    traj = larmor_trajectory(theta0=1.0, phi0=0.0, omega=1.0, times=np.linspace(0, 2 * np.pi, 200))
+    ax.plot(traj[:, 0], traj[:, 1], traj[:, 2], color="C1", lw=2, label="Larmor precession")
+    # the path R_y(phi) R_x(theta)|0>: first rotate about x, then about y
+    k0 = np.array([1, 0], dtype=complex)
+    th_path = np.linspace(0, 2.0, 60)
+    p1 = np.array([bloch_vector(G.Rx(t) @ k0) for t in th_path])
+    ph_path = np.linspace(0, 1.2, 60)
+    p2 = np.array([bloch_vector(G.Ry(p) @ G.Rx(2.0) @ k0) for p in ph_path])
+    ax.plot(p1[:, 0], p1[:, 1], p1[:, 2], color="C2", lw=2, label=r"$R_x(\theta)|0\rangle$")
+    ax.plot(p2[:, 0], p2[:, 1], p2[:, 2], color="C4", lw=2, label=r"$R_y(\phi)R_x(2)|0\rangle$")
+    ax.legend(loc="lower left", fontsize=7)
+    ax.set_title("(b) rotations and precession", fontsize=9)
+    out = savefig(fig, 2, "bloch_sphere")
+    print("figure written to", out)
@@ -0,0 +1,81 @@
+"""Chapter 2, Sec. 2.5.5: decoherence seen through the density matrix.
+
+Two elementary noise channels on one qubit, written as Kraus maps
+rho -> sum_k E_k rho E_k^dagger (Eq. 2.kraus): dephasing, which shrinks the
+off-diagonal elements of rho by a factor (1 - 2p) and leaves the populations
+untouched (Eq. 2.dephasing), and amplitude damping, which drives the qubit
+towards |0> (Eq. 2.damping).  Fig. 2.decoherence shows the Bloch vector of
+|+> under repeated application of each channel and the accompanying loss of
+purity.  Runs in a second.
+"""
+import sys, pathlib
+root = pathlib.Path(__file__).resolve().parents[1]
+sys.path[:0] = [str(p) for p in root.iterdir() if p.is_dir()]
+
+import numpy as np
+import gates as G
+from density_matrix import DensityMatrix
+from bookplot import savefig, plt
+
+
+def kraus_apply(rho, kraus):
+    """rho -> sum_k E_k rho E_k^dagger."""
+    return sum(E @ rho @ E.conj().T for E in kraus)
+
+
+def dephasing(p):
+    """With probability p a Z error hits the qubit."""
+    return [np.sqrt(1 - p) * G.I, np.sqrt(p) * G.Z]
+
+
+def amplitude_damping(gamma):
+    """|1> decays to |0> with probability gamma (spontaneous emission)."""
+    E0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]], dtype=complex)
+    E1 = np.array([[0, np.sqrt(gamma)], [0, 0]], dtype=complex)
+    return [E0, E1]
+
+
+def trajectory(rho0, kraus, steps):
+    rho = rho0.copy()
+    out = [DensityMatrix(rho)]
+    for _ in range(steps):
+        rho = kraus_apply(rho, kraus)
+        out.append(DensityMatrix(rho))
+    return out
+
+
+if __name__ == "__main__":
+    np.set_printoptions(precision=4, suppress=True)
+    plus = np.array([1, 1]) / np.sqrt(2)
+    rho0 = np.outer(plus, plus.conj())
+    # Kraus operators satisfy sum_k E_k^dagger E_k = I
+    for kraus in (dephasing(0.1), amplitude_damping(0.1)):
+        assert np.allclose(sum(E.conj().T @ E for E in kraus), G.I)
+    # one dephasing step: off-diagonals scaled by 1 - 2p, Eq. (2.dephasing)
+    p = 0.1
+    rho1 = kraus_apply(rho0, dephasing(p))
+    assert np.allclose(rho1, [[0.5, 0.5 * (1 - 2 * p)], [0.5 * (1 - 2 * p), 0.5]])
+    print("rho(|+>) after one dephasing step with p = 0.1:\n", rho1.real)
+    steps = 30
+    deph = trajectory(rho0, dephasing(0.05), steps)
+    damp = trajectory(rho0, amplitude_damping(0.08), steps)
+    print(f"{'step':>4s} {'purity (dephasing)':>20s} {'purity (damping)':>18s} {'<Z> (damping)':>14s}")
+    for k in (0, 5, 10, 20, 30):
+        print(f"{k:4d} {deph[k].purity():20.4f} {damp[k].purity():18.4f} {damp[k].bloch_vector()[2]:14.4f}")
+    # the fully dephased state is the maximally mixed state, entropy 1 bit
+    final = trajectory(rho0, dephasing(0.5), 1)[-1]
+    assert np.allclose(final.rho, np.eye(2) / 2) and abs(final.von_neumann_entropy() - 1) < 1e-12
+
+    fig, axes = plt.subplots(1, 2, figsize=(8, 3.2))
+    t = np.arange(steps + 1)
+    for ax, traj, title in ((axes[0], deph, "dephasing, $p=0.05$ per step"),
+                            (axes[1], damp, r"amplitude damping, $\gamma=0.08$ per step")):
+        r = np.array([d.bloch_vector() for d in traj])
+        ax.plot(t, r[:, 0], label=r"$\langle X\rangle$")
+        ax.plot(t, r[:, 2], label=r"$\langle Z\rangle$")
+        ax.plot(t, [d.purity() for d in traj], "k--", label=r"purity $\mathrm{Tr}\rho^2$")
+        ax.plot(t, [d.von_neumann_entropy() for d in traj], "k:", label=r"entropy $S(\rho)$")
+        ax.set_xlabel("channel applications"); ax.set_title(title, fontsize=9)
+    axes[0].set_ylim(-0.05, 1.05); axes[1].set_ylim(-0.05, 1.05)
+    axes[0].legend(fontsize=8, loc="center right")
+    print("figure written to", savefig(fig, 2, "decoherence"))
@@ -0,0 +1,90 @@
+"""Chapter 2, Sec. 2.2: the worked linear-algebra examples of the text.
+
+Inner and outer products in Dirac notation (Eqs. 2.innerex-2.outer),
+projection operators and their action on a superposition (Eq. 2.project),
+Pauli algebra (Eqs. 2.pauliprod-2.paulicomm), tensor products of states and
+matrices (Eqs. 2.tensorstate-2.tensormatrix) and the spectral decomposition
+of a Hermitian 2x2 matrix (Eq. 2.spectral).  Runs in well under a second.
+"""
+import sys, pathlib
+root = pathlib.Path(__file__).resolve().parents[1]          # BookPrograms/
+sys.path[:0] = [str(p) for p in root.iterdir() if p.is_dir()]
+
+import numpy as np
+import gates as G
+
+
+def ket(*amplitudes):
+    return np.array(amplitudes, dtype=complex)
+
+
+def bra(v):
+    return v.conj()
+
+
+def inner(u, v):
+    """<u|v> = sum_i u_i^* v_i."""
+    return np.vdot(u, v)
+
+
+def outer(u, v):
+    """|u><v|, the matrix u_i v_j^*."""
+    return np.outer(u, v.conj())
+
+
+def spectral_decomposition(A):
+    """Eigenvalues lambda_i and projectors P_i = |i><i| with A = sum_i lambda_i P_i."""
+    lam, V = np.linalg.eigh(A)
+    projectors = [outer(V[:, i], V[:, i]) for i in range(len(lam))]
+    return lam, projectors
+
+
+if __name__ == "__main__":
+    np.set_printoptions(precision=4, suppress=True)
+    # --- Eq. (2.innerex): a two-component complex vector and its norm
+    x = ket(1 - 1j, 2 + 1j)
+    print("<x|x> =", inner(x, x).real, " normalised x =", x / np.sqrt(inner(x, x).real))
+    # --- Eq. (2.innerasym): <x|y>^* = <y|x>
+    x = ket(-1, 2j, 1); y = ket(1, 0, 1j)
+    print("<x|y> =", inner(x, y), "  <y|x> =", inner(y, x))
+    assert inner(x, y) == np.conj(inner(y, x))
+    # --- computational basis, identity as a sum of projectors, Eq. (2.completeness)
+    k0, k1 = ket(1, 0), ket(0, 1)
+    P, Q = outer(k0, k0), outer(k1, k1)
+    assert np.allclose(P + Q, np.eye(2)) and np.allclose(P @ P, P) and np.allclose(P @ Q, 0)
+    print("P = |0><0| =\n", P.real, "\nQ = |1><1| =\n", Q.real)
+    # --- projecting a superposition, Eq. (2.project)
+    alpha, beta = 0.6, 0.8j
+    psi = alpha * k0 + beta * k1
+    print("P|psi> =", P @ psi, "  Q|psi> =", Q @ psi)
+    # --- the density matrix of the superposition, Eq. (2.rhopure2)
+    rho = outer(psi, psi)
+    print("rho = |psi><psi| =\n", rho, "\ntrace =", np.trace(rho).real)
+    # --- Pauli algebra, Eqs. (2.pauliprod)-(2.paulicomm)
+    for a, b, c in (("X", "Y", "Z"), ("Y", "Z", "X"), ("Z", "X", "Y")):
+        assert np.allclose(G.commutator(G.PAULIS[a], G.PAULIS[b]), 2j * G.PAULIS[c])
+        assert np.allclose(G.PAULIS[a] @ G.PAULIS[b], 1j * G.PAULIS[c])
+    print("[X,Y] = 2iZ, XY = iZ, and cyclic permutations: verified")
+    # --- Hadamard on the basis states, Eq. (2.hadamardaction)
+    print("H|0> =", G.H @ k0, "  H|1> =", G.H @ k1)
+    # --- tensor products, Eqs. (2.tensorstate)-(2.tensormatrix)
+    for a, ka in ((0, k0), (1, k1)):
+        for b, kb in ((0, k0), (1, k1)):
+            print(f"|{a}{b}> =", np.kron(ka, kb).real)
+    print("|000> =", np.kron(np.kron(k0, k0), k0).real)
+    A = np.array([[1, 2], [3, 4]]); B = np.array([[0, 1], [1, 0]])
+    print("A (x) B =\n", np.kron(A, B))
+    # --- spectral decomposition of a Hermitian matrix, Eq. (2.spectral)
+    Aherm = np.array([[2, 1 - 1j], [1 + 1j, 3]])
+    lam, projs = spectral_decomposition(Aherm)
+    rebuilt = sum(l * Pj for l, Pj in zip(lam, projs))
+    assert np.allclose(rebuilt, Aherm)
+    print("eigenvalues of A:", lam)
+    for l, Pj in zip(lam, projs):
+        print(f"lambda = {l:.4f}, projector\n", Pj)
+    # a function of the operator through its spectral decomposition
+    expA = sum(np.exp(l) * Pj for l, Pj in zip(lam, projs))
+    from math import factorial
+    series = sum(np.linalg.matrix_power(Aherm, k) / factorial(k) for k in range(40))
+    assert np.allclose(expA, series)
+    print("exp(A) from the spectral decomposition agrees with the power series")