#!/usr/bin/env python3
"""Supplement 1 checker.

Canonical affine readouts, a rank-losing integer layer, and exact transfer
between rotated copies of the twelve-node rule.

Python 3.10+, standard library only. No floating-point acceptance test.
SageMath is not called. A Smith factorization is accepted only after this
program checks U W V = D and det U, det V in {1, -1}.
The canonical sine is named and is not evaluated as a decimal.
The degree formula for Q(cos(2π/m)) is cited, not reproved.
"""
from __future__ import annotations

import json
from math import gcd
from pathlib import Path

ROOT = Path(__file__).resolve().parent


def fail(label: str) -> None:
    raise SystemExit(label)


def need(ok: bool, label: str) -> None:
    if not ok:
        fail(label)


class Q:
    __slots__ = ("n", "d")

    def __init__(self, n: int, d: int = 1):
        if isinstance(n, Q):
            n, d = n.n, n.d
        if d == 0:
            raise ZeroDivisionError("zero denominator")
        if d < 0:
            n, d = -n, -d
        g = gcd(n, d)
        self.n = n // g
        self.d = d // g

    def __eq__(self, other: object) -> bool:
        return isinstance(other, Q) and self.n == other.n and self.d == other.d

    def __hash__(self) -> int:
        return hash((self.n, self.d))

    def __add__(self, other: "Q") -> "Q":
        other = other if isinstance(other, Q) else Q(other)
        return Q(self.n * other.d + other.n * self.d, self.d * other.d)

    def __radd__(self, other: int) -> "Q":
        return self + Q(other)

    def __neg__(self) -> "Q":
        return Q(-self.n, self.d)

    def __sub__(self, other: "Q | int") -> "Q":
        return self + (-(other if isinstance(other, Q) else Q(other)))

    def __rsub__(self, other: int) -> "Q":
        return Q(other) - self

    def __mul__(self, other: "Q | int") -> "Q":
        other = other if isinstance(other, Q) else Q(other)
        return Q(self.n * other.n, self.d * other.d)

    def __rmul__(self, other: int) -> "Q":
        return self * other

    def __truediv__(self, other: "Q | int") -> "Q":
        other = other if isinstance(other, Q) else Q(other)
        return Q(self.n * other.d, self.d * other.n)

    def text(self) -> str:
        return f"{self.n}/{self.d}"


def qparse(text: str) -> Q:
    if not isinstance(text, str) or "/" not in text:
        fail("rational text")
    num, den = text.split("/")
    if not _canon_int(num) or not _canon_int(den):
        fail("rational text")
    value = Q(int(num), int(den))
    need(value.text() == text, "reduced coefficients")
    return value


def _canon_int(text: str) -> bool:
    if text == "0":
        return True
    if text.startswith("-"):
        body = text[1:]
    else:
        body = text
    return bool(body) and body.isdigit() and not (len(body) > 1 and body.startswith("0"))


class G:
    """a + b*phi, phi**2 = phi + 1."""

    __slots__ = ("a", "b")

    def __init__(self, a: Q | int = 0, b: Q | int = 0):
        self.a = a if isinstance(a, Q) else Q(int(a))
        self.b = b if isinstance(b, Q) else Q(int(b))

    def __eq__(self, other: object) -> bool:
        return isinstance(other, G) and self.a == other.a and self.b == other.b

    def __add__(self, other: "G") -> "G":
        return G(self.a + other.a, self.b + other.b)

    def __neg__(self) -> "G":
        return G(-self.a, -self.b)

    def __sub__(self, other: "G") -> "G":
        return self + (-other)

    def __mul__(self, other: "G") -> "G":
        return G(
            self.a * other.a + self.b * other.b,
            self.a * other.b + self.b * other.a + self.b * other.b,
        )

    def __truediv__(self, other: "G") -> "G":
        norm = other.a * other.a + other.a * other.b - other.b * other.b
        if norm == Q(0):
            raise ZeroDivisionError("zero in Q(phi)")
        return self * G((other.a + other.b) / norm, -other.b / norm)

    def scale(self, k: Q) -> "G":
        return G(self.a * k, self.b * k)


