Math primer

The minimum quantum math to read and write CTF solvers, in numpy. Math is written in plain ASCII: |0>, <a|b>, U^dag, sqrt(2), alpha, theta.

Complex numbers

Every amplitude in this CTF is a complex number. Four operations matter: modulus, phase, the exponential form, and conjugation.

z = 0.6 + 0.8j
abs(z)        # modulus |z|            -> 1.0
np.angle(z)   # phase arg(z), radians  -> 0.9273
np.conj(z)    # conjugate z*           -> 0.6-0.8j
abs(z)**2     # |z|^2                  -> 1.0

Run:

>>> z = 0.6 + 0.8j
>>> z, abs(z), np.angle(z), np.conj(z)
((0.6+0.8j), 1.0, 0.9272952180016123, (0.6-0.8j))

e^(i*theta) is a unit-modulus complex number that lands at angle theta on the unit circle: abs(np.exp(1j*theta)) == 1 always. This is the building block of every phase gate.

Global phase does not change any measurement outcome — multiply the whole state vector by e^(i*phi) for any real phi and every |amplitude|^2 is unchanged:

psi  = np.array([1, 1], dtype=complex) / np.sqrt(2)
psi2 = np.exp(1j * 0.77) * psi
np.abs(psi)**2, np.abs(psi2)**2
# -> [0.5 0.5]   [0.5 0.5]   (identical)

Kets and bras

|0> and |1> are column vectors:

ket0 = np.array([[1], [0]], dtype=complex)   # |0>
ket1 = np.array([[0], [1]], dtype=complex)   # |1>

A bra is the conjugate transpose of a ket: <0| = (|0>)^dag, a row vector. In numpy that's .conj().T:

>>> bra0 = ket0.conj().T
>>> bra0
array([[1.-0.j, 0.-0.j]])

The inner product <a|b> is a bra-ket contraction — a single complex number, the overlap between two states. |<a|b>|^2 is the probability of measuring outcome a given the state was prepared as b (with a from an orthonormal measurement basis). Normalisation is <psi|psi> = 1:

psi = (ket0 + 1j*ket1) / np.sqrt(2)          # the |i> state
overlap = (psi.conj().T @ psi).item()
# -> (0.9999999999999998+0j)   normalisation check

phi = (ket0 - 1j*ket1) / np.sqrt(2)          # the |-i> state
(phi.conj().T @ psi).item()
# -> 0j   orthogonal to psi

Normalising an arbitrary vector — divide by its norm, verify amplitudes-squared sum to 1:

>>> raw = np.array([[3], [4j]], dtype=complex)
>>> norm = np.sqrt((raw.conj().T @ raw).item().real)
>>> norm
5.0
>>> normed = raw / norm
>>> normed.ravel(), np.sum(np.abs(normed)**2)
(array([0.6+0.j , 0. +0.8j]), 1.0000000000000002)

The single-qubit state

A general qubit state is |psi> = alpha|0> + beta|1> with |alpha|^2 + |beta|^2 = 1. Four states recur constantly in this CTF:

|+> and |-> give identical Z-basis probabilities (both 50/50) but are physically distinct — the difference is a relative phase on the |1> component, and it shows up the moment you measure in the X basis instead:

H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
plus  = (zero + one) / np.sqrt(2)
minus = (zero - one) / np.sqrt(2)

np.abs(plus)**2, np.abs(minus)**2
# -> [0.5 0.5]   [0.5 0.5]      identical in Z basis

H @ plus,  H @ minus
# -> [1.+0.j 0.+0.j]   [0.+0.j 1.+0.j]
#    H|+> = |0> with certainty; H|-> = |1> with certainty

Measurement probabilities from amplitudes directly, and how to get probabilities in a different basis: project onto that basis's states and square.

alpha, beta = 0.6, 0.8j
psi = alpha*zero + beta*one

# Z basis: read off the amplitudes
abs(alpha)**2, abs(beta)**2
# -> 0.36  0.6400000000000001

# X basis: project onto |+>, |-> first
c_plus  = np.vdot(plus, psi)    # <+|psi>
c_minus = np.vdot(minus, psi)   # <-|psi>
abs(c_plus)**2, abs(c_minus)**2
# -> 0.4999999999999999  0.4999999999999999

Expectation value of an observable: <psi|Z|psi>, which for Z reduces to P(0) - P(1):

>>> Z = np.array([[1, 0], [0, -1]], dtype=complex)
>>> np.vdot(psi, Z @ psi).real
-0.28000000000000014
>>> abs(alpha)**2 - abs(beta)**2
-0.28000000000000014

Operators: unitary, Hermitian, eigenbasis

Gates are unitary matrices (U^dag U = I) — they preserve the norm of the state and are reversible, which is why every gate has an inverse gate. Observables are Hermitian matrices (M = M^dag) — their eigenvalues are the possible measurement outcomes and their eigenvectors are the states that give that outcome with certainty.

X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
H = (1/np.sqrt(2)) * np.array([[1, 1], [1, -1]], dtype=complex)
I = np.eye(2, dtype=complex)

