"""Exact, target-free test of the paper's Lie-closure forward program.

No empirical constants, stochastic search, floating-point rank or fitted weights.
Uses integer Gaussian matrices and Fraction Gaussian elimination.
This is a candidate amplitude lift of existing seed/loop data, not a claim
that this amplitude lift is already uniquely derived in the original Lean file.
"""
from fractions import Fraction
from functools import reduce
from itertools import product
from math import gcd
from pathlib import Path
import json


def mat(n, real=None, imag=None):
    return (n, tuple(real or [0] * (n*n)), tuple(imag or [0] * (n*n)))


def identity(n):
    return mat(n, [int(i == j) for i in range(n) for j in range(n)])


def add(a, b, factor=1):
    n, ar, ai = a
    _, br, bi = b
    return mat(n, [x+factor*y for x,y in zip(ar,br)], [x+factor*y for x,y in zip(ai,bi)])


def scale(a, k):
    return mat(a[0], [k*x for x in a[1]], [k*x for x in a[2]])


def times_i(a):
    return mat(a[0], [-x for x in a[2]], a[1])


def multiply(a, b):
    n, ar, ai = a
    _, br, bi = b
    rr, ii = [0]*(n*n), [0]*(n*n)
    for i in range(n):
        for k in range(n):
            x, y = ar[i*n+k], ai[i*n+k]
            if x or y:
                for j in range(n):
                    u, v = br[k*n+j], bi[k*n+j]
                    rr[i*n+j] += x*u-y*v
                    ii[i*n+j] += x*v+y*u
    return mat(n, rr, ii)


def commutator(a,b):
    return add(multiply(a,b),multiply(b,a),-1)