def gtext(value: G) -> list[str]:
    return [value.a.text(), value.b.text()]


ZERO, ONE, PHI = G(0), G(1), G(0, 1)
LABELS = ("d2", "d1", "d", "b", "b1", "b2")
# z = A + s, with A in Q(phi). s = sin(pi/1296000) is not evaluated.
# phi**2 = phi + 1, so -s - phi**2 = (-1 - phi) + (-1)*s.
NODES = (
    (G(-1, -1), -1),
    (G(-1, 0), -1),
    (G(0, 0), -1),
    (G(0, 0), 1),
    (G(1, 0), 1),
    (G(1, 1), 1),
)


def reading(alpha: G, beta: G, node: tuple[G, int]) -> tuple[G, G]:
    base, sign = node
    return alpha + beta * base, beta.scale(Q(sign))


def all_readings(alpha: G, beta: G) -> list[tuple[G, G]]:
    return [reading(alpha, beta, node) for node in NODES]


def pack(pair: tuple[G, G]) -> list[str]:
    real, coeff = pair
    return [real.a.text(), real.b.text(), coeff.a.text(), coeff.b.text()]


def unpack(values: list[str]) -> tuple[G, G]:
    need(isinstance(values, list) and len(values) == 4, "four coefficients")
    a0, a1, b0, b1 = (qparse(item) for item in values)
    return G(a0, a1), G(b0, b1)


def recover(readings: list[tuple[G, G]]) -> tuple[G, G]:
    left, right = readings[0], readings[-1]
    alpha = (left[0] + right[0]).scale(Q(1, 2))
    beta = (right[1] - left[1]).scale(Q(1, 2))
    return alpha, beta


def record_for(alpha: G, beta: G) -> dict:
    rows = all_readings(alpha, beta)
    return {
        "schema": "CANONICAL_AFFINE_READOUT/1",
        "sine": "sin(pi/1296000)",
        "frame": "normalized horizontal marks",
        "labels": list(LABELS),
        "field": "affine",
        "readings": [{"label": label, "coefficients": pack(row)} for label, row in zip(LABELS, rows)],
    }


def load_record(payload: dict) -> tuple[G, G]:
    need(isinstance(payload, dict), "record object")
    schema = payload.get("schema")
    if schema == "GEOMETRY_LEDGER_RECORD/1.0":
        fail("extinction ledger is a different format")
    if schema == "BYTE_LEDGER/23":
        fail("fixed 23-step byte ledger is a different format")
    if schema != "CANONICAL_AFFINE_READOUT/1":
        fail("schema")
    need(payload.get("sine") == "sin(pi/1296000)", "canonical sine binding")
    need(payload.get("frame") == "normalized horizontal marks", "frame")
    need(payload.get("labels") == list(LABELS), "label order")
    need(payload.get("field") == "affine", "field")
    rows = payload.get("readings")
    need(isinstance(rows, list) and len(rows) == 6, "six readings")
    parsed = []
    for label, row in zip(LABELS, rows):
        need(isinstance(row, dict) and row.get("label") == label, "label order")
        parsed.append(unpack(row.get("coefficients")))
    alpha, beta = recover(parsed)
    need(all_readings(alpha, beta) == parsed, "recalculated readings disagree")
    return alpha, beta


