"""Independent GRH9 presieved middle pieces and full M=11 sieve cover.
256-bit Arb; exact rational grid; no author module imported.
"""
import os
for k in ['OPENBLAS_NUM_THREADS','OMP_NUM_THREADS','MKL_NUM_THREADS','NUMEXPR_NUM_THREADS']:os.environ[k]='1'
from pathlib import Path
from fractions import Fraction
import math,json
from flint import arb,arb_poly,ctx
ctx.prec=256;ctx.threads=1
OUT=Path(__file__).resolve().parent/'resume1'
def primes(n):
    s=bytearray(b'\x01')*(n+1);s[:2]=b'\x00\x00'
    for p in range(2,math.isqrt(n)+1):
        if s[p]:s[p*p::p]=b'\x00'*len(s[p*p::p])
    return [p for p in range(2,n+1) if s[p]]
PR=primes(250000)
def D(l):return 1/(10*l*l)+arb(4)/(15*l**3)
def up(v):return arb(v.upper())
BHI=arb('.26149721284765');BLO=arb('.26149721284764')
CT=1+arb(1)/36260
def middle(w,Ltop,Llow,delta):
    ll=arb(Llow) if Llow else arb('4e18').log();lt=arb(Ltop);d=arb(delta)
    rempr=[p for p in PR if p>w];wp=rempr[0]
    km=int(float(lt/arb(wp).log()));assert km*arb(wp).log()<lt<(km+1)*arb(wp).log()
    Sp=up(lt.log()+BHI+D(lt)-sum((arb(1)/p for p in PR if p<=w),arb(0)))
    P=[arb(0),Sp];A=[arb(0),arb(0)]
    for k in range(2,km+1):
        cutoff=250000 if k<=12 else 5000
        pp=[p for p in rempr if p<=cutoff]
        val=[arb(p)**(-k) for p in pp]
        cut=arb(cutoff)
        P.append(up(sum(val,arb(0))+cut**(1-k)/(k-1)))
        A.append(up(k*(sum((arb(p).log()*v for p,v in zip(pp,val)),arb(0))+cut**(1-k)*(cut.log()/(k-1)+1/arb(k-1)**2))))
    R=[arb(1)]
    for n in range(1,km+1):R.append(up(sum((P[k]*R[n-k] for k in range(1,n+1)),arb(0))/n))
    phi=sum((CT*R[k-1]+sum((A[a]*R[k-a] for a in range(2,k+1)),arb(0)) for k in range(10,km+1)),arb(0))
    lhs=ll*(-d*ll).exp()+phi/(1-d)
    q=math.prod(p for p in PR if 3<=p<=w);totient=math.prod(p-1 for p in PR if 3<=p<=w)
    g=arb(1)
    for p in PR:
        if 3<=p<=w:g*=arb(p-2)/(p-1)
    assert q*q<ll.exp()-2 and d*ll>=1 and 0<d<1 and ll<=lt
    rhs=g*(1-4*(-ll).exp())-arb('.64')*totient*ll**2*(-ll/2).exp()*arb('1.01')-ll*(-ll).exp()-arb('1e-12')
    result={'w':w,'q':q,'phi':totient,'kmax':km,'g':g.str(28),'S':Sp.str(28),'lhs_majorant':lhs.str(28),'rhs':rhs.str(28),'q2_over_Yminus2':(q*q/(ll.exp()-2)).str(25)}
    assert lhs<rhs
    if w==11:
        uncond=135*(arb(1)/480-1/(840*(arb('4e18')-2).log()))*(1-arb('1e-18'))-arb('1e-18')
        assert lhs<uncond;result['unconditional_rhs']=uncond.str(28)
    print('MIDDLE',json.dumps(result),flush=True)
    return result

