#!/usr/bin/env python3
"""Verify the saved exceptional certificates using only Python's standard library.

The global E7 substitution matrices and all 100 support-subset factorizations
are reconstructed, not merely read as asserted bounds.  This checks finite
algebra, not the geometric support bounds or enumerative gluing theorems.
"""
from __future__ import annotations
import ast
import json
from collections import Counter
from fractions import Fraction as Q
from itertools import combinations
from math import comb, gcd
from pathlib import Path

ROOT=Path(__file__).resolve().parent

def require(condition, message):
    if not condition:
        raise AssertionError(message)

def determinant(rows):
    a=[[Q(c) for c in row] for row in rows]
    n=len(a)
    if not n:return Q(1)
    require(all(len(r)==n for r in a),'nonsquare determinant')
    ans=Q(1)
    for i in range(n):
        j=next((j for j in range(i,n) if a[j][i]),None)
        if j is None:return Q(0)
        if j!=i:a[i],a[j]=a[j],a[i];ans=-ans
        p=a[i][i];ans*=p
        for j in range(i+1,n):
            t=a[j][i]/p
            if t:
                for k in range(i+1,n):a[j][k]-=t*a[i][k]
    return ans

def verify_kernel(M, v, rank, ri, ci, minor):
    v=list(map(Q,v));M=[list(map(Q,r)) for r in M]
    require(all(sum(a*b for a,b in zip(row,v))==0 for row in M),'wrong kernel')
    require(any(v) and rank==len(v)-1,'wrong nullity certificate')
    require(len(ri)==len(ci)==rank,'wrong maximal minor size')
    d=determinant([[M[i][j] for j in ci] for i in ri])
    require(d==Q(minor) and abs(d)==1,'non-unit or incorrect minor')

# E6: exact kernels and integral normalizations.
e6=json.loads((ROOT/'e6_recognition_certificates.json').read_text())['packets']
require(len(e6)==8,'missing E6 cases')
for p in e6:
    verify_kernel(p['constraint_rows'],p['kernel'],p['rank'],p['independent_rows'],p['minor_columns'],p['minor'])
    require(Q(p['kernel'][p['normalization_column']])==1,'E6 normalization')
    require(all(Q(c).denominator==1 and Q(c)>=0 for c in p['kernel']),'E6 coefficient integrality')
require(max(map(Q,next(p for p in e6 if p['name']=='W12')['kernel']))==2,'W12 coefficient two')

# E7: eighteen boundary coefficient spaces.
e7b=json.loads((ROOT/'e7_boundary_certificates.json').read_text())['boundary_cases']
require([r['boundary'] for r in e7b]==list(range(18)),'missing E7 boundary cases')
for p in e7b:
    verify_kernel(p['rows'],[1]*len(p['support']),p['rank'],p['independent_rows'],p['minor_columns'],p['minor'])

# Integral source marking and twelve source-torus normalizations.
cal=json.loads((ROOT/'e7_calibration.json').read_text());C=cal['C'];N=len(C)
require(N==28 and determinant(C)==Q(cal['determinant']) and abs(cal['determinant'])==1,'full marking')
require(len(cal['source_exchanges'])==10,'source exchanges missing')
source_terms=0
for j,s in enumerate(cal['source_exchanges']):
    terms=s['exchange_terms'];source_terms+=len(terms);L=s['least_exponent']
    require(len(terms)==len(s['normalized_support']),'normalized support length')
    for (e,coeff),alpha in zip(terms,s['normalized_support']):
        require(Q(coeff)!=0,'zero exchange term')
        require(all(sum(C[k][i]*alpha[k] for k in range(N))==e[i]-L[i] for i in range(N)),'normalized source character')
        require(all(a>=0 for a in alpha[:20]),'not least in mutable partial order')
    H=[L[i]-s['normalizer_exponents'][i]-(i==j) for i in range(N)]
    require(H==s['first_mutant_leading_exponent'],'associated normalizer omitted')
