verify-six-point-obstruction.py

#!/usr/bin/env python3
"""Exact six-point audit of the proposed braid-arrangement filtration bridge.

Requires SymPy.  Every calculation is over Q or Q[x,y,z].
"""

from itertools import permutations

import sympy as sp


def inversion_sign(word):
    inversions = sum(
        word[i] > word[j]
        for i in range(len(word))
        for j in range(i + 1, len(word))
    )
    return -1 if inversions % 2 else 1


def main():
    # Five-point positive control.  With (s12,s13,s23)=(2,3,-5), the hard
    # invariants are (s1P,s2P,s3P)=(-5,3,2).  The natural diagonal turns the
    # degree-one row into the unique sum-zero chamber function.
    G5 = sp.Matrix([[2, -3]])
    R5 = sp.Matrix([[4, 9]])
    D5 = sp.diag(sp.Rational(1, 2), sp.Rational(1, 3))
    assert G5 * D5 == sp.Matrix([[1, -1]])
    assert G5.col_join(R5).rank() == 2

    # A valid generic five-hard-particle specialization.  The six pairwise
    # invariants sum to zero, so P=-(p1+p2+p3+p4) is massless.
    s12, s13, s14, s23, s24 = map(sp.Rational, (2, 3, 5, 7, 11))
    s34 = -(s12 + s13 + s14 + s23 + s24)
    sij = {
        (1, 2): s12,
        (1, 3): s13,
        (1, 4): s14,
        (2, 3): s23,
        (2, 4): s24,
        (3, 4): s34,
    }

    def s_pair(a, b):
        return sij[tuple(sorted((a, b)))]

    p = {a: -sum(s_pair(a, b) for b in range(1, 5) if b != a) for a in range(1, 5)}

    def u(a, b, c):
        return 2 * s_pair(b, c) + 2 * p[c] + p[b]

    def v(a, b, c):
        return u(a, b, c) + 2 * s_pair(a, c)

    # Dong--Stieberger (5.11) and the two degree-one rows in (5.12), in
    # column order 234,243,324,342,423,432.  The irrelevant factor 1/2 in
    # R_6 is retained.
    G = sp.Matrix(
        [
            [p[4], -p[3], 0, 0, -p[3], p[2]],
            [0, 0, p[4], -p[2], p[3], -p[2]],
        ]
    )
    R = sp.Rational(1, 2) * sp.Matrix(
        [
            [p[4] ** 2, -p[3] * u(2, 3, 4), 0, 0, -p[3] * v(2, 3, 4), -p[2] ** 2],
            [0, 0, p[4] ** 2, -p[2] * u(3, 2, 4), -p[3] ** 2, -p[2] * v(3, 2, 4)],
            [-p[4] * u(2, 4, 3), p[3] ** 2, -p[4] * v(2, 4, 3), -p[2] ** 2, 0, 0],
        ]
    )
    F2 = G.col_join(R)

    chambers = list(permutations((2, 3, 4)))
    pairs = ((2, 3), (2, 4), (3, 4))
    one = sp.ones(6, 1)
    heavisides = [
        sp.Matrix([int(word.index(a) < word.index(b)) for word in chambers])
        for a, b in pairs
    ]
    Q = sp.Matrix.hstack(one, *heavisides)

    assert s34 == -28
    assert [p[i] for i in range(1, 5)] == [-10, -20, 18, 12]
    assert F2.rank() == 5
    kernel = F2.nullspace()
    assert len(kernel) == 1
    d = kernel[0]
    assert list(d) == [
        sp.Rational(10, 13),
        sp.Rational(20, 39),
        sp.Rational(-10, 7),
        sp.Rational(6, 7),
        sp.Rational(-10, 9),
        sp.Rational(1),
    ]
    diagonal_defect = G * sp.diag(*d) * Q
    assert list(diagonal_defect[:, 0]) == [0, 0]
    assert list(diagonal_defect[:, 1]) == [20, -20]
    assert diagonal_defect != sp.zeros(2, 4)

    # The type-A2 Varchenko matrix: entry (sigma,tau) is the product of the
    # weights of the hyperplanes separating the two chambers.
    x, y, z = sp.symbols("x y z")
    weights = dict(zip(pairs, (x, y, z)))

    def side(word, pair):
        a, b = pair
        return word.index(a) < word.index(b)

    V = sp.Matrix(
        6,
        6,
        lambda i, j: sp.prod(
            weights[pair]
            for pair in pairs
            if side(chambers[i], pair) != side(chambers[j], pair)
        ),
    )

    # Already the first EYM row would have to annihilate 1,h23,h24,h34.
    first_row_equations = [sp.expand(e) for e in list(G[0, :] * V * Q)]
    assert list(sp.groebner(first_row_equations, x, y, z, order="lex")) == [1]

    # A transparent elimination certificate for the same contradiction.
    # Treat xy,x,y,1 as four linear unknowns for fixed z.
    M = sp.Matrix(
        [
            [-22 * z - 13, -18 * z - 13, -13 * z - 18, -13 * z - 22],
            [-10 * z - 10, -10, -3 * z - 18, -3 * z - 12],
            [-19 * z - 10, 6 - 9 * z, -19 * z - 9, -3 * z - 3],
            [6 - 28 * z, 6 - 18 * z, -19 * z, 6 - 19 * z],
        ]
    )
    det_expected = -64 * (z - 1) * (z + 1) * (189 * z**2 - 608 * z + 189)
    assert sp.factor(M.det() - det_expected) == 0
    adj_column = M.adjugate()[:, 3]
    segre_defect = sp.factor(adj_column[0] * adj_column[3] - adj_column[1] * adj_column[2])
    quartic = 7551 * z**4 + 29650 * z**3 - 70178 * z**2 + 29650 * z + 7551
    assert sp.factor(segre_defect + 64 * (z + 1) ** 2 * quartic) == 0
    assert sp.gcd(quartic, (z - 1) * (189 * z**2 - 608 * z + 189)) == 1
    special_solution = sp.linsolve(
        (M.subs(z, -1)[:, :3], -M.subs(z, -1)[:, 3]),
        sp.symbols("A B C"),
    )
    assert special_solution == sp.FiniteSet((sp.Rational(-1), sp.Rational(9, 5), sp.Rational(-9, 5)))

    # Check the only canonical one-dimensional character on either side.
    sign = sp.diag(*(inversion_sign(word) for word in chambers))
    identity = sp.eye(6)
    twist_certificates = {}
    for left_name, left in (("plain", identity), ("sign", sign)):
        for right_name, right in (("plain", identity), ("sign", sign)):
            equations = [sp.expand(e) for e in list(G * left * V * right * Q)]
            basis = list(sp.groebner(equations, x, y, z, order="lex"))
            twist_certificates[(left_name, right_name)] = basis
            assert basis == [1]

    determinant_factors = [
        p[1], p[2], p[3], p[4], s12 + s34, s13 + s24, s14 + s23
    ]
    assert all(value != 0 for value in determinant_factors)

    print("five-point positive control G5*D5 =", (G5 * D5).tolist())
    print("hard-particle invariants:", [s12, s13, s14, s23, s24, s34])
    print("s_iP:", [p[i] for i in range(1, 5)])
    print("G6 =", G.tolist())
    print("R6 =", R.tolist())
    print("rank([G6;R6]) =", F2.rank())
    print("right kernel =", list(d))
    print("diagonal Heaviside defect G6*diag(d)*[1,h23,h24,h34] =", diagonal_defect.tolist())
    print("Varchenko first-row ideal Groebner basis =", list(sp.groebner(first_row_equations, x, y, z, order="lex")))
    print("linearized determinant =", sp.factor(M.det()))
    print("Segre defect =", segre_defect)
    print("z=-1 linear solution (xy,x,y) =", special_solution)
    print("sign-twist Groebner certificates =", twist_certificates)
    print("six-point reconstruction factors =", determinant_factors)
    print("all exact six-point obstruction checks passed")


if __name__ == "__main__":
    main()