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_DENSE = (
    (0,0, 0,0),
    (0,1, 0,0),
    (0,0,-1,0),
    (0,0, 0,0),
)
B_DENSE = (
    (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 addvec(v,w):
    return tuple(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 trace_eta(H):
    return sum(ETA[a]*H[a][a] for a in range(4))

def transverse(v,H):
    return all(
        sum(v[a]*H[a][b] for a in range(4)) == 0
        for b in range(4)
    )

t0 = (1,0,0,0)
xx = (0,1,0,0)

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))

invariants = tuple(
    dot(addvec(KUP[0],KUP[j]),
        addvec(KUP[0],KUP[j]))
    for j in (1,2,3)
)
assert invariants == (-4,2,2)
assert sum(invariants) == 0
assert dot(t0,t0) == -1

for k,H in zip(KUP[1:],(A_DENSE,B_DENSE,B_DENSE)):
    assert all(H[a][b] == H[b][a]
               for a in range(4) for b in range(4))
    assert trace_eta(H) == 0
    assert transverse(k,H)
    assert transverse(t0,H)

GT_DENSE = gauge(KUP[0],t0)
GX_DENSE = gauge(KUP[0],xx)
assert trace_eta(GT_DENSE) == -2

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
    )

K = tuple(sparse_vector(lower(k)) for k in KUP)
A = sparse_tensor(A_DENSE)
B = sparse_tensor(B_DENSE)
G_TIME = sparse_tensor(GT_DENSE)
G_X = sparse_tensor(GX_DENSE)

assert K == (
    (((0,),-1),((3,), 1)),
    (((0,),-1),((3,),-1)),
    (((0,), 1),((1,),-1)),
    (((0,), 1),((1,), 1)),
)
assert G_TIME == (
    ((0,0),2),((0,3),-1),((3,0),-1)
)

# Matching slots are zero-based in the code.
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) == 945
assert len(set(P_LIST)) == 945
assert len(LABELS) == 15120

def contract_sparse(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_value(I,P,J,Kdata,Hdata):
    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 = (
            Kdata[sigma[0]], Kdata[sigma[1]],
            Hdata[sigma[0]], Hdata[sigma[1]],
            Hdata[sigma[2]], Hdata[sigma[3]],
        )
        # Two Fourier derivatives contribute i^2 = -1.
        total -= contract_sparse(P,factors)
    return total

J_PLUS = (1,1,2,2)
J_MINUS = (2,2,1,1)
H_TIME = (G_TIME,A,B,B)
H_BLIND = (G_X,A,B,B)

P0 = ((0,1),(2,3),(4,6),(5,8),(7,9))
LABEL_A = ((1,1,2,2),P0)
LABEL_B = ((2,2,1,1),P0)

assert LABEL_A in LABELS
assert LABEL_B in LABELS

# Target arithmetic from the proved restricted Einstein cubic.
# This is not an independent implementation of Einstein vertices.
def v_abb(p0,q0,r0):
    assert p0+q0+r0 == 0
    return -2*(p0*q0+p0*r0+q0*r0)

W = (
    -v_abb(2,-1,-1),
     v_abb(1,0,-1),
     v_abb(1,0,-1),
)
assert W == (-6,2,2)

V4 = -sum(W)
X = tuple(-w for w in W)
M_VALUE = Fraction(-2*W[0]+W[1]+W[2],3)

assert V4 == 2
assert X == (6,-2,-2)
assert M_VALUE == Fraction(16,3)
assert Fraction(2,3)*V4-X[1]-X[2] == M_VALUE

assert V4-sum(X) == 0
assert 0*V4-sum(0*x for x in X) == 0
assert 2*V4-sum(2*x for x in X) == 0

def reconstruct_and_check():
    plus = {}
    blind = {}

    for label in LABELS:
        I,P = label
        plus[label] = raw_value(I,P,J_PLUS,K,H_TIME)
        blind[label] = raw_value(I,P,J_PLUS,K,H_BLIND)
        assert blind[label] == 0

    columns = {}
    for I,P in LABELS:
        Ibar = tuple(3-i for i in I)
        minus = raw_value(I,P,J_MINUS,K,H_TIME)
        assert minus == plus[(Ibar,P)]

        col = (blind[(I,P)],plus[(I,P)],minus)
        columns[(I,P)] = col

        # Complete factorization N = B_piv T.
        t1 = Fraction(col[1],8)
        t2 = Fraction(col[2],8)
        assert (0,8*t1,8*t2) == col

    ca = columns[LABEL_A]
    cb = columns[LABEL_B]
    assert ca == (0,8,0)
    assert cb == (0,0,8)

    # This minor uses ORIGINAL columns, not an arbitrary factor.
    original_minor = ca[1]*cb[2]-cb[1]*ca[2]
    assert original_minor == 64

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

    N_rows = tuple(
        tuple(columns[label][i] for label in LABELS)
        for i in range(3)
    )
    assert all(value == 0 for value in N_rows[0])

    # Selected permutation regressions; the general identity
    # is proved analytically in the manuscript.
    for label in (LABEL_A,LABEL_B):
        I,P = label
        expected = plus[label]
        for pi in PERMS:
            Jp = tuple(J_PLUS[i] for i in pi)
            Kp = tuple(K[i] for i in pi)
            Hp = tuple(H_TIME[i] for i in pi)
            assert raw_value(I,P,Jp,Kp,Hp) == expected

    # Rank lower bound: original minor 64.
    # Rank upper bound: zero first row.
    # Augmented rank: the target is the printed column fit.
    return {
        "shape": (3,len(LABELS)),
        "labels": LABELS,
        "N_rows": N_rows,
        "target": target,
        "original_minor": original_minor,
        "rank": 2,
        "augmented_rank": 2,
        "nonzero_alpha": (
            (LABEL_A,Fraction(2,3)),
            (LABEL_B,Fraction(2,3)),
        ),
        "left_kernel_basis": ((1,0,0),),
    }

if __name__ == "__main__":
    certificate = reconstruct_and_check()
    print("shape =",certificate["shape"])
    print("target =",certificate["target"])
    print("original minor =",certificate["original_minor"])
    print("rank pair =",(
        certificate["rank"],
        certificate["augmented_rank"],
    ))
    print("status = SEP-finite")