selected=[]
for row in cal['rows']:
    for j,power,hd,minimum,defect in row['choices']:
        beta=row['beta']+[0]*18;mu=row['label']+[0]*8;s=cal['source_exchanges'][j]
        require(sum(beta)==1,'mixed direction mislabeled as single contact')
        got_hd=sum(x*y for x,y in zip(s['first_mutant_leading_exponent'],beta))
        got_min=min(sum(x*y for x,y in zip(alpha,mu)) for alpha in s['normalized_support'])
        got_defect=power*got_hd+sum(beta[i] for i in range(10) if i!=j)
        require((got_hd,got_min,got_defect)==(hd,minimum,defect)==(-1,-1,-1),'coefficient calibration')
        selected.append(row['target'])
require(sorted(selected)==[2,3,4,6,7,9,12,13,16,17,19,20],'wrong twelve directions')

# Small exact Laurent polynomial implementation for the Freudenthal certificate.
n=20;zero=(0,)*n

def add(*polys):
    ans={}
    for p in polys:
        for e,c in p.items():ans[e]=ans.get(e,Q(0))+c
    return {e:c for e,c in ans.items() if c}

def scale(p,c):return {e:v*c for e,v in p.items() if v*c}

def multiply(p,q):
    ans={}
    for a,c in p.items():
        for b,d in q.items():
            e=tuple(x+y for x,y in zip(a,b));ans[e]=ans.get(e,Q(0))+c*d
    return {e:c for e,c in ans.items() if c}

def power(p,k):
    if k<0:
        require(len(p)==1,'negative power of a nonmonomial')
        e,c=next(iter(p.items()));return {tuple(k*x for x in e):c**k}
    ans={zero:Q(1)}
    while k:
        if k%2:ans=multiply(ans,p)
        k//=2
        if k:p=multiply(p,p)
    return ans

def parse(text):
    def visit(node):
        if isinstance(node,ast.Constant):return {zero:Q(node.value)}
        if isinstance(node,ast.Name):
            require(node.id.startswith('v') and node.id[1:].isdigit(),'bad variable')
            i=int(node.id[1:])-1;require(0<=i<n,'bad index')
            return {tuple(int(j==i) for j in range(n)):Q(1)}
        if isinstance(node,ast.UnaryOp) and isinstance(node.op,ast.USub):return scale(visit(node.operand),-1)
        if isinstance(node,ast.BinOp):
            a=visit(node.left)
            if isinstance(node.op,ast.Pow):return power(a,int(ast.literal_eval(node.right)))
            b=visit(node.right)
            if isinstance(node.op,ast.Add):return add(a,b)
            if isinstance(node.op,ast.Sub):return add(a,scale(b,-1))
            if isinstance(node.op,ast.Mult):return multiply(a,b)
            if isinstance(node.op,ast.Div):return multiply(a,power(b,-1))
        raise ValueError(f'Unsupported expression: {ast.dump(node)}')
    return visit(ast.parse(text,mode='eval').body)

def normalize(p):
    ell=tuple(min(e[i] for e in p) for i in range(n))
    return ell,{tuple(e[i]-ell[i] for i in range(n)):c for e,c in p.items()}

def terms(data):return {tuple(e):Q(c) for e,c in data if Q(c)}
def polykey(p):return tuple(sorted(p.items()))
def factor_key(p):
    _,p=normalize(p)
    # Saved factors have integral coefficients. Normalize sign and content.
    nums=[int(c) for c in p.values()];require(all(c.denominator==1 for c in p.values()),'nonintegral factor')
    content=0
    for c in nums:content=gcd(content,abs(c))
    if p[max(p)]<0:content=-content
    return polykey(scale(p,Q(1,content)))

D=json.loads((ROOT/'e7_divisor_certificates.json').read_text())
func=[parse(f) for f in D['functions']]
require(len(func)==D['nfunctions']==24,'wrong E7 function set')
normalized=[normalize(f) for f in func]
expected_subsets=[]
for ell,p in normalized:
    its=list(p.items())
    for k in range(2,len(its)+1):
        for S in combinations(its,k):expected_subsets.append(polykey(normalize(dict(S))[1]))
require(len(expected_subsets)==D['subset_checks']==100,'wrong subset count')
require(Counter(expected_subsets)==Counter(polykey(terms(f['polynomial_terms'])) for f in D['factorizations']),'support-subset factorization coverage')
keys={factor_key(terms(c['factor_terms'])) for c in D['cases']}
require(len(keys)==D['nfactors']==80,'factor count')
for f in D['factorizations']:
    product={zero:Q(f['coefficient'])}
    for term in f['factors']:
        q=terms(term['terms']);require(factor_key(q) in keys,'unlisted subset factor')
        product=multiply(product,power(q,term['power']))
    require(product==terms(f['polynomial_terms']),'factorization identity')

