open-standards/wai/quantum-cross-simulator/tools/corpus.py

#!/usr/bin/env python3
"""The corpus: every circuit, its parameters, why it is here and what it must give.

Standard library only. Every circuit comes out of this file the same on every machine: random
circuits draw from SHA-256 in counter mode, not from a library's generator, and angles are
written exactly (see `exact.parse_angle`). `corpus/manifest.json` pins each circuit's digest, so
a change to any circuit shows as a change to the manifest.

A step is a gate {"g": name, "q": [qubits], "a": [angles]}, or {"measure": q, "bit": c},
{"reset": q}, or {"if": [c, v], "then": [gates]}. Gate names are those of OpenQASM 3's standard
library, and `exact.prims` says what each means.

An entry's `expect` says what it is held to:

  exact_state    the integer state vector (exact.py) of the circuit
  known_state    a state stated in closed form, which exact.py also reproduces
  distribution   the exact distribution of its classical bits (exact.py walks every branch)
  stabilizer     the stabilizer group, from the integer tableau (exact.py) and every engine
  amplitudes     chosen amplitudes in closed form, for circuits too wide for any state vector
  pairwise       no exact reference at this width: the engines against each other
  expectation    an expectation value under noise, in closed form
  qec            a noisy Clifford circuit's detector error model and detection rates
"""
import functools
import hashlib
import json
import os
import sys
from fractions import Fraction

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import exact  # noqa: E402


class Stream:
    """Integers from SHA-256 in counter mode: the same on every machine and every version."""

    def __init__(self, label):
        self.label = label.encode()
        self.counter = 0
        self.pool = b""

    def _more(self):
        self.pool += hashlib.sha256(self.label + b"/" + str(self.counter).encode()).digest()
        self.counter += 1

    def below(self, n):
        """Uniform in 0..n-1, by rejection."""
        bits = max(n - 1, 1).bit_length()
        size = (bits + 7) // 8
        mask = (1 << bits) - 1
        while True:
            while len(self.pool) < size:
                self._more()
            v = int.from_bytes(self.pool[:size], "big") & mask
            self.pool = self.pool[size:]
            if v < n:
                return v

    def choice(self, xs):
        return xs[self.below(len(xs))]


def g(name, *qubits, a=None):
    s = {"g": name, "q": list(qubits)}
    if a:
        s["a"] = list(a)
    return s


def t(m, k):
    """m/2^k of a full turn."""
    return "t:%d/%d" % (m, k)


# ------------------------------------------------------------------ programs as text


def qasm_angle(a):
    kind, v = exact.parse_angle(a)
    if kind == "r":
        return a[2:]
    f = v * 2
    n, d = f.numerator, f.denominator
    if n == 0:
        return "0"
    sign = "-" if n < 0 else ""
    n = abs(n)
    if n == 1 and d == 1:
        return sign + "pi"
    if d == 1:
        return "%spi*%d" % (sign, n)
    if n == 1:
        return "%spi/%d" % (sign, d)
    return "%spi*%d/%d" % (sign, n, d)


# The standard library's names, and the modifiers for the rest.
_PLAIN = {"id", "x", "y", "z", "h", "s", "sdg", "t", "tdg", "sx", "p", "rz", "rx", "ry", "cx", "cy", "cz", "ch", "cp", "crx", "cry", "crz", "cu", "swap", "cswap", "ccx"}
_DEFS = {
    "rzz": "gate rzz(theta) a, b { cx a, b; rz(theta) b; cx a, b; }",
    "rxx": "gate rxx(theta) a, b { h a; h b; cx a, b; rz(theta) b; cx a, b; h a; h b; }",
    "ryy": "gate ryy(theta) a, b { rx(pi/2) a; rx(pi/2) b; cx a, b; rz(theta) b; cx a, b; rx(-pi/2) a; rx(-pi/2) b; }",
}


def _gate_text(s):
    name, qs, args = s["g"], s["q"], s.get("a", [])
    params = "(" + ", ".join(qasm_angle(x) for x in args) + ")" if args else ""
    q = ", ".join("q[%d]" % x for x in qs)
    if name in _PLAIN or name in _DEFS:
        return "%s%s %s;" % (name, params, q)
    if name == "u":
        return "U%s %s;" % (params, q)
    if name == "sxdg":
        return "inv @ sx %s;" % q
    if name == "gphase":
        if qs:
            return "ctrl(%d) @ gphase%s %s;" % (len(qs), params, q)
        return "gphase%s;" % params
    if name == "negcx":
        return "negctrl @ x %s;" % q
    form = exact.controlled_form(name)
    if form:
        base, k = form
        return "ctrl(%d) @ %s%s %s;" % (k, base, params, q)
    raise ValueError("no text for gate %s" % name)


def qasm_of(entry):
    """The entry's program as OpenQASM 3: the text the engine reads and the text published."""
    out = ["OPENQASM 3.0;", 'include "stdgates.inc";']
    names = set()
    for s in entry["steps"]:
        for x in s.get("then", [s]):
            if "g" in x:
                names.add(x["g"])
    for name in sorted(names):
        if name in _DEFS:
            out.append(_DEFS[name])
    out.append("qubit[%d] q;" % entry["n"])
    if entry.get("bits"):
        out.append("bit[%d] c;" % entry["bits"])
    for s in entry["steps"]:
        if "g" in s:
            out.append(_gate_text(s))
        elif "measure" in s:
            out.append("c[%d] = measure q[%d];" % (s["bit"], s["measure"]))
        elif "reset" in s:
            out.append("reset q[%d];" % s["reset"])
        elif "if" in s:
            body = " ".join(_gate_text(x) for x in s["then"])
            out.append("if (c[%d] == %d) { %s }" % (s["if"][0], s["if"][1], body))
    return "\n".join(out) + "\n"


