# SPDX-License-Identifier: MIT

# quantum — simulate a quantum circuit, and reason about one.
#
# Two halves that share nothing but the subject.
#
# SIMULATION. A state of n qubits is a `list<complex>` of 2**n amplitudes, and
# qubit q is bit q of an amplitude's index, so qubit 0 is the lowest bit. Every
# gate is a value script that takes a state and returns the next one, as a
# Python function would, and nothing is written to the store:
#
#     import quantum
#     psi = quantum.zero(2)
#     psi = quantum.h(psi, 0)
#     psi = quantum.cnot(psi, 0, 1)
#     print(quantum.probabilities(psi))
#
# A Clifford circuit can instead be run on a stabilizer tableau, whose size
# grows with n² rather than 2**n. A run has no randomness, so a measurement
# takes its uniform value as an argument, and two generators supply a sequence
# of them from a seed you choose: `lcg_next` from a 31-bit seed, for runs you
# replay, and `chacha_uniform` from a 256-bit key, the one whose bits may become
# part of a secret. The same seed gives the same answers every time.
#
# REASONING. You supply a circuit as facts: its gates (`quantum.Gate`) and the
# pairs of qubits your chip can act on together (`quantum.Coupler`). The
# standing rules below conclude which gates must stay in order, which pairs
# commute although they share a qubit, and which two-qubit gates name a pair
# the chip cannot couple. Each conclusion answers `why`.
#
# DEVICES. A device is facts too: its physical qubits (`quantum.PhysicalQubit`)
# and the operations it runs natively on each qubit and each ordered pair
# (`quantum.Native`), with what each was measured at and when. dkepy's
# `dke.quantum.device` writes them from a calibration snapshot. A circuit names
# the device it was compiled for (`quantum.Circuit`) and its gates name the
# circuit. The rules then conclude which gates the device cannot run as
# written, and why; `stale_gates` asks which ran on a calibration older than
# you allow, at a moment you name.
#
# REVERSIBLE CIRCUITS. A circuit of `x`, `cnot`, `ccx`, `swap` and `cswap`, with
# the signs `z`, `cz` and `ccz`, takes every basis state to one basis state under
# a sign of +1 or -1, so it is checked exactly, on every input, with integers
# alone. `rev_same` asks whether two such circuits are equal under one sign
# shared by every input, and `rev_meets` whether a circuit computes a function
# you specify: inputs unchanged, outputs as specified, every other qubit zero
# again, the sign as specified.
#
# RESULTS. A result computed outside the store (a compilation, a comparison,
# the bits a run measured) is recorded as a `quantum.Result` naming the circuit
# it is about, the circuit it came from and the device it was computed for,
# each with the digest its producer read. While each of those is still what the
# store holds, `holds` concludes the result's outcome, and `why` on it is the
# answer to what produced it. When one has moved, a field names which.

import cmath
import math
import struct
import vector

# ── State vectors ────────────────────────────────────────────────────────────
#
# Each gate is its matrix, applied by one call to `vector.apply`: a 2×2 matrix
# along one qubit, a 4×4 along two. The first qubit named is the high bit of a
# row or column, so `cnot`'s matrix reads |control, target⟩ in the order
# 00, 01, 10, 11. One call is one operation, however many amplitudes it
# touches.

def zero(n: int) -> list<complex>:
    # |0…0⟩ on n qubits.
    psi = [0j] * (2 ** n)
    psi[0] = 1 + 0j
    return psi


def phase(theta: real) -> complex:
    # e**(i·theta), which is cos(theta) + i·sin(theta).
    return cmath.exp(1j * theta)


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


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


def y(psi: list<complex>, q: int) -> list<complex>:
    return vector.apply(psi, [[0, -1j], [1j, 0]], [q])


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


def s(psi: list<complex>, q: int) -> list<complex>:
    return vector.apply(psi, [[1, 0], [0, 1j]], [q])


def t(psi: list<complex>, q: int) -> list<complex>:
    w = (1 + 1j) * math.sqrt(0.5)
    return vector.apply(psi, [[1, 0], [0, w]], [q])


def rx(psi: list<complex>, q: int, theta: real) -> list<complex>:
    w = phase(theta / 2)
    c = w.real
    m = -1j * w.imag
    return vector.apply(psi, [[c, m], [m, c]], [q])


def ry(psi: list<complex>, q: int, theta: real) -> list<complex>:
    w = phase(theta / 2)
    c = w.real
    sn = w.imag
    return vector.apply(psi, [[c, -sn], [sn, c]], [q])


def rz(psi: list<complex>, q: int, theta: real) -> list<complex>:
    w = phase(theta / 2)
    return vector.apply(psi, [[w.conjugate(), 0], [0, w]], [q])


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])


def cz(psi: list<complex>, a: int, b: int) -> list<complex>:
    return vector.apply(psi, [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]],
                        [a, b])


# The gates a device speaks besides those above: `sx`, the square root of `x`,
# and `ecr`, the two-qubit gate of IBM's Eagle devices; and `swap`, which a
# compiler inserts to route. Their matrices are Qiskit's, with the first qubit
# named as the high bit, as everywhere here.

def sx(psi: list<complex>, q: int) -> list<complex>:
    return vector.apply(psi, [[0.5 + 0.5j, 0.5 - 0.5j], [0.5 - 0.5j, 0.5 + 0.5j]], [q])


def ecr(psi: list<complex>, a: int, b: int) -> list<complex>:
    s = math.sqrt(0.5)
    w = 1j * s
    return vector.apply(psi, [[0, 0, s, w], [0, 0, w, s], [s, -w, 0, 0], [-w, s, 0, 0]],
                        [a, b])


def swap(psi: list<complex>, a: int, b: int) -> list<complex>:
    return vector.apply(psi, [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]],
                        [a, b])


# The three-qubit gates of a reversible circuit, with their matrices as above
# (the first qubit named is the high bit of a row or column): `ccx`, the
# Toffoli gate, flips `tg` when both `a` and `b` are 1; `ccz` negates the
# amplitudes whose three qubits are all 1; `cswap`, the Fredkin gate, exchanges
# `a` and `b` when `c` is 1.

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


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


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


def probabilities(psi: list<complex>) -> list<real>:
    return [a.real * a.real + a.imag * a.imag for a in psi]


def weights(psi: list<complex>, q: int) -> tuple<real, real>:
    # The probability of qubit q being 0 and of it being 1: the squared parts
    # of the amplitudes whose qubit q is each, summed exactly (fsum) in one
    # call.
    p0, p1 = vector.sumsq(psi, [q])
    return (p0, p1)


def collapsed(psi: list<complex>, q: int, bit: int, norm: real) -> list<complex>:
    # The diagonal projector on qubit q with `norm` at the outcome `bit`,
    # applied in one call: the amplitudes whose qubit q is `bit`, each times
    # `norm`, and a zero where it is not.
    z = 0.0
    if bit == 0:
        return vector.apply(psi, [[norm, z], [z, z]], [q])
    return vector.apply(psi, [[z, z], [z, norm]], [q])


def measure(psi: list<complex>, q: int, u: real) -> tuple<int, list<complex>>:
    # Measure qubit q. The outcome is 0 when u falls below the probability of
    # 0, and 1 otherwise, so a u drawn uniformly from [0, 1) samples it. The
    # state returned is the one that outcome leaves, renormalised.
    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))


