#!/usr/bin/env python3
"""Exact finite certificates for the explicitly scoped OMEGA joining theorem.

Python 3.10+, standard library only. All acceptance checks use integers/Fraction.
The general theorem and its hypotheses are proved in JOINING_THEOREM.md.
This is not an end-to-end Lean proof, a maximum-fiber calculation, or a novelty claim.
Graph generation and polynomial-inverse method reproduce the supplied Paper IX
construction, rather than claiming to be an independent third-party reproduction.
"""
from __future__ import annotations
import argparse
import json
from collections import Counter
from fractions import Fraction as F
from hashlib import sha256
from itertools import combinations, permutations, product
from math import gcd, isqrt
from pathlib import Path
from random import Random

Q = 339949771344778
S = 10**15
D = 25*S
q = F(Q, S)
GRAPH_SHA256 = '40bc6a1d52efab3d75cf343e65ac06b314f9c9c375fc02bf6d45bf6bb21f9529'


def require(ok, message):
    if not ok:
        raise ArithmeticError(message)


def canonical(obj):
    return (json.dumps(obj, sort_keys=True, separators=(',', ':'))+'\n').encode()


def pmul(x, y):
    a,b=x; c,d=y
    return (a*c+b*d, a*d+b*c+b*d)


def padd(x, y):
    return (x[0]+y[0], x[1]+y[1])


def pscale(a, x):
    return (a*x[0], a*x[1])


def psum(xs):
    out=(0,0)
    for x in xs:
        out=padd(out,x)
    return out


def graph():
    vertices=set()
    for i in range(4):
        for sign in (-1,1):
            v=[(0,0)]*4; v[i]=(2*sign,0); vertices.add(tuple(v))
    for signs in product((-1,1), repeat=4):
        vertices.add(tuple((s,0) for s in signs))
    even=[p for p in permutations(range(4))
          if sum(p[i]>p[j] for i in range(4) for j in range(i+1,4))%2==0]
    for a,b,c in product((-1,1), repeat=3):
        v=((0,0),(a,0),(0,b),(-c,c))
        vertices.update(tuple(v[i] for i in p) for p in even)
    vertices=sorted(vertices)
    neighbors=[[] for _ in vertices]
    for i,j in combinations(range(len(vertices)),2):
        dist=psum(pmul((a-c,b-d),(a-c,b-d))
                  for (a,b),(c,d) in zip(vertices[i],vertices[j]))
        if dist==(8,-4):
            neighbors[i].append(j);neighbors[j].append(i)
    fixture={'coordinate_rule':'vertex coordinate is (a+b*phi)/2 for pair [a,b]',
             'neighbors':neighbors,'schema':'OMEGA_600CELL_GRAPH_FIXTURE/1',
             'vertex_order':'Python lexicographic order of tuples of four integer pairs',
             'vertices':vertices}
    require(sha256(canonical(fixture)).hexdigest()==GRAPH_SHA256,'graph identity')
    require(len(vertices)==120 and all(len(r)==12 for r in neighbors),'graph dimensions')
    return vertices,neighbors


def apply_M(v, neighbors):
    return [13*v[i]+sum(v[j] for j in row) for i,row in enumerate(neighbors)]


def apply_matrix(N,v):
    return [sum(a*b for a,b in zip(row,v)) for row in N]


def inverse_certificate(neighbors):
    factors=[[-r,1] for r in (10,11,13,16,25)]+[[211,-32,1],[205,-30,1]]
    poly=[1]
    for f in factors:
        out=[0]*(len(poly)+len(f)-1)
        for i,a in enumerate(poly):
            for j,b in enumerate(f):out[i+j]+=a*b
        poly=out
    n=120
    def mm(X):
        return [[13*X[i][k]+sum(X[j][k] for j in neighbors[i])
                 for k in range(n)] for i in range(n)]
    N=[[poly[-1]*int(i==j) for j in range(n)] for i in range(n)]
    for k in range(len(poly)-2,0,-1):
        N=mm(N)
        for i in range(n):N[i][i]+=poly[k]
    den=-poly[0];common=den
    for row in N:
        for a in row:common=gcd(common,a)
    h=den//common;N=[[a//common for a in row]for row in N]
    prod=mm(N)
    require(all(prod[i][j]==h*int(i==j)for i in range(n)for j in range(n)),'MN=hI')
    require(all(N[i][j]==N[j][i]for i in range(n)for j in range(n)),'symmetric inverse')
    require(h==4948372000,'inverse denominator')
    norm=max(sum(abs(a) for a in row)for row in N)
    require(norm==801320608,'inverse infinity-norm numerator')
    return N,h,norm


