verify-global-meridian-index.py

#!/usr/bin/env python3
"""Exact symbolic audit for the global meridian-index research note.

Python 3.10+, standard library only. Polynomial identities are universal,
not random tests. A short seed q-series check is a normalization sanity check,
not a proof of the sourced all-orders index recurrence or degree formula.
No manifolds, slopes, triangulations, sign assignments, or skein inputs are
enumerated. The one-block cone example is an abstract lattice model only.
"""

from fractions import Fraction as Q
from collections import defaultdict


class P:
    """Sparse commutative Laurent polynomials over Q, named indeterminates."""

    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, P) else P({(): Q(x)})

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

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

    __radd__ = __add__

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

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

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

    def __mul__(self, other):
        out = defaultdict(Q)
        for a, c in self.terms.items():
            for b, d in P.of(other).terms.items():
                exponents = defaultdict(int, a)
                for name, exponent in b:
                    exponents[name] += exponent
                m = tuple(sorted((n, e) for n, e in exponents.items() if e))
                out[m] += c*d
        return P(out)

    __rmul__ = __mul__

    def __pow__(self, n):
        if n < 0:
            assert len(self.terms) == 1
            (m, c), = self.terms.items()
            return P({tuple((name, exponent*n) for name, exponent in m): c**n})
        out = P.of(1)
        for _ in range(n):
            out = out*self
        return out

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

    def __repr__(self):
        return repr(self.terms)


def variables(names):
    return tuple(P.var(name) for name in names.split())


B2 = ((0, -1), (1, 0))
B3 = ((0, 1, -1), (-1, 0, 1), (1, -1, 0))
DS = ((0, 1, 1), (1, 0, 1), (1, 1, 0))
B6 = tuple(tuple(B3[i % 3][j % 3] for j in range(6)) for i in range(6))


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


def pair(a, matrix, b):
    return sum(a[i]*matrix[i][j]*b[j]
               for i in range(len(a)) for j in range(len(b)))