# ── Keeping a state ──────────────────────────────────────────────────────────
#
# A state is kept in the store as one value: its amplitudes packed, in order,
# into one blob, sixteen bytes each (`struct`'s `D`, a complex of two doubles).
# `save` writes it, replacing any state kept under that name, `load` reads it
# back bit for bit, and `drop` removes it, each one call however many
# amplitudes the state has.

class State:
    # A state vector kept in the store, by a name you choose.
    amplitudes: blob

def save(name: string, psi: list<complex>):
    s = State(name)
    update(s.amplitudes, struct.pack(f"<{len(psi)}D", *psi), "quantum")


def load(name: string) -> list<complex>:
    # The state saved under `name`, or an empty list when none is.
    kept = State(name).amplitudes
    match kept:
        case active_claim:
            match kept.value:
                case blob as b:
                    return list(struct.unpack(f"<{len(b) // 16}D", b))
                case default:
                    return []
        case empty:
            return []


def drop(name: string):
    forget(State(name).amplitudes)


def packed(name: string) -> blob:
    # The bytes kept under `name`, or an empty blob when none are.
    kept = State(name).amplitudes
    match kept:
        case active_claim:
            match kept.value:
                case blob as b:
                    return b
                case default:
                    return blob("")
        case empty:
            return blob("")


def keep(name: string, b: blob):
    # Keeps bytes `packed` answered, as `save` keeps a state.
    update(State(name).amplitudes, b, "quantum")


# ── Stabilizer tableau ───────────────────────────────────────────────────────
#
# Aaronson and Gottesman's tableau for n qubits: 2n + 1 rows of 2n + 1 bits.
# Rows 0..n-1 are destabilizers, rows n..2n-1 stabilizers, and row 2n is
# scratch. In each row, column j is the X bit of qubit j, column n + j its Z
# bit, and column 2n the sign. Bits are ints, and combine by `^` and `&`.

def tableau(n: int) -> list<list<int>>:
    # |0…0⟩ on n qubits.
    rows = [[0] * (2 * n + 1) for i in range(2 * n + 1)]
    for i in range(2 * n):
        rows[i][i] = 1
    return rows


def pauli_phase(x1: int, z1: int, x2: int, z2: int) -> int:
    # The power of i that multiplying Pauli (x1, z1) into (x2, z2) contributes.
    if x1 == 0 and z1 == 0:
        return 0
    if x1 == 1 and z1 == 1:
        return z2 - x2
    if x1 == 1:
        return z2 * (2 * x2 - 1)
    return x2 * (1 - 2 * z2)


def row_product(hrow: list<int>, irow: list<int>, n: int) -> list<int>:
    # Row h multiplied by row i, with the sign worked out mod 4.
    total = 2 * hrow[2 * n] + 2 * irow[2 * n]
    out = list(hrow)
    for j in range(n):
        total = total + pauli_phase(irow[j], irow[n + j], hrow[j], hrow[n + j])
        out[j] = hrow[j] ^ irow[j]
        out[n + j] = hrow[n + j] ^ irow[n + j]
    out[2 * n] = 0 if total % 4 == 0 else 1
    return out


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


def tab_s(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][n + a] = za ^ xa
    return out


def tab_cnot(tab: list<list<int>>, a: int, b: 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]
        xb = tab[i][b]
        zb = tab[i][n + b]
        out[i][2 * n] = tab[i][2 * n] ^ (xa & zb & (xb ^ za ^ 1))
        out[i][b] = xb ^ xa
        out[i][n + a] = za ^ zb
    return out


def tab_measure(tab: list<list<int>>, a: int, u: real) -> tuple<int, list<list<int>>>:
    # Measure qubit a. When the outcome is random, it is 1 exactly when
    # u >= 0.5; when it is determined, u is not read.
    n = (len(tab) - 1) // 2
    out = [list(row) for row in tab]
    p = -1
    for i in range(n, 2 * n):
        if tab[i][a] == 1:
            p = i
            break
    if p >= 0:
        for i in range(2 * n):
            if i != p and out[i][a] == 1:
                out[i] = row_product(out[i], out[p], n)
        out[p - n] = list(out[p])
        row = [0] * (2 * n + 1)
        row[n + a] = 1
        bit = 1 if u >= 0.5 else 0
        row[2 * n] = bit
        out[p] = row
        return (bit, out)
    scratch = [0] * (2 * n + 1)
    for i in range(n):
        if out[i][a] == 1:
            scratch = row_product(scratch, out[i + n], n)
    return (scratch[2 * n], out)


# ── Seeded uniform values ────────────────────────────────────────────────────
#
# Park and Miller's minimal standard generator (multiplier 48271, modulus
# 2**31 - 1). A seed is an int from 1 to 2**31 - 2; the product never leaves
# int64, since 48271 * (2**31 - 1) is below 2**47.
#
# Entropy enters as a seed, here or as the key of the generator below, and
# nowhere else: a run draws none of its own. Draw the seed outside the run, for example `secrets.randbelow(2**31 - 2)
# + 1` in Python, and pass it in as an argument. The run then stays replayable,
# and the seed is an argument of the call, so you can record it.

def lcg_next(state: int) -> int:
    return (48271 * state) % 2147483647


def lcg_uniform(state: int) -> real:
    # The state as a value in (0, 1).
    return state / 2147483647


# ── Uniform values from a 256-bit key ────────────────────────────────────────
#
# The generator above has 2**31 - 2 seeds, so anyone can try every one of
# them, and a value it made can stand for no secret. It is for runs you want
# to replay. ChaCha20 (RFC 8439) takes a key of 256 bits instead, as 8 words
# from 0 to 0xFFFFFFFF, and without the key its values cannot be told from
# random ones. It is the generator to use when a measured bit is to become
# part of a secret.
#
# The key is the seed, drawn outside the run as the seed above is, for example
# `list(struct.unpack("<8I", secrets.token_bytes(32)))` in Python. It is an
# argument like any other: it travels to the service in the program you send,
# and the service sees it. A secret that is to stay on your own machine is
# drawn on your own machine.
#
# Position n of the stream is value n % 8 of block n // 8, with a nonce of
# zero, made from that block's words 2·(n % 8) and 2·(n % 8) + 1 as an exact
# 53-bit value in [0, 1). Positions run from 0 to 2**35 - 1, and the same key
# and position give the same value every time.

def chacha_quarter(w: list<int>, a: int, b: int, c: int, d: int) -> list<int>:
    # ChaCha20's quarter round on words a, b, c and d of w, in place.
    w[a] = (w[a] + w[b]) & 0xFFFFFFFF
    x = w[d] ^ w[a]
    w[d] = ((x << 16) | (x >> 16)) & 0xFFFFFFFF
    w[c] = (w[c] + w[d]) & 0xFFFFFFFF
    x = w[b] ^ w[c]
    w[b] = ((x << 12) | (x >> 20)) & 0xFFFFFFFF
    w[a] = (w[a] + w[b]) & 0xFFFFFFFF
    x = w[d] ^ w[a]
    w[d] = ((x << 8) | (x >> 24)) & 0xFFFFFFFF
    w[c] = (w[c] + w[d]) & 0xFFFFFFFF
    x = w[b] ^ w[c]
    w[b] = ((x << 7) | (x >> 25)) & 0xFFFFFFFF
    return w