np.allclose(X.conj().T @ X, I)   # -> True   (X is unitary)
np.allclose(X, X.conj().T)       # -> True   (X is also Hermitian)

All four of X, Y, Z, H checked unitary; X, Y, Z checked Hermitian (H is unitary but the Pauli matrices are the ones used as observables, so they're the Hermitian check that matters).

Eigenvectors/eigenvalues of Z, X, Y — this is where |+>, |->, |i>, |-i> come from: they're the eigenbases of X and Y.

>>> np.linalg.eigh(Z)
(array([-1.,  1.]),
 array([[0.+0.j, 1.+0.j],
        [1.+0.j, 0.+0.j]]))          # eigenvectors are |1>, |0>

>>> np.linalg.eigh(X)
(array([-1.,  1.]),
 array([[-0.70710678+0.j,  0.70710678+0.j],
        [ 0.70710678+0.j,  0.70710678+0.j]]))   # |->, |+>  (up to sign)

>>> np.linalg.eigh(Y)
(array([-1.,  1.]),
 array([[-0.70710678-0.j        , -0.70710678+0.j        ],
        [ 0.        +0.70710678j,  0.        -0.70710678j]]))  # |-i>, |i>  (up to phase)

Norm preservation under a unitary — this is the operational meaning of "reversible":

>>> psi = np.array([0.6, 0.8j], dtype=complex); psi /= np.linalg.norm(psi)
>>> np.linalg.norm(psi), np.linalg.norm(H @ psi)
(1.0, 0.9999999999999999)

Tensor products and the endianness trap

An n-qubit state lives in a 2^n-dimensional space, built with the tensor (Kronecker) product: np.kron(a, b). Textbook convention: if you write |q0>|q1> with qubit 0 listed first, the vector is np.kron(state_q0, state_q1).

>>> np.kron(zero, one)   # textbook |0>|1>
array([0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j])
# index 0=|00>, 1=|01>, 2=|10>, 3=|11>

Concrete check — flip only qubit 0, then only qubit 1, and read off both the Qiskit label and the raw statevector:

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector

qc = QuantumCircuit(2); qc.x(0)          # flip qubit 0 only
Statevector.from_instruction(qc).data
# -> [0.+0.j 1.+0.j 0.+0.j 0.+0.j]
Statevector.from_instruction(qc).probabilities_dict()
# -> {'01': 1.0}          <- label reads q1=0, q0=1, qubit 0 is the rightmost char

qc2 = QuantumCircuit(2); qc2.x(1)        # flip qubit 1 only
Statevector.from_instruction(qc2).data
# -> [0.+0.j 0.+0.j 1.+0.j 0.+0.j]
Statevector.from_instruction(qc2).probabilities_dict()
# -> {'10': 1.0}

# to reproduce the first statevector by hand, kron in REVERSE index order:
np.kron(zero, one)   # kron(q1=|0>, q0=|1>)
# -> [0.+0.j 1.+0.j 0.+0.j 0.+0.j]   matches Statevector for X on qubit 0

All three outputs above were run and match exactly. When a challenge gives you a raw statevector or a bitstring and asks which qubit did what, write out the index arithmetic before trusting your intuition.

Entanglement

An entangled state is one that cannot be written as np.kron(a, b) for any single-qubit states a, b. The Bell state is the standard example:

>>> bell = (np.kron(zero, zero) + np.kron(one, one)) / np.sqrt(2)
>>> bell
array([0.70710678+0.j, 0.        +0.j, 0.        +0.j, 0.70710678+0.j])

Algebraic argument: if bell = np.kron(a, b) with a=(a0,a1), b=(b0,b1), then component-wise a0*b0=0.707, a0*b1=0, a1*b0=0, a1*b1=0.707. From a0*b0 != 0: both a0 != 0 and b0 != 0. From a1*b1 != 0: both a1 != 0 and b1 != 0. But then a0*b1 is a product of two nonzero numbers, so it can't be 0 — contradicts a0*b1 = 0. No product state reproduces the Bell state.

Numeric confirmation — the best fidelity any product state can reach against the Bell state, searched over 200k random product states, tops out at 0.5, never 1:

>>> best = 0.0
>>> for _ in range(200000):
...     a = rng.normal(size=2) + 1j*rng.normal(size=2); a /= np.linalg.norm(a)
...     b = rng.normal(size=2) + 1j*rng.normal(size=2); b /= np.linalg.norm(b)
...     fid = abs(np.vdot(np.kron(a, b), bell))**2
...     best = max(best, fid)
>>> best
0.5

Measurement correlation: on the Bell state, Z-measuring both qubits always gives matching bits, even though each individual qubit's outcome is 50/50 random:

>>> probs = np.abs(bell)**2
>>> dict(zip(["00","01","10","11"], probs.round(3)))
{'00': 0.5, '01': 0.0, '10': 0.0, '11': 0.5}

This perfect correlation across incompatible measurement bases (used with rotated measurement angles) is what lets a CHSH strategy exceed the classical bound of 0.75 win probability and approach cos^2(pi/8) = 0.8536 — see probability section.

Rotations and the Bloch sphere

Rx, Ry, Rz rotate the Bloch vector by angle theta about the X, Y, Z axis. The matrices use theta/2, not theta:

def Ry(theta):
    return np.array([[np.cos(theta/2), -np.sin(theta/2)],
                      [np.sin(theta/2),  np.cos(theta/2)]], dtype=complex)

>>> Ry(np.pi) @ zero
array([0.+0.j, 1.+0.j])          # correct: |0> -> |1>

>>> # buggy version, dropped the /2:
>>> Ry_buggy = lambda t: np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])
>>> Ry_buggy(np.pi) @ zero
array([-1.+0.j,  0.+0.j])        # wrong: rotated by 2*pi, landed back on |0>

