#!/usr/bin/env python3
"""Dependency-free exact verification for the Clear Onyx continuity snapshot."""

from dataclasses import dataclass
from fractions import Fraction as F
from itertools import combinations, permutations, product


@dataclass(frozen=True)
class Quad:
    """a + b*sqrt(d), with exact rational coefficients."""

    a: F = F(0)
    b: F = F(0)
    d: int = 3

    def __add__(self, other):
        other = q(other, self.d)
        return Quad(self.a + other.a, self.b + other.b, self.d)

    __radd__ = __add__

    def __neg__(self):
        return Quad(-self.a, -self.b, self.d)

    def __sub__(self, other):
        return self + (-q(other, self.d))

    def __rsub__(self, other):
        return q(other, self.d) - self

    def __mul__(self, other):
        other = q(other, self.d)
        return Quad(
            self.a * other.a + self.d * self.b * other.b,
            self.a * other.b + self.b * other.a,
            self.d,
        )

    __rmul__ = __mul__

    def inverse(self):
        den = self.a * self.a - self.d * self.b * self.b
        if den == 0:
            raise ZeroDivisionError
        return Quad(self.a / den, -self.b / den, self.d)

    def __truediv__(self, other):
        return self * q(other, self.d).inverse()

    def __bool__(self):
        return self.a != 0 or self.b != 0


def q(value, d=3):
    if isinstance(value, Quad):
        if value.d != d:
            raise TypeError("quadratic fields do not match")
        return value
    return Quad(F(value), F(0), d)


def matrix(rows, d=3):
    return tuple(tuple(q(x, d) for x in row) for row in rows)


def eye(n, d=3):
    return tuple(tuple(q(1 if i == j else 0, d) for j in range(n)) for i in range(n))


def madd(a, b):
    return tuple(tuple(a[i][j] + b[i][j] for j in range(len(a[0]))) for i in range(len(a)))


def mscale(k, a):
    k = q(k, a[0][0].d)
    return tuple(tuple(k * x for x in row) for row in a)


def mmul(a, b):
    return tuple(
        tuple(sum((a[i][k] * b[k][j] for k in range(len(b))), q(0, a[0][0].d))
              for j in range(len(b[0])))
        for i in range(len(a))
    )


def mpow(a, n):
    out = eye(len(a), a[0][0].d)
    for _ in range(n):
        out = mmul(out, a)
    return out


def rank(rows):
    a = [list(row) for row in rows]
    if not a:
        return 0
    m, n = len(a), len(a[0])
    r = 0
    for c in range(n):
        pivot = next((i for i in range(r, m) if a[i][c]), None)
        if pivot is None:
            continue
        a[r], a[pivot] = a[pivot], a[r]
        inv = a[r][c].inverse() if isinstance(a[r][c], Quad) else F(1, 1) / a[r][c]
        a[r] = [inv * x for x in a[r]]
        for i in range(m):
            if i != r and a[i][c]:
                factor = a[i][c]
                a[i] = [x - factor * y for x, y in zip(a[i], a[r])]
        r += 1
        if r == m:
            break
    return r


def parity(p):
    inversions = sum(p[i] > p[j] for i in range(len(p)) for j in range(i + 1, len(p)))
    return -1 if inversions % 2 else 1


def root_and_escapement():
    d = 5
    I = eye(2, d)
    Q = matrix(((1, 1), (1, 0)), d)
    C = matrix(((0, -1), (1, 0)), d)
    Omega = mscale(-1, C)
    G = madd(mmul(C, Q), mscale(-1, mmul(Q, C)))
    assert mpow(C, 2) == mscale(-1, I)
    assert mpow(C, 4) == I
    assert mpow(Q, 2) == madd(Q, I)
    assert G == matrix(((-2, 1), (1, 2)), d)
    assert mpow(G, 2) == mscale(5, I)

    flattened = [[M[i][j] for M in (I, Q, Omega, G)] for i in range(2) for j in range(2)]
    assert rank(flattened) == 4

    inv_sqrt5 = Quad(F(0), F(1, 5), d)
    H = mscale(inv_sqrt5, G)
    S = mmul(H, C)
    Qinv = matrix(((0, 1), (1, -1)), d)
    Lax = mscale(inv_sqrt5, madd(Q, mscale(2, I)))
    assert mpow(H, 2) == I
    assert mpow(S, 2) == I
    assert mmul(H, S) == C
    assert mmul(S, H) == mscale(-1, C)
    assert mmul(mmul(H, C), H) == mscale(-1, C)
    assert mmul(mmul(H, Q), H) == mscale(-1, Qinv)
    assert mpow(Lax, 2) == mpow(Q, 2)

    for r0 in range(4):
        for n in range(257):
            r, carries = r0, 0
            for _ in range(n):
                if r == 3:
                    carries += 1
                r = (r + 1) % 4
            assert r == (r0 + n) % 4
            assert carries == (r0 + n) // 4


def su3_constants():
    d = 3
    rt3 = Quad(F(0), F(1), d)
    seeds = {
        (0, 1, 2): q(1, d),
        (0, 3, 6): q(F(1, 2), d),
        (1, 3, 5): q(F(1, 2), d),
        (1, 4, 6): q(F(1, 2), d),
        (2, 3, 4): q(F(1, 2), d),
        (0, 4, 5): q(F(-1, 2), d),
        (2, 5, 6): q(F(-1, 2), d),
        (3, 4, 7): rt3 / 2,
        (5, 6, 7): rt3 / 2,
    }
    f = {}
    for base, value in seeds.items():
        for p in permutations(range(3)):
            key = tuple(base[i] for i in p)
            f[key] = q(parity(p), d) * value
    return f


