"""Independent audit: exact binary grid endpoints, Arb cell convolutions,
and a separate transcription of the Lemma T envelope. All writes are local.
Run single-threaded; no certificate module is imported.
"""
import os
for key in ['OPENBLAS_NUM_THREADS','OMP_NUM_THREADS','MKL_NUM_THREADS','NUMEXPR_NUM_THREADS']:
    os.environ[key]='1'
import json, math, time
from pathlib import Path
from flint import arb, arb_poly, ctx
ctx.prec=256
ctx.threads=1
HERE=Path(__file__).resolve().parent
FILES=['cert_qrhB3_K14_350_51046.json', 'cert_qrhB3p_K14_51046_159011874857757.json']
GRIDS=[json.loads((HERE/'frozen/second'/f).read_text()) for f in FILES]
WMAX=max(39,max(row[2]-1 for g in GRIDS for row in g['rows']))
H=200
J={(0,w):arb(1) for w in range(WMAX+1)}
def save_ball(x):
    # str includes a radius, not a rounded claim about a directed endpoint.
    return x.str(24)
def cells():
    print('Arb cell convolution H=',H,'WMAX=',WMAX,flush=True)
    c=arb_poly([arb(1)/(H+i) for i in range((WMAX-1)*H+1)])
    p=arb_poly([1])
    for m in range(1,WMAX):
        # All discarded coefficients have index too large for any later sum.
        nmax=(WMAX-m)*H
        p=(p*c).truncate(nmax+1)
        coeff=p.coeffs()
        prefix=arb(0)
        for n in range(nmax+1):
            prefix+=coeff[n]
            if n%H==0:
                w=m+n//H
                J[m,w]=prefix
        print('cells',m,'J_m(39)=',save_ball(J.get((m,39),arb(0))),flush=True)

def D(l): return arb(1)/(10*l*l)+arb(4)/(15*l**3)
U=2*(-arb.const_euler()).exp()*arb('0.6601618158')
CT=1+arb(1)/36260
def cofactors(l,M,crude=False):
    z=l.exp(); eps=2*D(l)+1/z; W=M-1
    mu=[arb(1)]
    rep=W*(2*eps).exp()/(z-1)
    assert eps>1/(z-1)
    for j in range(1,W):
        nj=sum((math.comb(j,i)*eps**i*J[j-i,W] for i in range(j+1)),arb(0))
        extra=arb(0) if j<2 else (rep if crude else mu[j-2]/(z-1))
        mu.append(nj/math.factorial(j)+extra)
    aa=[arb(0),arb(0)]+[a*(l/z**a+l/((a-1)*z**(a-1))+1/((a-1)**2*z**(a-1))) for a in range(2,M)]
    return mu,aa,eps

def bad(L,M,K,presieve=False,crude=False,delta=arb('0.05')):
    l=L/M; z=l.exp()
    assert l>arb(10372).log()
    if presieve: assert l>arb(8000000000).log()
    mu,aa,eps=cofactors(l,M,crude)
    cp=arb(9)/32*(1+arb(480)/(840*l)) if presieve else CT
    X=sum((cp*mu[k-1] for k in range(K+1,M)),arb(0))/((1-delta)*L)
    Y=sum((sum((aa[a]*mu[k-a] for a in range(2,k+1)),arb(0)) for k in range(K+1,M)),arb(0))/((1-delta)*L)
    Y+=(-delta*L).exp()+(4*(-L).exp() if presieve else 0)
    return X,Y,eps