def positive_definiteness_certificate(neighbors):
    n=120
    # M-9I = 4I+A. Bareiss pivots are leading principal determinants;
    # every pivot is required positive, so no row swaps are used.
    B=[[4 if i==j else int(j in neighbors[i]) for j in range(n)]for i in range(n)]
    previous=1;pivots=[]
    for k in range(n):
        pivot=B[k][k];require(pivot>0,'Sylvester leading minor');pivots.append(pivot)
        for i in range(k+1,n):
            for j in range(k+1,n):
                val,rem=divmod(B[i][j]*pivot-B[i][k]*B[k][j],previous)
                require(rem==0,'Bareiss exact division');B[i][j]=val
            B[i][k]=0
        previous=pivot
    return pivots



def lambda_enclosure():
    """Rational enclosure using phi and Machin's identity for pi."""
    scale=10**40; root=isqrt(5*scale*scale)
    require(root*root < 5*scale*scale < (root+1)**2, 'sqrt(5) enclosure')
    pl=F(scale+root,2*scale);pu=F(scale+root+1,2*scale)
    def atan_bounds(n,terms):
        a=sum((F((-1)**k,(2*k+1)*n**(2*k+1)) for k in range(terms)),F(0))
        b=a+F((-1)**terms,(2*terms+1)*n**(2*terms+1))
        return min(a,b),max(a,b)
    al,au=atan_bounds(5,48);bl,bu=atan_bounds(239,16)
    pil=16*al-4*bu;piu=16*au-4*bl
    lo=pl/(pl+piu);hi=pu/(pu+pil)
    require(q<lo<hi<F(Q+1,S), 'q != ideal lambda, exact enclosure')
    t=2*S
    require(2*Q+1<=lo*t and hi*t<2*Q+2, 'large-state ideal/rational floor difference')
    return {'input_scalar':str(t),'rational_floor':str(2*Q),
            'ideal_floor_certified':str(2*Q+1),'q_strictly_less_than_lambda':True,
            'method':'rational sqrt and alternating arctangent bounds via Machin identity'}


def check_d(d):
    if len(d)!=120 or any(type(x)is not int for x in d):
        raise ValueError('120 Python integers required')


def forward(d,neighbors):
    check_d(d)
    pairs=[divmod(Q*z,D)for z in apply_M(d,neighbors)]
    return [a for a,b in pairs],[b for a,b in pairs]


def decode_full(y,r,N,h):
    check_d(y);check_d(r)
    if any(not 0<=x<D for x in r):raise ValueError('remainder bounds')
    mixed=[]
    for yi,ri in zip(y,r):
        z,rem=divmod(D*yi+ri,Q)
        if rem:raise ValueError('mixed-coordinate divisibility')
        mixed.append(z)
    d=[]
    for z in apply_matrix(N,mixed):
        a,rem=divmod(z,h)
        if rem:raise ValueError('lattice divisibility')
        d.append(a)
    return d


