verify-actual-incidence.py

#!/usr/bin/env python3
"""Exact algebra for the actual-incidence two-angle vacuum theorem.

Standard library only. No topological configurations or slopes are searched.
Default: verify universal identities and a fractional-lift calibration.
Optional: --certificate FILE checks supplied native incidence/angle data.
The optional checker does NOT certify topological realization of matrices.
"""
from collections import defaultdict
from fractions import Fraction as Q
from pathlib import Path
import argparse
import json
import math


class Poly:
    """A small exact commutative polynomial ring with rational coefficients."""
    def __init__(self, terms=None):
        self.terms = {tuple(m): Q(c) for m, c in (terms or {}).items() if c}

    @staticmethod
    def of(x):
        return x if isinstance(x, Poly) else Poly({(): Q(x)})

    @staticmethod
    def var(name):
        return Poly({(name,): 1})

    def __add__(self, other):
        out = defaultdict(Q, self.terms)
        for m, c in Poly.of(other).terms.items():
            out[m] += c
        return Poly(out)

    __radd__ = __add__

    def __neg__(self):
        return Poly({m: -c for m, c in self.terms.items()})

    def __sub__(self, other):
        return self + (-Poly.of(other))

    def __rsub__(self, other):
        return Poly.of(other) - self

    def __mul__(self, other):
        out = defaultdict(Q)
        for m, c in self.terms.items():
            for n, d in Poly.of(other).terms.items():
                out[tuple(sorted(m+n))] += c*d
        return Poly(out)

    __rmul__ = __mul__

    def __truediv__(self, other):
        return self * (Q(1)/Q(other))

    def __pow__(self, n):
        assert isinstance(n, int) and n >= 0
        out = Poly.of(1)
        for _ in range(n):
            out = out*self
        return out

    def __eq__(self, other):
        return self.terms == Poly.of(other).terms


def dot(a, b):
    assert len(a) == len(b)
    return sum(x*y for x, y in zip(a, b))


def add(a, b):
    assert len(a) == len(b)
    return tuple(x+y for x, y in zip(a, b))


def scale(c, a):
    return tuple(c*x for x in a)


def collapse_six(a):
    assert len(a) % 6 == 0
    return tuple(a[j+k]+a[j+k+3] for j in range(0, len(a), 6)
                 for k in range(3))


def P(a):
    assert len(a) % 3 == 0
    return tuple(v for j in range(0, len(a), 3)
                 for v in (a[j]-a[j+1], a[j+2]-a[j+1]))


def J(a):
    assert len(a) % 3 == 0
    return tuple(v for j in range(0, len(a), 3)
                 for v in (a[j+1]-a[j+2], a[j+2]-a[j], a[j]-a[j+1]))


def omega3(a, b):
    return dot(a, J(b))


def omega2(a, b):
    return sum(a[j+1]*b[j]-a[j]*b[j+1] for j in range(0, len(a), 2))


def prime_sum(a):
    return sum(a[1::3])


def delta(a):
    return sum(a[j]*a[j+1]+a[j+1]*a[j+2]+a[j+2]*a[j]
               for j in range(0, len(a), 3))