def sieve():
    W,H=10,400
    cp=arb_poly([(1+arb(1)/(H+i)).log() for i in range(9*H+1)])
    poly=arb_poly([1]);J=[arb(1)]
    for m in range(1,W):
        poly=(poly*cp).truncate((W-m)*H+1);J.append(sum(poly.coeffs(),arb(0)))
    table=[p for p in PR if p<=10372]
    acc=arb(0);errors={}
    for p in table:
        acc+=arb(1)/p;errors[p]=up(acc-arb(p).log().log()-BLO)
    run=up(D(arb(10372).log()));suffix={}
    for p in reversed(table):
        run=max(run,errors[p],arb(0),key=lambda x:float(x))
        # Every value is an exact upper endpoint, so order is exact at 256 bits.
        assert run>=errors[p];suffix[p]=run
    U=2*(-arb.const_euler()).exp()*arb('.6601617536');good=arb('.146')*U
    rows=[];first=Fraction(73);top=Fraction(7700000)
    while first<top:
        last=min(top,first*Fraction(21,20));l1=arb(first.numerator)/first.denominator;l2=arb(last.numerator)/last.denominator
        l=l1/11;z=l.exp()
        if z>=10372:sz=up(D(l));pk=None
        else:
            pk=max(p for p in table if arb(p)<=arb(z.lower()))
            sz=suffix[pk]
            assert arb(pk)<=z and all(p<=pk or arb(p)>z for p in table)
        eps=sz+D(l)+1/z;mu=[arb(1)]
        for j in range(1,W):
            nj=sum((math.comb(j,i)*eps**i*J[j-i] for i in range(j+1)),arb(0))
            mu.append(nj/math.factorial(j)+(mu[j-2]/(z-1) if j>=2 else 0))
        def A(a):return a*(l/z**a+l/((a-1)*z**(a-1))+1/((a-1)**2*z**(a-1)))
        phi=CT*mu[9]+sum((A(a)*mu[10-a] for a in range(2,11)),arb(0))
        vals=[]
        for ds in ['.30','.25','.20','.15','.10','.05','.02','.01']:
            d=arb(ds)
            if d*l1<2:continue
            v=l2*phi/((1-d)*good)+(l1**2*(-d*l1).exp()+5*l2**2*(-l1).exp())/good
            vals.append((up(v),ds))
        best,ds=min(vals,key=lambda pair:float(pair[0]));assert best<1
        rows.append({'L1':str(first),'L2':str(last),'delta':ds,'pk':pk,'S':sz.str(25),'eps':eps.str(25),'Phi':phi.str(25),'upper':best.str(25)})
        first=last
    assert len(rows)==238 and rows[0]['L1']=='73' and rows[-1]['L2']=='7700000'
    assert all(a['L2']==b['L1'] for a,b in zip(rows,rows[1:]))
    idx=max(range(len(rows)),key=lambda i:float(arb(rows[i]['upper'])))
    assert arb(rows[idx]['upper'])<arb('.000265')
    assert all(arb(row['upper'])<arb('.000265') for row in rows)
    threshold=arb('15.85').exp();assert threshold<7700000
    print('SIEVE FIRST',json.dumps(rows[0]),flush=True)
    print('SIEVE WORST',idx,json.dumps(rows[idx]),flush=True)
    print('THRESHOLD',threshold.str(28),'overlap',(arb(7700000)-threshold).str(28),flush=True)
    return {'rows':rows,'worst_index':idx,'threshold':threshold.str(28),'J9':J[9].str(28)}

def qG():
    L=arb('4e18').log();x=arb('4e18');a=x**(arb(1)/4)
    val=arb('.165')+arb('12.683')/L+arb('254.980')/L**2+arb('2607.854')/L**3+arb('11605.056')/L**4+arb('1.314')*L/a+arb('.092')*L.log()/a+arb('60.883')/a+arb('8.250')*L.log()/(a*L)+arb('939.260')/(a*L)-arb('237.934')/(x.sqrt()*L)
    print('BS qG at 4e18',val.str(28));assert val<arb('.640')
    return val.str(28)
if __name__=='__main__':
    result={'qG':qG(),'middle':[middle(11,'62',None,'.20'),middle(19,'70','62','.15'),middle(23,'73','65','.15')],'sieve':sieve()}
    (OUT/'audit_grh9.json').write_text(json.dumps(result,indent=2)+'\n')