"""Exact Laurent-divisibility certificates for the eight E6 coefficient spaces.

No coefficient specialization is used.  Independent coefficient characters are
units factored out by the prescribed grading, as in the manuscript.
Requires SymPy. Run: python e6_recognition.py
"""
from pathlib import Path
import json
import sympy as sp
s,u,v=sp.symbols('s u v'); p=sp.symbols('p1:8'); variables=(s,u,v)+p
P=sp.prod(p); T=sp.prod(p[4:]); Q=p[3]*T; R=p[1]*p[2]*p[5]*p[6]
L=p[2]*T; H=p[0]*p[1]*p[3]*p[6]
E=sp.expand((1+R)*(1+Q)+u*v*L*(1+H))
monE=sp.Poly(E,*variables).monoms()
poly_terms=lambda f:sp.Poly(sp.expand(f),*variables).monoms()
packets=[
 ('f',E,s*u,E),
 ('ubar',1+Q,u,1+Q),
 ('vbar',1+Q+s+s*Q,v*Q,1+Q+s),
 ('k',1+s,P,1+s),
 ('g',E+s+s*R,s*v*T,E+s*(1+R)),
 ('W3',E+s+s*Q,s*u*p[2]*p[5]*p[6],E+s*(1+Q)),
 ('W9',E+s*(1+R+Q+R*Q+u*v*L),s*u*P,E+s*(1+R+Q+R*Q+u*v*L)),
 ('W12',E+s*(1+Q+R+u*v*L)+s**2,s*v*p[1]*p[2]*p[3]*p[4]*p[5]*p[6]**2,E+s*(2+Q+R+u*v*L)+s**2)
]
exchanges=[E,1+Q,1+s+Q]+[1+s]*7
report=[]
for name,cap,den,want in packets:
 mons=poly_terms(cap); denom=sp.Poly(den,*variables).monoms()[0]
 constraints={}
 for i,exchange in enumerate(exchanges):
  spectators=variables[:i]+variables[i+1:]
  values=sorted(set(m[i]-denom[i] for m in mons if m[i]<denom[i]))
  for degree in values:
   divisor=sp.Poly(exchange**(-degree),*spectators)
   for j,m in enumerate(mons):
    if m[i]-denom[i]!=degree:continue
    coefficient=sp.prod(variables[k]**m[k] for k in range(len(variables)) if k!=i)
    rem=sp.Poly(sp.reduced(coefficient,[divisor.as_expr()],*spectators)[1],*spectators)
    for exp,c in rem.terms():
     if c:
      key=(i,degree,exp)
      constraints.setdefault(key,[sp.Rational(0)]*len(mons))[j]+=c
 M=sp.Matrix(list(constraints.values()))
 rank=M.rank(); ns=M.nullspace()
 
 if len(ns)!=1:
  wantvec=sp.Matrix([sp.Poly(want,*variables).coeff_monomial(m) for m in mons])
  for key,row in constraints.items():
   val=sum(row[j]*wantvec[j] for j in range(len(mons)))
   if val: print('BAD',key,'VALUE',val,'row',row)
  raise AssertionError((name,M.shape,rank))
 const_index=mons.index((0,)*10)
 sol=ns[0]/ns[0][const_index]
 expected=sp.Poly(want,*variables)
 assert all(sol[j]==expected.coeff_monomial(m) for j,m in enumerate(mons))
 assert all(c.q==1 and c>=0 for c in sol)
 # An exact independent row subset and maximal minor certify rank >= n-1.
 _,pivot_rows=M.T.rref(); I=list(pivot_rows)
 M0=M[I,:]; _,pivot_cols=M0.rref()
 minor=M0[:,list(pivot_cols)].det()
 assert len(I)==len(mons)-1 and minor!=0
 report.append({'name':name,'variables':[str(x) for x in variables],
  'numerator_monomials':[list(m) for m in mons], 'denominator_exponents':list(denom),
  'constraint_rows':[[str(c) for c in row] for row in M.tolist()],
  'rank':rank,'kernel':[str(c) for c in sol],'normalization_column':const_index,
  'independent_rows':I,'minor_columns':list(pivot_cols),'minor':str(minor)})
 print(name,'support',len(mons),'equations',M.rows,'rank',rank,'minor',minor,flush=True)
out=Path(__file__).with_name('e6_recognition_certificates.json')
out.write_text(json.dumps({'convention':'Spectator polynomial divisibility; coefficients remain independent units.','packets':report},indent=2))
print('PASS: eight exact one-dimensional spaces; normalized kernels match all displayed formulas.')