verify-ck-pbw.py

#!/usr/bin/env python3
"""Independent exact CK audit in ORDERED PBW coordinates.

Python 3 standard library only. No float, CAS, source import, triangulation
enumeration, or random sampling. Coefficients: Z[i,t^+-1], i^2=-1,
t=A^(1/2), v0=i*t^2. This is the universal CK coefficient ring.

Unlike the ATTACK verifier, multiplication reorders ordinary monomials;
Weyl ordering is a conversion at the interface, not the multiplication law.
The fixed group of twelve oriented tetrahedral symmetries is checked in
full as a finite algebraic certificate, not a geometric configuration search.
"""
from collections import defaultdict
from itertools import permutations

LABELS = ('z', 'p', 'd', 'y', 'P', 'D')
ROWS = (('z', 'p', 'D'), ('p', 'd', 'y'),
        ('d', 'z', 'P'), ('y', 'P', 'D'))
NAMES = tuple((name, face) for face, row in enumerate(ROWS) for name in row)
N = len(NAMES)
B = [[0] * N for _ in range(N)]
for face in range(4):
    a, b, c = range(3 * face, 3 * face + 3)
    for j, k in ((a, b), (b, c), (c, a)):
        B[j][k], B[k][j] = 1, -1


def reordering(a, b):
    """t exponent when x^a x^b is sorted in increasing generator order."""
    return 2 * sum(a[j] * b[k] * B[j][k]
                   for j in range(N) for k in range(j))


def omega(a, b):
    return sum(a[j] * B[j][k] * b[k]
               for j in range(N) for k in range(N))


class Poly(dict):
    # key = (ordered exponent vector, t power, i parity); value = integer.
    def __add__(self, other):
        out = defaultdict(int, self)
        for key, value in other.items():
            out[key] += value
        return Poly({key: value for key, value in out.items() if value})

    def __neg__(self):
        return Poly({key: -value for key, value in self.items()})

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

    def __mul__(self, other):
        out = Poly()
        for (a, ta, ia), ca in self.items():
            for (b, tb, ib), cb in other.items():
                out += term(tuple(x + y for x, y in zip(a, b)),
                            ta + tb + reordering(a, b), ia + ib, ca * cb)
        return out

    def inverse(self):
        assert len(self) == 1, 'Only Laurent monomial units are inverted.'
        (a, t, i), c = next(iter(self.items()))
        assert abs(c) == 1
        b = tuple(-x for x in a)
        result = term(b, -t - reordering(a, b), -i, c)
        assert self * result == ONE == result * self
        return result

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