def check_dual(case, rhs, H):
    A=case['rows'];lam={int(i):Q(c) for i,c in case['inequalities']};mu={int(i):Q(c) for i,c in case['equations']}
    require(all(v>=0 for v in lam.values()),'negative inequality multiplier')
    coeff=[sum(lam[j]*A[j][i] for j in lam)+sum(mu[j]*H[j][i] for j in mu) for i in range(n)]
    require(all(c<=-1 for c in coeff),'dual does not bound sum of weights')
    bound=-sum(lam[j]*rhs[j] for j in lam)
    require(bound==Q(case['upper_bound']) and bound in (0,1),'wrong bound')
    return bound

z=D['zero_case'];require(set(map(tuple,z['rows']))==set(e for f in func for e in f),'zero-face rows')
require(check_dual(z,[0]*len(z['rows']),[])==0,'zero-face bound')
counts=Counter();max_den=max_ineq=max_eq=0
for case in D['cases']:
    q=terms(case['factor_terms']);require(parse(case['factor'])==q,'factor expression')
    p=case['pivot']-1
    require(max(e[p] for e in q)==1 and min(e[p] for e in q)==0,'factor not linear')
    A={tuple(e[i]-(i==p) for i in range(n)):c for e,c in q.items() if e[p]==1}
    B={e:c for e,c in q.items() if e[p]==0}
    require(len(A)==1 and B,'pivot must have a nonzero monomial leading coefficient')
    ae,ac=next(iter(A.items()))
    require(all(ae[i]==0 or min(e[i] for e in B)==0 for i in range(n)),'linear factor has a common monomial divisor')
    alpha=max(q)
    H=[tuple(e[i]-alpha[i] for i in range(n)) for e in sorted(q,reverse=True)[1:]]
    require(H==list(map(tuple,case['homogeneity'])),'wrong homogeneity equations')
    rows=set()
    for ell,P in normalized:
        degree=max(e[p] for e in P);R={}
        for k in range(degree+1):
            a={tuple(e[i]-k*(i==p) for i in range(n)):c for e,c in P.items() if e[p]==k}
            if not a:continue
            pre=multiply(a,power(A,degree-k))
            for j in range(k+1):
                R[j]=add(R.get(j,{}),scale(multiply(pre,power(B,k-j)),Q(comb(k,j)*(-1)**(k-j))))
        for j,Rj in R.items():
            for ex in Rj:
                row=tuple(ell[i]+ex[i]-degree*ae[i]+j*alpha[i] for i in range(n))
                rows.add((row,-j))
    require(rows==set(zip(map(tuple,case['rows']),case['rhs'])),'incorrect substitution inequalities')
    bound=check_dual(case,case['rhs'],H);counts[int(bound)]+=1
    values=[Q(x) for _,x in case['inequalities']+case['equations']]
    den=max([x.denominator for x in values]+[1]);require(den==case['denominator_max'],'denominator bound')
    max_den=max(max_den,den);max_ineq=max(max_ineq,len(case['inequalities']));max_eq=max(max_eq,len(case['equations']))
require(counts=={0:70,1:10},'incorrect outcome distribution')
report={
 'status':'PASS',
 'E6':{'coefficient_systems':8,'maximal_minors':'all +/-1','W12_coefficient':2},
 'E7':{'full_lattice_rank':28,'determinant':cal['determinant'],'source_exchange_terms':source_terms,
       'source_torus_calibrated_directions':sorted(selected),'mixed_contact_directions_retained':[1,5,8,10,11,14,15,18],
       'boundary_systems':18,'support_subset_factorizations':100,'irreducible_factor_cases':80,
       'bound_zero':70,'bound_one':10,'max_denominator':max_den,'max_inequality_rows_used':max_ineq,'max_homogeneity_rows_used':max_eq},
 'scope':'Exact finite algebra from the saved formulas. No verification of enumerative geometry, full source support bounds, or the eight mixed-contact incidence arguments.'}
(ROOT/'verification_report.json').write_text(json.dumps(report,indent=2)+'\n')
print(json.dumps(report,indent=2))