"""Exact verification of the new general-label upper estimates.

Only the Python standard library is used. No optimizer or floating-point
number enters an asserted inequality. Decimal is used solely to display
an already certified rational upper bound and to propose rational
brackets; each bracket is checked by exact polynomial substitution.
"""
from fractions import Fraction as F
from decimal import Decimal, localcontext
from pathlib import Path
import json

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

def dec(x):
    return Decimal(x.numerator)/Decimal(x.denominator)

def display(x,digits=15):
    with localcontext() as c:
        c.prec=digits+10
        return str(dec(x))

def count_criterion(r,A,mean,var,delta):
    """r*(A*r-mean-sqrt(var*delta)) > delta, exactly."""
    assert r>0 and var>=0 and delta>0
    residual=A*r-mean-delta/r
    assert residual>0
    assert residual*residual>var*delta

def cleanup_constants(ell,j):
    k=2**(ell-2);S=3*k-2*ell+4;p=k//2;n=p-j+1
    assert ell>=5 and k//4<j<=p
    eta=1 if j==p else 2;tau=2*j-p
    A=F(j*j)+(eta-F(5,2))*j+tau
    cs=[1,k]+[2**(ell-t+1)-2 for t in range(3,ell)]+[1]
    ls=[-2*(j-1),eta]+[1 if 2**(t-2)>k-2*j else 2 for t in range(3,ell)]+[0]
    assert sum(cs)==S
    total=sum(c*l for c,l in zip(cs,ls))
    mean=F(total,S)
    var=F(sum(c*l*l for c,l in zip(cs,ls)))-F(total*total,S)
    chi=F(4*(S-k),k*S)+F(2*(j-1),p*n)
    assert chi==F(2,n)-F(4,S)
    return k,S,p,n,A,mean,var,chi

def verify_cleanup_row(row):
    ell=row['ell'];j=row['q'];theta=F(row['theta']);delta=F(row['delta'])
    k,S,p,n,A,mean,var,chi=cleanup_constants(ell,j)
    B=F(2*S,n)/theta
    r=F(2,S)-chi/(4*B)-B*delta
    assert theta>0 and r>2*theta/S
    count_criterion(r,A,mean,var,delta)
    G=F(1,2)-F(1,2*S)-delta/2
    return G,delta

def verify_uniform():
    # The only nonuniform checks in the proof.
    for ell in [5,6,7]:
        k=2**(ell-2);S=3*k-2*ell+4;delta=F(1,324*k)
        if k==8:
            theta=F(1,2);r=F(11,144);A=F(14);mean=F(5,9);var=F(52)
        elif k==16:
            theta=F(2,3);r=F(37,1080)
            A=F(9*k*k,64)+F(13*k,16)+F(5,2)
            mean=2-F(3*k,4*S)-F(12,S);var=F(9*k*k,16)+4*S-32
        else:
            theta=F(4,9);r=F(17*k*k-S*S,9*S*k*k)
            A=F(9*k*k,64)+F(13*k,16)+F(5,2)
            mean=2-F(3*k,4*S)-F(12,S);var=F(9*k*k,16)+4*S-32
        B=F(16*S,k)/theta
        assert r==F(2,S)-F(16,k)/(4*B)-B*delta
        assert r>2*theta/S
        count_criterion(r,A,mean,var,delta)
    # Uniform k>=64 proof: r>=8/(27k), sqrt(V delta)<=k/144.
    assert 64**2-36*64-768>0
    assert 2*64-36>0  # derivative of the last polynomial remains positive
    assert F(5*64,144)-F(95,54)==F(25,54)
    assert F(25,54)>F(1,48)
    assert F(2,324)/F(8,27)==F(1,48)
    # Room r>2h is equivalent to S<3k, true for every ell>=5.
    print('PASS uniform gamma improvement 1/(648k), every ell>=5')

def path_constants(ell,b=None):
    k=2**(ell-2);S=3*k-2*ell+4
    if b is None:b=k//2
    T=S-k+b;L=b+1;a=1-F(2,L)
    A=a*a/2-F(1,T);B=F(2,T)-F(4,L);d0=F(k-b,S*T)
    return k,S,b,T,L,a,A,B,d0

def path_bracket(ell):
    k,S,b,T,L,a,A,B,d0=path_constants(ell)
    assert ell>=10 and A>0
    with localcontext() as ctx:
        ctx.prec=100
        u=(-dec(B)+(dec(B*B+4*A*d0)).sqrt())/(2*dec(A))
        denominator=10**70
        low=F(int(u*denominator),denominator)
        high=low+F(1,denominator)
    def poly(t):return A*t*t+B*t-d0
    def d(t):return (1-t)**2/T-F(1,S)
    def h(t):return a*a*t*t/2-F(4,L)*t
    assert 0<low<high<F(1,10)
    assert poly(low)<0<poly(high)
    assert 2*A*low+B>0
    assert a*a*low-F(4,L)>0
    delta=min(d(high),h(low))
    assert delta>0
    G=F(1,2)-F(1,2*S)-delta/2
    return G,delta,low,high