def bracket_vector(f, a, b, n, d=3):
    return tuple(f.get((a, b, c), q(0, d)) for c in range(n))


def su3_and_standard_model():
    d = 3
    zero = q(0, d)
    f = su3_constants()
    pairs8 = list(combinations(range(8), 2))
    beta = [[bracket_vector(f, a, b, 8, d)[row] for a, b in pairs8] for row in range(8)]
    assert rank(beta) == 8
    assert len(pairs8) - rank(beta) == 20

    for a, b, c, out in product(range(8), repeat=4):
        total = zero
        for e in range(8):
            total += (
                f.get((b, c, e), zero) * f.get((a, e, out), zero)
                + f.get((c, a, e), zero) * f.get((b, e, out), zero)
                + f.get((a, b, e), zero) * f.get((c, e, out), zero)
            )
        assert total == zero

    ad = []
    for a in range(8):
        ad.append(tuple(tuple(f.get((a, col, row), zero) for col in range(8)) for row in range(8)))
    for a in range(8):
        for b in range(8):
            killing = sum((mmul(ad[a], ad[b])[i][i] for i in range(8)), zero)
            assert killing == q(-3 if a == b else 0, d)

    def weak(a, b, c):
        if len({a, b, c}) < 3:
            return zero
        return q(parity((a, b, c)), d)

    pairs12 = list(combinations(range(12), 2))
    full = []
    for out in range(12):
        row = []
        for a, b in pairs12:
            if a < 8 and b < 8 and out < 8:
                row.append(f.get((a, b, out), zero))
            elif 8 <= a < 11 and 8 <= b < 11 and 8 <= out < 11:
                row.append(weak(a - 8, b - 8, out - 8))
            else:
                row.append(zero)
        full.append(row)
    assert rank(full) == 11
    assert len(pairs12) - rank(full) == 55
    assert all(x == zero for x in full[11])


def anomalies():
    fields = [
        # name, color dimension, weak dimension, Y, cubic color sign
        ("Q_L", 3, 2, F(1, 6), 1),
        ("u_c", 3, 1, F(-2, 3), -1),
        ("d_c", 3, 1, F(1, 3), -1),
        ("L_L", 1, 2, F(-1, 2), 0),
        ("e_c", 1, 1, F(1), 0),
    ]
    su3_cubic = sum(weak_dim * cubic for _, _, weak_dim, _, cubic in fields)
    su3_sq_u1 = sum(weak_dim * y * F(1, 2) for _, color_dim, weak_dim, y, _ in fields if color_dim == 3)
    su2_sq_u1 = sum(color_dim * y * F(1, 2) for _, color_dim, weak_dim, y, _ in fields if weak_dim == 2)
    grav_sq_u1 = sum(color_dim * weak_dim * y for _, color_dim, weak_dim, y, _ in fields)
    u1_cubic = sum(color_dim * weak_dim * y ** 3 for _, color_dim, weak_dim, y, _ in fields)
    doublets = sum(color_dim for _, color_dim, weak_dim, _, _ in fields if weak_dim == 2)
    assert (su3_cubic, su3_sq_u1, su2_sq_u1, grav_sq_u1, u1_cubic) == (0, 0, 0, 0, 0)
    assert doublets == 4 and doublets % 2 == 0


def invariant_forms():
    # Symmetric B on su(2)+u(1), variables indexed by i<=j.
    variables = [(i, j) for i in range(4) for j in range(i, 4)]
    index = {pair: k for k, pair in enumerate(variables)}

    def eps(i, j, k):
        if max(i, j, k) >= 3 or len({i, j, k}) < 3:
            return F(0)
        return F(parity((i, j, k)))

    def add_B(row, i, j, coefficient):
        pair = (i, j) if i <= j else (j, i)
        row[index[pair]] += coefficient

    equations = []
    for x, z, y in product(range(4), repeat=3):
        row = [F(0) for _ in variables]
        for out in range(4):
            add_B(row, out, y, eps(x, z, out))
            add_B(row, z, out, eps(x, y, out))
        if any(row):
            equations.append(row)
    assert rank(equations) == 8
    su2_scale = [F(1) if pair in ((0, 0), (1, 1), (2, 2)) else F(0) for pair in variables]
    u1_scale = [F(1) if pair == (3, 3) else F(0) for pair in variables]
    for row in equations:
        assert sum(x * y for x, y in zip(row, su2_scale)) == 0
        assert sum(x * y for x, y in zip(row, u1_scale)) == 0
    assert len(variables) - rank(equations) == 2


def bell_chsh():
    values = []
    for a0, a1, b0, b1 in product((-1, 1), repeat=4):
        s = a0 * b0 + a0 * b1 + a1 * b0 - a1 * b1
        values.append(s)
    assert max(abs(x) for x in values) == 2
    quantum = Quad(F(0), F(2), 2)  # 2*sqrt(2)
    assert quantum * quantum == q(8, 2)


def main():
    checks = [
        ("root algebra, mirrors, lift, and escapement", root_and_escapement),
        ("su(3) bracket and 12-to-66-to-11 count", su3_and_standard_model),
        ("one-generation anomaly cancellation", anomalies),
        ("su(2)+u(1) invariant-form obstruction", invariant_forms),
        ("Bell-CHSH local and quantum exact bounds", bell_chsh),
    ]
    for label, fn in checks:
        fn()
        print(f"PASS  {label}")
    print(f"PASS  all {len(checks)} exact verification blocks")


if __name__ == "__main__":
    main()