def normal_class(raw_sum, longitude):
    """Exact X_lambda. No saturation or rounding is allowed."""
    lam = P(raw_sum)
    if any(Q(v).denominator != 1 or int(v) % 2 for v in lam):
        raise ValueError("lambda must belong to the actual lattice AND 2L")
    U = add(longitude, scale(Q(1, 2), raw_sum))
    maxima = tuple(max(U[j:j+3]) for j in range(0, len(U), 3))
    X = tuple(maxima[j//3]-U[j] for j in range(len(U)))
    assert all(Q(v).denominator == 1 and v >= 0 for v in X)
    return U, maxima, X


def universal_checks():
    checks = {}

    def check(name, assertion):
        if not assertion:
            raise AssertionError(name)
        checks[name] = "pass"

    a, b, c, d, e, f = (Poly.var(n) for n in "abcdef")
    check("collapse preserves the skew pairing",
          omega2(P((a, b, c)), P((d, e, f)))
          == omega3((a, b, c), (d, e, f)))
    check("tetrahedron vector is radical", J((1, 1, 1)) == (0, 0, 0))
    check("skew sign", omega3((a, b, c), (d, e, f))
          == -omega3((d, e, f), (a, b, c)))

    E, n, S, rho, T, dd = (Poly.var(v) for v in
                           ("edge_sum", "meridian", "prime_sum", "rho",
                            "max_sum", "double_arcs"))
    char_degree = 2*E+2*n-S
    minimum_sum = rho+S/2-T
    native_degree = 2*rho-char_degree+2*n+2*dd-2*minimum_sum
    normal_degree = 2*(dd-E+T)
    check("native longitude degree equals twice normal Euler degree",
          native_degree == normal_degree)
    check("omitting Weyl commutator leaves a genuine meridian error",
          native_degree-2*n-normal_degree == -2*n
          and -2*n != 0)
    check("omitting the native longitude scalar is not harmless",
          native_degree-2*rho-normal_degree == -2*rho
          and -2*rho != 0)
    unit_degree = -char_degree+2*dd-2*(S/2-T)
    check("unit retains its filling correction",
          unit_degree == 2*(dd-E+T-n))
    shift = Poly.var("dependency_shift")
    check("Euler expression is independent of an edge dependency",
          (E+shift)-(T+shift) == E-T)
    check("complete angle pairing gives minus formal Euler",
          T-E == -(E-T))
    check("all-plus dependency exponent is zero",
          2*E-2*E == 0)

    # This is a local coordinate calibration, NOT a triangulation.
    # Odd raw coordinates give even reduced lambda and half-integral U.
    raw = (1, 3, -1)
    lng = (0, 1, -1)
    U, maxima, X = normal_class(raw, lng)
    check("half-integral lift gives an integral normal class",
          U == (Q(1, 2), Q(5, 2), Q(-3, 2)) and X == (2, 0, 4))
    check("normal-class map recovers the exact even lattice row",
          P(X) == scale(-1, add(P(lng), scale(Q(1, 2), P(raw)))))
    vroot = P(U)
    seed_triple = (-vroot[0], 0, -vroot[1])
    check("maximum/minimum conversion retains fractional translations",
          min(seed_triple) == U[1]-maxima[0]
          and X == tuple(t-min(seed_triple) for t in seed_triple))
    _, _, shifted_X = normal_class(add(raw, (3, 3, 3)), lng)
    check("tetrahedron translation leaves X unchanged", X == shifted_X)
    try:
        normal_class((1, 0, 0), (0, 0, 0))
    except ValueError:
        rejected_odd = True
    else:
        rejected_odd = False
    check("odd root row is not silently admitted", rejected_odd)

    p, q, u, v = (Poly.var(name) for name in ("p", "q", "u", "v"))
    r2 = p*p+q*q
    determinant = p*p-q*(-q)
    check("deformation has exact unit symplectic dual coordinate",
          determinant == r2)
    check("deformation gives meridian angular holonomy two",
          2*determinant == 2*r2)
    check("coordinatewise Cauchy bound is an exact sum of squares",
          (p*v-q*u)**2+(p*u+q*v)**2 == r2*(u*u+v*v))
    check("source-native meridian in the old witness forbids complete angles",
          (1-a-b)+a+3*b == 1+2*b)

    return {
        "scope": "universal exact algebra; no topological instances searched",
        "checks": checks,
        "fractional_lift_calibration": {
            "raw_row": raw, "longitude_row": lng,
            "U": [str(x) for x in U], "X": [int(x) for x in X],
            "double_arc_value": str(delta(X)),
            "topological_realization_asserted": False,
        },
        "proved_formula": "D_l(lambda)=2*(delta(X_lambda)-chi_top(X_lambda))",
        "bound": "p^2+q^2 > max_j((JR_a)_j^2+(JR_d)_j^2)/(theta0_j)^2",
    }


def certificate_check(data):
    """Checks supplied necessary data; this is not a triangulation recognizer."""
    N = data["tetrahedra"]
    assert isinstance(N, int) and N > 0
    edges6 = [tuple(r) for r in data["edges_six_root"]]
    assert len(edges6) == N, "torus-ended ideal triangulation has N edge classes"
    assert all(len(r) == 6*N and all(v in (0, 1) for v in r) for r in edges6)
    assert all(sum(r[j] for r in edges6) == 1 for j in range(6*N)), \
        "every bare edge must occur exactly once in the full edge list"
    edges = [collapse_six(r) for r in edges6]
    basis6 = [data["basis_a_six_root"], data["basis_d_six_root"]]
    basis6 += data.get("other_cusp_basis_six_root", [])
    assert len(basis6) % 2 == 0
    assert all(len(r) == 6*N and all(isinstance(x, int) for x in r)
               for r in basis6)
    basis = [collapse_six(r) for r in basis6]
    for e in edges:
        assert all(omega3(e, f) == 0 for f in edges+basis)
    for i, arow in enumerate(basis):
        for j, brow in enumerate(basis):
            expected = 2 if i % 2 == 0 and j == i+1 else \
                (-2 if i % 2 == 1 and j == i-1 else 0)
            assert omega3(arow, brow) == expected
    theta0 = tuple(Q(x) for x in data["theta_complete"])
    assert len(theta0) == 3*N and all(t > 0 for t in theta0)
    assert all(sum(theta0[j:j+3]) == 1 for j in range(0, 3*N, 3))
    assert all(dot(theta0, e) == 2 for e in edges)
    assert all(dot(theta0, b) == 0 for b in basis)
    p, q = data["slope"]
    assert isinstance(p, int) and isinstance(q, int) and math.gcd(p, q) == 1
    r2 = p*p+q*q
    arow, drow = basis[:2]
    u, v = J(arow), J(drow)
    Rs = add(scale(p, arow), scale(q, drow))
    C2 = max((x*x+y*y)/(t*t) for x, y, t in zip(u, v, theta0))
    constructed = tuple(t+Q(p*y-q*x, r2) for t, x, y in zip(theta0, u, v))
    theta_s = tuple(Q(x) for x in data.get("theta_filling", constructed))
    assert len(theta_s) == 3*N
    affine_ok = (all(sum(theta_s[j:j+3]) == 1 for j in range(0, 3*N, 3))
                 and all(dot(theta_s, e) == 2 for e in edges)
                 and dot(theta_s, Rs) == 2)
    assert affine_ok, "supplied filling angles violate exact native affine equations"
    strict = all(t > 0 for t in theta_s)
    if r2 > C2:
        assert all(t > 0 for t in constructed), "contradiction to proved bound"
    return {
        "scope": "exact necessary incidence and angle checks, conditional on actual native input",
        "topological_realization_verified": False,
        "complete_angular_certificate": True,
        "filling_angular_certificate": strict,
        "two_angle_theorem_applicable_if_input_is_actual": strict,
        "C_squared": str(C2), "slope_squared_norm": r2,
        "outside_uniform_bound": r2 > C2,
        "constructed_filling_angles": [str(x) for x in constructed],
        "angles_used": [str(x) for x in theta_s],
        "warning": "failure of positivity is not target or vacuum vanishing",
    }


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--certificate", type=Path)
    parser.add_argument("--output", type=Path)
    args = parser.parse_args()
    result = universal_checks()
    if args.certificate:
        result["supplied_certificate"] = certificate_check(
            json.loads(args.certificate.read_text()))
    payload = json.dumps(result, indent=2, sort_keys=True)+"\n"
    if args.output:
        args.output.write_text(payload)
    print(payload, end="")


if __name__ == "__main__":
    main()