def ceil_div(a,b):
    return -((-a)//b)


def decode_mod8(y,s,neighbors,N,h,max_nodes=200000):
    check_d(y);check_d(s)
    if any(not 0<=x<8 for x in s):raise ValueError('residue bounds')
    den=2*Q*h
    centers=[D*sum(a*(2*b+1)for a,b in zip(row,y))for row in N]
    radii=[D*sum(abs(a)for a in row)for row in N]
    ds=[]
    for center,radius,residue in zip(centers,radii,s):
        lo=ceil_div(center-radius,den);hi=(center+radius)//den
        first=lo+(residue-lo)%8
        choices=tuple(range(first,hi+1,8))
        if not choices:raise ValueError('empty coordinate domain')
        require(len(choices)<=2,'coordinate candidate count')
        ds.append(choices)
    low=[ceil_div(D*a,Q)for a in y]
    high=[ceil_div(D*(a+1),Q)-1 for a in y]
    rows=[[(i,13)]+[(j,1)for j in near]for i,near in enumerate(neighbors)]
    stats={'nodes':0,'sweeps':0,'initial_ambiguous':sum(len(v)>1 for v in ds)}
    def propagate(domains):
        change=True
        while change:
            change=False;stats['sweeps']+=1
            for i,row in enumerate(rows):
                mn=sum(w*domains[j][0]for j,w in row)
                mx=sum(w*domains[j][-1]for j,w in row)
                if mn>high[i]or mx<low[i]:return False
                for j,w in row:
                    if len(domains[j])==1:continue
                    otherlo=mn-w*domains[j][0];otherhi=mx-w*domains[j][-1]
                    keep=tuple(a for a in domains[j]
                               if otherlo+w*a<=high[i] and otherhi+w*a>=low[i])
                    if not keep:return False
                    if keep!=domains[j]:domains[j]=keep;change=True
        return True
    def search(domains):
        stats['nodes']+=1
        if stats['nodes']>max_nodes:raise RuntimeError('UNRESOLVED: decoder resource limit')
        if not propagate(domains):return None
        undecided=[i for i,v in enumerate(domains)if len(v)>1]
        if not undecided:
            d=[v[0]for v in domains]
            require(forward(d,neighbors)[0]==y,'decoded forward verification')
            return d
        i=max(undecided,key=lambda j:sum(len(domains[k])>1 for k in neighbors[j]))
        for a in sorted(domains[i],key=lambda a:(abs(a*den-centers[i]),a)):
            branch=domains[:];branch[i]=(a,)
            ans=search(branch)
            if ans is not None:return ans
        return None
    ans=search(ds)
    if ans is None:raise ValueError('no valid predecessor')
    require([x%8 for x in ans]==s,'decoded residues')
    return ans,stats


def run():
    vertices,neighbors=graph()
    ideal_floor_counterexample=lambda_enclosure()
    N,h,norm=inverse_certificate(neighbors)
    pivots=positive_definiteness_certificate(neighbors)
    kappa=F(D*norm,Q*h)
    require(F(11)<kappa<F(12),'inverse bound')
    require(0<Q<S and 73*Q<D<74*Q and 80*Q>D,'quantizer separation')
    # Actual rows of the coordinate matrix C are half these integer pairs.
    coord_checks=0
    for i in range(120):
        for a in range(4):
            require(psum(vertices[j][a]for j in neighbors[i])==
                    pmul((0,6),vertices[i][a]),'A C = 6 phi C')
            coord_checks+=1
    for a in range(4):
        require(psum(v[a]for v in vertices)==(0,0),'centered coordinate columns')
        for b in range(4):
            require(psum(pmul(v[a],v[b])for v in vertices)==
                    (120*int(a==b),0),'C transpose C = 30 I')
    # Restriction-of-scalars bridge: 2C=C0+phi*C1, with integer C0,C1.
    C0=[[v[j][0] for j in range(4)]for v in vertices]
    C1=[[v[j][1] for j in range(4)]for v in vertices]
    for a in range(4):
        u=[v[a]for v in C0];w=[v[a]for v in C1]
        require(sum(u)==sum(w)==0,'centered integer geometry modes')
        require(apply_M(u,neighbors)==[13*x+6*y for x,y in zip(u,w)],'M C0=13 C0+6 C1')
        require(apply_M(w,neighbors)==[6*x+19*y for x,y in zip(u,w)],'M C1=6 C0+19 C1')
        for j in range(4):
            require(sum(v[a]*v[j]for v in C0)==72*int(a==j),'C0 Gram')
            require(sum(v[a]*v[j]for v in C1)==48*int(a==j),'C1 Gram')
            require(sum(u[a]*v[j]for u,v in zip(C0,C1))==-24*int(a==j),'mixed Gram')
    def integer_geom_encode(t,x,z):
        return [t+sum(a*b+c*e for a,b,c,e in zip(u,x,v,z))for u,v in zip(C0,C1)]
    def integer_geom_decode(d):
        t=F(sum(d),120)
        x=[F(sum((2*C0[i][j]+C1[i][j])*d[i]for i in range(120)),120)for j in range(4)]
        z=[F(sum((C0[i][j]+3*C1[i][j])*d[i]for i in range(120)),120)for j in range(4)]
        return t,x,z
    integer_geometry_records=[]
    rng0=Random(7072026)
    for _ in range(12):
        t=rng0.randrange(-20,21)
        x=[rng0.randrange(-50,51)for _ in range(4)]
        z=[rng0.randrange(-50,51)for _ in range(4)]
        d=integer_geom_encode(t,x,z)
        require(integer_geom_decode(d)==(t,x,z),'integer geometric inverse')
        require(apply_M(d,neighbors)==integer_geom_encode(25*t,[13*a+6*b for a,b in zip(x,z)],
                                                        [6*a+19*b for a,b in zip(x,z)]),
                'integer geometric matrix intertwiner')
        y,r=forward(d,neighbors)
        back,_=decode_mod8(y,[a%8 for a in d],neighbors,N,h)
        require(integer_geom_decode(back)==(t,x,z),'geometry-to-quantizer-to-recovery composite')
        integer_geometry_records.append({'t':t,'x':x,'z':z,'d':d,'y':y,'witness_mod8':[a%8 for a in d]})
    # Exact 5-dimensional geometric-plus-scalar intertwiner in Q(phi).
    rng=Random(20261107)
    geometric_cases=0
    for _ in range(20):
        t=(F(rng.randrange(-8,9),3),F(rng.randrange(-3,4),7))
        w=[(F(rng.randrange(-8,9),5),F(rng.randrange(-3,4),11))for _ in range(4)]
        signal=[padd(t,pscale(F(1,2),psum(pmul(v[a],w[a])for a in range(4))))
                for v in vertices]
        left=[padd((1,0),pscale(q/F(25),padd(pscale(13,padd(signal[i],(-1,0))),
                psum(padd(signal[j],(-1,0))for j in neighbors[i]))))for i in range(120)]
        tnew=padd((1,0),pscale(q,padd(t,(-1,0))))
        wnew=[pscale(q/F(25),pmul((13,6),x))for x in w]
        right=[padd(tnew,pscale(F(1,2),psum(pmul(v[a],wnew[a])for a in range(4))))for v in vertices]
        require(left==right,'geometric-scalar intertwiner')
        require(pscale(F(1,120),psum(signal))==t,'scalar recovery by mean')
        recovered=[pscale(F(1,60),psum(pmul(v[a],x)for v,x in zip(vertices,signal)))for a in range(4)]
        require(recovered==w,'four-coordinate recovery by Gram inverse')
        geometric_cases+=1
    scalar_values=list(range(-128,129))+[-10**100-17,10**100+17]
    for t in scalar_values:
        y,r=forward([t]*120,neighbors)
        ty,rho=divmod(Q*t,S)
        require(y==[ty]*120 and r==[25*rho]*120,'scalar diagonal specialization')
        require((S*ty+rho)//Q==t,'scalar full record inverse')
    samples=[[0]*120,[1]*120,[-1]*120,list(range(120)),
             [(-1)**i*(10**100+i)for i in range(120)]]
    samples += [[rng.randrange(-100000,100001)for _ in range(120)]for _ in range(15)]
    # Include actual collisions of T to rule out accidental assumptions of injectivity.
    require(forward([0]*120,neighbors)[0]==forward([1]*120,neighbors)[0], 'T collision')
    max_nodes=0;nodes=0
    for d in samples:
        y,r=forward(d,neighbors);s=[x%8 for x in d]
        a=decode_full(y,r,N,h);c,stats=decode_mod8(y,s,neighbors,N,h)
        require(a==c==d,'full versus mod8 witness equivalence')
        require(forward(c,neighbors)[1]==r,'mod8 to full remainder conversion')
        require([x%8 for x in a]==s,'full remainder to mod8 conversion')
        lhs=q*F(sum(d),120)
        require(lhs==F(sum(y),120)+F(sum(r),120*D),'exact mean ledger')
        max_nodes=max(max_nodes,stats['nodes']);nodes+=stats['nodes']
    # Quantization does not commute with the ordinary mean floor.
    z=[100]+[0]*119;y,_=forward(z,neighbors)
    mean_counterexample={'input_rule':'d=100*e_0','mean_of_T':str(F(sum(y),120)),
                         'floor_q_mean':(Q*sum(z))//(S*120)}
    require(F(sum(y),120)!=(Q*sum(z))//(S*120),'mean noncommutation control')
    history_count=0;history_steps=0
    # Exact error decomposition for gamma in Q; the written proof is for any real gamma.
    for seed in [samples[2],samples[3],samples[5]]:
        states=[seed];rem=[];witnesses=[]
        # q and a different rational parameter validate the coefficient-defect term.
        ideal={gamma:[F(x)for x in seed]for gamma in (q,F(1,3))}
        error={gamma:[F(0)for _ in seed]for gamma in ideal}
        Bpow=list(seed);ledger=[0]*120
        for k in range(6):
            d=states[-1];y,r=forward(d,neighbors)
            states.append(y);rem.append(r);witnesses.append([x%8 for x in d])
            for gamma in ideal:
                ideal[gamma]=[gamma*x/25 for x in apply_M(ideal[gamma],neighbors)]
                error[gamma]=[gamma*a/25+(gamma-q)*b/25+F(ri,D)
                              for a,b,ri in zip(apply_M(error[gamma],neighbors),apply_M(d,neighbors),r)]
                require(ideal[gamma]==[F(yi)+ei for yi,ei in zip(y,error[gamma])], 'ideal quantized defect identity')
            Bpow=[Q*x for x in apply_M(Bpow,neighbors)]
            ledger=[Q*x+D**k*ri for x,ri in zip(apply_M(ledger,neighbors),r)]
            require(Bpow==[D**(k+1)*yi+li for yi,li in zip(y,ledger)],'matrix telescoping identity')
        x=states[-1]
        for k in range(5,-1,-1):
            x,_=decode_mod8(x,witnesses[k],neighbors,N,h)
            require(x==states[k],'finite history reconstruction')
            history_steps+=1
        history_count+=1
    # Difference between preserving a coded object and asserting all real objects are coded.
    messages=[b'',b'\0',b'\0\0',b'OMEGA',b' {"scale":"3/7","frame":"fixed"} ']
    for message in messages:
        integer=int.from_bytes(b'\x01'+message,'big')
        d=[integer]+[0]*119;y,r=forward(d,neighbors)
        decoded=decode_full(y,r,N,h)[0]
        raw=decoded.to_bytes((decoded.bit_length()+7)//8,'big')
        require(raw[0]==1 and raw[1:]==message,'finite record adapter')
    return {
      'schema':'OMEGA_EXPLICIT_JOIN_CHECK/1','status':'PASS','date':'2026-11-07',
      'graph_sha256':GRAPH_SHA256,'source_sha256':sha256(Path(__file__).read_bytes()).hexdigest(),
      'exact_arithmetic':'integers, Fraction, and integer/rational pairs in Q(phi)',
      'inverse_certificate':{'h':h,'infinity_norm_numerator':norm,'C_inverse_infinity_norm':str(kappa)},
      'positive_definiteness':{'matrix':'M-9I','positive_leading_minors':len(pivots),
                              'minor_list_sha256':sha256(canonical(pivots)).hexdigest()},
      'geometry':{'AC_equals_6phiC_entries':coord_checks,'centered_columns':4,'gram_entries':16,
                  'gram_identity':'C^T C=30 I_4','joint_geometric_scalar_cases':geometric_cases,
                  'integer_basis_Gram_values':{'C0_C0':72,'C0_C1':-24,'C1_C1':48},
                  'integer_geometric_compositions_checked':len(integer_geometry_records),
                  'integer_geometric_mixer':[[13,6],[6,19]],
                  'integer_geometric_example':integer_geometry_records[0]},
      'scalar_diagonal_cases':len(scalar_values),'witness_conversion_cases':len(samples),
      'decoder_search_nodes_total':nodes,'decoder_search_nodes_max':max_nodes,
      'mean_noncommutation_counterexample':mean_counterexample,
      'history_cases':history_count,'history_steps_recovered':history_steps,
      'rational_parameter_cross_checks':[str(q),'1/3'],
      'ideal_vs_rational_counterexample':ideal_floor_counterexample,
      'finite_message_adapter_cases':len(messages),
      'full_geometric_proof_in_lean':False,'exact_F_max_computed':False,
      'changed_operator':False,'new_historical_priority_claim':False,
      'source_scope':'Reuses provided graph construction; fresh executions by this assistant, not independent third-party review.'
    }

if __name__=='__main__':
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output',type=Path,default=Path(__file__).with_name('verification_results.json'))
    args=parser.parse_args()
    result=run();data=json.dumps(result,indent=2,sort_keys=True)+'\n'
    args.output.write_text(data);print(data,end='')