verify-clean-excess.py

#!/usr/bin/env python3
"""Exact algebra audit for the clean-excess meridian research note.

Python 3.10+, standard library only. Sparse Laurent/Poisson arithmetic over
Q(i), exact Jacobian rank, and a formal-series recurrence. The all-orders
extinction and general clean-intersection theorem are proved in the note,
not inferred from truncation. No manifold, slope, triangulation, or skein
configuration is enumerated.
"""

from dataclasses import dataclass
from fractions import Fraction as F


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

    @staticmethod
    def of(x):
        return x if isinstance(x, G) else G(x)

    def __add__(self, y):
        y = G.of(y)
        return G(self.re + y.re, self.im + y.im)

    __radd__ = __add__

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

    def __sub__(self, y):
        return self + (-G.of(y))

    def __mul__(self, y):
        y = G.of(y)
        return G(self.re*y.re - self.im*y.im,
                 self.re*y.im + self.im*y.re)

    __rmul__ = __mul__

    def __truediv__(self, y):
        y = G.of(y)
        norm = y.re*y.re + y.im*y.im
        assert norm
        return self * G(y.re/norm, -y.im/norm)

    def __pow__(self, n):
        if n < 0:
            return (G(1)/self)**(-n)
        out = G(1)
        while n:
            if n & 1:
                out = out*self
            self = self*self
            n //= 2
        return out

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


ZERO = (0, 0, 0, 0)
OMEGA = ((0, -1, 0, 0), (1, 0, 0, 0),
         (0, 0, 0, -1), (0, 0, 1, 0))


def pair(a, b):
    return sum(a[j]*OMEGA[j][k]*b[k] for j in range(4) for k in range(4))


class L:
    def __init__(self, terms=None):
        self.terms = {tuple(a): G.of(c) for a, c in (terms or {}).items() if c}

    @staticmethod
    def of(x):
        return x if isinstance(x, L) else L({ZERO: G.of(x)})

    @staticmethod
    def mono(a, c=1):
        return L({tuple(a): G.of(c)})

    def __add__(self, other):
        other = L.of(other)
        out = dict(self.terms)
        for a, c in other.terms.items():
            out[a] = out.get(a, G()) + c
        return L(out)

    __radd__ = __add__

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

    def __sub__(self, y):
        return self + (-L.of(y))

    def __rsub__(self, y):
        return L.of(y) - self

    def __mul__(self, other):
        other = L.of(other)
        out = L()
        for a, c in self.terms.items():
            for b, d in other.terms.items():
                out += L.mono(tuple(x+y for x, y in zip(a, b)), c*d)
        return out

    __rmul__ = __mul__

    def __pow__(self, n):
        if n < 0:
            assert len(self.terms) == 1
            (a, c), = self.terms.items()
            return L.mono(tuple(n*x for x in a), c**n)
        out = L.of(1)
        while n:
            if n & 1:
                out = out*self
            self = self*self
            n //= 2
        return out

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

    def logdiff(self, j):
        return L({a: c*a[j] for a, c in self.terms.items()})

    def bracket(self, other):
        other = L.of(other)
        out = L()
        for a, c in self.terms.items():
            for b, d in other.terms.items():
                out += L.mono(tuple(x+y for x, y in zip(a, b)), c*d*pair(a, b))
        return out

    def fold(self):
        """Classical substitution y=x^-1, v=i*x*w; no f1 reduction."""
        out = L()
        for (a, b, c, d), coeff in self.terms.items():
            out += L.mono((a-c+d, b+d, 0, 0), coeff*(G(0, 1)**d))
        return out

    def at(self, p):
        out = G()
        for a, coeff in self.terms.items():
            term = coeff
            for value, exponent in zip(p, a):
                term = term*(G.of(value)**exponent)
            out += term
        return out


def rank(matrix):
    a = [[G.of(x) for x in row] for row in matrix]
    r = 0
    for j in range(len(a[0])):
        pivot = next((i for i in range(r, len(a)) if a[i][j]), None)
        if pivot is None:
            continue
        a[r], a[pivot] = a[pivot], a[r]
        leading = a[r][j]
        a[r] = [x/leading for x in a[r]]
        for i in range(len(a)):
            if i != r:
                coeff = a[i][j]
                a[i] = [x - coeff*y for x, y in zip(a[i], a[r])]
        r += 1
    return r


COUNT = 0


def check(name, condition):
    global COUNT
    if not condition:
        raise AssertionError(name)
    COUNT += 1
    print("PASS:", name)


def conv(a, b, n):
    return [sum((a[j]*b[k-j] for j in range(k+1)
                 if j < len(a) and k-j < len(b)), F(0)) for k in range(n)]


