Quantum computing, read from a simulator's source

Co-authored by Claude · Version 1.0.1, changed in rev 170

A rectangular tank of water seen from above, with two sets of ripples spreading from a pair of dippers and crossing, and five people leaning in around its rim.
Ripples from two sources crossing in a tank of water, the everyday picture of interference that the essay's simulator reproduces in numbers.

A quantum computer is a machine that obeys a short list of mathematical rules. Building one is among the hardest engineering problems anyone has attempted. Writing the rules down is much easier: they fit in a few dozen lines of code. The quantum library for DKE, the Deterministic Knowledge Engine, is those lines, written in DKE Python, and this essay teaches the subject by reading them.

The library is a simulator. There are no qubits anywhere in it, no cryostat and no laser. It runs on an ordinary processor, doing ordinary floating-point arithmetic on a list of numbers, and every result printed below was computed that way. A simulator works out what an ideal quantum computer would do, and it pays for that in memory that doubles with every qubit. It gives no quantum speedup of any kind. What it gives instead is close to the opposite of a real machine's character: you can look at everything, repeat anything, and get the same bits every time.

A simulator is also where the subject began. At a conference at MIT in 1981, published the next year as "Simulating Physics with Computers", Richard Feynman pointed out that a classical computer simulating quantum physics pays the cost this essay meets in section 7, and proposed computers built from quantum parts instead. Paul Benioff had described a quantum Turing machine the year before, and Yuri Manin made Feynman's suggestion independently. Two algorithms then gave the field its reason to try: Peter Shor's for factoring (1994), which on a large enough error-corrected machine would break RSA encryption, and Lov Grover's for search (1996).

That makes it a good place to learn. Each idea below appears as a line of source you can read and a result you can reproduce. Every program runs as shown once the library's source file has been compiled into a store, and every output is copied from a run.

1. The state is a list of numbers

Three classical bits hold one of eight patterns, 000 through 111. To describe three qubits you need one complex number for each of those eight patterns. That number is called an amplitude, and the list of them is the state. The library starts every computation from this:

def zero(n: int) -> list<complex>:
    psi = [0j] * (2 ** n)
    psi[0] = 1 + 0j
    return psi

A position in the list stands for a bit pattern. The library's convention is that qubit q is bit q of the index, so qubit 0 is the lowest bit. For two qubits:

Index   Qubit 1   Qubit 0   Written as
  0        0         0        |00⟩
  1        0         1        |01⟩
  2        1         0        |10⟩
  3        1         1        |11⟩

The bracket notation looks exotic and is only a label for a position in this list. zero(2) is amplitude 1 at position 0 and amplitude 0 everywhere else, which is the state written in the first row.

Two rules govern the list. The squared magnitudes of the amplitudes always add up to 1. And when the qubits are measured, pattern i appears with probability equal to the squared magnitude of amplitude i. Almost everything else in the subject concerns how the list changes between the start and that measurement.

Notice how the length grows: 2 qubits take 4 numbers, 10 take 1,024, and 30 take more than a billion. Hold on to that. It is why quantum computers are interesting and why simulators hit a wall.

2. A gate is a small matrix applied in pairs

A gate replaces the list with a new one. Here is the Hadamard gate, the most used gate in the subject, exactly as the library writes it:

def h(psi: list<complex>, q: int) -> list<complex>:
    s = math.sqrt(0.5)
    return vector.apply(psi, [[s, s], [s, -s]], [q])

Every gate in the library has this shape: a matrix, and one call to vector.apply. The published language specification defines that function as a short Python function. For one qubit it comes down to this:

from fractions import Fraction

def exact(a, x, b, y):
    return float(Fraction(a) * Fraction(x) + Fraction(b) * Fraction(y))

def apply_one(xs, m, q):
    out = list(xs)
    b = 1 << q
    for i in range(len(xs)):
        if i & b == 0:
            a0, a1 = xs[i], xs[i | b]
            out[i]     = exact(m[0][0], a0, m[0][1], a1)
            out[i | b] = exact(m[1][0], a0, m[1][1], a1)
    return out