def sieve(left,right,M,lam):
    # arb(float) represents the exact binary value, consistently throughout.
    L1,L2,lam=arb(left),arb(right),arb(lam)
    assert L1>=300 and L2>=L1 and lam>0 and lam<=100 and M>=2
    h1=L1/8+lam.log()-3*L1.log(); h2=L2/8+lam.log()-3*L2.log()
    assert h1>=arb(25).log() and h1<L1 and h2<L2 and h2>=h1
    l=L1/M; z=l.exp()
    s1=(h1-arb(1155).log())/l; s2=(h2-arb(1155).log())/(L2/M)
    assert s1>2 and s2<=4 and s2>=s1
    assert z>694 and z>2 and L1*(M-1)/M>arb(10372).log()
    g=arb.const_euler(); pi=arb.pi(); e=arb(1).exp()
    F=2*g.exp()*(s1-1).log()/s1-arb('0.00477')*154*e**2*(-s1).exp()
    v1=1-arb('1.002')*M/(z-1); v2=1-1/(2*l*l)
    assert F>0 and v1>0 and v2>0
    main=U*M/L2**2*v1*v2*F
    a,b=arb('0.298'),arb('4.358'); const=(2*pi*e).log()
    nq=(h2-const)/pi+a*h2+b
    def lp(x): return x+(1+(-x).exp()).log()
    t=lp(L2/8); it=(-t).exp()
    Z1=7*nq+((h2-const)*t+t*t/2)/pi+(a*h2+b)*(1-it)+a*(1-(t+1)*it)+(4/pi)*(h2+t-const+1)
    Z2=arb('6.12')*nq+2/pi*(h2-const+1)+a*h2+a/2+b
    T1=arb(16)/7*(1+(-L1/16).exp()/(L1/2))
    T2=(L2+1)*(arb('0.65')*(-L1/4).exp()+(((L2+1)/arb(2).log()).log()+1)*(-7*L1/16).exp())
    T3=arb('0.5')+arb('1.5')*(-7*L1/16).exp()
    T4=(1+(-L1/16).exp())**(arb(7)/8)*(Z1+2*(a*h2+a/2+b)*(-lp(L1/16)).exp()+2*a/e)
    T5=(arb('4.67')*Z2+L2+1+arb('3.1'))*(-7*L1/16).exp()
    ce=T1+T2+T3+T4+T5
    rem=arb('0.45')*lam/L1**3*(ce/L1*(1+arb(32)/(7*L1))+3*(-3*L1/8).exp()+(2*h2/arb(3).log()+arb('1.3841')*L2/L2.log())*(-7*L1/8).exp())
    return main,rem,(-L1).exp(),s1,s2,F

def js():
    L0=arb('2736.86'); L3=arb('32.7').exp(); delta=arb('0.01')
    print('U_low',save_ball(U),'L3',save_ball(L3),flush=True)
    out={}
    mu,aa,eps=cofactors(L0/40,40,True)
    for K in [18,19,20]:
        phi=sum((CT*mu[k-1]+sum((aa[a]*mu[k-a] for a in range(2,k+1)),arb(0)) for k in range(K+1,40)),arb(0))
        first=L3*phi/((1-delta)*5*U)
        second=(L0**2*(-delta*L0).exp()+L3**2*(-L0).exp())/(5*U)
        out[K]={key:save_ball(val) for key,val in [('Phi',phi),('first',first),('second',second),('lhs',first+second)]}
        print('S3',K,json.dumps(out[K]),flush=True)
        if K==19: assert first+second<arb('0.226172')
    out['eps']=save_ball(eps)
    out['inv_z_minus_1']=save_ball(1/((L0/40).exp()-1))
    return out

def grids():
    out=[]
    for p,g in enumerate(GRIDS):
        rows=g['rows']; assert len(rows)==[501,449][p]
        assert rows[0][0]==g['Lstart'] and rows[-1][1]==g['Lend']
        assert all(a[1]==b[0] for a,b in zip(rows,rows[1:]))
        worst=min(range(len(rows)),key=lambda i:rows[i][4])
        selected=sorted({0,1,len(rows)//4,len(rows)//2,3*len(rows)//4,len(rows)-1,worst})
        summary={'nrows':len(rows),'original_worst_index':worst,'samples':[]}
        minrel=None; minlz=None
        for idx,row in enumerate(rows):
            left,right,M,lam,*_=row
            main,rem,one,s1,s2,F=sieve(left,right,M,lam)
            X,Y,eps=bad(arb(left),M,14,bool(p))
            margin=main-rem-one-X-Y
            assert margin>0
            if p: assert main>X # necessary direction for Pi cancellation
            rel=margin/main
            if minrel is None or rel.lower()<minrel.lower(): minrel=rel; minidx=idx
            lz=arb(left)/M
            if minlz is None or lz<minlz: minlz=lz
            if idx in selected:
                sample={'index':idx,'row':row,'main':save_ball(main),'rem':save_ball(rem),'X':save_ball(X),'Y':save_ball(Y),'margin':save_ball(margin),'relative_margin':save_ball(rel),'s1':save_ball(s1),'s2':save_ball(s2),'log_z1':save_ball(lz)}
                summary['samples'].append(sample)
                print('GRID',p,json.dumps(sample),flush=True)
        summary.update({'min_relative_margin':save_ball(minrel),'min_index':minidx,'min_log_z1':save_ball(minlz),'min_M':min(r[2] for r in rows),'max_M':max(r[2] for r in rows)})
        print('GRID SUMMARY',p,json.dumps({k:v for k,v in summary.items() if k!='samples'}),flush=True)
        out.append(summary)
    return out

if __name__=='__main__':
    start=time.time(); cells(); result={'JS':js(),'grids':grids()}
    result['elapsed_seconds']=time.time()-start
    (HERE/'audit_v6.json').write_text(json.dumps(result,indent=2)+'\n')
    print('ALL ASSERTIONS PASSED',result['elapsed_seconds'],flush=True)