#!/usr/bin/env python3
"""Exact rational checks for the finite-dimensional seven-label Tiles proof.

No third-party package is used.  The script verifies coefficient flatness,
absence of complex rank-one seed matrices by exhaustive linear branches,
six exceptional contractions, the displayed integer minors, and the
pairwise and triple skew-evaluation minors.
"""
from fractions import Fraction as Q
from itertools import product, combinations, combinations_with_replacement

def rref(A):
    A=[[Q(x) for x in row] for row in A]
    m=len(A); n=len(A[0]) if m else 0
    piv=[]; i=0
    for j in range(n):
        p=next((k for k in range(i,m) if A[k][j]),None)
        if p is None:continue
        A[i],A[p]=A[p],A[i]
        z=A[i][j];A[i]=[x/z for x in A[i]]
        for k in range(m):
            if k!=i and A[k][j]:
                z=A[k][j];A[k]=[x-z*y for x,y in zip(A[k],A[i])]
        piv.append(j);i+=1
        if i==m:break
    return A,piv

def rank(A): return len(rref(A)[1]) if A else 0
def det(A):
    A=[[Q(x) for x in row] for row in A]; n=len(A);v=Q(1)
    for j in range(n):
        p=next((i for i in range(j,n) if A[i][j]),None)
        if p is None:return Q(0)
        if p!=j:A[j],A[p]=A[p],A[j];v=-v
        z=A[j][j];v*=z
        for i in range(j+1,n):
            t=A[i][j]/z
            for k in range(j+1,n):A[i][k]-=t*A[j][k]
    return v

def transpose(A): return [list(x) for x in zip(*A)]
def mul(A,B): return [[sum(x*y for x,y in zip(row,col)) for col in zip(*B)] for row in A]

# Integer basis matrices, in parameter order a,b,c,d.
S=[[[1,1,0],[0,-2,0],[0,0,0]],
   [[0,0,1],[0,-2,1],[0,0,0]],
   [[0,0,0],[0,-2,0],[0,1,1]],
   [[0,0,0],[1,-2,0],[1,0,0]]]
L=[[[S[h][rho][b] for b in range(3)] for h in range(4)] for rho in range(3)]
flat=[[S[h][rho][b] for rho in range(3) for b in range(3)] for h in range(4)]
G=mul(flat,transpose(flat))
aug=[G[i]+[int(i==j) for j in range(4)] for i in range(4)]
QG=[row[4:] for row in rref(aug)[0]]
P=mul(mul(transpose(flat),QG),flat)
assert rank(P)==4
assert all(P[3*r+b][3*s+c]==P[3*r+c][3*s+b]
           for r,s,b,c in product(range(3),repeat=4))
assert mul(P,P)==P

# Each product equation chooses its vanishing u-factor or v-factor.
uf=[(1,0,0),(1,-1,0),(0,1,-1),(0,0,1)]
vf=[(1,-1,0),(0,0,1),(1,0,0),(0,1,-1)]
for choices in product((0,1),repeat=4):
    uc=[uf[i] for i in range(4) if choices[i]==0]
    vc=[vf[i] for i in range(4) if choices[i]==1]
    for stopper_on_u in (False,True):
        uu=uc+([(1,1,1)] if stopper_on_u else [])
        vv=vc+([] if stopper_on_u else [(1,1,1)])
        assert rank(uu)==3 or rank(vv)==3

special=[(1,-1,0),(0,0,1),(0,1,-1),(1,0,0),(1,1,1),(-2,1,-2)]
for lam in special:
    contraction=[[sum(S[h][rho][b]*lam[b] for b in range(3))
                  for rho in range(3)] for h in range(4)]
    assert rank(contraction)==2

sym3=list(combinations_with_replacement(range(3),2))
eval_cols=[]
for rho in range(3):
    for i,j in sym3:
        H=[[0]*3 for _ in range(3)];H[i][j]=H[j][i]=1
        image=mul(L[rho],H)
        eval_cols.append([x for row in image for x in row])
E=transpose(eval_cols)
chosen=[0,1,2,4,5,6,7,8,9,10,13,14]
eval_minor=det([[row[j] for j in chosen] for row in E])
assert eval_minor==20

pairs=list(combinations(range(3),2))
ord_rows=[]
for i,j in combinations_with_replacement(range(4),2):
    row=[]
    for rho,sigma in pairs:
        for b,c in pairs:
            row.append(L[rho][j][b]*L[sigma][i][c]-L[sigma][j][b]*L[rho][i][c]
                      +L[rho][i][b]*L[sigma][j][c]-L[sigma][i][b]*L[rho][j][c])
    ord_rows.append(row)
ordinary_minor=det(ord_rows[:9])
assert ordinary_minor==448

print('Positive flat Gram: rational projection rank 4, partial transpose fixed.')
print('Complex rank-one coefficient matrices: all 32 factor branches excluded.')
print('All six displayed exceptional contractions have exact rank 2.')
print('Symmetric-evaluation minor =',eval_minor,'; ordinary coefficient minor =',ordinary_minor)

# Pairwise skew-evaluation minors in the manuscript.
skew_matrices = []
for i, j in pairs:
    matrix = [[0] * 3 for _ in range(3)]
    matrix[i][j] = 1
    matrix[j][i] = -1
    skew_matrices.append(matrix)
for (rho, sigma), selected, expected in [
    ((0, 1), (0, 1, 2, 3, 4, 5), -4),
    ((0, 2), (0, 2, 3, 6, 7, 10), -1),
    ((1, 2), (0, 2, 3, 6, 7, 10), -4),
]:
    columns = [[x for row in mul(L[a], matrix) for x in row]
               for a in (rho, sigma) for matrix in skew_matrices]
    evaluation = transpose(columns)
    value = det([evaluation[i] for i in selected])
    assert value == expected, (rho, sigma, value, expected)
    print('Skew-evaluation pair', (rho, sigma), 'minor =', value)

# Triple skew-evaluation minor (injectivity of so(3)^3 -> Mat_{4,3}).
columns = [[x for row in mul(L[a], matrix) for x in row]
           for a in (0, 1, 2) for matrix in skew_matrices]
evaluation = transpose(columns)
value = det([evaluation[i] for i in (0, 1, 2, 3, 4, 5, 6, 7, 10)])
assert value == -4, value
print('Skew-evaluation triple minor =', value)