Ry(pi/2) should put |0> exactly on the equator, i.e. at |+>:

>>> Ry(np.pi/2) @ zero
array([0.70710678+0.j, 0.70710678+0.j])   # == |+>

Measuring at a Bloch angle theta away from Z: apply Ry(-theta) first, then measure Z. Example — measure |+> at theta = pi/3:

>>> rotated = Ry(-np.pi/3) @ plus
>>> rotated
array([0.96592583+0.j, 0.25881905+0.j])
>>> abs(rotated[0])**2, abs(rotated[1])**2
(0.9330127018922194, 0.06698729810778065)

This rotate-then-measure-Z pattern is how CHSH measurement angles and any "measure at angle theta" challenge get implemented on a simulator or real backend — there is no native "measure along an arbitrary axis" instruction.

Number theory for Shor

Shor's algorithm factors N by finding the multiplicative order r of a random a coprime to N — the smallest r > 0 with a^r mod N == 1 — then using r to build candidate factors. The quantum part (order-finding via QPE) is the hard part; everything downstream is classical Euclid.

import math
N, a = 15, 7
math.gcd(a, N)   # -> 1   (must be coprime; a factor found trivially otherwise)

[pow(a, k, N) for k in range(8)]
# -> [1, 7, 4, 13, 1, 7, 4, 13]     period 4 -> r=4

r = next(k for k in range(1, N) if pow(a, k, N) == 1)
# -> r = 4

If r is even, a^(r/2) is a nontrivial square root of 1 mod N (as long as it isn't -1 mod N), and N divides (a^(r/2)-1)(a^(r/2)+1) without dividing either factor — so gcd(a^(r/2) +/- 1, N) peels off a nontrivial factor of N:

>>> half = pow(a, r // 2, N)          # a^(r/2) mod N
>>> half
4
>>> math.gcd(half - 1, N), math.gcd(half + 1, N)
(3, 5)
>>> 3 * 5
15

Second example, N=21, a=2:

>>> N, a = 21, 2
>>> r = next(k for k in range(1, N) if pow(a, k, N) == 1)
>>> r
6
>>> half = pow(a, r // 2, N); half
8
>>> math.gcd(half - 1, N), math.gcd(half + 1, N)
(7, 3)

Fermat/Euler context in one line: for gcd(a,N)=1, Euler's theorem gives a^phi(N) == 1 mod N, and the order r always divides phi(N) (checked: phi(15)=8, and r=4 divides 8). This is why QPE, which estimates the phase k/r to enough bits to recover r by continued fractions, works at all — r is bounded by N, so a phase register of O(log N) qubits suffices.

Probability: CHSH bounds and sampling noise

CHSH win probability per trial: classical local-hidden-variable strategies cap out at 0.75; quantum strategies using a Bell pair and angles spaced 22.5 degrees apart reach cos^2(pi/8):

>>> import math
>>> math.cos(math.pi/8)**2
0.8535533905932737

In expectation-value form, the CHSH statistic S = E(a,b) - E(a,b') + E(a',b) + E(a',b') relates to win probability by S = 8p - 4 for this encoding, giving the classical bound |S| <= 2 and the Tsirelson bound |S| <= 2*sqrt(2):

>>> for p in (0.75, 0.8535533905932737):
...     print(p, 8*p - 4)
0.75 2.0
0.8535533905932737 2.8284271247461903      # == 2*sqrt(2)

Binomial sampling noise: with n shots, a measured win rate has standard error sqrt(p(1-p)/n). To tell p=0.75 apart from p=0.8536 at roughly 95%/99.7% confidence (z=2 / z=3), solving z*(sigma1+sigma2)/sqrt(n) <= |p1-p2| for n:

>>> def shots_needed(p1, p2, z):
...     s1, s2 = math.sqrt(p1*(1-p1)), math.sqrt(p2*(1-p2))
...     return ((z*(s1+s2)) / abs(p1-p2))**2
>>> math.ceil(shots_needed(0.75, 0.8536, z=2.0))
231
>>> math.ceil(shots_needed(0.75, 0.8536, z=3.0))
520

Empirical check — simulate 2000 shots at the true quantum win rate and see how many standard errors it lands from the classical bound:

>>> samples = rng.binomial(1, 0.8535533905932737, size=2000)
>>> phat = samples.mean(); phat
0.854
>>> se = math.sqrt(phat*(1-phat)/2000); se
0.007889...
>>> (phat - 0.75) / se
13.17

Notation crib table