#!/usr/bin/env python3
"""Exact local audit of the accompanying note, using the standard library.

No numerical sampling, manifold census, or approximation is used. The
universal rank/geometric/quantum conclusions are proved in the note, not
inferred from these checks. Graph tests exercise loops and multiple edges.
"""
from dataclasses import dataclass
from fractions import Fraction as F
from itertools import combinations


@dataclass(frozen=True)
class G:
    re: F = F(0)
    im: F = F(0)

    def __post_init__(self):
        object.__setattr__(self, "re", F(self.re))
        object.__setattr__(self, "im", F(self.im))

    def __add__(self, other):
        o = gauss(other)
        return G(self.re + o.re, self.im + o.im)

    __radd__ = __add__

    def __neg__(self):
        return G(-self.re, -self.im)

    def __sub__(self, other):
        return self + (-gauss(other))

    def __rsub__(self, other):
        return gauss(other) + (-self)

    def __mul__(self, other):
        o = gauss(other)
        return G(self.re*o.re-self.im*o.im, self.re*o.im+self.im*o.re)

    __rmul__ = __mul__

    def conj(self):
        return G(self.re, -self.im)

    def inv(self):
        n = self.re**2 + self.im**2
        if not n:
            raise ZeroDivisionError
        return G(self.re/n, -self.im/n)

    def __truediv__(self, other):
        return self*gauss(other).inv()

    def __rtruediv__(self, other):
        return gauss(other)*self.inv()

    def __pow__(self, n):
        if n < 0:
            return self.inv()**(-n)
        out = G(1)
        for _ in range(n):
            out = out*self
        return out

    def __bool__(self):
        return bool(self.re or self.im)


def gauss(x):
    return x if isinstance(x, G) else G(x)


O, I = G(1), G(0, 1)
checks = []


def check(name, left, right=True):
    assert left == right, (name, left, right)
    checks.append(name)
    print("PASS:", name)


def rank(matrix):
    m = [[gauss(x) for x in row] for row in matrix]
    r = 0
    for c in range(len(m[0])):
        p = next((j for j in range(r, len(m)) if m[j][c]), None)
        if p is None:
            continue
        m[r], m[p] = m[p], m[r]
        pivot = m[r][c]
        m[r] = [x/pivot for x in m[r]]
        for j in range(len(m)):
            if j != r:
                factor = m[j][c]
                m[j] = [x-factor*y for x, y in zip(m[j], m[r])]
        r += 1
        if r == len(m):
            break
    return r


def prod(values):
    out = O
    for x in values:
        out *= x
    return out


# Exact Gaussian-rational point on the universal folded-shape family.
x1, w1 = G(F(3, 4), F(1, 4)), G(F(-3, 5), F(-4, 5))
x2, w2 = 1/x1, I*x1*w1
Z1, Z2 = x1**2, x2**2
D1, D2 = 1/(Z1-1), 1/(Z2-1)
check("first tetrahedron equation", x1**-2+w1**2-1, G())
check("second tetrahedron equation", x2**-2+w2**2-1, G())
check("first moment equation", x1*x2, O)
check("second moment equation", w2/(x1*w1), I)
check("mixed orientation is strict", Z1.im > 0 and Z2.im < 0)
check("no zero or forbidden shapes", all(z and z != O for z in (Z1, Z2)))
check("folded tangent identity", D1+D2, G(-1))
for z, d in ((Z1, D1), (Z2, D2)):
    check("exact imaginary-part identity at " + str(z),
          d.im, -z.im/((z.re-1)**2+z.im**2))

P, Q = [[1, 1], [-1, 0]], [[0, 0], [-1, 1]]
PQ = [[sum(P[i][k]*Q[j][k] for k in range(2)) for j in range(2)] for i in range(2)]
QP = [[sum(Q[i][k]*P[j][k] for k in range(2)) for j in range(2)] for i in range(2)]
check("isotropic moment rows", PQ, QP)
check("full moment row rank", rank([P[i]+Q[i] for i in range(2)]), 2)
restricted = [[gauss(P[i][j])+Q[i][j]*[D1, D2][j] for j in range(2)] for i in range(2)]
check("mixed tangent rank drops", rank(restricted), 1)
jac = [[-2*x1**-2, 0, 2*w1**2, 0],
       [0, -2*x2**-2, 0, 2*w2**2], P[0]+Q[0], P[1]+Q[1]]
