# Lemma (shell 2): for all real xi >= 1e9,  xi_2 = sqrt(2 xi)(1 - (log2/2)/t + v/t^2) with 0.47 <= v <= 0.55.
import sympy as sp, pickle, math
from flint import arb as A, ctx
ctx.prec = 160
d = pickle.load(open('lemma.pkl','rb'))
s, vv = sp.symbols('s vv')
L = A(2).log()
def load(name, vval):
    e = eval(d[name], sp.__dict__ | {'Symbol': sp.Symbol})
    P = sp.Poly(sp.expand(e.subs(vv, vval)), s)
    cs = []
    for (k,), c in P.terms():
        c = sp.nsimplify(c)
        # c is polynomial in log(2) with rational coeffs
        q = sp.Poly(c, sp.log(2))
        val = A(0)
        for (j,), r in q.terms(): val += A(int(sp.Rational(r).p))/A(int(sp.Rational(r).q))*L**j
        cs.append((k, val))
    return cs
def ev(cs, S):
    return sum((c*S**k for k, c in cs), A(0))
s0 = 1/math.log(1e9)
up = load('Pup', sp.Rational(55,100)); lo = load('Plo', sp.Rational(47,100))
N = 4000
mup = mlo = 1e9
for i in range(N):
    a = s0*i/N; b = s0*(i+1)/N
    S = A(a).union(A(b)); Sb = A(b)
    w = -L/2*S + A("0.47")*S*S
    # upper: Pup/s^2 >= e^{-1/s}/(2 s^2)
    rhs_u = (-1/Sb).exp()/(2*Sb*Sb)
    mu_ = ev(up, S) - rhs_u
    # lower: Plo/s^2 >= delta (1 + l s + 2 w s)/s^2, delta = (sqrt2 (1+w) e^{-t/2} + e^{-t}/2)/2 <= (sqrt2 e^{-1/(2s)} + e^{-1/s}/2)/2
    delta_over_s2 = (A(2).sqrt()*(-1/(2*Sb)).exp() + (-1/Sb).exp()/2)/2/(Sb*Sb)
    ml = ev(lo, S) - delta_over_s2*(1 + L*S + 2*w*S)
    mup = min(mup, float(mu_.lower())); mlo = min(mlo, float(ml.lower()))
    assert mu_ > 0 and ml > 0, (i, mu_, ml)
print("Lemma shell2 verified on", N, "cells of s in (0, 1/log 1e9]; min margins:", mup, mlo)