Listing: Listing

from fractions import Fraction
from itertools import permutations, product
from math import prod

ETA = (-1, 1, 1, 1)
PERMS = tuple(permutations(range(4)))

KUP = (
    ( 1, 0, 0,  1),
    ( 1, 0, 0, -1),
    (-1,-1, 0,  0),
    (-1, 1, 0,  0),
)

A = (
    (0,0, 0,0),
    (0,1, 0,0),
    (0,0,-1,0),
    (0,0, 0,0),
)
B = (
    (0,0, 0,0),
    (0,0, 0,0),
    (0,0,-1,0),
    (0,0, 0,1),
)

def lower(v):
    return tuple(ETA[a] * v[a] for a in range(4))

def dot(v, w):
    return sum(ETA[a] * v[a] * w[a] for a in range(4))

def gauge(k, xi):
    kl, xl = lower(k), lower(xi)
    return tuple(
        tuple(kl[a]*xl[b] + kl[b]*xl[a]
              for b in range(4))
        for a in range(4)
    )

def sparse_vector(v):
    return tuple(((a,), v[a])
                 for a in range(4) if v[a] != 0)

def sparse_tensor(H):
    return tuple(
        ((a,b), H[a][b])
        for a in range(4) for b in range(4)
        if H[a][b] != 0
    )

def matchings(slots):
    slots = tuple(slots)
    if not slots:
        yield ()
        return
    a = slots[0]
    for j in range(1, len(slots)):
        b = slots[j]
        rest = slots[1:j] + slots[j+1:]
        for tail in matchings(rest):
            yield ((a,b),) + tail

P_LIST = tuple(matchings(range(10)))
I_LIST = tuple(product((1,2), repeat=4))
LABELS = tuple((I,P) for I in I_LIST for P in P_LIST)

assert len(P_LIST) == len(set(P_LIST)) == 945
assert len(LABELS) == 15120

K = tuple(sparse_vector(lower(k)) for k in KUP)
GT = sparse_tensor(gauge(KUP[0], (1,0,0,0)))
GX = sparse_tensor(gauge(KUP[0], (0,1,0,0)))
AS = sparse_tensor(A)
BS = sparse_tensor(B)

H_TIME = (GT, AS, BS, BS)
H_BLIND = (GX, AS, BS, BS)
J_PLUS = (1,1,2,2)
J_MINUS = (2,2,1,1)

assert all(any(k) and dot(k,k) == 0 for k in KUP)
assert all(sum(k[a] for k in KUP) == 0
           for a in range(4))

def contract(P, factors):
    total = 0
    for chosen in product(*factors):
        indices = tuple(
            a
            for index_tuple, value in chosen
            for a in index_tuple
        )
        assert len(indices) == 10
        if any(indices[u] != indices[v] for u,v in P):
            continue
        value = prod(value for index_tuple, value in chosen)
        value *= prod(ETA[indices[u]] for u,v in P)
        total += value
    return total

def raw(I, P, J, H):
    if sorted(I) != sorted(J):
        return 0
    total = 0
    for sigma in PERMS:
        if any(I[r] != J[sigma[r]] for r in range(4)):
            continue
        factors = (
            K[sigma[0]], K[sigma[1]],
            H[sigma[0]], H[sigma[1]],
            H[sigma[2]], H[sigma[3]],
        )
        total -= contract(P, factors)
    return total

P0 = ((0,1),(2,3),(4,6),(5,8),(7,9))
LA = ((1,1,2,2), P0)
LB = ((2,2,1,1), P0)

def v_abb(p0, q0, r0):
    assert p0 + q0 + r0 == 0
    return -2*(p0*q0 + p0*r0 + q0*r0)

# Uses the proved restricted Einstein cubic.
W = (-v_abb(2,-1,-1),
      v_abb(1,0,-1),
      v_abb(1,0,-1))
assert W == (-6,2,2)

target_value = Fraction(-2*W[0] + W[1] + W[2], 3)
assert target_value == Fraction(16,3)
target = (Fraction(0), target_value, target_value)

def reconstruct():
    plus = {}
    blind = {}
    for label in LABELS:
        I, P = label
        plus[label] = raw(I, P, J_PLUS, H_TIME)
        blind[label] = raw(I, P, J_PLUS, H_BLIND)
        assert blind[label] == 0
        assert abs(plus[label]) <= 768

    columns = {}
    for I, P in LABELS:
        Ibar = tuple(3-i for i in I)
        minus = raw(I, P, J_MINUS, H_TIME)
        assert minus == plus[(Ibar,P)]
        assert abs(minus) <= 768
        columns[(I,P)] = (0, plus[(I,P)], minus)

    ca, cb = columns[LA], columns[LB]
    assert ca == (0,8,0)
    assert cb == (0,0,8)

    # An actual minor of original columns.
    minor = ca[1]*cb[2] - cb[1]*ca[2]
    assert minor == 64

    fit = tuple(
        Fraction(2,3)*ca[i] + Fraction(2,3)*cb[i]
        for i in range(3)
    )
    assert fit == target

    # delta(0),delta(1),delta(2),delta(3) = 0,+,-,0.
    row_map = (0,1,2,0)
    target4 = tuple(target[i] for i in row_map)
    ca4 = tuple(ca[i] for i in row_map)
    cb4 = tuple(cb[i] for i in row_map)
    assert target4 == (
        Fraction(0), Fraction(16,3),
        Fraction(16,3), Fraction(0)
    )
    assert tuple(
        Fraction(2,3)*ca4[i] + Fraction(2,3)*cb4[i]
        for i in range(4)
    ) == target4

    # Analytic upper bound: every first entry is zero.
    # Analytic lower bound: original minor 64.
    # Augmented rank: the explicit fit above.
    return columns, minor

if __name__ == "__main__":
    columns, minor = reconstruct()
    print("shape =", (3, len(columns)))
    print("target =", target)
    print("original minor =", minor)
    print("status = SEP-finite; declared Grade S model")