check("joint Jacobian has rank three", rank(jac), 3)

# Symbolic Laurent identity after x2=x1^-1, w2=i*x1*w1.
# Polynomials in x1,w1, integer coefficients, represented by exponent dicts.
f1 = {(-2, 0): 1, (0, 2): 1, (0, 0): -1}
f2_substituted = {(2, 0): 1, (2, 2): -1, (0, 0): -1}
minus_xsq_f1 = {(a+2, b): -c for (a, b), c in f1.items()}
check("universal folded syzygy f2=-x1^2 f1", f2_substituted, minus_xsq_f1)

# Local PP branch factors: exact fourth roots of unity, no Weyl changes.
check("outgoing triangle branch factor", I**2, G(-1))
check("incoming triangle branch factor", (-I)**2, G(-1))
check("face bigon branch factor", I*(-I), O)
scalars = [I, I, -I, -I]
check("Euler tetrahedron half-face product", prod(scalars), O)
# Faces are indexed by the opposite vertex. Edge ab uses the other two faces.
edge_sign = {}
for a, b in combinations(range(4), 2):
    faces = [f for f in range(4) if f not in (a, b)]
    edge_sign[(a, b)] = prod(scalars[f] for f in faces)
for a, b in ((0, 1), (0, 2), (0, 3)):
    c, d = [v for v in range(4) if v not in (a, b)]
    check("opposite-root signs " + str((a, b)), edge_sign[(a, b)], edge_sign[(c, d)])
for v in range(4):
    check("vertex central branch factor " + str(v),
          prod(edge_sign[tuple(sorted((v, w)))] for w in range(4) if w != v), G(-1))
check("Lagrangian squares unchanged", all(e**2 == O for e in edge_sign.values()))


def euler_orient(n, edges):
    """Orient closed walks in an even multigraph; endpoint occurrences retain identity."""
    incident = [[] for _ in range(n)]
    for e, (a, b) in enumerate(edges):
        incident[a].append((e, 0))
        incident[b].append((e, 1))
    unused = set(range(len(edges)))
    oriented = {}
    while unused:
        seed = min(unused)
        start = edges[seed][0]
        v = start
        while True:
            choice = next(((e, h) for e, h in incident[v] if e in unused), None)
            if choice is None:
                assert v == start
                break
            e, h = choice
            unused.remove(e)
            oriented[e] = (h, 1-h)
            v = edges[e][1-h]
    return oriented


def audit_graph(label, n, edges):
    orientation = euler_orient(n, edges)
    outgoing, incoming = [0]*n, [0]*n
    half = {}
    for e, (a, b) in enumerate(edges):
        h, other = orientation[e]
        vertices = (a, b)
        outgoing[vertices[h]] += 1
        incoming[vertices[other]] += 1
        half[(e, h)] = I
        half[(e, other)] = -I
        assert half[(e, 0)]*half[(e, 1)] == O
    check(label + ": every vertex has two arrivals/departures",
          outgoing == [2]*n and incoming == [2]*n)
    for v in range(n):
        occurrences = [half[(e, h)] for e, pair in enumerate(edges)
                       for h in (0, 1) if pair[h] == v]
        check(label + ": half-face product at " + str(v), prod(occurrences), O)


audit_graph("two self-glued face pairs", 1, [(0, 0), (0, 0)])
audit_graph("four parallel faces", 2, [(0, 1)]*4)
audit_graph("loops with parallel faces", 2, [(0, 0), (1, 1), (0, 1), (0, 1)])
check("original PP common scalar takes c_B=-i to 1", I*(-I), O)
check("shape factor for original-to-default vertex", I**3*G(-1), I)

print("\nAll exact checks passed.")
print("The rank lemma, signed geometric lift, and quotient conclusions are proved in the note.")
print("No test count is used as a measure of mathematical significance.")