def chacha_block(key: list<int>, counter: int, nonce: list<int>) -> list<int>:
    # One ChaCha20 block, as its 16 words, for a key of 8 words, a block
    # counter and a nonce of 3 words.
    if len(key) != 8:
        raise ValueError("a ChaCha20 key is 8 words, each from 0 to 0xFFFFFFFF")
    for k in key:
        if k < 0 or k > 0xFFFFFFFF:
            raise ValueError("a ChaCha20 key is 8 words, each from 0 to 0xFFFFFFFF")
    if counter < 0 or counter > 0xFFFFFFFF:
        raise ValueError("a ChaCha20 block counter is from 0 to 0xFFFFFFFF")
    if len(nonce) != 3:
        raise ValueError("a ChaCha20 nonce is 3 words, each from 0 to 0xFFFFFFFF")
    for k in nonce:
        if k < 0 or k > 0xFFFFFFFF:
            raise ValueError("a ChaCha20 nonce is 3 words, each from 0 to 0xFFFFFFFF")
    start = [0x61707865, 0x3320646E, 0x79622D32, 0x6B206574] + key + [counter] + nonce
    w = start.copy()
    for _ in range(10):
        w = chacha_quarter(w, 0, 4, 8, 12)
        w = chacha_quarter(w, 1, 5, 9, 13)
        w = chacha_quarter(w, 2, 6, 10, 14)
        w = chacha_quarter(w, 3, 7, 11, 15)
        w = chacha_quarter(w, 0, 5, 10, 15)
        w = chacha_quarter(w, 1, 6, 11, 12)
        w = chacha_quarter(w, 2, 7, 8, 13)
        w = chacha_quarter(w, 3, 4, 9, 14)
    return [(w[i] + start[i]) & 0xFFFFFFFF for i in range(16)]


def chacha_unit(block: list<int>, i: int) -> real:
    # Value i of a block, from its words 2·i and 2·i + 1, as an exact 53-bit
    # value in [0, 1): 27 bits of the first, 26 of the second.
    return ((block[2 * i] >> 5) * 67108864 + (block[2 * i + 1] >> 6)) / 9007199254740992


