# Exact symbolic part of the proof: the shell-2 main terms (deficit + shell energy) and P1.
# All cancellations at orders s^0, s^1 are verified exactly here; the remaining polynomials are
# exported (with exact rational/log2/sqrt2 coefficients) for interval evaluation.
import sympy as sp, pickle
s,v,mu,kap = sp.symbols('s v mu kap')
l = sp.log(2); r2=sp.sqrt(2)
def Qser(sig,n):   # z e^z E1(z) = sum_{k<n} (-1)^k k! z^-k + (-1)^n n! kap z^-n, kap in [0,1]
    return sum((-1)**k*sp.factorial(k)*sig**k for k in range(n)) + (-1)**n*sp.factorial(n)*kap*sig**n
NQ=8
# ---- P1 = t^2 (U J(1/2) - 2 + 2/t),  U J(1/2) = Q(t/2) + 1, 1/z = 2s
P1 = sp.expand((Qser(2*s,NQ) - 1 + 2*s)/s**2)
assert sp.simplify(P1.subs(s,0)-8)==0
# ---- shell 2 
w   = -l/2*s + v*s**2
lam = w + mu*s**2                    # lam = log(1+w)
D   = 1 + l*s + 2*lam*s              # L2 = t D / 2
c   = r2*(2+l/2)
# H = U E1(L2) = sqrt2 * Q(L2) / ((1+w) D),  sigma = 1/L2 = 2s/D
# Num = sqrt2 * sum a_k (2s)^k D^(n-k) - (sqrt2 - c s)(1+w) D^(n+1) ;  H - sqrt2 + c s = Num / ((1+w) D^(n+1))
a = [(-1)**k*sp.factorial(k) for k in range(NQ)] + [(-1)**NQ*sp.factorial(NQ)*kap]
Num = sp.expand(r2*sum(a[k]*(2*s)**k*D**(NQ-k) for k in range(NQ+1)) - (r2 - c*s)*(1+w)*D**(NQ+1))
p = sp.Poly(Num, s)
c0 = sp.simplify(p.coeff_monomial(1)); c1 = sp.simplify(p.coeff_monomial(s))
print("Num coeff s^0:", c0, "  s^1:", c1)
assert c0==0 and c1==0
NumS2 = sp.expand(sp.cancel(Num/s**2))
Den  = sp.expand((1+w)*D**(NQ+1))
# Mcomb = -P2 - sqrt2 v ,  P2 = NumS2/Den
# leading constants:
P2_0 = sp.simplify((NumS2/Den).subs(s,0))
print("P2(0) =", P2_0, "=", sp.N(P2_0))
print("Mcomb(0) =", sp.simplify(-P2_0 - r2*v), "=", sp.N(-P2_0-r2*v))
print("P1(0) =", P1.subs(s,0))
print("c0 = P1(0)+Mcomb(0) =", sp.nsimplify(sp.simplify(8 - P2_0 - r2*v)), sp.N(8-P2_0-r2*v))
pickle.dump({'P1':sp.srepr(P1),'NumS2':sp.srepr(NumS2),'Den':sp.srepr(Den),'NQ':NQ}, open('sym.pkl','wb'))
# ---- Lemma (shell 2) polynomials
vv=sp.symbols('vv')
ww = -l/2*s + vv*s**2
Pup = sp.expand((1+ww)**2*(1+l*s) + 2*ww*s*(1+ww) - 1)       # need Pup >= exp(-1/s)/2
Plo = sp.expand(1 - (1+ww)**2*(1+l*s+2*ww*s))                # need Plo >= delta*(1+l s+2 w s)
for name,P in [('Pup',Pup),('Plo',Plo)]:
    q=sp.Poly(P,s); print(name, [sp.nsimplify(q.coeff_monomial(s**k)) for k in range(3)])
pickle.dump({'Pup':sp.srepr(sp.expand(sp.cancel(Pup/s**2))),'Plo':sp.srepr(sp.expand(sp.cancel(Plo/s**2)))},open('lemma.pkl','wb'))
print("asymptotic v* =", sp.N((l+sp.Rational(3,4)*l**2)/2))