def primitive(a):
    d = reduce(gcd, (abs(x) for x in a[1]+a[2]), 0)
    if not d:
        return a
    return mat(a[0], [x//d for x in a[1]], [x//d for x in a[2]])


def skew_coords(a):
    n, r, s = a
    assert all(r[i*n+j] == -r[j*n+i] and s[i*n+j] == s[j*n+i]
               for i in range(n) for j in range(n))
    return [s[i*n+i] for i in range(n)] + [v for i in range(n) for j in range(i+1,n)
                                          for v in (r[i*n+j],s[i*n+j])]


class Span:
    def __init__(self):
        self.rows = []

    def add(self, coords):
        row = list(map(Fraction, coords))
        for pivot, old in self.rows:
            v = row[pivot]
            if v:
                row = [x-v*y for x,y in zip(row,old)]
        p = next((i for i,x in enumerate(row) if x), None)
        if p is None:
            return False
        v = row[p]
        row = [x/v for x in row]
        # Keep echelon rows sorted: this handles newly introduced earlier pivots.
        self.rows.append((p,row))
        self.rows.sort(key=lambda z:z[0])
        return True


def closure(generators):
    span, basis, words = Span(), [], []
    def accept(a, word):
        a = primitive(a)
        if span.add(skew_coords(a)):
            basis.append(a)
            words.append(word)
    for k,g in enumerate(generators):
        accept(g, f'g{k}')
    i = 0
    while i < len(basis):
        for j,g in enumerate(generators):
            accept(commutator(basis[i],g), f'[{words[i]},g{j}]')
        i += 1
    # Independently verify closure under ALL pairwise brackets, not just originals.
    for a in basis:
        for b in basis:
            assert not span.add(skew_coords(commutator(a,b)))
    return basis, words


def edge(n, i, j, sign):
    # Positive energy |psi_i-sign*psi_j|^2: rank-one signed edge Laplacian.
    a = [0]*(n*n)
    a[i*n+i] = a[j*n+j] = 1
    a[i*n+j] = a[j*n+i] = -sign
    return mat(n, a)


def traceless(a):
    n,r,s = a
    tr = sum(r[i*n+i] for i in range(n))
    assert not any(s)
    return add(scale(a,n),scale(identity(n),tr),-1)


def kron(a,b):
    n, ar, ai = a
    m, br, bi = b
    r, s = [], []
    for i in range(n):
        for k in range(m):
            for j in range(n):
                for l in range(m):
                    x,y,u,v = ar[i*n+j],ai[i*n+j],br[k*m+l],bi[k*m+l]
                    r.append(x*u-y*v)
                    s.append(x*v+y*u)
    return mat(n*m,r,s)


def dump(a):
    n,r,s = a
    return {'real':[list(r[i*n:(i+1)*n]) for i in range(n)],
            'imag':[list(s[i*n:(i+1)*n]) for i in range(n)]}


def determinant(rows):
    a = [list(map(Fraction,row)) for row in rows]
    det = Fraction(1)
    for k in range(len(a)):
        p = next((i for i in range(k,len(a)) if a[i][k]),None)
        if p is None: return 0
        if p != k:
            a[k],a[p]=a[p],a[k]
            det=-det
        v=a[k][k]
        det*=v
        for i in range(k+1,len(a)):
            ratio=a[i][k]/v
            for j in range(k+1,len(a)): a[i][j]-=ratio*a[k][j]
    assert det.denominator==1
    return int(det)


def main():
    rows = []
    for signs in product((-1,1), repeat=3):
        edges = [edge(3,i,j,s) for (i,j),s in zip(((0,1),(1,2),(2,0)),signs)]
        gs = [times_i(traceless(e)) for e in edges]
        bs, words = closure(gs)
        parity = signs[0]*signs[1]*signs[2]
        assert len(bs) == (8 if parity == -1 else 4)
        rows.append({'edge_signs':signs,'loop_holonomy':parity,'traceless_lie_dimension':len(bs)})
    color_generators = [times_i(traceless(edge(3,i,j,-1))) for i,j in ((0,1),(1,2),(2,0))]
    color, color_words = closure(color_generators)
    assert len(color) == 8
    g0,g1,g2=color_generators
    raw_words=[g0,g1,g2,commutator(g0,g1),commutator(g0,g2),commutator(g1,g2),
               commutator(commutator(g0,g1),g0),commutator(commutator(g0,g1),g2)]
    # Drop the third imaginary diagonal: trace zero determines it from the first two.
    raw_coordinates=[skew_coords(a)[:2]+skew_coords(a)[3:] for a in raw_words]
    det=determinant(raw_coordinates)
    assert det != 0
    # Explicit inclusion of every standard su(3) direction provides a basis check.
    span = Span()
    for b in color:
        span.add(skew_coords(b))
    for i,j in ((0,1),(1,2),(0,2)):
        r,s=[0]*9,[0]*9
        r[i*3+j],r[j*3+i]=1,-1
        s[i*3+j]=s[j*3+i]=1
        assert not span.add(skew_coords(mat(3,r)))
        assert not span.add(skew_coords(mat(3,imag=s)))
    for diagonal in ((1,-1,0),(0,1,-1)):
        s=[0]*9
        for i,x in enumerate(diagonal): s[i*3+i]=x
        assert not span.add(skew_coords(mat(3,imag=s)))

    X = mat(2,[0,1,1,0])
    Z = mat(2,[1,0,0,-1])
    weak_generators = [times_i(X), times_i(Z)]
    weak,_ = closure(weak_generators)
    assert len(weak) == 3
    one_weak_edge,_ = closure([times_i(traceless(edge(2,0,1,-1)))])
    assert len(one_weak_edge) == 1 # contrast observable really is required here

    # Local action on the existing ternary x binary factorization.
    factors = [kron(g,identity(2)) for g in color_generators] + [kron(identity(3),g) for g in weak_generators]
    local,_ = closure(factors)
    with_phase,_ = closure(factors+[times_i(identity(6))])
    assert len(local) == 11 and len(with_phase) == 12
    # Test the importance of factor locality, rather than hiding it.
    mixed = times_i(kron(mat(3,[1,0,0,0,-1,0,0,0,0]),Z))
    enlarged,_ = closure(factors+[mixed])
    assert len(enlarged) == 35

    # Uniform all-negative triangle: spectrum {1,1,4}, exact characteristic identity.
    triangle = add(add(edge(3,0,1,-1),edge(3,1,2,-1)),edge(3,2,0,-1))
    assert multiply(add(triangle,identity(3),-1),add(triangle,scale(identity(3),4),-1)) == mat(3)

    # Same seed: the six-step lift of a three-step closure has deck mirror T^3.
    # CRT: k -> (k mod 2, k mod 3); +3 flips Role and leaves Axis fixed.
    role_mirror=[(k+3)%6 for k in range(6)]
    closure_inverse=[(-k)%6 for k in range(6)]
    assert all(role_mirror[role_mirror[k]]==k and role_mirror[k]!=k for k in range(6))
    assert all(role_mirror[k]%2!=k%2 and role_mirror[k]%3==k%3 for k in range(6))
    L6=mat(6,[2*int(i==j)-int((i-j)%6 in (1,5)) for i in range(6) for j in range(6)])
    M=mat(6,[int(j==role_mirror[i]) for i in range(6) for j in range(6)])
    T=mat(6,[int(j==(i+1)%6) for i in range(6) for j in range(6)])
    assert multiply(multiply(T,T),T)==M
    assert multiply(M,M)==identity(6)
    assert commutator(L6,M)==mat(6)
    mirror_blocks=[]
    for parity in (1,-1):
        r=[]
        for i in range(3):
            for j in range(3):
                r.append(L6[1][i*6+j]+parity*L6[1][i*6+j+3])
        block=mat(3,r)
        # Verify the entire six-dimensional invariant embedding, not just upper rows.
        for i in range(6):
            for j in range(3):
                lhs=L6[1][i*6+j]+parity*L6[1][i*6+j+3]
                rhs=(1 if i<3 else parity)*block[1][(i%3)*3+j]
                assert lhs==rhs
        mirror_blocks.append(block)
    balanced=add(add(edge(3,0,1,1),edge(3,1,2,1)),edge(3,2,0,1))
    D=mat(3,[1,0,0,0,-1,0,0,0,1])
    assert mirror_blocks[0]==balanced
    assert multiply(multiply(D,mirror_blocks[1]),D)==triangle
    assert multiply(balanced,add(balanced,scale(identity(3),3),-1))==mat(3)
    uniform_only,_=closure([times_i(traceless(triangle))])
    assert len(uniform_only)==1
    assert any(commutator(g,triangle)!=mat(3) for g in color_generators)

    # Kernel of (A,B,z) -> z*(A tensor B) in the defining 3 x 2 representation.
    # Centers A=exp(2*pi*i*a/3), B=(-1)^b, z=exp(pi*i*k/3).
    center_kernel=[(a,b,k) for a in range(3) for b in range(2) for k in range(6)
                   if (2*a+3*b+k)%6==0]
    cyclic_kernel=[(j%3,j%2,j%6) for j in range(6)]
    assert set(center_kernel)==set(cyclic_kernel)
    # This yields representation congruences, not chirality or hypercharge uniquely.
    charge_residues={name:next(q for q in range(6) if (2*t+3*p+q)%6==0)
                     for name,t,p in [('color_doublet',1,1),('color_singlet',1,0),
                                       ('sterile_doublet',0,1),('sterile_singlet',0,0)]}
    assert charge_residues=={'color_doublet':1,'color_singlet':4,'sterile_doublet':3,'sterile_singlet':0}

    result = {
        'status':'Exact conditional construction; no empirical fitting; not a full seed-only derivation',
        'construction':'L_ij=(e_i-s_ij*e_j)(e_i-s_ij*e_j)^T; G_ij=i*(3*L_ij-2*I)',
        'triangle_sign_census':rows,
        'color':{'dimension':len(color),'algebra':'su(3)','generators':[dump(g) for g in color_generators],
                 'basis':[dump(b) for b in color],'basis_commutator_words':color_words,
                 'raw_coordinate_certificate':raw_coordinates,'raw_certificate_determinant':det},
        'weak':{'dimension':len(weak),'algebra':'su(2)','generators':'iX, iZ','one_edge_only_dimension':len(one_weak_edge)},
        'factor_local':{'dimension_without_center':len(local),'with_one_center':len(with_phase),
                        'algebra':'su(3) + su(2) + u(1)', 'one_mixed_generator_dimension':len(enlarged)},
        'negative_triangle_spectrum':[1,1,4],
        'same_seed_six_cycle_lift':{
            'common_clock':'T: k -> k+1 mod 6; T^3=M, T^6=I',
            'CRT_role_mirror':role_mirror,'closure_negation':closure_inverse,
            'operators_are_distinct':role_mirror!=closure_inverse,
            'mirror_even_block':dump(mirror_blocks[0]),'mirror_odd_block':dump(mirror_blocks[1]),
            'mirror_even_spectrum':[0,3,3],'mirror_odd_spectrum':[1,1,4],
            'bridge':'odd block is switching-equivalent to the inversion triangle',
            'scope':'collective parity modes across the two roles, not separate spatial and temporal coordinate sets'},
        'gauge_symmetry_guard':{
            'uniform_generator_alone_dimension':len(uniform_only),
            'all_color_generators_commute_with_uniform_energy':False,
            'uniform_energy_commutant_u3_dimension':5,
            'scope':'A dynamical/operator Lie algebra is not yet a local gauge symmetry. Gauge-covariant link dynamics remains to derive.'},
        'tensor_representation_center':{'kernel':center_kernel,'cyclic_generator':[1,1,1],
            'quotient':'(SU(3) x SU(2) x U(1))/Z6 for this representation',
            'hypercharge_integer_congruence':'2*color_triality + 3*weak_parity + q = 0 mod 6',
            'charge_residues':charge_residues,
            'scope':'Congruences only; integer winding shifts q -> q+6m remain, chirality not derived'},
        'assumptions_to_derive':['complex amplitude space / nontrivial sign representation of inversion',
            'rank-one edge mismatch energy and commutator closure',
            'independent availability of edge generators, not just their uniform sum',
            'binary contrast observable in addition to binary exchange',
            'factor-local gauge action', 'one central phase generator and its physical role'],
        'checks':['integer/Fraction arithmetic only','all 8 triangle signings','all pairwise bracket closure',
            'containment of a standard su(3) basis','mixed-generator countercheck','exact triangle polynomial'],
    }
    dest=Path(__file__).with_name('seed-lie-results.json')
    dest.write_text(json.dumps(result,indent=2)+'\n',encoding='utf8')
    print(json.dumps({k:v for k,v in result.items() if k!='color'},indent=2))
    print('color dimension:',len(color),'basis words:',color_words)
    print('Saved',dest)


if __name__=='__main__':
    main()