def main():
    print("Clean-excess meridian obstruction: exact algebra audit")
    print("Coefficient field for finite checks: Q(i); no floating point.")
    x, w, y, v = [L.mono(tuple(int(j == k) for j in range(4))) for k in range(4)]
    ii = L.of(G(0, 1))
    f1, f2 = x**-2+w**2-1, y**-2+v**2-1
    U, V = x*y, x**-1*w**-1*v
    g1, g2 = U-1, V-ii
    s = [f1, f2, g1, g2]
    a = [x**2, L.of(1), x**2*(U+1)*U**-2, -x**2*w**2*(V+ii)]
    edge_rows = [(1, 0, 1, 0), (-1, -1, 0, 1)]
    check("the ambient logarithmic Poisson form is nondegenerate", rank(OMEGA) == 4)
    check("the two fixed moment rows are independent and isotropic",
          rank(edge_rows) == 2 and pair(*edge_rows) == 0)
    check("the two tetrahedron symbols Poisson-commute", f1.bracket(f2) == 0)
    check("the two moment symbols Poisson-commute", g1.bracket(g2) == 0)
    check("the displayed excess syzygy is an exact Laurent identity",
          sum((aa*ss for aa, ss in zip(a, s)), L()) == 0)
    check("moment reduction makes the second tetrahedron redundant",
          f2.fold() == -x**2*f1)
    check("moment reduction kills both normalized moment symbols",
          g1.fold() == 0 and g2.fold() == 0)
    p = [G(F(5, 3)), G(F(4, 5)), G(F(3, 5)), G(0, F(4, 3))]
    check("the fixed witness is a nonzero root point", all(p) and not any(t.at(p) for t in s))
    jacobian = [[t.logdiff(j).at(p) for j in range(4)] for t in s]
    check("the full logarithmic Jacobian has rank three", rank(jacobian) == 3)
    check("the tetrahedron and moment blocks each have rank two",
          rank(jacobian[:2]) == rank(jacobian[2:]) == 2)
    beta_scalar = F(1, 2)*(a[0].bracket(f1)+a[1].bracket(f2)
                               +g1.bracket(a[2])+g2.bracket(a[3]))
    check("the first Bockstein scalar is exactly -2*x^2*w^2",
          beta_scalar == -2*x**2*w**2)
    vf = lambda b: a[0]*b.bracket(f1)+a[1]*b.bracket(f2)
    # Four coordinate identities verify equality of the derivations; this is
    # a finite coordinate audit, not a sampling of monomial exponents.
    defects = []
    for b in [x, w, y, v]:
        mixed = F(1, 2)*(a[0]*b.bracket(f1)+a[1]*b.bracket(f2)
                        +a[2]*g1.bracket(b)+a[3]*g2.bracket(b))
        boundary = F(1, 2)*sum((ss*b.bracket(aa) for aa, ss in zip(a, s)), L())
        defects.append(mixed-vf(b)-boundary)
    check("the mixed-side principal symbol equals the tangent field modulo the ideal",
          all(t == 0 for t in defects))
    check("the tangent field kills both tetrahedron symbols", vf(f1) == vf(f2) == 0)
    check("the tangent field preserves the moment ideal on the shape curve",
          vf(g1).fold() == 0 and vf(g2).fold() == 2*ii*x**2*f1)
    check("the derivative coefficient in log(x) is exactly -2*x^2*w^2",
          (vf(x)*x**-1).fold() == -2*x**2*w**2)
    check("the derivative coefficient does not vanish at the root point",
          bool((-2*x**2*w**2).at(p)))
    check("the alpha quantum shift contributes -2*x^2", (-a[2]).fold() == -2*x**2)
    check("the beta quantum shift contributes -2*x^2*w^2", (-ii*a[3]).fold() == -2*x**2*w**2)

    # Exact rational audit of the formal ODE at x=5/3. The theorem's proof
    # establishes all orders; this truncation audits the recurrence only.
    degree = 16
    x0 = F(5, 3)
    x_series = [x0, F(1)]
    square = conv(x_series, x_series, 3)
    square[0] -= 1
    coeff_c = [-2*t for t in square]
    coeff_v = conv(x_series, coeff_c, 4)
    solution = [F(0)]*(degree+2)
    for n in range(degree+1):
        residual = F(int(n == 0))
        residual -= sum((coeff_v[j]*(n-j+1)*solution[n-j+1]
                         for j in range(1, min(n, len(coeff_v)-1)+1)), F(0))
        residual -= sum((coeff_c[j]*solution[n-j]
                         for j in range(min(n, len(coeff_c)-1)+1)), F(0))
        solution[n+1] = residual/(coeff_v[0]*(n+1))
    derivative = [(j+1)*solution[j+1] for j in range(degree+1)]
    residual = [u+v for u, v in zip(conv(coeff_v, derivative, degree+1),
                                   conv(coeff_c, solution, degree+1))]
    check("the formal unit preimage recurrence is exact through degree 16",
          residual == [F(1)]+[F(0)]*degree)
    print("First four nonconstant unit-preimage coefficients:",
          ", ".join(str(c) for c in solution[1:5]))

    # For the cyclic peripheral model the obstruction to resonance is exactly
    # the integer equation c+2*n=2. Its two fixed cases are a divisibility
    # check, not a finite scan of exponents.
    check("the cyclic c=1 model has no integer resonant exponent", (2-1) % 2 != 0)
    check("the cyclic c=2 model has the unique resonant exponent zero", (2-2)//2 == 0)
    print(f"All {COUNT} exact checks passed.")
    print("Written proofs, not this truncation, establish all-orders extinction,")
    print("formal ODE surjectivity, the global Fourier functional, and all-exponent resonance.")
    print("No actual triangulation realization or PR 7.1 counterexample is certified.")


if __name__ == "__main__":
    main()