def chacha_uniform(key: list<int>, n: int) -> real:
    # The value at position n of the key's stream.
    if n < 0 or n > 34359738367:
        raise ValueError("a ChaCha20 stream position is from 0 to 2**35 - 1")
    return chacha_unit(chacha_block(key, n // 8, [0, 0, 0]), n % 8)


# ── A circuit in steps ───────────────────────────────────────────────────────
#
# `advance` runs a list of steps on a kept state and keeps the state it ends
# with, so a circuit longer than one run is a series of calls that pass nothing
# but the steps. A step is a gate's name, the qubits it acts on and an angle,
# which only the rotations read, or `("measure", [q], 0.0)`. Measurements take
# their uniform values from `lcg_uniform(u)`, then `u = lcg_next(u)`, and
# `advance` answers the generator state it ended with and the bits it measured.
# `advance_chacha` takes them from a 256-bit key instead, at a stream position
# it answers in the generator state's place.
#
# Every gate is applied by one line, with its matrix from `unitary`, so what a
# run returns does not grow with the number of steps it takes. `unitary` writes
# each matrix as the gate functions above apply it, each real number as a
# complex with a zero imaginary part, which multiplies exactly as the real
# number does.

def unitary(name: string, theta: real) -> list<list<complex>>:
    # The matrix of gate `name`. `theta` is read by `rx`, `ry` and `rz` only.
    if name == "h":
        s = math.sqrt(0.5)
        return [[complex(s, 0.0), complex(s, 0.0)], [complex(s, 0.0), complex(-s, 0.0)]]
    if name == "x":
        return [[0j, 1 + 0j], [1 + 0j, 0j]]
    if name == "y":
        return [[0j, -1j], [1j, 0j]]
    if name == "z":
        return [[1 + 0j, 0j], [0j, -1 + 0j]]
    if name == "s":
        return [[1 + 0j, 0j], [0j, 1j]]
    if name == "t":
        return [[1 + 0j, 0j], [0j, (1 + 1j) * math.sqrt(0.5)]]
    if name == "rx":
        w = phase(theta / 2)
        m = -1j * w.imag
        return [[complex(w.real, 0.0), m], [m, complex(w.real, 0.0)]]
    if name == "ry":
        w = phase(theta / 2)
        return [[complex(w.real, 0.0), complex(-w.imag, 0.0)],
                [complex(w.imag, 0.0), complex(w.real, 0.0)]]
    if name == "rz":
        w = phase(theta / 2)
        return [[w.conjugate(), 0j], [0j, w]]
    if name == "cnot":
        return [[1 + 0j, 0j, 0j, 0j], [0j, 1 + 0j, 0j, 0j],
                [0j, 0j, 0j, 1 + 0j], [0j, 0j, 1 + 0j, 0j]]
    if name == "cz":
        return [[1 + 0j, 0j, 0j, 0j], [0j, 1 + 0j, 0j, 0j],
                [0j, 0j, 1 + 0j, 0j], [0j, 0j, 0j, -1 + 0j]]
    if name == "sx":
        return [[0.5 + 0.5j, 0.5 - 0.5j], [0.5 - 0.5j, 0.5 + 0.5j]]
    if name == "ecr":
        s = math.sqrt(0.5)
        r = complex(s, 0.0)
        w = 1j * s
        return [[0j, 0j, r, w], [0j, 0j, w, r], [r, -w, 0j, 0j], [-w, r, 0j, 0j]]
    if name == "swap":
        return [[1 + 0j, 0j, 0j, 0j], [0j, 0j, 1 + 0j, 0j],
                [0j, 1 + 0j, 0j, 0j], [0j, 0j, 0j, 1 + 0j]]
    if name == "ccx":
        return [[1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 0j, 1 + 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j]]
    if name == "ccz":
        return [[1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 0j, -1 + 0j]]
    if name == "cswap":
        return [[1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 1 + 0j, 0j, 0j],
                [0j, 0j, 0j, 0j, 0j, 0j, 0j, 1 + 0j]]
    raise ValueError(f"no gate named {name}")


def advance(name: string, steps: list<tuple<string, list<int>, real>>, u: int) -> tuple<int, list<int>>:
    psi = load(name)
    bits = []
    for step in steps:
        if step[0] == "measure":
            bit, psi = measure(psi, step[1][0], lcg_uniform(u))
            bits.append(bit)
            u = lcg_next(u)
        else:
            psi = vector.apply(psi, unitary(step[0], step[2]), step[1])
    save(name, psi)
    return (u, bits)


def advance_chacha(name: string, steps: list<tuple<string, list<int>, real>>, key: list<int>, n: int) -> tuple<int, list<int>>:
    # `advance` with its uniform values from `chacha_uniform(key, n)`, then
    # `n = n + 1`. It answers the next position and the bits it measured.
    if n < 0 or n > 34359738367:
        raise ValueError("a ChaCha20 stream position is from 0 to 2**35 - 1")
    psi = load(name)
    bits = []
    held = -1
    words = []
    for step in steps:
        if step[0] == "measure":
            if n > 34359738367:
                raise ValueError("a ChaCha20 stream position is from 0 to 2**35 - 1")
            if n // 8 != held:
                held = n // 8
                words = chacha_block(key, held, [0, 0, 0])
            bit, psi = measure(psi, step[1][0], chacha_unit(words, n % 8))
            bits.append(bit)
            n = n + 1
        else:
            psi = vector.apply(psi, unitary(step[0], step[2]), step[1])
    save(name, psi)
    return (n, bits)


# The two halves of `measure` on a kept state, each reached by name: `weigh`
# answers a qubit's weights, and `collapse` collapses the state to one outcome.
# A caller that keeps one state in parts, across several stores, weighs every
# part before it collapses any.

def weigh(name: string, q: int) -> tuple<real, real>:
    return weights(load(name), q)


def collapse(name: string, q: int, bit: int, norm: real):
    save(name, collapsed(load(name), q, bit, norm))


# ── Equivalence ──────────────────────────────────────────────────────────────
#
# Whether two circuits on n qubits are the same unitary up to a global phase,
# which is what a compiler must preserve. Two ways, one exact.
#
# A Clifford circuit — every gate one `tab_gate` names — has a tableau, and two
# tableaux built from |0…0⟩ are equal, sign bits included, exactly when their
# circuits are equal up to a global phase: the rows are what the circuit makes
# of each X and Z, and those fix it.
#
# Any circuit can be compared by `overlap`: run one circuit, then the other
# backwards, on |j⟩, and read amplitude j. Equal circuits answer the same
# e**(i·φ) for every j, so the caller asks every j and compares the answers
# with one another. One phase for all of them is the point: a phase that may
# differ from column to column accepts an `rz` dropped or of the wrong sign.

def adjoint(m: list<list<complex>>) -> list<list<complex>>:
    # The conjugate transpose: the inverse of a unitary matrix, exactly.
    return [[m[c][r].conjugate() for c in range(len(m))] for r in range(len(m))]


def overlap(n: int, steps: list<tuple<string, list<int>, real>>, back: list<tuple<string, list<int>, real>>, j: int) -> complex:
    # Amplitude j after `steps`, then `back` undone, run on |j⟩. `back` is
    # written in the order it runs, and is undone last step first, each step by
    # its matrix's adjoint.
    psi = [0j] * (2 ** n)
    psi[j] = 1 + 0j
    for step in steps:
        psi = vector.apply(psi, unitary(step[0], step[2]), step[1])
    for i in range(len(back)):
        step = back[len(back) - 1 - i]
        psi = vector.apply(psi, adjoint(unitary(step[0], step[2])), step[1])
    return psi[j]


def tab_gate(tab: list<list<int>>, name: string, q: list<int>) -> list<list<int>>:
    # A Clifford gate on the tableau, made of `tab_h`, `tab_s` and `tab_cnot`,
    # each innermost call first. A tableau has no global phase, so each is its
    # gate up to one: `z` is `s s`, `x` is `h s s h`, `y` is `z` then `x`, `sx`
    # is `h s h`, `cz` is `cnot` between two `h` on the second qubit, and `ecr`
    # is the definition Qiskit's exporter writes: `s` on the first qubit, `sx`
    # on the second, `cnot`, then `x` on the first. An `rz` by a quarter turn
    # is `s`, but whether an angle is exactly one is the caller's to know.
    a = q[0]
    if name == "h":
        return tab_h(tab, a)
    if name == "s":
        return tab_s(tab, a)
    if name == "z":
        return tab_s(tab_s(tab, a), a)
    if name == "x":
        return tab_h(tab_s(tab_s(tab_h(tab, a), a), a), a)
    if name == "y":
        return tab_h(tab_s(tab_s(tab_h(tab_s(tab_s(tab, a), a), a), a), a), a)
    if name == "sx":
        return tab_h(tab_s(tab_h(tab, a), a), a)
    if name == "cnot":
        return tab_cnot(tab, a, q[1])
    if name == "cz":
        return tab_h(tab_cnot(tab_h(tab, q[1]), a, q[1]), q[1])
    if name == "swap":
        return tab_cnot(tab_cnot(tab_cnot(tab, a, q[1]), q[1], a), a, q[1])
    if name == "ecr":
        b = q[1]
        t1 = tab_h(tab_s(tab_h(tab_s(tab, a), b), b), b)
        return tab_h(tab_s(tab_s(tab_h(tab_cnot(t1, a, b), a), a), a), a)
    raise ValueError(f"no Clifford gate named {name}")


def clifford(n: int, steps: list<tuple<string, list<int>, real>>) -> list<list<int>>:
    # The tableau of a Clifford circuit on n qubits, from |0…0⟩.
    tab = tableau(n)
    for step in steps:
        tab = tab_gate(tab, step[0], step[1])
    return tab


# ── Reversible circuits ──────────────────────────────────────────────────────
#
# A circuit of `x`, `cnot`, `ccx` and `swap`, `cswap`, with the signs `z`, `cz`
# and `ccz`, never makes a superposition: every basis state goes to one basis
# state, with a sign of +1 or -1. So it can be checked exactly, on every input,
# with integers alone, however many qubits a state vector could not hold.
#
# Qubit q of an input is bit q of the input's index, as everywhere here. The
# inputs are run 64 at a time: a wire is one int whose bit k is that qubit in
# input 64·p + k, so one bitwise operation applies a gate to 64 inputs. A sign
# is one int the same way, its bit k set where input k's sign is -1.
#
# Two circuits are EQUAL when, on every input, they give the same output bits
# under one sign shared by every input: a circuit that differs by a sign on
# some inputs only is not equal, and a circuit that differs by -1 on all of
# them is, as a global phase is.
#
# A circuit MEETS a specification when, on every input x of its input qubits,
# the input qubits are unchanged, the output qubits hold the value the
# specification gives for x, every other qubit is 0 again, and the sign is the
# one specified. The specification's values are a list, one per input, that the
# caller computes with a function of its own; the circuit starts with every
# qubit but the input qubits at 0.
#
# A circuit may MEASURE. `("mx", [q], 0.0)` measures qubit q in the X basis and
# leaves it at 0. Its outcome is not drawn from anything: it is one bit of the
# case being checked, as a seed is an argument. Under outcome 1 it negates the
# sign on every input where q held 1, and under outcome 0 it changes no sign.
# A step named `x_if`, `cnot_if`, `ccx_if`, `z_if`, `cz_if`, `ccz_if`,
# `swap_if` or `cswap_if` is the gate it names applied only where an outcome
# was 1. It takes the measured qubit FIRST and then the gate's own qubits in
# the gate's order: `("cz_if", [t, a, b], 0.0)` is `cz` on a and b where the
# latest `mx` of t before it gave 1, and nothing where it gave 0. A step whose
# qubit no earlier `mx` measured is refused.
#
# Outcomes are numbered by the `mx` steps in step order: outcome bit j is the
# j-th `mx`. A specification is held over EVERY assignment of the outcomes,
# so a circuit with M measurements is checked on 2**M times as many cases, and
# a case is an input under one assignment: case i is input i mod 2**k under the
# outcomes i div 2**k, with k the number of input qubits, and the sign stays
# exact. A circuit with no `mx` has one assignment, and is checked as it
# always was. Two circuits are not compared when either measures: `rev_same`
# looks at every input of all n qubits, a measured qubit that starts at 1 is
# reset and not undone, and a comparison there means nothing. A circuit that
# measures is held to a specification, where every qubit but the inputs
# starts at 0, and `rev_same` refuses it.
#
# Both answer (holds, cases checked, first failing case, why). The reason is
# "" when it holds, and otherwise "output" or "sign" for two circuits, and
# "input" (an input qubit changed), "ancilla" (a qubit that should be 0 is
# not), "output" or "sign" for a specification. A specification check of more
# than 2**62 cases is refused.

def rev_code(n: int, steps: list<tuple<string, list<int>, real>>) -> list<tuple<int, int, int, int>>:
    # The steps of a reversible circuit on n qubits as (operation, a, b, c)
    # for `rev_apply`: operations 0 `x`, 1 `cnot`, 2 `ccx`, 3 `z`, 4 `cz`,
    # 5 `ccz`, 6 `swap`, 7 `cswap`, with the qubits in the order the gate
    # function takes them, and 8 `mx` on a, its number among the `mx` steps in b. A step `<gate>_if` is its gate's
    # operation plus 16·(1 + j), j the number of the latest `mx` of its first
    # qubit, and the qubits that follow are the gate's own. A gate that makes
    # a superposition, a qubit outside the circuit, a qubit named twice in
    # one step (the measured one included) and a measured qubit no earlier
    # `mx` measured are refused.
    base = ["x", "cnot", "ccx", "z", "cz", "ccz", "swap", "cswap"]
    arity = [1, 2, 3, 1, 2, 3, 2, 3]
    last = [-1] * n
    seen = 0
    out = []
    for step in steps:
        name = step[0]
        q = step[1]
        op = -1
        k = 0
        cond = False
        for g in range(8):
            if name == base[g]:
                op = g
                k = arity[g]
            elif name == base[g] + "_if":
                op = g
                k = arity[g] + 1
                cond = True
        if name == "mx":
            op = 8
            k = 1
        if op < 0:
            raise ValueError(f"no reversible gate named {name}")
        if len(q) != k:
            raise ValueError(f"{name} takes {k} qubit(s)")
        for i in range(k):
            if q[i] < 0 or q[i] >= n:
                raise ValueError(f"{name} names qubit {q[i]}, outside a circuit of {n}")
            for j in range(i):
                if q[i] == q[j]:
                    raise ValueError(f"{name} names qubit {q[i]} twice")
        s = 0
        if cond:
            if last[q[0]] < 0:
                raise ValueError(f"{name} names qubit {q[0]}, which no earlier mx measured")
            op = op + 16 * (1 + last[q[0]])
            s = 1
        if name == "mx":
            out.append((8, q[0], seen, 0))
            last[q[0]] = seen
            seen = seen + 1
        elif k - s == 1:
            out.append((op, q[s], 0, 0))
        elif k - s == 2:
            out.append((op, q[s], q[s + 1], 0))
        else:
            out.append((op, q[s], q[s + 1], q[s + 2]))
    return out


def rev_measures(code: list<tuple<int, int, int, int>>) -> int:
    # The number of `mx` steps in a circuit `rev_code` read.
    m = 0
    for op, a, b, c in code:
        if op == 8:
            m = m + 1
    return m


def rev_apply(w: list<int>, code: list<tuple<int, int, int, int>>) -> int:
    # The circuit `rev_code` read, run on the wires `w` (one int per qubit,
    # 64 inputs each) and left in them; the answer is the sign of the 64
    # inputs, bit k set where input k's is -1. Every outcome is 0.
    return rev_apply_under(w, code, 0)


def rev_apply_under(w: list<int>, code: list<tuple<int, int, int, int>>, outcomes: int) -> int:
    # `rev_apply` with the outcomes of the circuit's `mx` steps given: outcome
    # j is bit j of `outcomes`, the same for all 64 inputs.
    sg = 0
    for step in code:
        op = step[0]
        a = step[1]
        b = step[2]
        c = step[3]
        run = True
        if op >= 16:
            run = (outcomes >> (op // 16 - 1)) & 1 == 1
            op = op % 16
        if not run:
            op = -1
        if op == 2:
            w[c] = w[c] ^ (w[a] & w[b])
        elif op == 1:
            w[b] = w[b] ^ w[a]
        elif op == 0:
            w[a] = ~w[a]
        elif op == 4:
            sg = sg ^ (w[a] & w[b])
        elif op == 3:
            sg = sg ^ w[a]
        elif op == 5:
            sg = sg ^ (w[a] & w[b] & w[c])
        elif op == 6:
            t = w[a]
            w[a] = w[b]
            w[b] = t
        elif op == 7:
            d = w[a] & (w[b] ^ w[c])
            w[b] = w[b] ^ d
            w[c] = w[c] ^ d
        elif op == 8:
            if (outcomes >> b) & 1 == 1:
                sg = sg ^ w[a]
            w[a] = 0
    return sg


def lowest_lane(m: int) -> int:
    # The lowest bit set in a nonzero 64-bit int, 0 to 63.
    lane = 0
    x = m
    if x & 4294967295 == 0:
        lane = lane + 32
        x = x >> 32
    if x & 65535 == 0:
        lane = lane + 16
        x = x >> 16
    if x & 255 == 0:
        lane = lane + 8
        x = x >> 8
    if x & 15 == 0:
        lane = lane + 4
        x = x >> 4
    if x & 3 == 0:
        lane = lane + 2
        x = x >> 2
    if x & 1 == 0:
        lane = lane + 1
    return lane


def rev_same(n: int, a: list<tuple<string, list<int>, real>>, b: list<tuple<string, list<int>, real>>) -> tuple<bool, int, int, string>:
    # Whether two reversible circuits on n qubits are equal, as above, on all
    # 2**n inputs: (holds, 2**n, the first input they differ at or -1, why).
    # A circuit that measures is refused: it is held to a specification.
    ca = rev_code(n, a)
    cb = rev_code(n, b)
    if rev_measures(ca) + rev_measures(cb) != 0:
        raise ValueError("rev_same takes no circuit that measures: hold it to a specification with rev_meets, where every qubit but the inputs starts at 0")
    pat = [-6148914691236517206, -3689348814741910324, -1085102592571150096,
           -71777214294589696, -281470681808896, -4294967296]
    total = 2 ** n
    valid = -1
    if total < 64:
        valid = (1 << total) - 1
    wa = [0] * n
    wb = [0] * n
    qs = list(range(n))
    shared = 0
    for p in range((total + 63) // 64):
        for q in qs:
            v = pat[q] if q < 6 else -((p >> (q - 6)) & 1)
            wa[q] = v
            wb[q] = v
        sa = rev_apply(wa, ca)
        sb = rev_apply(wb, cb)
        diff = 0
        for q in qs:
            diff = diff | (wa[q] ^ wb[q])
        diff = diff & valid
        sd = (sa ^ sb) & valid
        if p == 0:
            shared = sd & 1
        bad_sign = (sd ^ -shared) & valid
        bad = diff | bad_sign
        if bad != 0:
            lane = lowest_lane(bad)
            why = "output" if (diff >> lane) & 1 == 1 else "sign"
            return (False, total, 64 * p + lane, why)
    return (True, total, -1, "")


def rev_meets(n: int, steps: list<tuple<string, list<int>, real>>, ins: list<int>, outs: list<int>, want: list<int>, signs: list<int>) -> tuple<bool, int, int, string>:
    # Whether a reversible circuit on n qubits meets a specification, as above,
    # on all 2**len(ins) inputs under every outcome of its measurements:
    # (holds, the cases, the first case it fails at or -1, why). Input x puts bit j of x on qubit ins[j]; `want[x]` is the
    # value the qubits `outs` hold at the end, bit j on qubit outs[j], and
    # `signs[x]` is 1 where the sign is -1 (an empty `signs` is +1 for every
    # input). The three lists of qubits share none, and `want` has one value
    # of the outputs' width for each input.
    code = rev_code(n, steps)
    k = len(ins)
    m = len(outs)
    if k + rev_measures(code) > 62:
        raise ValueError(f"a check of 2**{k + rev_measures(code)} cases does not fit")
    total = 2 ** k
    outcomes = 2 ** rev_measures(code)
    if len(want) != total:
        raise ValueError(f"the specification has {len(want)} values and {k} input qubits take {total}")
    if len(signs) != 0 and len(signs) != total:
        raise ValueError(f"the signs are {len(signs)} and there are {total} inputs")
    role = [0] * n
    for q in ins:
        if q < 0 or q >= n or role[q] != 0:
            raise ValueError(f"qubit {q} is outside the circuit or named twice")
        role[q] = 1
    for q in outs:
        if q < 0 or q >= n or role[q] != 0:
            raise ValueError(f"qubit {q} is outside the circuit or named twice")
        role[q] = 2
    for x in want:
        if x < 0 or x >= 2 ** m:
            raise ValueError(f"a specified value {x} does not fit {m} output qubit(s)")
    pat = [-6148914691236517206, -3689348814741910324, -1085102592571150096,
           -71777214294589696, -281470681808896, -4294967296]
    valid = -1
    if total < 64:
        valid = (1 << total) - 1
    w = [0] * n
    qs = list(range(n))
    ks = list(range(k))
    ms = list(range(m))
    lane_ids = list(range(64))
    lane_bit = [1 << i for i in range(63)] + [-9223372036854775807 - 1]
    wantw = [0] * m
    for oc in range(outcomes):
        for p in range((total + 63) // 64):
            for q in qs:
                w[q] = 0
            for j in ks:
                w[ins[j]] = pat[j] if j < 6 else -((p >> (j - 6)) & 1)
            sg = rev_apply_under(w, code, oc)
            bad_in = 0
            for j in ks:
                was = pat[j] if j < 6 else -((p >> (j - 6)) & 1)
                bad_in = bad_in | (w[ins[j]] ^ was)
            bad_anc = 0
            for q in qs:
                if role[q] == 0:
                    bad_anc = bad_anc | w[q]
            signw = 0
            for j in ms:
                wantw[j] = 0
            for ln in lane_ids:
                x = 64 * p + ln
                if x < total:
                    v = want[x]
                    for j in ms:
                        wantw[j] = wantw[j] | (((v >> j) & 1) * lane_bit[ln])
                    if len(signs) != 0 and signs[x] != 0:
                        signw = signw | lane_bit[ln]
            bad_out = 0
            for j in ms:
                bad_out = bad_out | (w[outs[j]] ^ wantw[j])
            bad_sign = sg ^ signw
            bad = (bad_in | bad_anc | bad_out | bad_sign) & valid
            if bad != 0:
                lane = lowest_lane(bad)
                why = "sign"
                if (bad_out >> lane) & 1 == 1:
                    why = "output"
                if (bad_anc >> lane) & 1 == 1:
                    why = "ancilla"
                if (bad_in >> lane) & 1 == 1:
                    why = "input"
                return (False, total * outcomes, oc * total + 64 * p + lane, why)
    return (True, total * outcomes, -1, "")


# ── A circuit as facts ───────────────────────────────────────────────────────

class Gate:
    # One gate: its operation (an `Op`), the qubit it acts on, the second qubit
    # of a two-qubit gate, and the step it runs at. Qubits are named, and the
    # rules take every name a gate uses as a qubit, recorded or not, a number
    # included. `circuit` names the `Circuit` it belongs to. Gates of different
    # circuits are never ordered against each other, though they name the same
    # qubit, and the gates that name no circuit are one circuit together.
    op: string
    a: string
    b: string
    step: int
    circuit: string

class Qubit:
    # A qubit, by the name gates and couplers use. Recording one is optional,
    # and no rule reads it: the rules compare the names gates use.
    label: string

class Coupler:
    # Two qubits the chip can act on together, in either order.
    a: string
    b: string

class Op:
    # An operation's properties: how many qubits it takes, and whether it is
    # diagonal. Diagonal operations commute with one another. The module
    # records its own operations when it is imported; record more the same way.
    qubits: int
    diagonal: bool

remember(Op.h.qubits, 1, "quantum")
remember(Op.h.diagonal, False, "quantum")
remember(Op.x.qubits, 1, "quantum")
remember(Op.x.diagonal, False, "quantum")
remember(Op.y.qubits, 1, "quantum")
remember(Op.y.diagonal, False, "quantum")
remember(Op.z.qubits, 1, "quantum")
remember(Op.z.diagonal, True, "quantum")
remember(Op.s.qubits, 1, "quantum")
remember(Op.s.diagonal, True, "quantum")
remember(Op.t.qubits, 1, "quantum")
remember(Op.t.diagonal, True, "quantum")
remember(Op.rx.qubits, 1, "quantum")
remember(Op.rx.diagonal, False, "quantum")
remember(Op.ry.qubits, 1, "quantum")
remember(Op.ry.diagonal, False, "quantum")
remember(Op.rz.qubits, 1, "quantum")
remember(Op.rz.diagonal, True, "quantum")
remember(Op.cnot.qubits, 2, "quantum")
remember(Op.cnot.diagonal, False, "quantum")
remember(Op.cz.qubits, 2, "quantum")
remember(Op.cz.diagonal, True, "quantum")
remember(Op.sx.qubits, 1, "quantum")
remember(Op.sx.diagonal, False, "quantum")
remember(Op.id.qubits, 1, "quantum")
remember(Op.id.diagonal, True, "quantum")
remember(Op.ecr.qubits, 2, "quantum")
remember(Op.ecr.diagonal, False, "quantum")
remember(Op.swap.qubits, 2, "quantum")
remember(Op.swap.diagonal, False, "quantum")
remember(Op.rzz.qubits, 2, "quantum")
remember(Op.rzz.diagonal, True, "quantum")

# ── A device as facts ────────────────────────────────────────────────────────

class Device:
    # A device a circuit can be compiled for: how many physical qubits it has,
    # when the snapshot its facts were read from was taken, and that snapshot's
    # digest, which names the calibration exactly.
    qubits: int
    snapshot_at: datetime
    digest: string

class PhysicalQubit:
    # One of a device's qubits. `index` is its number as text, which is how a
    # gate of a circuit compiled for the device names it. `t1` and `t2` are its
    # relaxation and dephasing times in microseconds, and `measured_at` is when
    # the latest of its values was measured.
    device: string
    index: string
    t1: real
    t2: real
    readout_error: real
    measured_at: datetime

class Native:
    # An operation a device runs natively on qubit `a`, or on the ordered pair
    # `a` then `b`. The pair is ordered because a device may run a two-qubit
    # operation one way round only. `error` is the measured probability that it
    # fails, `length` its duration in nanoseconds, and `measured_at` when the
    # error was measured.
    device: string
    op: string
    a: string
    b: string
    error: real
    length: real
    measured_at: datetime

class Circuit:
    # The device a circuit was compiled for, and the digest of the text it was
    # written from. Its gates name it in `circuit`. The digest is the writer's
    # to keep true: the store compares digests and never computes one.
    device: string
    digest: string

class Result:
    # A result computed outside the store, and what produced it. `circuit` is
    # the circuit it is about and `source` the circuit that one was derived
    # from, if any; `device` is the device it was computed for, if any. Each
    # carries the digest its producer read. `producer` is what produced it,
    # with its version, `settings` how it was asked, `seed` the seed it drew
    # from, and `outcome` the result itself, as text. Those four are a recipe
    # and are not judged: no rule can know whether a producer is still current.
    circuit: string
    digest: string
    source: string
    source_digest: string
    device: string
    device_digest: string
    producer: string
    settings: string
    seed: int
    outcome: string

# ── Standing rules ───────────────────────────────────────────────────────────

@rule
def gate_diagonal():
    for g in Gate, o in Op:
        if g.op == o:
            g.diagonal = o.diagonal

@rule
def track_of_circuit():
    # The circuit a gate is ordered within: the one it names.
    for g in Gate:
        if g.circuit != "":
            g.track = g.circuit

@rule
def track_of_none():
    # The gates that name no circuit are ordered among themselves. Only a gate
    # with a step can be ordered at all, and that premise names the gates the
    # absence is asked of.
    for g in Gate:
        if g.step == g.step and absent(g.circuit):
            g.track = ""

# g runs before h on a qubit both act on, in one circuit: one rule for each
# pairing of g's qubits with h's. A qubit is compared by its name, never bound
# as a subject, since a subject's name cannot be a number and a compiled
# circuit names its qubits by number.

@rule
def shares_earlier_aa():
    for g in Gate, h in Gate:
        if g.a == h.a and g.track == h.track and g.step < h.step:
            RunsBefore.g.h = True

@rule
def shares_earlier_ab():
    for g in Gate, h in Gate:
        if g.a == h.b and g.track == h.track and g.step < h.step:
            RunsBefore.g.h = True

@rule
def shares_earlier_ba():
    for g in Gate, h in Gate:
        if g.b == h.a and g.track == h.track and g.step < h.step:
            RunsBefore.g.h = True

@rule
def shares_earlier_bb():
    for g in Gate, h in Gate:
        if g.b == h.b and g.track == h.track and g.step < h.step:
            RunsBefore.g.h = True

@rule
def ordered_first():
    # The order of g and h matters: they share a qubit and g is not diagonal.
    for g in Gate, h in Gate:
        if RunsBefore.g.h == True and g.diagonal == False:
            MustPrecede.g.h = True

@rule
def ordered_second():
    for g in Gate, h in Gate:
        if RunsBefore.g.h == True and h.diagonal == False:
            MustPrecede.g.h = True

@rule
def commutes():
    # g and h share a qubit and both are diagonal, so either may run first.
    for g in Gate, h in Gate:
        if RunsBefore.g.h == True and g.diagonal == True and h.diagonal == True:
            CommutesWith.g.h = True

@rule
def uncoupled():
    # A two-qubit gate whose qubits no coupler joins, in either order. A gate
    # on a device is judged by the device's own pairs instead, in the device
    # rules below, so couplers say nothing about it.
    for g in Gate, o in Op:
        if g.op == o and o.qubits == 2 and absent(g.device) and none(c in Coupler: c.a == g.a and c.b == g.b) and none(c in Coupler: c.a == g.b and c.b == g.a):
            g.uncoupled = True

@constraint
def every_pair_coupled():
    for g in Gate:
        if g.uncoupled == True:
            violated

# ── Device rules ─────────────────────────────────────────────────────────────
#
# A gate of a circuit compiled for a device runs on it as one of the device's
# native operations, or is refused for one of four reasons, each its own field
# so `why` names the one that holds: `off_device`, a qubit the device does not
# have; `not_native`, an operation the device runs nowhere; `not_here`, an
# operation the device runs, on qubits it has, but not on these; and
# `reversed`, a pair the device runs the operation on the other way round only.
# A gate of no circuit, or of a circuit that names no device, is not judged at
# all. A circuit naming a device the store holds no facts for has every gate
# refused, as `off_device` and `not_native`: nothing is known to run there,
# which is also how a misspelt device name shows itself.
#
# What a reason needs to exist — the device runs the operation somewhere, has
# the first qubit, has the second — is concluded once per gate, as
# `op_on_device`, `first_on_device` and `second_on_device`, and each reason's
# rule reads that field.

@rule
def gate_device():
    # A gate runs on the device its circuit was compiled for.
    for g in Gate, c in Circuit:
        if g.circuit == c:
            g.device = c.device

@rule
def gate_width():
    # How many qubits the gate's operation takes. An operation with no `Op`
    # fact has no width, and only the reasons that need none can hold for it.
    for g in Gate, o in Op:
        if g.op == o:
            g.width = o.qubits

@rule
def op_on_device():
    # The device runs the gate's operation on some qubit or pair.
    for g in Gate, n in Native:
        if n.op == g.op and n.device == g.device:
            g.op_on_device = True

@rule
def first_on_device():
    # The device has the gate's first qubit.
    for g in Gate, p in PhysicalQubit:
        if p.index == g.a and p.device == g.device:
            g.first_on_device = True

@rule
def second_on_device():
    # The device has the gate's second qubit.
    for g in Gate, p in PhysicalQubit:
        if p.index == g.b and p.device == g.device:
            g.second_on_device = True

@rule
def runs_as_one():
    # A one-qubit gate runs as the device's native operation on its qubit.
    for g in Gate, n in Native:
        if g.width == 1 and n.a == g.a and n.op == g.op and n.device == g.device:
            RunsAs.g.n = True

@rule
def runs_as_pair():
    # A two-qubit gate runs as the device's native operation on its ordered pair.
    for g in Gate, n in Native:
        if n.a == g.a and n.b == g.b and n.op == g.op and n.device == g.device:
            RunsAs.g.n = True

@rule
def off_device_first():
    # The gate names a qubit its device does not have.
    for g in Gate, d in Device:
        if g.device == d and none(p in PhysicalQubit: p.index == g.a and p.device == g.device):
            g.off_device = True

@rule
def off_device_second():
    for g in Gate, d in Device:
        if g.device == d and g.width == 2 and none(p in PhysicalQubit: p.index == g.b and p.device == g.device):
            g.off_device = True

@rule
def not_native():
    # The device runs this operation on no qubit at all.
    for g in Gate, d in Device:
        if g.device == d and none(n in Native: n.op == g.op and n.device == g.device):
            g.not_native = True

@rule
def not_here_one():
    # The device runs the operation and has the qubit, but not the one on the other.
    for g in Gate:
        if g.width == 1 and g.op_on_device == True and g.first_on_device == True and none(n in Native: n.a == g.a and n.op == g.op and n.device == g.device):
            g.not_here = True

@rule
def not_here_pair():
    # The device runs the operation and has both qubits, but runs it on this
    # pair in neither order.
    for g in Gate:
        if g.width == 2 and g.op_on_device == True and g.first_on_device == True and g.second_on_device == True and none(n in Native: n.a == g.a and n.b == g.b and n.op == g.op and n.device == g.device) and none(n in Native: n.a == g.b and n.b == g.a and n.op == g.op and n.device == g.device):
            g.not_here = True

@rule
def reversed_pair():
    # The device runs the operation on this pair the other way round only.
    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

@rule
def illegal_off_device():
    for g in Gate:
        if g.off_device == True:
            g.illegal = True

@rule
def illegal_not_native():
    for g in Gate:
        if g.not_native == True:
            g.illegal = True

@rule
def illegal_not_here():
    for g in Gate:
        if g.not_here == True:
            g.illegal = True

@rule
def illegal_reversed():
    for g in Gate:
        if g.reversed == True:
            g.illegal = True

@constraint
def every_gate_legal():
    for g in Gate:
        if g.illegal == True:
            violated

# ── Result rules ─────────────────────────────────────────────────────────────
#
# A result holds while every circuit and device it names is, in the store, what
# its producer read: each named digest equals the store's own. Each part is
# concluded current on its own, a part the result does not name counting as
# current, and `holds` joins the three, so its proof lists every digest the
# result was checked against. A result that does not hold has a reason, each
# its own field so `why` names the one that holds:
#
#     unrecorded       the result omits what it needs to be checked: a circuit,
#                      an outcome, or the digest of a circuit or device it names
#     circuit_unknown  the circuit it is about has no digest in the store
#     circuit_changed  that circuit's digest differs from the one it read
#     source_unknown   likewise for the circuit it came from
#     source_changed
#     device_unknown   the device it names has no digest in the store
#     recalibrated     that device's digest differs from the one it read
#
# `holds` also asks that no reason holds. That matters only when a digest holds
# two values, as two sources that disagree leave it: a result then holds
# through neither, and the reason names the other value.

@rule
def circuit_current():
    for r in Result, c in Circuit:
        if r.circuit == c and r.digest == c.digest:
            r.circuit_current = True

@rule
def source_current():
    for r in Result, c in Circuit:
        if r.source == c and r.source_digest == c.digest:
            r.source_current = True

@rule
def source_none():
    # A result derived from no circuit has no source to check.
    for r in Result:
        if r.outcome == r.outcome and absent(r.source):
            r.source_current = True

@rule
def device_current():
    for r in Result, d in Device:
        if r.device == d and r.device_digest == d.digest:
            r.device_current = True

@rule
def device_none():
    # A result computed for no device has no calibration to check.
    for r in Result:
        if r.outcome == r.outcome and absent(r.device):
            r.device_current = True

@rule
def holds():
    for r in Result:
        if r.circuit_current == True and r.source_current == True and r.device_current == True and absent(r.stopped):
            r.holds = r.outcome

@rule
def unrecorded_circuit():
    for r in Result:
        if r.outcome == r.outcome and absent(r.circuit):
            r.unrecorded = True

@rule
def unrecorded_outcome():
    for r in Result:
        if r.circuit == r.circuit and absent(r.outcome):
            r.unrecorded = True

@rule
def unrecorded_digest():
    for r in Result:
        if r.circuit == r.circuit and absent(r.digest):
            r.unrecorded = True

@rule
def unrecorded_source_digest():
    for r in Result:
        if r.source == r.source and absent(r.source_digest):
            r.unrecorded = True

@rule
def unrecorded_device_digest():
    for r in Result:
        if r.device == r.device and absent(r.device_digest):
            r.unrecorded = True

@rule
def circuit_unknown():
    for r in Result:
        if r.circuit == r.circuit and none(c in Circuit: r.circuit == c and c.digest == c.digest):
            r.circuit_unknown = True

@rule
def circuit_changed():
    for r in Result, c in Circuit:
        if r.circuit == c and r.digest != c.digest:
            r.circuit_changed = True

@rule
def source_unknown():
    for r in Result:
        if r.source == r.source and none(c in Circuit: r.source == c and c.digest == c.digest):
            r.source_unknown = True

@rule
def source_changed():
    for r in Result, c in Circuit:
        if r.source == c and r.source_digest != c.digest:
            r.source_changed = True

@rule
def device_unknown():
    for r in Result:
        if r.device == r.device and none(d in Device: r.device == d and d.digest == d.digest):
            r.device_unknown = True

@rule
def recalibrated():
    for r in Result, d in Device:
        if r.device == d and r.device_digest != d.digest:
            r.recalibrated = True

@rule
def stopped_unrecorded():
    for r in Result:
        if r.unrecorded == True:
            r.stopped = True

@rule
def stopped_circuit_unknown():
    for r in Result:
        if r.circuit_unknown == True:
            r.stopped = True

@rule
def stopped_circuit_changed():
    for r in Result:
        if r.circuit_changed == True:
            r.stopped = True

@rule
def stopped_source_unknown():
    for r in Result:
        if r.source_unknown == True:
            r.stopped = True

@rule
def stopped_source_changed():
    for r in Result:
        if r.source_changed == True:
            r.stopped = True

@rule
def stopped_device_unknown():
    for r in Result:
        if r.device_unknown == True:
            r.stopped = True

@rule
def stopped_recalibrated():
    for r in Result:
        if r.recalibrated == True:
            r.stopped = True

# ── Queries ──────────────────────────────────────────────────────────────────

def after(g: string) -> list<string>:
    # Every gate that must run after gate `g`.
    try:
        later = list_attributes(MustPrecede.g)
        return sorted(later)
    except refuse:
        # No gate of this store must run after another yet.
        return []


def commuting(g: string) -> list<string>:
    # Every later gate that shares a qubit with `g` and may run before it.
    try:
        later = list_attributes(CommutesWith.g)
        return sorted(later)
    except refuse:
        # No two gates of this store commute yet.
        return []


def uncoupled_gates() -> list<string>:
    # Every two-qubit gate the chip cannot run as written.
    gs = subjects(Gate.uncoupled, True)
    return sorted(gs)


def illegal_gates(c: string) -> list<string>:
    # Every gate of circuit `c` its device cannot run as written, as
    # "<gate> <reason>", once for each of the four reasons that holds.
    of_c = subjects(Gate.circuit, c)
    mine = sorted(of_c)
    out = []
    for reason in ["off_device", "not_native", "not_here", "reversed"]:
        flagged = subjects(Gate.reason, True)
        for g in flagged:
            if g in mine:
                out = out + [g + " " + reason]
    return sorted(out)


def stale_gates(c: string, at: datetime, max_age: duration) -> list<string>:
    # Every gate of circuit `c` whose native operation was not calibrated as of
    # `at`: measured after `at`, more than `max_age` before it, or never. The
    # store keeps no clock, so the moment is yours to name, and the same
    # question gives the same answer every time. A gate the device cannot run
    # at all is in `illegal_gates`, not here.
    out = []
    gs = subjects(Gate.circuit, c)
    for g in gs:
        try:
            ns = list_attributes(RunsAs.g)
            for n in ns:
                m = current(Native.n.measured_at)
                match m:
                    case active_claim:
                        match m.value:
                            case datetime as when:
                                if when > at or when + max_age < at:
                                    out = out + [g]
                            case default:
                                pass
                    case empty:
                        out = out + [g]
        except refuse:
            # No gate of this store runs natively yet, so there is no
            # calibration to judge.
            pass
    return sorted(out)


def results_not_holding() -> list<string>:
    # Every result that does not hold, as "<result> <reason>", once for each
    # reason that holds.
    out = []
    for reason in ["unrecorded", "circuit_unknown", "circuit_changed", "source_unknown",
                   "source_changed", "device_unknown", "recalibrated"]:
        flagged = subjects(Result.reason, True)
        for r in flagged:
            out = out + [r + " " + reason]
    return sorted(out)


def results_of(c: string) -> list<string>:
    # Every result that read circuit `c`, as the circuit it is about or as the
    # one it came from: the results a change to `c` stops.
    about = subjects(Result.circuit, c)
    came = subjects(Result.source, c)
    out = sorted(about)
    for r in came:
        if r not in out:
            out = out + [r]
    return sorted(out)


def produced(r: string) -> list<string>:
    # What result `r` records about what produced it, as "<field> <value>", one
    # line for each field it records, in this order.
    out = []
    for f in ["circuit", "digest", "source", "source_digest", "device", "device_digest",
              "producer", "settings", "seed", "outcome"]:
        v = current(Result.r.f)
        match v:
            case active_claim:
                out = out + [f + " " + str(v.value)]
            case empty:
                pass
    return out