"""Directed-ball transcription of the lemma statements of an earlier version of the paper, written before the certificate code was read.

No certificate code is imported.  Exact rational/decimal inputs, Arb at 256 bits.
The four initial T witnesses use our own parameters and the grid described in prose.
"""
import json
from pathlib import Path
from flint import arb, ctx

ctx.prec = 256
ctx.threads = 1
A = arb
PI = A.pi()
GAMMA = A.const_euler()
ULOW = 2 * (-GAMMA).exp() * A("0.6601618158")

def dump(v):
    if isinstance(v, arb):
        return v.str(35)
    if isinstance(v, dict):
        return {k: dump(x) for k, x in v.items()}
    if isinstance(v, (tuple, list)):
        return [dump(x) for x in v]
    return v

def save(name, v):
    Path(name).write_text(json.dumps(dump(v), indent=2) + "\n")

def primes(n):
    b = bytearray(b"\x01") * (n+1)
    b[:2] = b"\0\0"
    for p in range(2, int(n**0.5)+1):
        if b[p]:
            b[p*p::p] = b"\0" * ((n-p*p)//p+1)
    return [p for p in range(2,n+1) if b[p]]

def mobius(n):
    v, k = n, 0
    p = 2
    while p*p <= v:
        if v % p == 0:
            v //= p
            k += 1
            if v % p == 0:
                return 0
        p += 1
    return (-1)**(k + (v > 1))

def constants():
    # B = gamma + sum_{m>=2} mu(m) log(zeta(m))/m.
    # log(zeta(s)) < 3*2^{-s}, s>=2, gives the explicit tails below.
    cache = {}
    def lz(s):
        if s not in cache:
            cache[s] = A(s).zeta().log()
        return cache[s]
    b = GAMMA + sum((mobius(m)*lz(m)/m for m in range(2,101)), A(0))
    b += A(0, 3*A(2)**(-100))
    # log C_tw = -sum_{k>=2} (2^k-2) P_odd(k)/k.
    ctlog = A(0)
    for k in range(2,101):
        j = max(3, (340+k-1)//k)
        ps = sum((mobius(m)*lz(k*m)/m for m in range(1,j+1)), A(0))
        ps += A(0, 3*A(2)**(-k*(j+1))/(1-A(2)**(-k)))
        ps -= A(2)**(-k)
        ctlog -= (A(2)**k-2)*ps/k
    # P_odd(k) <= 4*3^{-k}; drop 1/k in the outer tail.
    tail = 12*(A(2)/3)**101
    ctlog = ctlog.union(ctlog-tail)
    twin = ctlog.exp()
    assert b < A("0.26149721284765")
    assert twin > A("0.6601618158")
    return {"B": b, "twin_constant": twin, "U_low": ULOW}

def correction(ell):
    return 1/(10*ell**2)+4/(15*ell**3)

def c2(t, ps, cutoff):
    finite = sum((-(1-t/p).log()-t/p for p in ps if p >= 13), A(0))
    tail = t*t/(2*(1-t/cutoff)*cutoff)
    return finite.union(finite+tail)

def c_lemma(L, K, t, ps, cutoff):
    sp = L.log()+A("0.26149721284765")+correction(L)-sum((1/A(p) for p in [2,3,5,7,11]), A(0))
    g0 = 135*(1/A(480)-1/(840*(A(4)*10**18-2).log()))
    rhs = (g0*(1-A("1e-18"))-A("1e-18")).log()
    lhs = L.log()+t*sp+c2(t,ps,cutoff)-(K+1)*t.log()
    return {"L":L,"K":K,"t":t,"Sprime":sp,"G0":g0,"lhs":lhs,"rhs":rhs,"slack":rhs-lhs}

def r_correction(ell, M, t):
    ez = (-ell).exp()
    return t*(correction((M-1)*ell)+correction(ell)+ez)+t*t*ez/(2*(1-t*ez)*(1-ez))

def s_lemma(K,t):
    L=A("32.7").exp()
    ell=L/40
    logbad=(1+A("1.2762")/ell).log()-ell.log()+t*A(39).log()+r_correction(ell,40,t)-K*t.log()
    loggood=(5*ULOW/L**2-(-L).exp()).log()
    # A uniform proof freezes decreasing factors at the LEFT endpoint.
    L0=A("2736.86")
    q0=r_correction(L0/40,40,t)
    logratio_uniform=(40*(L+A("51.048"))/(5*ULOW-L0**2*(-L0).exp())).log()+t*A(39).log()-K*t.log()+q0
    return {"K":K,"t":t,"L":L,"logbad":logbad,"loggood":loggood,"slack":loggood-logbad,"left_correction":q0,"uniform_log_slack":-logratio_uniform}

def cE(L1,L2,lam):
    # EXPLICIT_EX A.4, taking the positive increasing terms at the top.
    w=L1/2
    logH=L2/8-3*L2.log()+lam.log()
    a,b=A("0.298"),A("4.358")
    c=(2*PI).log()+1
    tau=L2/8+(1+(-L2/8).exp()).log()
    invT=1/(1+(L1/16).exp())
    n1=(logH-c)/PI+a*logH+b
    z1=7*n1+((logH-c)*tau+tau*tau/2)/PI+a*logH+b+a+4*(logH+tau-c+1)/PI
    z1+=2*(a*logH+a/2+b)*invT+2*a/A(1).exp()
    z2=A("6.12")*n1+2*(logH-c+1)/PI+a*logH+a/2+b
    e8=(-w/8).exp()
    eb=(-7*w/8).exp()
    logth=w+(1+e8).log()
    ramp=A(16)/7*(1+e8/w)
    pp=logth*(A("0.65")*(-w/2).exp()+((logth/A(2).log()).log()+1)*eb)
    delta=A(1)/2+A("1.5")*eb
    zeros=(1+e8)**(A(7)/8)*z1
    fixed=A("4.67")*z2*eb
    trivial=(logth+A("3.1"))*eb
    return {"total":ramp+pp+delta+zeros+fixed+trivial,"logHmax":logH,"ramp":ramp,"pp":pp,"delta":delta,"zeros":zeros,"fixed":fixed,"trivial":trivial}

def t_lemma(L1,L2,M,lam,K,t):
    logH1=L1/8-3*L1.log()+lam.log()
    logH2=L2/8-3*L2.log()+lam.log()
    s1=M*(logH1-A(1155).log())/L1
    s2=M*(logH2-A(1155).log())/L2
    ell=L1/M
    ez=(-ell).exp()
    z1=1/ez
    F=2*GAMMA.exp()*(s1-1).log()/s1-A("0.00477")*154*A(2).exp()*(-s1).exp()
    main=ULOW*M/L2**2*(1-A("1.002")*M*ez/(1-ez))*(1-1/(2*ell**2))*F
    ce=cE(L1,L2,lam)
    rem=A("0.45")*lam/L1**3*(ce["total"]/L1*(1+A(32)/(7*L1))+3*(-3*L1/8).exp()+(2*logH2/A(3).log()+A("1.3841")*L2/L2.log())*(-7*L1/8).exp())
    bad=A(0) if K >= M-1 else (1+A("1.2762")/ell)/ell*(-K*t.log()+t*A(M-1).log()+r_correction(ell,M,t)).exp()
    loss=(-L1).exp()
    net=main-rem-bad-loss
    assert s1 > 2 and s2 <= 4
    assert z1 >= 694 and t < z1 and logH1 > A(25).log()
    assert L1*(M-1)/M > A(10372).log()
    assert F > 0 and 1-A("1.002")*M*ez/(1-ez) > 0
    assert net > 0
    return {"L1":L1,"L2":L2,"M":M,"lam":lam,"K":K,"t":t,"s1":s1,"s2":s2,"log_z1":ell,"cE":ce,"main":main,"rem":rem,"bad":bad,"unit_loss":loss,"net":net,"relative_margin":net/main}

def blind():
    out={"provenance":"Computed from the paper's statements and the cited sources, before the certificate code was read."}
    out["constants"]=constants()
    ps=primes(1000000)
    out["C"]=[c_lemma(A("2736.86"),26,A(193)/50,ps,1000000),c_lemma(A(300),18,A(79)/20,ps,1000000),c_lemma(A(73),13,A(103)/25,ps,1000000)]
    out["C_K25_trials"]=[c_lemma(A("2736.86"),25,A(t),ps,1000000) for t in ["3.6","3.7","3.75","3.8","3.9"]]
    # Strong convexity on 1<t<13 certifies failure everywhere, not just on a grid.
    t=A("3.7")
    row=c_lemma(A("2736.86"),25,t,ps,1000000)
    deriv=row["Sprime"]+sum((1/(A(p)-t)-1/A(p) for p in ps if p>=13),A(0))-26/t
    deriv=deriv.union(deriv+t/(1000000*(1-t/1000000)))
    global_lower=row["lhs"]-deriv**2/(2*(A(26)/169))
    assert global_lower > row["rhs"]
    out["C_K25_global"]={"derivative_at_3.7":deriv,"strong_convexity":A(26)/169,"global_lhs_lower":global_lower,"rhs":row["rhs"]}
    out["S"]=[s_lemma(31,A(423)/50),s_lemma(32,A(873)/100)]
    out["T_own_witnesses"]=[]
    for j,M,lam,t in [(0,45,A(100),A(24)/5),(231,26,A(1),A(28)/5),(514,22,A(1)/10,A(59)/10)]:
        L1=300*(A(101)/100)**j
        L2=L1*A(101)/100
        out["T_own_witnesses"].append(t_lemma(L1,L2,M,lam,18,t))
    L1=51046*(A(21)/20)**448
    out["T_own_witnesses"].append(t_lemma(L1,A(159011874857757),20,A(1)/10,19,A(1)))
    out["Q_sqrt_coefficient"]={"low":2/A(2).log(),"prime_powers":A(4).log()*(1/A(300)+8/A(300)**2)+(2+A(12)/300)*(-A(300)/6).exp()}
    save("blind_results.json",out)
    print(json.dumps(dump(out),indent=2))

if __name__ == "__main__":
    blind()