In words: walk through the list in pairs of positions that differ only in bit q, and multiply each pair by the 2×2 matrix. Each new amplitude is its exact value rounded once, which exact computes with fractions for real numbers; a complex amplitude is rounded the same way, one part at a time. Which positions pair up depends on the qubit, as a list holding 1, 2, 3 and 4 shows:

import vector
import math
s = math.sqrt(0.5)
print(vector.apply([1.0, 2.0, 3.0, 4.0], [[s, s], [s, -s]], [0]))
print(vector.apply([1.0, 2.0, 3.0, 4.0], [[s, s], [s, -s]], [1]))

This prints:

[2.121320343559643, -0.7071067811865476, 4.949747468305833, -0.7071067811865476]
[2.8284271247461903, 4.242640687119286, -1.4142135623730951, -1.4142135623730951]

On qubit 0 the pairs are (1, 2) and (3, 4). On qubit 1 they are (1, 3) and (2, 4).

This is the first thing worth noticing about quantum computing. A gate on one qubit out of fifty is described by four numbers, yet it can change every amplitude in a list of about 10¹⁵. A physical machine does that by acting on one physical object. A simulator has to do all the arithmetic.

The matrices are not arbitrary. Each must be unitary, which means its conjugate transpose is its inverse. That single condition does two jobs: it keeps the squared magnitudes adding up to 1, and it makes every gate reversible. The library's adjoint function builds that inverse. Measurement, as we will see, is the one step that cannot be undone.

The gates differ only in what they do to a pair of amplitudes (a0, a1):

GateEffect on the pairIn words
x(a1, a0)Flips the bit, the quantum NOT
z(a0, −a1)Flips the sign of the "1" half
y(−i·a1, i·a0)Both of the above, with a factor of i
h((a0 + a1)·√½, (a0 − a1)·√½)Sums and differences
s(a0, i·a1)A quarter turn of phase
t(a0, w·a1), with w an eighth of a turnAn eighth turn of phase
sxsee the source√x, a square root of x: applied twice, it gives x
rx, ry, rzsee the sourceContinuous rotations by an angle you choose

3. Interference is what makes this different from chance

It is tempting to read the Hadamard gate as a coin flip, since one application turns a definite 0 into a 50/50 result. A real coin flipped twice is still 50/50. Apply Hadamard twice and watch what happens:

import quantum
one = quantum.h(quantum.zero(1), 0)
two = quantum.h(one, 0)
print(one)
print(two)
print(quantum.probabilities(two))

This prints:

[(0.7071067811865476+0j), (0.7071067811865476+0j)]
[(1.0000000000000002+0j), 0j]
[1.0000000000000004, 0.0]

The second application returns the qubit to a certain 0. (Set aside the stray digits at the end of those numbers until section 9.) Follow the arithmetic for position 1. It receives s·s from one amplitude and (−s)·s from the other, and these add to exactly zero. Two ways of reaching the result 1 cancelled.

This is interference, and it is the central fact of the subject. Probabilities are never negative, so they can only pile up. Amplitudes have signs, and more generally phases, so they can cancel. "Superposition" is often described as a qubit being 0 and 1 at the same time, and that picture hides the part that matters: a superposition is a list of amplitudes that can later be made to cancel or reinforce. Every quantum algorithm is a choreography of those cancellations.

The sign carries real information even when no probability shows it:

import quantum
plus = quantum.h(quantum.zero(1), 0)
minus = quantum.z(plus, 0)
print(quantum.probabilities(plus))
print(quantum.probabilities(minus))
print(quantum.probabilities(quantum.h(plus, 0)))
print(quantum.probabilities(quantum.h(minus, 0)))

This prints:

[0.5000000000000001, 0.5000000000000001]
[0.5000000000000001, 0.5000000000000001]
[1.0000000000000004, 0.0]
[0.0, 1.0000000000000004]

plus and minus give identical statistics if measured straight away. One more Hadamard turns the hidden sign into a certain, opposite answer. That difference between amplitudes is called relative phase. It need not be a plain sign, which is why amplitudes are complex and why the library has s, t and rz:

import math
import quantum
psi = quantum.h(quantum.zero(1), 0)
tilted = quantum.rz(psi, 0, math.pi / 2)
print(tilted)
print(quantum.probabilities(tilted))

This prints:

[(0.5000000000000001-0.5j), (0.5000000000000001+0.5j)]
[0.5000000000000001, 0.5000000000000001]

There is one kind of phase that truly means nothing. Multiply the entire list by the same unit complex number, and no measurement can ever detect it, because every later interference is unchanged. That is global phase. Two lists that differ only by it are the same physical state. The distinction will matter again in section 9, when we ask whether two circuits are the same.

4. Entanglement is why the whole list is needed

A two-qubit gate is the same idea with a 4×4 matrix along two axes. The first qubit named is the high bit of the matrix's rows and columns, so cnot's matrix reads the pattern (control, target) in the order 00, 01, 10, 11, and swaps the last two:

def cnot(psi: list<complex>, c: int, tg: int) -> list<complex>:
    return vector.apply(psi, [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]],
                        [c, tg])

CNOT flips the target exactly when the control is 1. Apply it after a Hadamard and you get the most famous state in the subject:

import quantum
psi = quantum.zero(2)
print(psi)
psi = quantum.h(psi, 0)
print(psi)
psi = quantum.cnot(psi, 0, 1)
print(psi)
print(quantum.probabilities(psi))

This prints:

[(1+0j), 0j, 0j, 0j]
[(0.7071067811865476+0j), (0.7071067811865476+0j), 0j, 0j]
[(0.7071067811865476+0j), 0j, 0j, (0.7071067811865476+0j)]
[0.5000000000000001, 0.0, 0.0, 0.5000000000000001]

This is a Bell pair: amplitude only on the patterns 00 and 11. Now ask what state qubit 0 is in on its own. Two independent qubits with amplitudes (a0, a1) and (b0, b1) would make the list (a0·b0, a1·b0, a0·b1, a1·b1). To match the Bell pair, a0·b0 and a1·b1 must both be nonzero, while a1·b0 must be zero. No numbers do that. The qubits have no separate states. The only description that exists is the joint one, and that is what entanglement means.

Here is the consequence for a simulator. If qubits were never entangled, you could store two numbers per qubit, and 100 qubits would take 200 numbers. Entanglement is what forces the full list of 2ⁿ amplitudes in the general case. One-qubit gates alone can never create it. It takes a gate like CNOT.

5. Measurement is where the list collapses

Measurement is short enough to read in full:

def measure(psi: list<complex>, q: int, u: real) -> tuple<int, list<complex>>:
    p0, p1 = weights(psi, q)
    bit = 0 if u < p0 else 1
    keep = p0 if bit == 0 else p1
    return (bit, collapsed(psi, q, bit, keep ** -0.5))

weights adds the squared magnitudes of the amplitudes where qubit q is 0, and separately of those where it is 1. That is the Born rule. The value u picks an outcome. collapsed then sets every amplitude inconsistent with that outcome to zero and rescales the rest so they add up to 1 again. The discarded amplitudes are gone, which is why this step cannot be undone.