def add(a, 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 T(u):
    return (-u[0], P.of(0), -u[1])


def delta(x):
    return x[0]*x[1]+x[1]*x[2]+x[2]*x[0]


def collapse(u):
    r, s, t = (u[j]+u[j+3] for j in range(3))
    return (r-s, t-s)


checks = []


def check(name, assertion):
    if not assertion:
        raise AssertionError(name)
    checks.append(name)
    print(f"PASS {len(checks):02d}: {name}")


check("Omega3 is skew", all(B3[i][j] == -B3[j][i]
                            for i in range(3) for j in range(3)))
check("both quadratic/skew positivity matrices are entrywise nonnegative",
      all(DS[i][j] + sign*B3[i][j] >= 0
          for i in range(3) for j in range(3) for sign in (-1, 1)))
check("tetrahedron translation lies in the skew radical",
      all(sum(row) == 0 for row in B3))

a, b, c, d = variables("a b c d")
check("T pulls the triple skew form back to the native root form",
      pair(T((a, b)), B3, T((c, d))) == pair((a, b), B2, (c, d)))
u6 = variables("u0 u1 u2 u3 u4 u5")
v6 = variables("v0 v1 v2 v3 v4 v5")
check("full six-root central elimination preserves every Weyl pairing",
      pair(u6, B6, v6) == pair(collapse(u6), B2, collapse(v6)))

r, s, t, gamma, alpha = variables("r s t gamma alpha")
beta = 1-gamma-alpha
check("exact central scalar exponent yields the affine filling angle equation",
      gamma*(r-s)+alpha*(t-s)+s == gamma*r+beta*s+alpha*t)

check("right multiplication by z^-2 has Weyl exponent -u_w",
      Q(1, 2)*pair((a, b), B2, (-2, 0)) == -b)
check("right multiplication by w^2 has Weyl exponent -u_z",
      Q(1, 2)*pair((a, b), B2, (0, 2)) == -a)
m, e = variables("m e")
check("q=A^4 translates the two seed recurrence factors without a sign swap",
      -2*m == 4*(-m)*Q(1, 2) and 2*e == 4*e*Q(1, 2))
check("right tetrahedron shifts preserve each root parity coset",
      (-2) % 2 == 0 and 2 % 2 == 0)

# Recurrence (3.1a) telescoping check for the source summand T_n.
# X=q^n, Y=q^m, Z=q^e, Qv=q. Negative factorial boundary terms are
# handled in the cited source, not inferred from this rational identity.
X, Y, Z, Qv = variables("X Y Z q")
den = Y*(1-Qv*X)*(1-Qv*X*Z)
lhs = (X**-1-1)*den + (1-Qv*X)
rhs = (X**-1-1)*den + (Qv**-1*X**-1-1)*Qv*X
check("tetrahedron-index first recurrence telescoper: cleared numerator identity",
      lhs == rhs)
check("tetrahedron-index second recurrence summand identity",
      X*Z+(1-X*Z) == 1)

u, v, minimum = variables("input_z input_w minimum")
theta = (gamma, beta, alpha)
vv = (a, b)
uu = (u, v)
xx = add(T(vv), scale(-minimum, (1, 1, 1)))
vartheta = lambda z: gamma*z[0]+alpha*z[1]
lam = add(scale(2, vv), scale(-1, uu))
raw_degree = -vartheta(lam)+Q(1, 2)*pair(lam, B2, uu)+2*delta(xx)-2*minimum
reduced_degree = 2*delta(xx)+2*dot(theta, xx)+vartheta(uu)+pair(xx, B3, T(uu))
check("all-input summand valuation identity including the Weyl phase",
      raw_degree == reduced_degree)

bvec = variables("b1 b2 b3")
fvec = variables("f1 f2 f3")
xvec = add(bvec, fvec)
degree = lambda x: 2*delta(x)+2*dot(theta, x)+pair(x, B3, T(uu))
rhs_growth = 2*delta(fvec)+2*dot(theta, fvec) + 2*pair(bvec, DS, fvec)-2*pair(bvec, B3, fvec)
residual = degree(xvec)-degree(bvec)-rhs_growth
check("growth identity residual is exactly the term killed by isotropy",
      residual == pair(fvec, B3, T(uu))+2*pair(bvec, B3, fvec))

mu = variables("mu_z mu_w")
lamb = variables("lambda_z lambda_w")
b_affine = T(scale(Q(1, 2), add(uu, lamb)))
f_homogeneous = T(scale(Q(1, 2), mu))
check("affine phase cancellation has precisely the homogeneous isotropy residual",
      pair(f_homogeneous, B3, T(uu))+2*pair(b_affine, B3, f_homogeneous)
      == Q(1, 2)*pair(lamb, B2, mu))

n = P.var("n")
check("null-ray degree scales quadratically plus linearly",
      2*delta(scale(n, fvec))+2*dot(theta, scale(n, fvec))
      == 2*n*n*delta(fvec)+2*n*dot(theta, fvec))
D, E, norm = variables("delta_value angle_value norm")
check("compactness coercivity estimate has a nonnegative exact remainder",
      (norm*norm*D+norm*E)-norm*(D+E) == norm*(norm-1)*D)

shift = variables("h_z h_w")
orbit = variables("orbit_z orbit_w")
check("moment reindex phase differs only by the isotropic orbit/moment pairing",
      Q(1, 2)*pair(orbit, B2, shift)+Q(1, 2)*pair(add(orbit, shift), B2, uu)
      -Q(1, 2)*pair(add(orbit, shift), B2, uu)
      == Q(1, 2)*pair(orbit, B2, shift))
check("a same-longitude square has no hidden Weyl scalar",
      pair(shift, B2, shift) == 0)

cone_x, cone_z = variables("cone_x cone_z")
weak_vector = (cone_x, cone_x, cone_z)
check("abstract weak-only cone has the stated quadratic degree",
      delta(weak_vector) == cone_x**2+2*cone_x*cone_z)
check("abstract weak-only moment forces alpha=1, not a strict angle",
      vartheta((0, 2)) == 2*alpha)


def index_series(magnetic, electric, max_half_degree):
    """Exact I_Delta coefficients in t=q^(1/2), up to the given degree.

    The summand base degree is quadratic in n. The stopping condition
    occurs only past its turning point, so no omitted summand contributes.
    """
    out = defaultdict(int)
    n = max(0, -electric)
    while True:
        base = n*(n+1)-(2*n+electric)*magnetic
        if n >= max(0, magnetic) and base > max_half_degree:
            break
        bound = (max_half_degree-base)//2
        if bound >= 0:
            inv = [0]*(bound+1)
            inv[0] = 1
            for length in (n, n+electric):
                for step in range(1, length+1):
                    for degree in range(step, bound+1):
                        inv[degree] += inv[degree-step]
            for degree, coefficient in enumerate(inv):
                out[base+2*degree] += (-1)**n*coefficient
        n += 1
    return {degree: coefficient for degree, coefficient in out.items() if coefficient}


precision = 16
i00 = index_series(0, 0, precision)
i10 = index_series(1, 0, precision)
i01 = index_series(0, 1, precision)
check("seed normalization has constant coefficient 1 and no negative powers",
      min(i00) == 0 and i00[0] == 1)
check("exact seed recurrence at (m,e)=(0,0) through q^8",
      all(i10.get(j, 0)+i01.get(j, 0)-i00.get(j, 0) == 0
          for j in set(i00) | set(i10) | set(i01)))
check("degree convention agrees with the nontrivial seed I(1,0)", min(i10) == 2)
print("I_Delta(0,0), indexed by q^(1/2) exponent:", sorted(i00.items()))
print(f"\n{len(checks)} exact checks passed.")
print("Convergence, all-relations descent, and nonzero unit are proved in the note.")
print("No assertion about the actual longitude difference or unrestricted fillings.")