# ------------------------------------------------------------------ circuit builders


def qft_steps(n):
    """The QFT with qubit 0 the lowest bit of the index: |x> to 2^(-n/2) Σ e^{2πi xy/2^n} |y>."""
    steps = []
    for j in range(n - 1, -1, -1):
        steps.append(g("h", j))
        for k in range(j - 1, -1, -1):
            steps.append(g("cp", k, j, a=[t(1, j - k + 1)]))
    for i in range(n // 2):
        steps.append(g("swap", i, n - 1 - i))
    return steps


def inverse_steps(steps):
    out = []
    dag = {"h": "h", "x": "x", "y": "y", "z": "z", "s": "sdg", "sdg": "s", "t": "tdg", "tdg": "t", "cx": "cx", "cy": "cy", "cz": "cz", "swap": "swap", "sx": "sxdg", "sxdg": "sx", "ccx": "ccx"}
    for s in reversed(steps):
        name = s["g"]
        if name in dag:
            out.append({"g": dag[name], "q": s["q"]})
        elif name in ("p", "rz", "rx", "ry", "cp", "crz", "crx", "cry", "rzz", "rxx", "ryy"):
            out.append({"g": name, "q": s["q"], "a": [exact._neg(s["a"][0])]})
        else:
            raise ValueError("no inverse for %s" % name)
    return out


def basis_prep(n, x):
    return [g("x", q) for q in range(n) if (x >> q) & 1]


def grover_steps(n, marked, iterations):
    steps = [g("h", q) for q in range(n)]
    all_q = list(range(n))
    cz = "c%dz" % (n - 1) if n > 3 else {2: "cz", 3: "ccz"}.get(n, "z")
    for _ in range(iterations):
        flip = [g("x", q) for q in range(n) if not (marked >> q) & 1]
        steps += flip + [{"g": cz, "q": all_q}] + flip
        steps += [g("h", q) for q in range(n)] + [g("x", q) for q in range(n)]
        steps += [{"g": cz, "q": all_q}]
        steps += [g("x", q) for q in range(n)] + [g("h", q) for q in range(n)]
    return steps


def qpe_steps(bits, phase_turns):
    """Phase estimation of the eigenvalue e^{2πi φ} of p(2πφ) on qubit `bits`, prepared in |1>.
    Counting qubits 0..bits-1; the outcome y estimates φ = y/2^bits, read as an integer."""
    num, den = phase_turns
    steps = [g("x", bits)] + [g("h", q) for q in range(bits)]
    for j in range(bits):
        # 2^j applications of the unitary, as one controlled phase.
        m = (num * (1 << j)) % den
        angle = "q:%d/%d" % (m, den) if den & (den - 1) else t(m, den.bit_length() - 1)
        steps.append(g("cp", j, bits, a=[angle]))
    steps += inverse_steps_qft(bits)
    return steps


def inverse_steps_qft(n):
    return inverse_steps(qft_steps(n))


def random_clifford_t(n, depth, seed, gateset=("h", "s", "t", "sdg", "tdg", "x")):
    r = Stream("clifford_t/%d/%d/%d" % (n, depth, seed))
    steps = []
    for _ in range(depth):
        for q in range(n):
            if r.below(3):
                steps.append(g(r.choice(gateset), q))
        order = list(range(n))
        for i in range(n - 1, 0, -1):
            j = r.below(i + 1)
            order[i], order[j] = order[j], order[i]
        for i in range(0, n - 1, 2):
            steps.append(g(r.choice(("cx", "cz")), order[i], order[i + 1]))
    return steps


def random_clifford(n, depth, seed):
    return random_clifford_t(n, depth, seed, gateset=("h", "s", "sdg", "x", "y", "z", "sx"))


def random_dyadic(n, depth, seed, k):
    r = Stream("dyadic/%d/%d/%d/%d" % (n, depth, seed, k))
    steps = []
    for layer in range(depth):
        for q in range(n):
            steps.append(g(r.choice(("rx", "ry", "rz")), q, a=[t(r.below(1 << k), k)]))
        for q in range(layer % 2, n - 1, 2):
            steps.append(g("cz", q, q + 1))
    return steps


def random_rotations(n, depth, seed):
    """U gates whose angles are decimal radians, no dyadic among them, on a brickwork of cx."""
    r = Stream("rotations/%d/%d/%d" % (n, depth, seed))

    def angle():
        return "r:%d.%017d" % (r.below(3), r.below(10**17))

    steps = []
    for layer in range(depth):
        for q in range(n):
            steps.append(g("u", q, a=[angle(), angle(), angle()]))
        for q in range(layer % 2, n - 1, 2):
            steps.append(g("cx", q, q + 1))
    return steps


def ghz_steps(n):
    return [g("h", 0)] + [g("cx", q, q + 1) for q in range(n - 1)]


def distinct_product(n):
    """A product state whose every amplitude differs: any exchange of qubits shows."""
    steps = []
    for q in range(n):
        steps.append(g("ry", q, a=[t(1, q + 4)]))
        steps.append(g("rz", q, a=[t(3, q + 6)]))
    return steps


def measure_all(n):
    return [{"measure": q, "bit": q} for q in range(n)]


# ------------------------------------------------------------------ the entries


def entries():
    out = []

    def add(id_, family, params, reason, n, steps, expect, tier="full", bits=0, engines=None, shots=0, **extra):
        e = {"id": id_, "family": family, "params": params, "reason": reason, "n": n, "bits": bits, "steps": steps, "expect": expect, "tier": tier, "shots": shots}
        if engines:
            e["engines"] = engines
        e.update(extra)
        out.append(e)

    # Bit order.
    add("bit_order_x_low", "bit_order", {"n": 5}, "An X on qubit 0 must read as index 1, the first thing that differs between engines' conventions.", 5, [g("x", 0)], "exact_state", tier="ci")
    add("bit_order_x_high", "bit_order", {"n": 5}, "The other end of the same convention: an X on the highest qubit is index 16.", 5, [g("x", 4)], "exact_state")
    add("bit_order_product", "bit_order", {"n": 6}, "Every amplitude differs, so any exchange of two qubits changes the state.", 6, distinct_product(6), "exact_state", tier="ci")
    add("bit_order_control", "bit_order", {"n": 4}, "A cx from qubit 0 to qubit 3: the control and target must not swap roles.", 4, [g("x", 0), g("cx", 0, 3)], "exact_state")
    add("bit_order_named_bits", "bit_order", {"n": 2}, "Classical bits named out of order: qubit 1 into bit 0. The counts must follow the bits, not the qubits.", 2, [g("x", 1), {"measure": 1, "bit": 0}, {"measure": 0, "bit": 1}], "distribution", bits=2, shots=4000, tier="ci")
    add("empty_register", "bit_order", {"n": 3}, "A program of no gates: the state is |000>.", 3, [], "exact_state")
    add("one_qubit", "bit_order", {"n": 1}, "The smallest register.", 1, [g("h", 0)], "exact_state")
    add("top_index_n20", "bit_order", {"n": 20}, "An index at bit 19: index arithmetic past 16 bits, with no entanglement to hide it.", 20, [g("x", 19), g("h", 18)], "pairwise", tier="ci")

    # GHZ.
    for n in (2, 3, 5, 8, 12):
        add("ghz_n%d" % n, "ghz", {"n": n}, "Known answer: two amplitudes of 1/sqrt(2), every other zero. It exposes any leakage out of those two.", n, ghz_steps(n), "known_state", tier="ci" if n == 3 else "full")
    add("ghz_counts_n8", "ghz", {"n": 8}, "Sampling: only the two outcomes may appear, equally often.", 8, ghz_steps(8) + measure_all(8), "distribution", bits=8, shots=100000)

    # QFT.
    for n, x in ((3, 5), (5, 19), (8, 173), (12, 2741), (14, 10923)):
        add("qft_n%d_x%d" % (n, x), "qft", {"n": n, "x": x}, "Closed form e^{2πi xy/2^n}/2^(n/2), and every controlled phase up to 2π/2^n: the rotations a fixed-point table holds least well.", n, basis_prep(n, x) + qft_steps(n), "known_state", tier="ci" if n == 5 else "full", qft_x=x)
    add("qft_superposed_n10", "qft", {"n": 10, "seed": 11}, "A QFT over a random Clifford+T state: its phases act on every amplitude, where from a basis state most cancel.", 10, random_clifford_t(10, 4, 11) + qft_steps(10), "exact_state")
    add("qft_counts_n8", "qft", {"n": 8, "x": 173}, "Sampling a uniform distribution of 256 outcomes.", 8, basis_prep(8, 173) + qft_steps(8) + measure_all(8), "distribution", bits=8, shots=200000)
    add("qft_wide_n20", "qft", {"n": 20, "x": 699051}, "Past any integer reference here: engines against each other, and chosen amplitudes against the closed form.", 20, basis_prep(20, 699051) + qft_steps(20), "pairwise", qft_x=699051)

    # Grover.
    for n, marked, it in ((3, 5, 2), (5, 19, 4), (8, 173, 12), (10, 683, 25)):
        add("grover_n%d" % n, "grover", {"n": n, "marked": marked, "iterations": it}, "Closed form: the marked state's probability is sin²((2k+1)θ) with sin θ = 2^(-n/2). Multi-controlled gates, and interference that rounding drains.", n, grover_steps(n, marked, it), "exact_state", tier="ci" if n == 3 else "full", grover=[marked, it])
    add("grover_counts_n5", "grover", {"n": 5, "marked": 19, "iterations": 4}, "Sampling a peaked distribution.", 5, grover_steps(5, 19, 4) + measure_all(5), "distribution", bits=5, shots=100000)

    # Phase estimation.
    for bits, m in ((4, 5), (8, 173), (12, 2741)):
        add("qpe_t%d_m%d" % (bits, m), "qpe", {"bits": bits, "phase": "%d/2^%d" % (m, bits)}, "A phase of exactly m/2^t: the outcome is m with probability 1, so every other amplitude must cancel to zero.", bits + 1, qpe_steps(bits, (m, 1 << bits)), "exact_state", tier="ci" if bits == 4 else "full", qpe=[m, bits])
    for bits in (6, 10):
        add("qpe_t%d_third" % bits, "qpe", {"bits": bits, "phase": "1/3"}, "A phase that is no dyadic: the outcome's law is the closed form sin²(π2^t δ)/(2^(2t) sin²(πδ)), and the bit-exact method must refuse it rather than round it.", bits + 1, qpe_steps(bits, (1, 3)), "exact_state", qpe_third=bits)

    # Random circuits.
    for n, depth, seed in ((4, 20, 1), (8, 40, 2), (12, 40, 3)):
        add("clifford_t_n%d_d%d" % (n, depth), "random_clifford_t", {"n": n, "depth": depth, "seed": seed}, "Universal, in the gates the bit-exact method holds as table constants.", n, random_clifford_t(n, depth, seed), "exact_state", tier="ci" if n == 4 else "full")
    for n, depth, seed, k in ((6, 20, 4, 16), (10, 30, 5, 32)):
        add("dyadic_n%d_k%d" % (n, depth), "random_dyadic", {"n": n, "depth": depth, "seed": seed, "k": k}, "Rotations by m/2^k of a turn: each is many table phases in fixed point, so its error is largest here.", n, random_dyadic(n, depth, seed, k), "exact_state", tier="ci" if n == 6 else "full")
    for n, depth, seed in ((6, 20, 6), (10, 30, 7), (12, 40, 8)):
        add("rotations_n%d" % n, "random_rotations", {"n": n, "depth": depth, "seed": seed}, "Arbitrary decimal angles: double-precision engines only, and the bit-exact method must refuse.", n, random_rotations(n, depth, seed), "exact_state", tier="ci" if n == 6 else "full")
    add("rotations_counts_n6", "random_rotations", {"n": 6, "depth": 20, "seed": 6}, "Sampling a distribution with no structure.", 6, random_rotations(6, 20, 6) + measure_all(6), "distribution", bits=6, shots=200000)
    add("rotations_wide_n20", "random_rotations", {"n": 20, "depth": 10, "seed": 9}, "No integer reference at this width: engines against each other, with fusion and threads in play.", 20, random_rotations(20, 10, 9), "pairwise")

    # Mirror circuits: known answer, any width.
    for n, depth, seed in ((10, 40, 21), (16, 30, 22), (20, 20, 23)):
        forward = random_dyadic(n, depth, seed, 16)
        add("mirror_n%d" % n, "mirror", {"n": n, "depth": depth, "seed": seed}, "U then its inverse: the state is |0…0> exactly, at a width no integer reference reaches, and what remains is the engine's accumulated rounding.", n, forward + inverse_steps(forward), "mirror", tier="ci" if n == 10 else "full")

    # Small angles.
    for k in (8, 16, 24, 26, 28, 30, 32):
        add("small_angle_k%d" % k, "small_angle", {"k": k}, "h, p(2π/2^k), h: P(1) = sin²(π/2^k). Within the phase table, and at its resolution floor as k nears 32.", 1, [g("h", 0), g("p", 0, a=[t(1, k)]), g("h", 0)], "exact_state", tier="ci" if k == 32 else "full", probe=1)
    for k in (33, 40, 47, 48, 52, 60):
        add("below_grid_k%d" % k, "small_angle", {"k": k}, "An angle below the table's last entry: the bit-exact method refuses it, or reads it as nothing when it is within its reader's tolerance of 0. Whether each engine holds the phase or drops it is recorded beside its distance.", 1, [g("h", 0), g("p", 0, a=[t(1, k)]), g("h", 0)], "exact_state", tier="ci" if k == 47 else "full", probe=1)
    add("tiny_radians", "small_angle", {"radians": "1e-12"}, "A small decimal angle, not on any dyadic grid.", 1, [g("h", 0), g("rz", 0, a=["r:1e-12"]), g("h", 0)], "exact_state", probe=1)
    add("controlled_small_angle_k30", "small_angle", {"k": 30}, "The same phase under a control.", 2, [g("h", 0), g("h", 1), g("cp", 0, 1, a=[t(1, 30)])], "exact_state", probe=2)
    # The phase is the result here, so these two are held to the state with its phase (`strict`),
    # not modulo one: a phase dropped, doubled or misread as a whole turn shows.
    add("repeated_phase_k19", "small_angle", {"k": 19, "count": 262144}, "p(2π/2^19) applied 2^18 times to |1> is a half turn, -|1>. The extension's section 5 states what fixed point leaves of it: modulus 1.0000189, within 1e-5 rad of the half turn.", 1, [g("x", 0)] + [g("p", 0, a=[t(1, 19)])] * 262144, "exact_state", phase_is_result=True, extension_figure=True, engines=["wai.exact", "wai.statevector", "aer.statevector", "ddsim.statevector", "qsim.statevector"])
    add("repeated_phase_k19_quarter", "small_angle", {"k": 19, "count": 131072}, "p(2π/2^19) applied 2^17 times to |1> is a quarter turn, i|1>: the half turn is its own negative, and this one is not, so a phase applied with the wrong sign shows.", 1, [g("x", 0)] + [g("p", 0, a=[t(1, 19)])] * 131072, "exact_state", phase_is_result=True, engines=["wai.exact", "wai.statevector", "aer.statevector", "ddsim.statevector", "qsim.statevector"])
    add("h_squared_chain", "small_angle", {"count": 4096}, "2048 pairs of H on a superposition: the identity, where fixed point's rounding of 1/sqrt(2) adds up.", 1, [g("h", 0)] * 4096, "exact_state", engines=["wai.exact", "wai.statevector", "aer.statevector", "ddsim.statevector", "qsim.statevector", "quimb.dense"])

    # Gate conventions: every gate of the library on a state that shows what it does.
    gate_cases = [
        ("id", 1, []), ("x", 1, []), ("y", 1, []), ("z", 1, []), ("h", 1, []), ("s", 1, []), ("sdg", 1, []), ("t", 1, []), ("tdg", 1, []), ("sx", 1, []), ("sxdg", 1, []),
        ("p", 1, [t(3, 7)]), ("rz", 1, [t(3, 7)]), ("rx", 1, [t(3, 7)]), ("ry", 1, [t(3, 7)]), ("u", 1, [t(3, 7), t(5, 8), t(7, 9)]),
        ("cx", 2, []), ("cy", 2, []), ("cz", 2, []), ("ch", 2, []), ("csx", 2, []), ("swap", 2, []),
        ("cp", 2, [t(3, 7)]), ("crx", 2, [t(3, 7)]), ("cry", 2, [t(3, 7)]), ("crz", 2, [t(3, 7)]), ("cu", 2, [t(3, 7), t(5, 8), t(7, 9), t(1, 5)]),
        ("rzz", 2, [t(3, 7)]), ("rxx", 2, [t(3, 7)]), ("ryy", 2, [t(3, 7)]), ("negcx", 2, []),
        ("ccx", 3, []), ("ccz", 3, []), ("cswap", 3, []), ("c2p", 3, [t(3, 7)]), ("c2rz", 3, [t(3, 7)]), ("c2ry", 3, [t(3, 7)]), ("c3x", 4, []), ("c3z", 4, []),
    ]
    for name, arity, args in gate_cases:
        qs = list(range(arity))
        add("gate_%s" % name, "gate_conventions", {"gate": name}, "What this gate means, from OpenQASM 3's standard library, applied to a state with no symmetry; engines that define it otherwise show here.", arity, distinct_product(arity) + [{"g": name, "q": qs, **({"a": args} if args else {})}], "exact_state", tier="ci" if name in ("sx", "rz", "cu", "c2rz", "negcx") else "full")
    for name, args in (("rz", [t(1, 0)]), ("rx", [t(1, 0)]), ("ry", [t(1, 0)]), ("p", [t(1, 0)]), ("rz", [t(2, 0)]), ("u", [t(1, 0), t(0, 0), t(0, 0)])):
        add("full_turn_%s_%s" % (name, args[0].replace(":", "").replace("/", "_")), "gate_conventions", {"gate": name, "angle": args[0]}, "A rotation by whole turns: the identity up to a sign, and the sign is the convention. Engines that lower it through a phase gate lose the sign; a global phase changes no outcome, and under a control it would.", 1, distinct_product(1) + [{"g": name, "q": [0], "a": args}], "exact_state")
    add("gate_gphase_controlled", "gate_conventions", {"gate": "ctrl @ gphase"}, "A global phase under a control is a relative phase, so no engine may drop it.", 2, distinct_product(2) + [g("gphase", 0, a=[t(3, 7)])], "exact_state")
    add("controlled_rz_full_turn", "gate_conventions", {"gate": "crz", "angle": "1 turn"}, "crz by a full turn is diag(1, 1, -1, -1): the sign rz loses globally is here relative.", 2, [g("h", 0), g("h", 1), g("crz", 0, 1, a=[t(1, 0)])], "exact_state", tier="ci")

    # Measurement in the middle, reset and feedback.
    alpha = t(1, 5)
    add("teleport", "dynamic", {"angle": alpha}, "Teleportation: Bell measurement, then corrections on classical bits. Qubit 2 must end as qubit 0 began, whatever the outcomes.", 3, [g("ry", 0, a=[alpha]), g("h", 1), g("cx", 1, 2), g("cx", 0, 1), g("h", 0), {"measure": 0, "bit": 0}, {"measure": 1, "bit": 1}, {"if": [1, 1], "then": [g("x", 2)]}, {"if": [0, 1], "then": [g("z", 2)]}, {"measure": 2, "bit": 2}], "distribution", bits=3, shots=200000, tier="ci")
    add("reset_after_superposition", "dynamic", {}, "A reset of a qubit in superposition leaves |0> whatever it read, and leaves its Bell partner mixed.", 2, [g("h", 0), g("cx", 0, 1), {"reset": 0}, {"measure": 0, "bit": 0}, {"measure": 1, "bit": 1}], "distribution", bits=2, shots=100000, tier="ci")
    add("measure_twice", "dynamic", {}, "A second measurement in the same basis repeats the first: 00 and 11 only.", 1, [g("h", 0), {"measure": 0, "bit": 0}, {"measure": 0, "bit": 1}], "distribution", bits=2, shots=100000)
    add("feedback_flip", "dynamic", {}, "Flip a measured 1 back: the second reading is always 0.", 1, [g("h", 0), {"measure": 0, "bit": 0}, {"if": [0, 1], "then": [g("x", 0)]}, {"measure": 0, "bit": 1}], "distribution", bits=2, shots=100000)
    add("ghz_by_fusion", "dynamic", {}, "A GHZ state made from two Bell pairs by a parity measurement and a correction: Clifford, with feedback, so a stabilizer engine must also get it.", 4, [g("h", 0), g("cx", 0, 1), g("h", 2), g("cx", 2, 3), g("cx", 1, 2), {"measure": 2, "bit": 0}, {"if": [0, 1], "then": [g("x", 3)]}, {"measure": 0, "bit": 1}, {"measure": 1, "bit": 2}, {"measure": 3, "bit": 3}], "distribution", bits=4, shots=100000)
    add("repeat_until_success", "dynamic", {"rounds": 3}, "Three rounds of a non-Clifford trial on a second qubit, which is measured and then flipped back to |0> by an x conditioned on what was read (no reset instruction): the branches weigh differently, and the law is exact.", 2, _rus_steps(3), "distribution", bits=4, shots=200000)

    # Stabilizer circuits at width.
    for n in (100, 1000, 5000):
        add("clifford_ghz_n%d" % n, "clifford_large", {"n": n}, "A GHZ state across the whole register: sparse generators, past any state vector.", n, ghz_steps(n), "stabilizer", tier="ci" if n == 100 else "full")
    for n in (65, 1000, 5000):
        add("clifford_cluster_n%d" % n, "clifford_large", {"n": n}, "A line cluster state: every generator X Z neighbours, with its signs flipped by a layer of Z.", n, [g("h", q) for q in range(n)] + [g("cz", q, q + 1) for q in range(n - 1)] + [g("z", q) for q in range(0, n, 3)], "stabilizer")
    for n, depth, seed in ((63, 12, 31), (64, 12, 32), (65, 12, 33), (127, 12, 34), (128, 12, 35), (129, 12, 36), (1000, 30, 37), (2000, 20, 38)):
        add("clifford_random_n%d" % n, "clifford_large", {"n": n, "depth": depth, "seed": seed}, "Random Clifford layers, with widths on and either side of a machine word: packing bugs live at those edges.", n, random_clifford(n, depth, seed), "stabilizer", tier="ci" if n == 65 else "full")
    add("clifford_random_n10_as_state", "clifford_large", {"n": 10, "depth": 30, "seed": 39}, "A Clifford circuit small enough for a state vector: the two formalisms must describe one state.", 10, random_clifford(10, 30, 39), "exact_state")
    add("clifford_counts_n40", "clifford_large", {"n": 40, "depth": 8, "seed": 40}, "Sampling past the provider receipt's 22-qubit comparison: outcomes must lie in the affine space the stabilizers fix.", 40, random_clifford(40, 8, 40) + measure_all(40), "stabilizer_counts", bits=40, shots=20000, tier="ci")

    # Noise, in closed form.
    for p, k in (("0.01", 20), ("0.1", 7)):
        add("noise_depol1_x_chain_p%s_k%d" % (p.replace(".", "_"), k), "noise", {"channel": "depolarizing_1q", "p": p, "gates": k}, "k X gates, each followed by a depolarizing channel: <Z> = (-1)^k (1 - 4p/3)^k exactly.", 1, [g("x", 0)] * k, "expectation", noise={"depolarizing_1q": float(p)}, observable="Z0", closed=["depol1", p, k], tier="ci" if k == 7 else "full")
    add("noise_depol2_cx_chain", "noise", {"channel": "depolarizing_2q", "p": "0.02", "gates": 9}, "k cx gates, each followed by a two-qubit depolarizing channel: <Z Z> = (1 - 16p/15)^k exactly.", 2, [g("cx", 0, 1)] * 9, "expectation", noise={"depolarizing_2q": 0.02}, observable="Z0 Z1", closed=["depol2", "0.02", 9])

    # Wide circuits, by methods that reach past a state vector.
    for n in (30, 50, 100):
        add("mps_ghz_n%d" % n, "wide_mps", {"n": n}, "Bond dimension 2 at any width: the amplitudes of the two outcomes, in closed form.", n, ghz_steps(n), "amplitudes", wide_kind="ghz", tier="ci" if n == 30 else "full")
    for n, x in ((30, 715827883), (40, 366503875925)):
        add("mps_qft_n%d" % n, "wide_mps", {"n": n, "x": x}, "The QFT of a basis state is a product state, reached through modest entanglement. Chosen amplitudes in closed form.", n, basis_prep(n, x) + qft_steps(n), "amplitudes", wide_kind="qft", qft_x=x)
    # The bond asked for is not the bond reached: at 256, past what an estimate at that bond admits.
    add("mps_qft_n30_bond256", "wide_mps", {"n": 30, "x": 715827883, "max_bond": 256}, "The 30-qubit QFT at a bond far above what it reaches: a job naming the method runs under a meter of the bond it holds, and stays exact.", 30, basis_prep(30, 715827883) + qft_steps(30), "amplitudes", wide_kind="qft", qft_x=715827883, max_bond=256, only_engines=["wai.mps"])
    add("tn_shallow_n26", "wide_tn", {"n": 26, "depth": 5, "seed": 41}, "A shallow rotation brickwork: amplitudes by contraction, against matrix product states and a state vector.", 26, random_rotations(26, 5, 41), "wide_pairwise")
    add("tn_mirror_n40", "wide_tn", {"n": 40, "depth": 3, "seed": 42}, "A shallow mirror circuit, too wide for a state vector here: the amplitude of |0…0> is exactly 1.", 40, random_dyadic(40, 3, 42, 16) + inverse_steps(random_dyadic(40, 3, 42, 16)), "wide_mirror", tier="ci")

    return out


def _rus_steps(rounds):
    steps = []
    for r in range(rounds):
        steps += [g("h", 1), g("t", 1), g("cx", 0, 1), g("h", 1), g("t", 1), {"measure": 1, "bit": r}, {"if": [r, 1], "then": [g("x", 1)]}, g("h", 0), g("t", 0)]
    steps.append({"measure": 0, "bit": rounds})
    return steps


# ------------------------------------------------------------------ what the bit-exact method must do

# The bit-exact method's reader (wai-quantum/src/qasm.rs, `dyadic`) takes an angle as 2πm/2^k for
# the smallest k ≤ 32 within 10^-14 of a turn, and refuses one it cannot (`[Q500]`). This is that
# rule applied to the corpus's angles as written, exactly: each angle is "exact" (m/2^k, k ≤ 32),
# "moved" (within the tolerance of such an angle, not on it) or "refused". A gate whose lowering
# halves an angle (crz, crx, cry, cu) or takes it under modifiers is given angles of k ≤ 30 here,
# which the corpus checks, so that the rule on the written angle is the rule on every angle the
# lowering reaches.
READER_TOLERANCE = Fraction(1, 10**14)
GRID = 32
HALVING = {"crz", "crx", "cry", "cu", "c2p", "c2rz", "c2ry", "gphase"}


def _turns(a):
    kind, v = exact.parse_angle(a)
    if kind == "t":
        return v
    return v / (2 * Fraction(exact._PI, 1 << (exact.PREC + exact.GUARD)))


@functools.lru_cache(maxsize=None)
def reader_angle(a):
    """("exact" | "moved" | "refused", how far the reader moves it, in turns)."""
    t = _turns(a) % 1
    for k in range(GRID + 1):
        x = t * (1 << k)
        m = round(x)
        if abs(x - m) <= READER_TOLERANCE * (1 << k):
            return ("exact" if x == m else "moved", abs(x - m) / (1 << k))
    return ("refused", None)


def reader_class(e):
    """What the bit-exact method must do with an entry: ("exact", 0), ("moved", the sum of how
    far its angles move, in turns) or ("refused", 0)."""
    moved = Fraction(0)
    worst = "exact"
    for s in e["steps"]:
        for x in s.get("then", [s]):
            for a in x.get("a", []):
                kind, by = reader_angle(a)
                if kind == "refused":
                    return ("refused", 0)
                if kind == "moved":
                    worst = "moved"
                    moved += by
                if x["g"] in HALVING and kind == "exact" and exact.parse_angle(a)[0] == "t" and int(a[2:].partition("/")[2]) > GRID - 2:
                    raise ValueError("%s: %s at %s halves past the grid; the rule above would not hold" % (e["id"], x["g"], a))
    return (worst, moved)


# Refusals by the engine under test that a run accepts, each with the words it must give. Any
# other refusal, crash or time-out of the engine under test fails the run. "must": the engine
# must refuse (accepting fails the run); "may": a refusal with these words is recorded as one,
# and an answer is checked like any other.
QEC_REFUSALS = [
    {"entry": "qec_refused_channel_with_no_independent_form", "engine": "wai.qec", "must": "refuse", "text": "no form as independent", "reason": "a Pauli channel no independent flips can give is refused, never approximated"},
    {"entry": "qec_refused_random_detector", "engine": "wai.qec", "must": "refuse", "text": "is not deterministic", "reason": "a detector that is random without noise is no detector"},
]


def expected_refusals():
    """QEC_REFUSALS, and every entry the bit-exact method must refuse by the reader's rule."""
    out = []
    for e in entries():
        runs_exact = e["expect"] in ("exact_state", "known_state", "mirror") and "wai.exact" in e.get("engines", ["wai.exact"])
        if runs_exact and reader_class(e)[0] == "refused":
            out.append({"entry": e["id"], "engine": "wai.exact", "must": "refuse", "text": "[Q500]", "reason": "an angle no grid of 2^32 holds within 10^-14 of a turn"})
    return out + QEC_REFUSALS


# ------------------------------------------------------------------ the error-correction circuits


def qec_entries():
    """The noisy Clifford circuits: id, family, what generates it, why it is here, its tier."""
    out = []

    def add(id_, code, d, r, p, reason, tier="full"):
        out.append({"id": id_, "family": "qec", "code": code, "distance": d, "rounds": r, "p": p, "reason": reason, "tier": tier, "expect": "qec"})

    for d in (3, 5, 7, 9, 11, 15):
        add("surface_z_d%d" % d, "surface_code:rotated_memory_z", d, d, 0.001, "The rotated surface-code memory experiment in Z: the circuit decoders are measured on. The model must match mechanism by mechanism, and every detector's rate must follow it.", tier="ci" if d == 3 else "full")
    for d in (3, 5, 7):
        add("surface_x_d%d" % d, "surface_code:rotated_memory_x", d, d, 0.001, "The same in the X basis: resets and measurements in X.")
    for d in (3, 5):
        add("unrotated_z_d%d" % d, "surface_code:unrotated_memory_z", d, d, 0.002, "The unrotated layout: another gate order, and detectors with other coordinates.")
    for d in (3, 9, 25):
        add("repetition_d%d" % d, "repetition_code:memory", d, d, 0.005, "A line of data qubits: the simplest model, with long detector chains at d = 25.")
    add("surface_z_d5_noiseless", "surface_code:rotated_memory_z", 5, 5, 0.0, "No noise at all: every detector and the observable are silent on every shot, and the model is empty.", tier="ci")
    add("color_d3", "color_code:memory_xyz", 3, 3, 0.002, "A colour code, whose syndrome cycle turns the axes with C_XYZ gates: the model must match mechanism by mechanism, and every detector's rate must follow it.")
    return out


HAND_WRITTEN_QEC = {
    "qec_inverted_and_tagged": ("A measurement target inverted, which flips its value and no detection event, and a tag on noise.", """\
R 0 1
X_ERROR[leak](0.01) 0
CX 0 1
M !0 1
DETECTOR rec[-1] rec[-2]
OBSERVABLE_INCLUDE(0) rec[-1]
"""),
    "qec_measure_flip_noise": ("Measurement flip noise in both bases, and a reset in X.", """\
R 0
RX 1
M(0.02) 0
MX(0.03) 1
DETECTOR rec[-2]
DETECTOR rec[-1]
"""),
    "qec_depolarize_merge": ("Two channels whose mechanisms have the same symptoms merge by the XOR rule, p1 + p2 - 2 p1 p2.", """\
R 0 1
DEPOLARIZE1(0.01) 0
X_ERROR(0.02) 0
DEPOLARIZE2(0.03) 0 1
M 0 1
DETECTOR rec[-2]
DETECTOR rec[-1]
"""),
    "qec_pauli_observable": ("An observable on Pauli targets rather than measurement records.", """\
R 0 1 2
H 0
CX 0 1 1 2
Z_ERROR(0.01) 1
OBSERVABLE_INCLUDE(0) X0 X1 X2
"""),
    "qec_correlated": ("A correlated error, with a Y among its targets: one mechanism with every symptom, and the targets flipping together.", """\
R 0 1 2
E(0.01) X0 Y1
E(0.02) X2
M 0 1 2
DETECTOR rec[-3]
DETECTOR rec[-2]
DETECTOR rec[-1]
"""),
    "qec_pauli_channel": ("A general Pauli channel: exclusive outcomes, written as independent flips at the rates that give the same law. The model must match mechanism by mechanism, and every detector's rate must follow it.", """\
R 0
RX 1
PAULI_CHANNEL_1(0.01, 0.02, 0.03) 0 1
M 0
MX 1
DETECTOR rec[-2]
DETECTOR rec[-1]
"""),
    "qec_refused_channel_with_no_independent_form": ("A Pauli channel whose exclusive outcomes no independent flips can give (Z's rate would be below zero): it must be refused with its line, never approximated in silence.", """\
R 0
PAULI_CHANNEL_1(0.1, 0.1, 0) 0
M 0
DETECTOR rec[-1]
"""),
    "qec_refused_random_detector": ("A detector that is random without noise is no detector: it must be refused.", """\
RX 0
M 0
DETECTOR rec[-1]
"""),
}


def manifest_entry(e, text):
    return {
        "id": e["id"],
        "family": e["family"],
        "params": e.get("params", {"distance": e.get("distance"), "rounds": e.get("rounds"), "p": e.get("p"), "code": e.get("code")}),
        "reason": e["reason"],
        "expect": e["expect"],
        "tier": e["tier"],
        "qubits": e.get("n"),
        "gates": sum(1 for s in e.get("steps", []) if "g" in s) if "steps" in e else None,
        "sha256": hashlib.sha256(text.encode()).hexdigest(),
    }


def build(tier=None):
    """Every circuit entry in the tier (or all), with its text."""
    es = entries()
    if tier == "ci":
        es = [e for e in es if e["tier"] == "ci"]
    for e in es:
        e["qasm"] = qasm_of(e)
    return es


COMMIT_LIMIT = 16 * 1024  # a program up to this many bytes is committed as a file


def main():
    root = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "corpus")
    os.makedirs(os.path.join(root, "qasm"), exist_ok=True)
    rows = []
    for e in build():
        rows.append(manifest_entry(e, e["qasm"]))
        if len(e["qasm"]) <= COMMIT_LIMIT:
            with open(os.path.join(root, "qasm", e["id"] + ".qasm"), "w") as f:
                f.write(e["qasm"])
    with open(os.path.join(root, "qec-hand-written.json"), "w") as f:
        json.dump({k: {"reason": v[0], "text": v[1]} for k, v in HAND_WRITTEN_QEC.items()}, f, indent=1)
        f.write("\n")
    # The error-correction circuits are the pinned reference generator's, committed as text
    # (corpus/qec/); a run checks that the generator still writes them.
    for e in qec_entries():
        text = open(os.path.join(root, "qec", e["id"] + ".circuit")).read()
        rows.append({"id": e["id"], "family": "qec", "params": {"code": e["code"], "distance": e["distance"], "rounds": e["rounds"], "p": e["p"]}, "reason": e["reason"], "expect": "qec", "tier": e["tier"], "sha256": hashlib.sha256(text.encode()).hexdigest()})
    with open(os.path.join(root, "expected-refusals.json"), "w") as f:
        json.dump({"about": "Refusals by the engine under test that a run accepts, each with the words it must give; tools/corpus.py writes this file from the reader's rule and its list, and a run checks it against them. Any other refusal, crash or time-out of the engine under test fails the run. must: refuse (accepting fails the run) or may (a refusal with these words is recorded as one).", "refusals": expected_refusals()}, f, indent=1, ensure_ascii=False)
        f.write("\n")
    with open(os.path.join(root, "manifest.json"), "w") as f:
        json.dump({"about": "Every circuit of the corpus: parameters, reason, what it is held to, and the SHA-256 of its OpenQASM text. tools/corpus.py writes this file; the run fails if the generator no longer gives these bytes.", "circuits": rows}, f, indent=1, ensure_ascii=False)
        f.write("\n")
    print("%d circuits" % len(rows))


if __name__ == "__main__":
    main()