def verify_path_uniform():
    # For ell>=10, k>=256, S>=14k/5; the latter holds at ell10
    # and persists when k doubles and ell increases by one.
    k=256;ell=10;S=3*k-2*ell+4
    assert S>=F(14,5)*k
    assert F(1,5)*k>=2*ell-4
    assert F(1,5)*k>2
    # Positive intersection: h changes sign before 1/12, while d(1/12)>0.
    L=129
    assert (L-2)**2>96*L
    assert 2*L-100>0
    assert F(1,3)-4*F(1,12)+2*F(1,12)**2==F(1,72)>0
    assert F(1,2)*F(127,129)**2-F(1,384)>0
    # Path-length optimum: these constants prove the derivative is negative
    # at every positive intersection for every real b in [k/2,k).
    assert F(23,28)>F(81,100)
    assert F(9,20)*F(127,129)**2>F(63,160)
    assert F(42,5)*F(12,256)==F(63,160)
    print('PASS path radical applicability and optimal b=k/2 for every ell>=10')

def pa(a,b):
    return [(a[i] if i<len(a) else F(0))+(b[i] if i<len(b) else F(0))
            for i in range(max(len(a),len(b)))]

def pm(a,b):
    out=[F(0)]*(len(a)+len(b)-1)
    for i,x in enumerate(a):
        for j,y in enumerate(b):out[i+j]+=x*y
    return out

def ps(c,a):return [c*x for x in a]
def pd(a):return [i*a[i] for i in range(1,len(a))]
def pe(a,x):
    out=F(0)
    for c in reversed(a):out=out*x+c
    return out

def budget_polynomials(ell,b=None):
    k,S,b,T,L,a,A,B,d0=path_constants(ell,b)
    q0=F(S,k*(S-k));cc=a*a/2;dd=F(4,L)
    m11=[q0,F(1,b)-F(1,k)]
    m12=[F(0),-F(1,b)]
    m22=[cc,F(1,b)-cc]
    v1=[F(1,S-k)];v2=[dd/2,-dd/2]
    det=pa(pm(m11,m22),ps(-1,pm(m12,m12)))
    num=pa(pa(pm(m22,pm(v1,v1)),ps(-2,pm(m12,pm(v1,v2)))),pm(m11,pm(v2,v2)))
    deriv=pa(pm(num,pd(det)),ps(-1,pm(pd(num),det)))
    while deriv and not deriv[-1]:deriv.pop()
    assert len(deriv)<=5
    assert pe(det,F(0))>0 and pe(det,F(1))==F(1,b*(S-k))
    return k,S,b,F(k,S*(S-k)),det,num,deriv

def budget_certificate(ell):
    k,S,b,constant,det,num,deriv=budget_polynomials(ell)
    # Concavity proves the sign bracket encloses the global maximum.
    assert pe(deriv,F(0))>0>pe(deriv,F(1))
    lo=F(0);hi=F(1)
    for _ in range(150):
        md=(lo+hi)/2
        if pe(deriv,md)>0:lo=md
        else:hi=md
    assert pe(deriv,lo)>0>pe(deriv,hi)
    lam=(lo+hi)/2
    assert pe(det,lam)>0
    delta=constant-pe(num,lam)/pe(det,lam)
    assert delta>0
    G=F(1,2)-F(1,2*S)-delta/2
    # Positive-semidefinite completed-square certificate reconstructed here.
    q0=F(S,k*(S-k));L=b+1;a=1-F(2,L);c=a*a/2;d=F(4,L)
    xx=q0+lam*(F(1,b)-F(1,k));xm=-lam/b;mm=lam/b+(1-lam)*c
    v1=F(1,S-k);v2=(1-lam)*d/2
    determinant=xx*mm-xm*xm
    assert xx>0 and determinant>0
    assert constant-delta==(mm*v1*v1-2*xm*v1*v2+xx*v2*v2)/determinant
    return G,delta,{'ell':ell,'b':b,'lambda_low':str(lo),'lambda_high':str(hi),
                   'lambda_certificate':str(lam),'delta':str(delta),
                   'stationary_quartic':[str(c) for c in deriv]}

def main():
    verify_uniform();verify_path_uniform()
    rows=json.loads(HERE.joinpath('cleanup-certificates.json').read_text())
    clean={}
    for row in rows:
        G,delta=verify_cleanup_row(row);clean[row['ell']]=(G,delta)
    print('PASS',len(rows),'exact optimized cleanup certificates')
    out=[];budget_rows=[]
    print('ell | paper | earlier-general | new upper | source')
    for ell in range(5,18):
        k=2**(ell-2);S=3*k-2*ell+4
        paper=F(1,2)-F(1,2*S)
        earlier=paper-F(1,256*k*k)
        G=paper-F(1,648*k);source='uniform cleanup'
        if ell in clean and clean[ell][0]<G:G=clean[ell][0];source='optimized cleanup'
        path=None
        if ell>=10:
            PG,PD,lo,hi=path_bracket(ell)
            if PG<G:G=PG;source='optimized paths'
            path={'u_low':str(lo),'u_high':str(hi),'delta_lower':str(PD)}
            BG,BD,brow=budget_certificate(ell);budget_rows.append(brow)
            assert BG<PG
            if BG<G:G=BG;source='full-budget paths'
        assert G<earlier<paper
        print(ell,display(paper),display(earlier),display(G),source,sep=' | ')
        out.append({'ell':ell,'upper_rational':str(G),'upper_display':display(G,30),
                    'source':source,'path_root':path})
    HERE.joinpath('verified-general-table.json').write_text(json.dumps(out,indent=2)+'\n')
    HERE.joinpath('path-budget-certificates.json').write_text(json.dumps(budget_rows,indent=2)+'\n')
    print('PASS',len(budget_rows),'full-budget 2D certificates and globally maximizing quartic brackets')
    print('PASS all displayed new general bounds improve the prior general estimate')

if __name__=='__main__':main()