The unusual part is u. A real measurement draws its randomness from nature. A DKE run may not draw any, because the language's contract is that the same source, the same store and the same arguments give a byte-identical result. So the outcome's randomness comes in as an argument, and the module supplies a seeded generator (Park and Miller's minimal standard) to produce it. Entropy enters in one place, the seed, which you draw outside the run and can record. These numbers are pseudo-random. They are good for repeatable tests and modest sampling, with a period of about 2.1 billion, and they are no substitute for entropy in anything that must be unpredictable.

A real device returns one bit pattern per run, so it runs the circuit many times, and each repetition is called a shot. Here are a thousand shots of the Bell pair, with the four counts in index order:

import quantum
psi = quantum.zero(2)
psi = quantum.h(psi, 0)
psi = quantum.cnot(psi, 0, 1)
counts = [0, 0, 0, 0]
u = 20261016
for shot in range(1000):
    a, after = quantum.measure(psi, 0, quantum.lcg_uniform(u))
    u = quantum.lcg_next(u)
    b, final = quantum.measure(after, 1, quantum.lcg_uniform(u))
    u = quantum.lcg_next(u)
    i = a + 2 * b
    counts[i] = counts[i] + 1
print(counts)

This prints:

[525, 0, 0, 475]

Neither qubit's result can be predicted, yet the two always agree. Replace the CNOT with a second Hadamard, on qubit 1, which leaves the qubits independent, and all four patterns appear about equally:

[253, 251, 272, 224]

This is also the place to be clear about what a simulator can do and nature forbids. probabilities(psi) reads the whole distribution exactly, with no shots at all, where a real device can only estimate it. b = a copies a state, which the no-cloning theorem says no physical process can do for an unknown state. adjoint runs a circuit backwards after the fact. Nothing decoheres. These privileges make a simulator a fine instrument for study, and a program that leans on them is a program about a quantum computer rather than one that could run on a quantum computer.

6. Two small algorithms

Teleportation moves a state from one qubit to another using a shared Bell pair and two classical bits. Qubit 0 holds a message, an arbitrary state made here with two rotations. Qubits 1 and 2 share a Bell pair. The sender entangles the message with their half and measures both of their qubits. The receiver repairs qubit 2 according to the two bits. This program runs all four possible outcomes by choosing u values that force each one:

import quantum
message = quantum.rz(quantum.ry(quantum.zero(1), 0, 0.7), 0, 0.3)
print(message)
for u0, u1 in [(0.25, 0.25), (0.25, 0.75), (0.75, 0.25), (0.75, 0.75)]:
    psi = quantum.zero(3)
    psi = quantum.ry(psi, 0, 0.7)
    psi = quantum.rz(psi, 0, 0.3)
    psi = quantum.h(psi, 1)
    psi = quantum.cnot(psi, 1, 2)
    psi = quantum.cnot(psi, 0, 1)
    psi = quantum.h(psi, 0)
    m0, psi = quantum.measure(psi, 0, u0)
    m1, psi = quantum.measure(psi, 1, u1)
    if m1 == 1:
        psi = quantum.x(psi, 2)
    if m0 == 1:
        psi = quantum.z(psi, 2)
    base = m0 + 2 * m1
    print([m0, m1, psi[base], psi[base + 4]])

This prints:

[(0.9288245698658071-0.1403781039045709j), (0.3390474346996321+0.051242007975434455j)]
[0, 0, (0.9288245698658072-0.14037810390457092j), (0.33904743469963217+0.05124200797543446j)]
[0, 1, (0.9288245698658072-0.14037810390457092j), (0.33904743469963217+0.05124200797543446j)]
[1, 0, (0.9288245698658072-0.14037810390457092j), (0.33904743469963217+0.05124200797543446j)]
[1, 1, (0.9288245698658072-0.14037810390457092j), (0.33904743469963217+0.05124200797543446j)]

The first line is the message. The other four are what qubit 2 holds after each outcome, the same state to within rounding, every time. Three facts come out of this. Nothing was copied, because measuring destroyed the original. Without the two classical bits the receiver cannot make the repair, so nothing travels faster than light. And the program branches on a measured bit with an ordinary if, which is what real hardware calls mid-circuit feed-forward.

Grover's search finds a marked item among N in about √N queries. With four items it needs one query, and it shows the choreography plainly:

import quantum
psi = quantum.zero(2)
psi = quantum.h(psi, 0)
psi = quantum.h(psi, 1)
print(psi)
psi = quantum.cz(psi, 0, 1)
print(psi)
psi = quantum.h(psi, 0)
psi = quantum.h(psi, 1)
psi = quantum.z(psi, 0)
psi = quantum.z(psi, 1)
psi = quantum.cz(psi, 0, 1)
psi = quantum.h(psi, 0)
psi = quantum.h(psi, 1)
print(psi)
print(quantum.probabilities(psi))

This prints:

[(0.5000000000000001+0j), (0.5000000000000001+0j), (0.5000000000000001+0j), (0.5000000000000001+0j)]
[(0.5000000000000001+0j), (0.5000000000000001+0j), (0.5000000000000001+0j), (-0.5000000000000001+0j)]
[0j, 0j, 0j, (1.0000000000000002+0j)]
[0.0, 0.0, 0.0, 1.0000000000000004]

The first line is every item at equal amplitude. The second is the query: cz marks item 3 by flipping its sign, a relative phase that no measurement would reveal. The remaining gates reflect every amplitude about the average, which makes the three unmarked ones cancel and the marked one reach certainty.

Eight items need more than the library's named gates, since marking one takes a three-qubit gate. That costs nothing extra, because vector.apply accepts a matrix of any matching size on any number of axes. Here it runs a three-qubit gate written as an 8×8 matrix, for three rounds:

import quantum
import vector
ccz = [[1 if r == c else 0 for c in range(8)] for r in range(8)]
ccz[7][7] = -1
psi = quantum.zero(3)
for q in range(3):
    psi = quantum.h(psi, q)
for iteration in range(3):
    psi = vector.apply(psi, ccz, [0, 1, 2])
    for q in range(3):
        psi = quantum.h(psi, q)
    for q in range(3):
        psi = quantum.x(psi, q)
    psi = vector.apply(psi, ccz, [0, 1, 2])
    for q in range(3):
        psi = quantum.x(psi, q)
    for q in range(3):
        psi = quantum.h(psi, q)
    print(quantum.probabilities(psi)[7])

This prints:

0.7812500000000007
0.9453125000000019
0.3300781250000011

The chance of finding the marked item climbs to about 95% and then, on the third round, falls to 33%. Grover's algorithm is a rotation, and it can overshoot. A quantum algorithm is not a search that keeps improving. It is interference aimed at one moment, and you have to stop at that moment.

Note what this does and does not show. Grover's speedup counts queries on a quantum computer. On this simulator, each query still walks the whole list in ordinary time.

7. The exponential wall, and a door in it

At 16 bytes per amplitude, which is how the library keeps a state, the wall is easy to compute:

QubitsAmplitudesMemory for one state
19524,2888 MiB
30about 1.07 billion16 GiB
40about 1.1 trillion16 TiB
50about 1.1 quadrillion16 PiB
300about 2×10⁹⁰More numbers than atoms in the observable universe

Each qubit doubles the cost. The largest published state-vector simulation we know of reached 50 qubits on the JUPITER supercomputer, and it needed compression to fit in about 2 petabytes. This wall is the reason a quantum computer of a few hundred good qubits could, in principle, do what no classical machine can.

The same wall decides arguments about who got there first. In 2019 Google reported that its 53-qubit Sycamore processor had sampled a random circuit in about 200 seconds, a task it estimated would take a supercomputer 10,000 years. IBM answered at once with an estimate of about two and a half days, holding the state on disk, and in 2022 Pan, Chen and Zhang solved the sampling problem in about 15 hours on 512 GPUs with tensor networks, which do not list amplitudes at all. The line between what a quantum machine can do and what a classical one can follows moves whenever someone finds a better way to simulate.

So it is a surprise that some large, heavily entangled circuits are easy. The Gottesman–Knill theorem says that circuits built only from h, s and cnot (the Clifford gates), with measurement, can be simulated efficiently on a classical computer. The library has this second simulator too, following Aaronson and Gottesman. They released their own program, CHP, with that paper in 2004, and Craig Gidney's Stim (2021) is, in the words of its title, "a fast stabilizer circuit simulator".

The idea is to stop listing amplitudes and instead record what leaves the state unchanged. The 0 state is the one that the Pauli Z operator leaves alone. The state a Hadamard makes from it is the one that X leaves alone. Clifford gates have the special property that they turn Pauli operators into other Pauli operators, so it is enough to keep a table of bits saying which ones, and that table is the tableau. Each row is a Pauli product: an X bit and a Z bit per qubit, and a sign.

import quantum
print(quantum.tableau(1))
print(quantum.tab_h(quantum.tableau(1), 0))

This prints:

[[1, 0, 0], [0, 1, 0], [0, 0, 0]]
[[0, 1, 0], [1, 0, 0], [0, 0, 0]]

In each row the columns are X bit, Z bit and sign. The middle row is the one that describes the state: at first it is Z, and after the Hadamard it is X. The first row keeps track of bookkeeping needed for measurement, and the last is scratch space. A Hadamard simply swaps the X and Z columns, and fixes up a sign:

def tab_h(tab: list<list<int>>, a: int) -> list<list<int>>:
    n = (len(tab) - 1) // 2
    out = [list(row) for row in tab]
    for i in range(2 * n):
        xa = tab[i][a]
        za = tab[i][n + a]
        out[i][2 * n] = tab[i][2 * n] ^ (xa & za)
        out[i][a] = za
        out[i][n + a] = xa
    return out

The table has 2n + 1 rows of 2n + 1 bits, so it grows with the square of n where the state vector doubles. Here is a Bell pair's big sibling, 50 qubits all entangled together, counting how many of the 50 measured bits came out 1:

import quantum
n = 50
tab = quantum.tableau(n)
tab = quantum.tab_h(tab, 0)
for q in range(1, n):
    tab = quantum.tab_cnot(tab, 0, q)
u = 1500000000
bits = []
for q in range(n):
    b, tab = quantum.tab_measure(tab, q, quantum.lcg_uniform(u))
    u = quantum.lcg_next(u)
    bits.append(b)
print(sum(bits))

This prints:

50

All 50 bits came out 1. Starting from the seed 20261016 instead, all 50 come out 0. As a state vector this needs 16 PiB. As a tableau it is 10,201 small integers, and it runs in seconds. The first measurement is random, and the other 49 are already determined, which the tableau knows without reading u at all.

This shows that entanglement alone does not make a quantum computation hard to simulate. What does? Try the T gate:

import quantum
tab = quantum.tab_gate(quantum.tableau(1), "t", [0])
print(tab)

The run is refused, with this message:

no Clifford gate named t

The refusal is an honest statement of the mathematics. T turns X into an equal combination of X and Y, which is not one Pauli product, so the table has no way to record it. And the Clifford gates plus T are enough to approximate any quantum computation. So the two simulators in this library stand on either side of a line that runs through the middle of the field: on one side are circuits a laptop can handle at any useful width, and on the other is the exponential wall.

8. How far the library goes

The DKE service gives each run a budget, and a program that goes over it is refused with a message saying so. Width and depth share that budget. One run holds a state vector of up to 19 qubits, with room at that width for one gate, and at 12 qubits more than 250 gates fit. The tableau goes much further: a short circuit fits past 150 qubits in one run.

Depth has an exact way out. save keeps a state in the store, at 16 bytes an amplitude, and advance runs a list of steps on a kept state. A circuit longer than one run becomes a series of calls that pass nothing but the steps:

import quantum
quantum.save("demo", quantum.zero(2))
u, bits = quantum.advance("demo", [("h", [0], 0.0), ("cnot", [0, 1], 0.0), ("measure", [0], 0.0)], 20261016)
print(u)
print(bits)
print(quantum.load("demo"))
quantum.drop("demo")

This prints the generator's next state, the bit measured, and the state that bit leaves:

914443951
[0]
[(1+0j), 0j, 0j, 0j]

A kept state reaches 20 qubits, and takes a step up to 19, one step a call at that width. One restriction to know is that advance's steps are gates and measurements only, so a circuit that branches on a measured bit, like teleportation, has to be written as a program.

Width has a costlier way out. A client can keep one state in parts across several stores. Some qubits are local to a part, and the rest say which part it is. Gates run on local qubits only, and before a gate on any other qubit the client exchanges bytes between the parts to bring that qubit local. Large distributed simulators use the same technique. The module supplies the pieces: packed and keep carry a state's bytes, and weigh and collapse split a measurement across parts. Each added qubit doubles the parts, so this buys a few qubits rather than dozens.

This is not a contest the library enters. Dedicated state-vector simulators handle about 30 qubits in 16 GB, and the largest run we know of reached 50. The library's value lies elsewhere, in the next two sections.

9. What "exact" buys, and what it does not

Now for those stray digits. A Hadamard twice printed 1.0000000000000002, because √½ has no exact binary floating-point form. Two rounded numbers multiplied and added don't quite return 1. Exact mathematics says the answer is 1. Every floating-point simulator has this. The library cannot avoid it either, and the word "exact" needs care.

What the library does give is reproducibility. Its arithmetic is correctly rounded and its runs are deterministic, so the same circuit gives the same bits on every host and every replay. The module is also checked against a separate implementation in CPython: every value is compared as text, with no tolerance, so a difference of one unit in the last place counts as a failure. The practical payoff is regression testing: if a simulated result changes, the circuit changed, because rounding noise cannot be the cause.

Reproducibility is not the same as telling whether two circuits are equivalent, meaning they do the same thing up to global phase. A compiler rewrites circuits and must preserve exactly that. The library answers it two ways, and the split between them is the Clifford line again.

For Clifford circuits a tableau is bits, so comparison involves no rounding at all. Here are three questions: is H Z H the same as X, is it the same as Z, and is H CZ H the same as CNOT?

import quantum
hzh = [("h", [0], 0.0), ("z", [0], 0.0), ("h", [0], 0.0)]
x = [("x", [0], 0.0)]
print(quantum.clifford(1, hzh) == quantum.clifford(1, x))
print(quantum.clifford(1, hzh) == quantum.clifford(1, [("z", [0], 0.0)]))
hczh = [("h", [1], 0.0), ("cz", [0, 1], 0.0), ("h", [1], 0.0)]
print(quantum.clifford(2, hczh) == quantum.clifford(2, [("cnot", [0, 1], 0.0)]))

This prints:

True
False
True

Any other circuit is compared by overlap: run one circuit forwards and the other backwards on a basis state, then read that state's amplitude. If the two circuits match, you get back where you started, up to a phase. The subtle part is that the phase must be the same one for every starting state. s and rz by a quarter turn are the same gate up to global phase, and rz by the wrong sign is not:

import math
import quantum
right = [("rz", [0], math.pi / 2)]
wrong = [("rz", [0], -math.pi / 2)]
s = [("s", [0], 0.0)]
print([quantum.overlap(1, right, s, j) for j in range(2)])
print([quantum.overlap(1, wrong, s, j) for j in range(2)])

This prints:

[(0.7071067811865476-0.7071067811865475j), (0.7071067811865475-0.7071067811865476j)]
[(0.7071067811865476+0.7071067811865475j), (-0.7071067811865475-0.7071067811865476j)]

In the first line the two answers share one phase. In the second, each answer has magnitude 1, so each looks fine alone, and only comparing them exposes the difference. That is global and relative phase from section 3, now doing practical work. This comparison does need a small tolerance, because of those stray digits.

Together these make a check on a compiler. dkepy's dke.quantum.qasm reads OpenQASM 2.0 into the step lists these functions take, so a compiled circuit can be compared with its source. A simulation needs only the qubits a circuit touches, so a circuit compiled for a large chip often fits.

10. A real machine, as facts

Everything so far describes an ideal machine that can run any gate on any qubits. Real devices are far less forgiving, and much of practical quantum computing is about the gap:

The second half of the library handles this with no amplitudes at all. A device and a circuit are recorded as facts in the store, and standing rules draw conclusions from them. Because it holds no state vector, the qubit ceiling does not bind it, and dkepy's dke.quantum.device writes a real device's calibration snapshot as these facts. One rule, read aloud, says that a gate is reversed if the device runs its operation on this pair the other way round and never this way:

def reversed_pair():
    for g in Gate, n in Native:
        if n.a == g.b and n.b == g.a and n.op == g.op and n.device == g.device and none(m in Native: m.a == g.a and m.b == g.b and m.op == g.op and m.device == g.device):
            g.reversed = True

Take the toy device from the DKE Library Tutorial, which writes its facts out: two qubits, sx on qubit 0, and cnot from 0 to 1 only, measured on stated dates. Record a three-gate circuit for it, h on 0, cnot from 1 to 0, and sx on 0, and then ask:

import quantum
print(quantum.illegal_gates("bell"))
print(quantum.after("bell_1"))
print(quantum.stale_gates("bell", datetime("2026-10-16T00:00:00Z"), duration("P7D")))

This prints:

['bell_1 not_native', 'bell_2 reversed']
['bell_2', 'bell_3']
['bell_3']

The device has no h at all, and it runs that cnot only the other way round. bell_1 must precede both later gates. And bell_3's calibration is more than seven days old on the date asked about. The store keeps no clock, so you name the moment, and the same question over the same facts gets the same answer.

The rules also know some physics about ordering. Diagonal gates, such as z, cz and rz, commute with one another, so a compiler may reorder them:

import quantum
update(quantum.Gate.d1.op, "z", "me")
update(quantum.Gate.d1.a, "q0", "me")
update(quantum.Gate.d1.step, 1, "me")
update(quantum.Gate.d2.op, "cz", "me")
update(quantum.Gate.d2.a, "q0", "me")
update(quantum.Gate.d2.b, "q1", "me")
update(quantum.Gate.d2.step, 2, "me")
update(quantum.Gate.d3.op, "h", "me")
update(quantum.Gate.d3.a, "q0", "me")
update(quantum.Gate.d3.step, 3, "me")
print(quantum.commuting("d1"))
print(quantum.after("d1"))
print(quantum.after("d2"))

This prints:

['d2']
['d3']
['d3']

This test is deliberately cautious. Two gates sharing a qubit are held in order unless both are diagonal, even where some other pair would in fact commute (two CNOTs with the same control, for instance). The rule never licenses a wrong reordering. It only declines to license some right ones.

Results get the same treatment. A Result records what produced an answer, and the digest of each circuit and device it read. While those digests still match the store, the result holds, and when a recalibrated device is loaded, results_not_holding() names the result and says why. That answers a question that matters in real experiments: which of my numbers are still about this chip?

Each conclusion also answers why, which lists the claims behind it:

import quantum
w = why(quantum.Gate.bell_2.reversed)
for c in w:
    print(str(c.path) + " = " + str(c.value))

Among the lines it prints are these:

quantum.Gate.bell_2.reversed = True
quantum.Gate.bell_2.a = 1
quantum.Gate.bell_2.b = 0
quantum.Native.toy_cx01.a = 0
quantum.Native.toy_cx01.b = 1
quantum.Gate.bell_2.op = cnot
quantum.Native.toy_cx01.op = cnot
quantum.Gate.bell_2.device = toy
quantum.Native.toy_cx01.device = toy

That is the conclusion and the eight facts the rule compares, and you can check the verdict by reading them: the gate wants 1 to 0, and the device offers 0 to 1. The full answer is longer, because it also names other fields held on the gate.

11. What the library leaves out

The list of absences is part of understanding the subject.

Closing

Quantum computing, at the level of rules, comes down to four things you have now read: zero builds a list of 2ⁿ complex numbers, vector.apply multiplies pairs of them by a small unitary matrix, probabilities squares them, and measure samples one outcome and discards the rest. Interference is what amplitudes can do that probabilities cannot. Entanglement is why the list cannot be split into smaller ones. The T gate is what takes you past the point where a tableau can follow.

Everything difficult about real quantum computers lies outside those lines: keeping a physical object in that list long enough to use it, on a chip where only some qubits can talk and every number drifts. A simulator cannot help with the physics. It can show you, exactly and repeatably, what the physics is supposed to do, and it can check whether a circuit still says what it said before the compiler touched it. Reading it is a good way to learn what all the effort is for.

Free to read

Ciaran Hughes, Joshua Isaacson, Anastasia Perry, Ranbel F. Sun and Jessica Turner, Quantum Computing for the Quantum Curious (Springer, Cham, 2021). Its chapters on qubits, quantum gates, entanglement and quantum algorithms teach the amplitudes and gates the essay works with. Read it free at Springer Link.

Read on