verify-empty-sector-obstruction.py

#!/usr/bin/env python3
"""Exact audit for the empty-sector Dehn-filling obstruction.

The script uses only the Python standard library.  It verifies:

1. tr(U G U G^{-1}) = 2-c^2 symbolically modulo ad-bc=1;
2. the boundary-character agreement between a nontrivial parabolic meridian
   and the meridional-filling solid-torus representation on an exact grid.

The topological and character-variety inputs are proved in the companion
research note.  The computation is a reproducibility check, not a finite
search for examples.
"""

from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction
from math import comb
from typing import Dict, Iterable, Tuple


Monomial = Tuple[int, int, int, int]  # exponents of a,b,c,d


def _reduce_ad(monomial: Monomial, coefficient: int) -> Dict[Monomial, int]:
    """Reduce (ad)^k to (bc+1)^k using ad-bc=1."""

    ea, eb, ec, ed = monomial
    k = min(ea, ed)
    reduced: Dict[Monomial, int] = {}
    for j in range(k + 1):
        key = (ea - k, eb + j, ec + j, ed - k)
        reduced[key] = reduced.get(key, 0) + coefficient * comb(k, j)
    return reduced


@dataclass(frozen=True)
class Poly:
    terms: Dict[Monomial, int]

    def __post_init__(self) -> None:
        normalized: Dict[Monomial, int] = {}
        for monomial, coefficient in self.terms.items():
            for key, value in _reduce_ad(monomial, coefficient).items():
                normalized[key] = normalized.get(key, 0) + value
        object.__setattr__(
            self, "terms", {key: value for key, value in normalized.items() if value}
        )

    def __add__(self, other: "Poly") -> "Poly":
        result = dict(self.terms)
        for key, value in other.terms.items():
            result[key] = result.get(key, 0) + value
        return Poly(result)

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

    def __sub__(self, other: "Poly") -> "Poly":
        return self + (-other)

    def __mul__(self, other: "Poly") -> "Poly":
        result: Dict[Monomial, int] = {}
        for left, left_coefficient in self.terms.items():
            for right, right_coefficient in other.terms.items():
                key = tuple(x + y for x, y in zip(left, right))
                result[key] = result.get(key, 0) + left_coefficient * right_coefficient
        return Poly(result)


ZERO = Poly({})
ONE = Poly({(0, 0, 0, 0): 1})
TWO = Poly({(0, 0, 0, 0): 2})
A = Poly({(1, 0, 0, 0): 1})
B = Poly({(0, 1, 0, 0): 1})
C = Poly({(0, 0, 1, 0): 1})
D = Poly({(0, 0, 0, 1): 1})

Matrix = Tuple[Tuple[Poly, Poly], Tuple[Poly, Poly]]


def matmul(left: Matrix, right: Matrix) -> Matrix:
    return (
        (
            left[0][0] * right[0][0] + left[0][1] * right[1][0],
            left[0][0] * right[0][1] + left[0][1] * right[1][1],
        ),
        (
            left[1][0] * right[0][0] + left[1][1] * right[1][0],
            left[1][0] * right[0][1] + left[1][1] * right[1][1],
        ),
    )


def symbolic_trace_check() -> None:
    u: Matrix = ((ONE, ONE), (ZERO, ONE))
    g: Matrix = ((A, B), (C, D))
    g_inverse: Matrix = ((D, -B), (-C, A))
    product = matmul(matmul(matmul(u, g), u), g_inverse)
    trace = product[0][0] + product[1][1]
    assert trace == TWO - C * C


NumericMatrix = Tuple[
    Tuple[Fraction, Fraction],
    Tuple[Fraction, Fraction],
]


def numeric_matmul(left: NumericMatrix, right: NumericMatrix) -> NumericMatrix:
    return (
        (
            left[0][0] * right[0][0] + left[0][1] * right[1][0],
            left[0][0] * right[0][1] + left[0][1] * right[1][1],
        ),
        (
            left[1][0] * right[0][0] + left[1][1] * right[1][0],
            left[1][0] * right[0][1] + left[1][1] * right[1][1],
        ),
    )


IDENTITY: NumericMatrix = (
    (Fraction(1), Fraction(0)),
    (Fraction(0), Fraction(1)),
)


def numeric_inverse(matrix: NumericMatrix) -> NumericMatrix:
    (a, b), (c, d) = matrix
    determinant = a * d - b * c
    assert determinant != 0
    return ((d / determinant, -b / determinant), (-c / determinant, a / determinant))


def numeric_power(matrix: NumericMatrix, exponent: int) -> NumericMatrix:
    if exponent < 0:
        return numeric_power(numeric_inverse(matrix), -exponent)
    result = IDENTITY
    factor = matrix
    power = exponent
    while power:
        if power & 1:
            result = numeric_matmul(result, factor)
        factor = numeric_matmul(factor, factor)
        power //= 2
    return result


def trace(matrix: NumericMatrix) -> Fraction:
    return matrix[0][0] + matrix[1][1]


def boundary_character_checks() -> int:
    checked = 0
    parameters: Iterable[Tuple[int, Fraction]] = (
        (1, Fraction(0)),
        (1, Fraction(3, 5)),
        (-1, Fraction(0)),
        (-1, Fraction(-7, 4)),
    )
    u: NumericMatrix = (
        (Fraction(1), Fraction(1)),
        (Fraction(0), Fraction(1)),
    )
    for epsilon, upper_right in parameters:
        longitude: NumericMatrix = (
            (Fraction(epsilon), upper_right),
            (Fraction(0), Fraction(epsilon)),
        )
        filled_meridian = IDENTITY
        filled_longitude: NumericMatrix = (
            (Fraction(epsilon), Fraction(0)),
            (Fraction(0), Fraction(epsilon)),
        )
        assert numeric_matmul(u, longitude) == numeric_matmul(longitude, u)
        for m in range(-4, 5):
            for n in range(-4, 5):
                cusp_word = numeric_matmul(numeric_power(u, m), numeric_power(longitude, n))
                filled_word = numeric_matmul(
                    numeric_power(filled_meridian, m),
                    numeric_power(filled_longitude, n),
                )
                assert trace(cusp_word) == trace(filled_word) == 2 * (epsilon ** n)
                checked += 1
    return checked


def main() -> None:
    symbolic_trace_check()
    checks = boundary_character_checks()
    print("symbolic quotient identity: tr(U G U G^-1) = 2 - c^2")
    print("determinant relation used exactly: ad - bc = 1")
    print(f"boundary-character identities checked exactly: {checks}")
    print("all exact checks passed")


if __name__ == "__main__":
    main()