"""Exact reconstruction of the Fock-gap manuscript's Appendix A.

Usage: python openai_laughlin_certificate.py PATH_TO_ROWS OUTPUT_JSON
Only the fixed rational certificate is checked; no analytic gap proof is checked.
"""
from fractions import Fraction as F
from functools import cache
from itertools import combinations
from math import comb, factorial as fac
from pathlib import Path
import hashlib, json, sys

if not __debug__:
    raise RuntimeError('Assertions must be enabled; do not use Python -O or -OO.')

rows_path, output_path = map(Path, sys.argv[1:])
rows=[]
for line in rows_path.read_text().splitlines():
    if not line.strip(): continue
    parts=line.split()
    if not line[0].isspace():
        t,ell=map(int,parts[:2]); rows.append((t,ell,{})); parts=parts[2:]
    for item in parts:
        labels,val=item.split(':'); p,j=map(int,labels)
        assert (p,j) not in rows[-1][2]
        rows[-1][2][p,j]=int(val)
assert len(rows)==7
P=10**7
eta=F(sum(ell**2 for t,ell,a in rows),P**2)
assert eta==F(93527408868499,10**14)

@cache
def U(c,z,T,p):
    if not (0<=z<=T and 0<=p<=T): return 0
    return sum((-1)**(z-h)*comb(z,h)*comb(T-z,p-h)*c**(p-h)
               for h in range(max(0,p+z-T),min(p,z)+1))

@cache
def V(r,D,T,p,j,k):
    if r>j+k: return 0
    return U(1,r,j+k,j)*U(1,D-r,T-r,p)

def two(exp):
    assert isinstance(exp,int)
    return F(2)**exp

floors=[]; margins=[]
for z in [2,*range(4,16)]:
    ez=F(0)
    for t,ell,a in rows:
        for T in range(max(z,t),16):
            x=sum(val*U(2,z,T,p) for (p,j),val in a.items() if p+j==T)
            ez+=F(2**z*comb(T,z),3**T*2**t*comb(T,t))*x*(x+2*ell*U(2,z,T,t))/P**2
    margin=F(3,2)*(3*z-1)*F(-1,2)**z-ez
    assert margin>0
    floors.append((10**6*margin).__floor__()); margins.append(str(margin))
assert floors==[1298313,211317,249,250,249,250,250,250,249,250,250,6470,250]

def mul(a,b):
    return [[sum((x*y for x,y in zip(row,col)),F(0)) for col in zip(*b)] for row in a]

expected=['0','0','00','00','+00','000','+000','+000','++000','+0000','+++000','++0000','++++000','+++0000','+++++000','++++0000','++++++000','+++++0000','+++++++000','++++++0000','++++++++000','+++++++0000','+++++++++000']
blocks=[]
for D in range(1,24):
    rr=list(range(1,D+1,2)); d=len(rr); u=[F(1,fac(r)*fac(D-r)) for r in rr]
    Y=[[F(0) for s in rr] for r in rr]
    for t,ell,a in rows:
        for (p,j),apj in a.items():
            for (q,l),aql in a.items():
                i=p+j-t; k=q+l-t; T=p+j+k
                assert i>=0 and k>=0
                if T<D: continue
                vv=[V(r,D,T,p,j,k) for r in rr]
                ww=[V(s,D,T,q,l,i) for s in rr]
                factor=F(-apj*aql*fac(i)*fac(k)*fac(t),P**2*fac(T-D))*two(1-2*T+p+q-t)
                for ri,r in enumerate(rr):
                    if vv[ri]:
                        for si,s in enumerate(rr):
                            if ww[si]: Y[ri][si]+=factor*vv[ri]*ww[si]*two((r+s)//2)
    assert Y==[list(x) for x in zip(*Y)]
    Z=[[F(0) for s in rr] for r in rr]
    for A in combinations(range(D+2),4):
        if sum(A)!=D+1: continue
        L=[]
        for r in rr:
            val=0
            for x,y in combinations(A,2):
                j,k=[v for v in A if v not in (x,y)]
                seq=[x,y,j,k]; inversions=sum(seq[ii]>seq[jj] for ii in range(4) for jj in range(ii+1,4))
                p=x+y-1
                val+=(-1)**inversions*(x-y)*fac(p)*fac(j)*fac(k)*V(r,D,D,p,j,k)
            L.append(two((1-2*D+r)//2)*val)
        f=1
        for x in A: f*=fac(x)
        for ri in range(d):
            for si in range(d): Z[ri][si]+=L[ri]*L[si]/f
    eps=F(3,10**6)
    B=[[(F(2*(r==1))+eps)*u[ri]*(ri==si)-u[ri]*Y[ri][si]*u[si] for si,s in enumerate(rr)] for ri,r in enumerate(rr)]
    M=mul(mul(Z,B),Z)
    assert M==[list(x) for x in zip(*M)]
    J=[[F(i==j) for j in range(d)] for i in range(d)]
    coeff=[]
    for b in range(1,d+1):
        X=mul(M,J); c=sum(X[i][i] for i in range(d))/b; coeff.append(c)
        J=[[c*(i==j)-X[i][j] for j in range(d)] for i in range(d)]
    assert all(c>=0 for c in coeff)
    signs=''.join('+' if c>0 else '0' for c in coeff)
    assert signs==expected[D-1],(D,signs,expected[D-1])
    blocks.append({'D':D,'dimension':d,'coefficient_signs':signs,'coefficients':[str(c) for c in coeff]})
    print('PASS block',D,signs,flush=True)

allowance=sum((D+1)**2//4 for D in range(1,24)); assert allowance==1222
tail=F(3,2)*F(1,2**17)*(900*F(4,3)+208*F(4,9)+12*F(20,27))
assert tail==F(61,4096)
gamma=1-eta-tail-allowance*F(3,10**6)
assert gamma==F(4616733319001,10**14)>F(1,25)
report={'success':True,'assertions_enabled':__debug__,'rows_sha256':hashlib.sha256(rows_path.read_bytes()).hexdigest(),'eta':str(eta),'three_body_margin_floors':floors,'three_body_margins':margins,'four_body_blocks':blocks,'allowance_factor':allowance,'tail_limit':str(tail),'gamma_star':str(gamma),'scope':'Reconstructed Appendix A rational arithmetic from supplied rows: 13 positive three-body margins, 23 four-body characteristic-polynomial sign patterns, and the final scalar margin. No finite-flux transfer, local stability argument, or Lean check.'}
output_path.write_text(json.dumps(report,indent=2)+'\n')
print('PASS all rational certificate checks',flush=True)