from decimal import Decimal as D, getcontext
import json, time
getcontext().prec = 80
PI = D('3.1415926535897932384626433832795028841971693993751058209749445923078164062862090')
def Phi(x):
    if x < 0: x = -x
    e2x = (2*x).exp()
    s = D(0)
    for n in range(1, 40):
        t = PI*n*n*e2x
        if t > D(300): break
        s += (2*PI*PI*n**4*(D(9)*x/2).exp() - 3*PI*n*n*(D(5)*x/2).exp())*(-t).exp()
    return 2*s
def W(v, h=D('0.004'), Y=D('3.5')):
    n = int(Y/h); n += n % 2
    tot = D(0); half = v/2
    for i in range(n+1):
        y = h*i
        f = (2*y)**2 * Phi(half+y) * Phi(half-y)
        w = 1 if i in (0,n) else (4 if i%2==1 else 2)
        tot += w*f
    return tot*h/3
t0=time.time()
us=[D(i)/D(10) for i in range(2,92,2)]   # u=0.2..9.0
hvals=[]
for u in us:
    hvals.append(W(u.sqrt()).ln())
print(" u      h(u)=ln W(sqrt u)      d2h (log-convexity: must be >=0 for Gaussian mixture)")
out=[]
first_neg=None
for i in range(1,len(us)-1):
    d2=hvals[i+1]-2*hvals[i]+hvals[i-1]
    out.append((float(us[i]),float(hvals[i]),float(d2)))
    mark=" <-- NEGATIVE: log-convexity FAILS" if d2<0 else ""
    if d2<0 and first_neg is None: first_neg=float(us[i])
    print(f"{float(us[i]):4.1f}  {float(hvals[i]):+12.5f}   {float(d2):+.6e}{mark}")
json.dump(out,open('gaussmix.json','w'))
print(f"\nfirst negative second difference at u = {first_neg}  (v = {first_neg**0.5:.3f})" if first_neg else "\nno violation found")
print(f"[{time.time()-t0:.0f}s]")