def term(a=None, t=0, i=0, c=1):
    a = (0,) * N if a is None else tuple(a)
    c *= (-1) ** ((i // 2) % 2)
    return Poly({(a, t, i % 2): c}) if c else Poly()


ONE = term()
HALF = {key: term(tuple(int(j == n) for j in range(N)))
        for n, key in enumerate(NAMES)}


def Weyl(a):
    # [x^a] = t^(-sum_{j<k} B[j,k] a[j] a[k]) x_0^a0 ... x_11^a11.
    return term(a, -sum(B[j][k] * a[j] * a[k]
                        for j in range(N) for k in range(j + 1, N)))


def unpack_weyl(m):
    assert len(m) == 1
    (a, t, i), c = next(iter(m.items()))
    correction = sum(B[j][k] * a[j] * a[k]
                     for j in range(N) for k in range(j + 1, N))
    return a, t + correction, i, c


def bracket(*monomials):
    data = [unpack_weyl(m) for m in monomials]
    a = tuple(sum(row[0][j] for row in data) for j in range(N))
    coefficient = term(t=sum(row[1] for row in data),
                       i=sum(row[2] for row in data))
    for row in data:
        coefficient = coefficient * term(c=row[3])
    return coefficient * Weyl(a)


def vector(m):
    return next(iter(m))[0]


SHAPE = {name: bracket(*(HALF[key] for key in NAMES if key[0] == name))
         for name in LABELS}
z, p, d, y, P, D = (SHAPE[name] for name in LABELS)
V0 = term(t=2, i=1)
C = V0.inverse()
TINV = term(t=-1)


def gamma(a, b):
    face, = [f for f, row in enumerate(ROWS) if a in row and b in row]
    return C * bracket(HALF[a, face], HALF[b, face])


def local_pair(labels):
    z0, p0, d0, y0, P0, D0 = labels
    vv = bracket(SHAPE[z0], SHAPE[p0], SHAPE[d0])
    vertex = gamma(z0, p0) * gamma(p0, d0) * gamma(d0, z0) - C ** 2
    raw = (bracket(gamma(z0, P0).inverse(), gamma(P0, d0))
           - bracket(gamma(d0, y0).inverse(), gamma(y0, p0)) * gamma(z0, p0)
           - gamma(p0, d0).inverse() * bracket(gamma(p0, D0).inverse(), gamma(D0, z0)))
    face, = [f for f, row in enumerate(ROWS) if z0 in row and d0 in row]
    unit = HALF[d0, face].inverse() * HALF[z0, face]
    lag = SHAPE[z0] ** 2 + SHAPE[d0] ** -2 - ONE
    correction = (vv - V0) * (vv.inverse() * SHAPE[z0] ** 2 - V0.inverse() * SHAPE[d0] ** -2)
    assert vertex == C ** 3 * (vv - V0)
    assert raw * unit == TINV * (-lag + correction)
    # Explicit reverse generator witnesses: lag in (vertex,raw), raw in (vv-V0,lag).
    assert lag == -(term(t=1) * raw * unit) + correction
    assert raw == TINV * (-lag + correction) * unit.inverse()
    return vv, raw, unit, lag


# Sanity-check the ordered engine against its defining relations and a
# Weyl product identity on the full ambient lattice basis and their negatives.
for j in range(N):
    for k in range(N):
        x, yy = HALF[NAMES[j]], HALF[NAMES[k]]
        assert x * yy == term(t=2 * B[j][k]) * yy * x
        for sa in (-1, 1):
            for sb in (-1, 1):
                a = tuple(sa * int(n == j) for n in range(N))
                b = tuple(sb * int(n == k) for n in range(N))
                assert Weyl(a) * Weyl(b) == term(t=omega(a, b)) * Weyl(tuple(x + y for x, y in zip(a, b)))
print('PASS 01: ordered PBW engine, all 144 defining commutators, signed Weyl conversion.')

v, raw, unit, f = local_pair(LABELS)
assert raw * unit == TINV - C * d.inverse() * p * z - C.inverse() * d.inverse() * p.inverse() * z
assert d.inverse() * p * z == TINV * v * d ** -2
assert d.inverse() * p.inverse() * z == TINV * v.inverse() * z ** 2
print('PASS 02: base vertex, raw trinomial, ordered right unit, exact correction, both ideal inclusions.')

EDGE_ENDS = {'z': (0, 1), 'p': (0, 2), 'd': (0, 3),
             'y': (2, 3), 'P': (1, 3), 'D': (1, 2)}
FROM_ENDS = {frozenset(value): key for key, value in EDGE_ENDS.items()}
symmetries = []
for perm in permutations(range(4)):
    if sum(perm[j] > perm[k] for j in range(4) for k in range(j + 1, 4)) % 2:
        continue
    mapping = {key: FROM_ENDS[frozenset(perm[j] for j in ends)]
               for key, ends in EDGE_ENDS.items()}
    index_map = []
    for name, face in NAMES:
        target_face, = [g for g, row in enumerate(ROWS)
                        if set(row) == {mapping[x] for x in ROWS[face]}]
        index_map.append(NAMES.index((mapping[name], target_face)))
    assert all(B[j][k] == B[index_map[j]][index_map[k]]
               for j in range(N) for k in range(N))
    local_pair(tuple(mapping[name] for name in LABELS))
    symmetries.append(mapping)
assert len(symmetries) == 12
print('PASS 03: all twelve oriented tetrahedral symmetries preserve the ambient form and both identities.')

fprime = p ** 2 + z ** -2 - ONE
assert f * z ** -2 + fprime == (v ** 2 - V0 ** 2) * bracket(z ** -2, d ** -2)
for labels in (('z', 'p', 'd'), ('p', 'd', 'z'), ('d', 'z', 'p')):
    aa, bb, cc = (SHAPE[name] for name in labels)
    ff = aa ** 2 + cc ** -2 - ONE
    gg = bb ** 2 + aa ** -2 - ONE
    assert ff * aa ** -2 + gg == (v ** 2 - V0 ** 2) * bracket(aa ** -2, cc ** -2)
print('PASS 04: cyclic Lagrangians have explicit right-ideal witnesses over the universal coefficient ring.')

base = ((0, 1, -1), (-1, 0, 1), (1, -1, 0))
form = [[omega(vector(SHAPE[a]), vector(SHAPE[b])) for b in LABELS] for a in LABELS]
assert form == [[base[j % 3][k % 3] for k in range(6)] for j in range(6)]
r, rr, rrr = y * z.inverse(), P * p.inverse(), D * d.inverse()
for ratio in (r, rr, rrr):
    assert all(ratio * shape == shape * ratio for shape in SHAPE.values())
vertices = [v, bracket(z, P, D), bracket(y, p, D), bracket(y, P, d)]
assert vertices[1] == v * rr * rrr
assert vertices[2] == v * r * rrr
assert vertices[3] == v * r * rr
q1, q2, q3 = rr * rrr - ONE, r * rrr - ONE, r * rr - ONE
assert r - rr == (q2 - q1) * rrr.inverse()
assert r - rrr == (q3 - q1) * rr.inverse()
assert r ** 2 - ONE == r * (q2 - q1) * rrr.inverse() + q3
print('PASS 05: shape skew form, central opposite ratios, division-free equal-ratio and square-one witnesses.')


def degree(a):
    pairs = [[j for j, key in enumerate(NAMES) if key[0] == name] for name in LABELS]
    return tuple(a[j] - a[k] for j, k in pairs)


def project(poly):
    return Poly({key: value for key, value in poly.items() if degree(key[0]) == (0,) * 6})


assert project(raw) == Poly() and project(raw * unit) == raw * unit
assert len({degree(a) for a, _, _ in raw}) == 1
for shape in SHAPE.values():
    for half in HALF.values():
        test = ONE + half + shape + half.inverse()
        assert project(shape * test) == shape * project(test)
        assert project(test * shape) == project(test) * shape
assert project(unit) == Poly() and project(unit.inverse()) == Poly()
assert project(unit * unit.inverse()) == ONE  # Deliberate non-algebra-map test.
print('PASS 06: raw degree, normalized degree zero, bimodule projection tests, explicit failure of multiplicativity.')

# The edge argument is universal in its incidence vectors. On each paired
# face, the three bigon vectors have opposite skew contributions, so pair
# cancellation works for arbitrary repeated visits to that face suspension.
for j in range(3):
    for k in range(3):
        assert base[j][k] + (-base[j][k]) == 0
e = bracket(z, p, y)  # One noncentral illustrative edge monomial, not a manifold.
assert e.inverse() + term(t=-4) == term(t=-4) * e.inverse() * (e + term(t=4))
print('PASS 07: universal paired-face commutator cancellation; exact redundant inverse-edge left-ideal witness.')

# Changed-sidedness countercheck: v is central in S but not in U.
half = gamma('z', 'D')  # In rho of the triangle algebra, not just an arbitrary half variable.
assert half * v != v * half
pairing = omega(vector(half), vector(v))
q = term(t=2 * pairing)
assert half * (v - V0) * half.inverse() - q * (v - V0) == (q - ONE) * V0
print('PASS 08: changing to two-sided closure creates a nonzero scalar generically; ambient centrality is forbidden.')

print('All independent ordered-PBW checks passed.')
print('SCOPE: exact finite algebraic certificates; the general lattice-intersection proof is in the paper.')
print('NOT VERIFIED: source trace theorems, external priority, or an original six-root peripheral named-action theorem.')