def symbolic_inverse() -> None:
    """The outer formulas invert p(z)=alpha+beta*z for every coefficient in Q."""
    # Each coefficient is a linear form in (alpha0, alpha1, beta0, beta1).
    eye = tuple(tuple(Q(1 if i == j else 0) for j in range(4)) for i in range(4))

    def add(u, v):
        return tuple(u[i] + v[i] for i in range(4))

    def scale(u, k: Q):
        return tuple(u[i] * k for i in range(4))

    def gadd(x, y):
        return (add(x[0], y[0]), add(x[1], y[1]))

    def gscale(x, k: Q):
        return (scale(x[0], k), scale(x[1], k))

    def gmul_const(x, const: G):
        # (a0+a1 phi)(c0+c1 phi) with a linear and c constant.
        a0, a1 = x
        return (
            add(scale(a0, const.a), scale(a1, const.b)),
            add(add(scale(a0, const.b), scale(a1, const.a)), scale(a1, const.b)),
        )

    alpha = (eye[0], eye[1])
    beta = (eye[2], eye[3])
    parsed = []
    for base, sign in NODES:
        real = gadd(alpha, gmul_const(beta, base))
        coeff = gscale(beta, Q(sign))
        parsed.append((real, coeff))
    left, right = parsed[0], parsed[-1]
    got_a = gscale(gadd(left[0], right[0]), Q(1, 2))
    got_b = gscale((add(right[1][0], scale(left[1][0], Q(-1))), add(right[1][1], scale(left[1][1], Q(-1)))), Q(1, 2))
    need(got_a == alpha and got_b == beta, "symbolic inverse")


def independence_arithmetic() -> dict:
    """Integer half of the independence argument. The degree formula is cited."""
    m = 1296000
    need(m == 2**7 * 3**4 * 5**3, "factor 1296000")
    phi_m = m // 2 * 2 // 3 * 4 // 5
    need(phi_m == 345600, "euler totient")
    degree = phi_m // 2
    need(degree == 172800 and degree > 2, "cited real degree")
    # cos(2π/5) = (sqrt(5)-1)/4, so sqrt(5) = 4*cos(2π/5)+1.
    # Checked here only as an identity in Q(phi): phi = (1+sqrt(5))/2.
    sqrt5 = G(-1, 2)  # -1 + 2 phi = sqrt(5)
    need(sqrt5 * sqrt5 == G(5), "sqrt(5) squared")
    return {
        "m": m,
        "euler_totient": phi_m,
        "cited_degree_of_Q_cos": degree,
        "classical_formula": "For m >= 3, [Q(cos(2*pi/m)):Q] = phi(m)/2. Not reproved here.",
        "conclusion": "sin(pi/1296000) is not in Q(phi). The set {1, phi, s, phi*s} is a Q-basis. s*s leaves that span.",
    }


def canonical_checks() -> dict:
    need(PHI * PHI == PHI + ONE, "golden relation")
    symbolic_inverse()
    alpha, beta = G(3), G(-2)
    record = record_for(alpha, beta)
    outers = (record["readings"][0]["coefficients"], record["readings"][-1]["coefficients"])
    need(outers == (
        ["5/1", "2/1", "2/1", "0/1"],
        ["1/1", "-2/1", "-2/1", "0/1"],
    ), "outer readings of 3-2z")
    got_a, got_b = load_record(json.loads(json.dumps(record)))
    need((got_a, got_b) == (alpha, beta), "recover 3 and -2")
    fraction_a, fraction_b = G(Q(1, 2), Q(3, 7)), G(Q(-4, 3), Q(5, 2))
    need(load_record(record_for(fraction_a, fraction_b)) == (fraction_a, fraction_b), "fractional coordinates")
    rejected = 0
    for bad in (
        {"schema": "GEOMETRY_LEDGER_RECORD/1.0"},
        {"schema": "BYTE_LEDGER/23", "horizon": 23},
        {"schema": "CANONICAL_AFFINE_READOUT/1", "sine": "1/2"},
        {**record, "readings": record["readings"][:-1]},
        {**record, "readings": [{**record["readings"][0], "coefficients": ["6/1", "2/1", "2/1", "0/1"]}] + record["readings"][1:]},
    ):
        try:
            load_record(bad)
        except SystemExit:
            rejected += 1
    need(rejected == 5, "rejected foreign or inconsistent records")
    path = ROOT / "canonical-readout.example.json"
    text = json.dumps(record, indent=2) + "\n"
    if not path.exists():
        path.write_text(text, encoding="utf-8")
    need(path.read_text(encoding="utf-8") == text, "hosted example")
    return {
        "example_field": "p(z)=3-2z",
        "recovered": [gtext(got_a), gtext(got_b)],
        "symbolic_inverse": True,
        "fractional_example": True,
        "rejected_records": rejected,
        "outer_readings": [list(outers[0]), list(outers[1])],
        "rational_independence": independence_arithmetic(),
    }


def egcd(a: int, b: int) -> tuple[int, int, int]:
    if b == 0:
        return (1, 0, a) if a >= 0 else (-1, 0, -a)
    _q, r = divmod(a, b)
    x, y, g = egcd(b, r)
    return y, x - _q * y, g


def det_int(matrix: list[list[int]]) -> int:
    n = len(matrix)
    need(n > 0 and all(len(row) == n for row in matrix), "square determinant")
    m = [row[:] for row in matrix]
    sign = 1
    prev = 1
    for k in range(n - 1):
        if m[k][k] == 0:
            swap = next((r for r in range(k + 1, n) if m[r][k] != 0), None)
            if swap is None:
                return 0
            m[k], m[swap] = m[swap], m[k]
            sign = -sign
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                numerator = m[k][k] * m[i][j] - m[i][k] * m[k][j]
                need(numerator % prev == 0, "bareiss divisibility")
                m[i][j] = numerator // prev
            m[i][k] = 0
        prev = m[k][k]
    return sign * m[-1][-1]


def matmul(a, b):
    p, q, r = len(a), len(b[0]), len(b)
    need(all(len(row) == r for row in a), "multiply")
    return [[sum(a[i][k] * b[k][j] for k in range(r)) for j in range(q)] for i in range(p)]


def smith(matrix: list[list[int]]):
    if not matrix or not matrix[0]:
        fail("empty matrix")
    a = [list(map(int, row)) for row in matrix]
    p, m = len(a), len(a[0])
    need(all(len(row) == m for row in a), "rectangular matrix")
    u = [[1 if i == j else 0 for j in range(p)] for i in range(p)]
    v = [[1 if i == j else 0 for j in range(m)] for i in range(m)]

    def swap_rows(i, j):
        if i != j:
            a[i], a[j] = a[j], a[i]
            u[i], u[j] = u[j], u[i]

    def swap_cols(i, j):
        if i == j:
            return
        for row in a:
            row[i], row[j] = row[j], row[i]
        for row in v:
            row[i], row[j] = row[j], row[i]

    def add_row(src, dst, factor):
        if factor == 0:
            return
        for j in range(m):
            a[dst][j] -= factor * a[src][j]
        for j in range(p):
            u[dst][j] -= factor * u[src][j]

    def add_col(src, dst, factor):
        if factor == 0:
            return
        for i in range(p):
            a[i][dst] -= factor * a[i][src]
        for i in range(m):
            v[i][dst] -= factor * v[i][src]

    def negate_row(i):
        for j in range(m):
            a[i][j] = -a[i][j]
        for j in range(p):
            u[i][j] = -u[i][j]

    t = 0
    while t < min(p, m):
        guard = 0
        while True:
            guard += 1
            need(guard < 10000, "smith did not terminate")
            pivot = None
            for i in range(t, p):
                for j in range(t, m):
                    if a[i][j] != 0 and (pivot is None or abs(a[i][j]) < abs(a[pivot[0]][pivot[1]])):
                        pivot = (i, j)
            if pivot is None:
                t = min(p, m)
                break
            swap_rows(t, pivot[0])
            swap_cols(t, pivot[1])
            if a[t][t] < 0:
                negate_row(t)
            progressed = False
            for i in range(t + 1, p):
                if a[i][t] != 0:
                    add_row(t, i, a[i][t] // a[t][t])
                    if a[i][t] != 0:
                        swap_rows(t, i)
                    progressed = True
                    break
            if progressed:
                continue
            for j in range(t + 1, m):
                if a[t][j] != 0:
                    add_col(t, j, a[t][j] // a[t][t])
                    if a[t][j] != 0:
                        swap_cols(t, j)
                    progressed = True
                    break
            if progressed:
                continue
            bad = None
            for i in range(t + 1, p):
                for j in range(t + 1, m):
                    if a[i][j] % a[t][t] != 0:
                        bad = j
                        break
                if bad is not None:
                    break
            if bad is None:
                break
            add_col(bad, t, -1)
        if t >= min(p, m):
            break
        t += 1
    rank = 0
    while rank < min(p, m) and a[rank][rank] != 0:
        rank += 1
    for i in range(rank - 1):
        guard = 0
        while a[i + 1][i + 1] % a[i][i] != 0:
            guard += 1
            need(guard < 10000, "divisibility repair")
            add_col(i + 1, i, -1)
            while a[i + 1][i] != 0:
                add_row(i, i + 1, a[i + 1][i] // a[i][i])
                if a[i + 1][i] != 0:
                    swap_rows(i, i + 1)
            for j in range(i + 1, m):
                while a[i][j] != 0:
                    add_col(i, j, a[i][j] // a[i][i])
                    if a[i][j] != 0:
                        swap_cols(i, j)
                        break
                    break
            if a[i][i] < 0:
                negate_row(i)
            if i + 1 < p and a[i + 1][i + 1] < 0:
                negate_row(i + 1)
    for i in range(rank):
        need(a[i][i] > 0, "positive diagonal")
        for j in range(m):
            if j != i:
                need(a[i][j] == 0, "diagonal row")
        for r in range(p):
            if r != i:
                need(a[r][i] == 0, "diagonal column")
    for i in range(rank - 1):
        need(a[i + 1][i + 1] % a[i][i] == 0, "divisibility chain")
    return u, a, v, rank


def certify_smith(original, u, d, v) -> None:
    need(matmul(matmul(u, original), v) == d, "UWV = D")
    need(det_int(u) in (1, -1) and det_int(v) in (1, -1), "unimodular")


def invert_unimodular(matrix):
    n = len(matrix)
    aug = [row[:] + [1 if i == j else 0 for j in range(n)] for i, row in enumerate(matrix)]
    for col in range(n):
        pivot = next((r for r in range(col, n) if aug[r][col] != 0), None)
        need(pivot is not None, "singular unimodular")
        aug[col], aug[pivot] = aug[pivot], aug[col]
        piv = aug[col][col]
        need(abs(piv) == 1, "unimodular pivot")
        if piv == -1:
            aug[col] = [-x for x in aug[col]]
        for r in range(n):
            if r != col and aug[r][col] != 0:
                factor = aug[r][col]
                aug[r] = [x - factor * y for x, y in zip(aug[r], aug[col])]
    return [row[n:] for row in aug]


def floor_div(numer: int, divisor: int) -> tuple[int, int]:
    need(isinstance(divisor, int) and divisor > 0, "positive divisor")
    quotient, remainder = divmod(numer, divisor)
    need(0 <= remainder < divisor, "remainder range")
    return quotient, remainder


def layer(w, b, s0, x):
    p, m = len(w), len(w[0])
    need(len(b) == p and len(x) == m, "layer shapes")
    linear = [sum(w[i][j] * x[j] for j in range(m)) + b[i] for i in range(p)]
    pairs = [floor_div(value, s0) for value in linear]
    y = [item[0] for item in pairs]
    rho = [item[1] for item in pairs]
    return y, rho


def witness_record(w, b, s0, x):
    u, d, v, rank = smith(w)
    certify_smith(w, u, d, v)
    v_inv = invert_unimodular(v)
    coords = [sum(v_inv[i][j] * x[j] for j in range(len(x))) for i in range(len(x))]
    k = coords[rank:]
    y, rho = layer(w, b, s0, x)
    z = [s0 * y[i] + rho[i] - b[i] for i in range(len(b))]
    uz = [sum(u[i][j] * z[j] for j in range(len(z))) for i in range(len(u))]
    need(uz[rank:] == [0] * (len(uz) - rank), "trailing image")
    head = []
    for i in range(rank):
        need(uz[i] % d[i][i] == 0, "diagonal divides")
        head.append(uz[i] // d[i][i])
    recovered = [sum(v[i][j] * value for j, value in enumerate(head + k)) for i in range(len(v))]
    need(recovered == list(x), "recover input")
    payload = {
        "schema": "INTEGER_LAYER_WITNESS/1",
        "s0": str(s0),
        "y": [str(item) for item in y],
        "rho": [str(item) for item in rho],
        "k": [str(item) for item in k],
    }
    encoded = json.dumps(payload, separators=(",", ":")).encode("ascii")
    return {
        "U": u, "D": d, "V": v, "rank": rank, "k": k, "y": y, "rho": rho,
        "recovered": recovered, "bytes": len(encoded), "payload": payload,
    }


def witness_checks() -> dict:
    w = [[1, 1]]
    b = [1]
    s0 = 5
    first = witness_record(w, b, s0, [4, 7])
    second = witness_record(w, b, s0, [5, 6])
    need(first["y"] == second["y"] == [2], "shared quotient")
    need(first["rho"] == second["rho"] == [2], "shared remainder")
    need(first["k"] != second["k"], "witness separates the collision")
    need(first["recovered"] == [4, 7] and second["recovered"] == [5, 6], "both inputs recover")
    without = json.dumps({"schema": "INTEGER_LAYER_WITNESS/1", "s0": "5", "y": ["2"], "rho": ["2"]}, separators=(",", ":"))
    need(b'"k"' not in without.encode("ascii"), "collision omits the witness")
    signed_cases = 0
    for x in ([-3, 1], [-8, 2], [0, -1], [-5, -5]):
        record = witness_record(w, b, s0, list(x))
        need(all(0 <= item < s0 for item in record["rho"]), "signed remainder")
        need(record["recovered"] == list(x), "signed recovery")
        signed_cases += 1
    square = [[2, 4], [6, 8]]
    u, d, v, rank = smith(square)
    certify_smith(square, u, d, v)
    need(rank == 2 and [d[i][i] for i in range(2)] == [2, 4], "example smith diagonal")
    path = ROOT / "witness-collision.example.json"
    text = json.dumps({"left": first["payload"], "right": second["payload"]}, indent=2) + "\n"
    if not path.exists():
        path.write_text(text, encoding="utf-8")
    need(json.loads(path.read_text(encoding="utf-8")) == {"left": first["payload"], "right": second["payload"]}, "hosted witness example")
    return {
        "collision_inputs": [[4, 7], [5, 6]],
        "shared_y": first["y"],
        "shared_rho": first["rho"],
        "witnesses": [first["k"], second["k"]],
        "serialized_bytes": [first["bytes"], second["bytes"]],
        "signed_recoveries": signed_cases,
        "smith_example_diagonal": [2, 4],
        "storage_note": "The serialized witness is sufficient. Its byte length is not claimed to be optimal.",
    }


def vertices():
    nodes = []
    for e in (-1, 1):
        for f in (-1, 1):
            nodes.append((ZERO, G(e), G(0, f)))
            nodes.append((G(e), G(0, f), ZERO))
            nodes.append((G(0, f), ZERO, G(e)))
    return nodes


def rotate(point):
    x, y, z = point
    return ((G(3) * x - G(4) * y) / G(5), (G(4) * x + G(3) * y) / G(5), z)


def features(point):
    x, y, z = point
    return [ONE, x, y, z, x * y, x * z, y * z, x * x - y * y, G(2) * z * z - x * x - y * y]


def gram(nodes):
    rows = [features(point) for point in nodes]
    out = []
    for i in range(9):
        out.append([])
        for j in range(9):
            total = ZERO
            for row in rows:
                total = total + row[i] * row[j]
            out[-1].append(total / G(12))
    return out


def continuous_gram():
    radius2 = PHI + G(2)
    radius4 = radius2 * radius2
    diag = (
        [ONE]
        + [radius2 / G(3)] * 3
        + [radius4 / G(15)] * 3
        + [G(4) * radius4 / G(15), G(4) * radius4 / G(5)]
    )
    return [[diag[i] if i == j else ZERO for j in range(9)] for i in range(9)]


def decode_matrix(nodes, gram_matrix):
    rows = [features(point) for point in nodes]
    out = []
    for mode in range(9):
        out.append([])
        inv = ONE / gram_matrix[mode][mode]
        for sample in range(12):
            out[-1].append(inv * rows[sample][mode] / G(12))
    return out


def transfer(target, source, gram_matrix):
    decoder = decode_matrix(source, gram_matrix)
    target_rows = [features(point) for point in target]
    out = []
    for row in target_rows:
        out.append([])
        for sample in range(12):
            total = ZERO
            for mode in range(9):
                total = total + row[mode] * decoder[mode][sample]
            out[-1].append(total)
    return out


def apply(matrix, vector):
    return [sum((matrix[i][j] * vector[j] for j in range(len(vector))), ZERO) for i in range(len(matrix))]


def transfer_checks() -> dict:
    base = vertices()
    turned = [rotate(point) for point in base]
    turned_twice = [rotate(point) for point in turned]
    radius2 = PHI + G(2)
    for grid in (base, turned, turned_twice):
        need(all(point[0] * point[0] + point[1] * point[1] + point[2] * point[2] == radius2 for point in grid), "same sphere")
    expected = continuous_gram()
    for grid in (base, turned, turned_twice):
        need(gram(grid) == expected, "discrete gram")
    eye = transfer(base, base, expected)
    round_trip = matmul_g(transfer(base, turned, expected), transfer(turned, base, expected))
    need(round_trip == eye, "round trip is the model projection")
    composed = matmul_g(transfer(turned_twice, turned, expected), transfer(turned, base, expected))
    need(composed == transfer(turned_twice, base, expected), "transfer composition")
    basis_ok = 0
    for mode in range(9):
        column = [features(point)[mode] for point in base]
        decoded = apply(decode_matrix(base, expected), column)
        need(decoded == [ONE if i == mode else ZERO for i in range(9)], "basis recovery")
        basis_ok += 1
    ones = [ONE] * 12
    need(apply(transfer(turned, base, expected), ones) == ones, "constants preserved")
    foreign = [G(i - 3, i) for i in range(12)]
    projected = apply(eye, foreign)
    residual = [foreign[i] - projected[i] for i in range(12)]
    need(any(item != ZERO for item in residual), "arbitrary array has a residual")
    need(apply(eye, projected) == projected, "projection is idempotent")
    return {
        "grids": 3,
        "nodes_each": 12,
        "harmonic_degree": 2,
        "basis_recoveries": basis_ok,
        "gram_entries_each": 81,
        "composition": True,
        "round_trip": True,
        "constants_preserved": True,
        "higher_frequency": "not enabled",
    }


def matmul_g(a, b):
    p, q = len(a), len(b[0])
    mid = len(b)
    out = []
    for i in range(p):
        out.append([])
        for j in range(q):
            total = ZERO
            for k in range(mid):
                total = total + a[i][k] * b[k][j]
            out[-1].append(total)
    return out


def main() -> None:
    report = {
        "status": "PASS",
        "release": "supplement-1",
        "groups": {
            "canonical_readouts": canonical_checks(),
            "integer_layer": witness_checks(),
            "grid_transfer": transfer_checks(),
        },
        "limits": [
            "Not a proof-assistant kernel certificate.",
            "The cyclotomic degree formula is cited, not reproved.",
            "The canonical sine is not evaluated.",
            "Witness byte length is not an optimality claim.",
            "Different subdivision frequencies are not enabled.",
            "Items 4, 5, and 6 of the directive are not in this release.",
        ],
    }
    print(json.dumps(report, indent=2, sort_keys=True))


if __name__ == "__main__":
    main()