# Copied 2026-10-08 from an earlier exploratory run.
# (sha256 e971720c8a81f4e3c8db55dde3d13f09623f2ccba96f33cf6043d6c436c6dfdd); unchanged below this header. Run from the repository root:
#   .venv/bin/python experiments/reproducibility/F/arb_F.py
# Rigorous (Arb / python-flint) enclosure of F = sup_tau f_*(tau), independent implementation.
# f_*(tau) for tau <= 0 (ideal star, r = 2sqrt2/3, s = 1/3):
#   f(tau) = (3/pi) int_0^{pi/3} A exp(-tau^2/(2q)) Phi(-a tau) dtheta,
#   sigma = s cos(theta), q = sigma^2 + r^2, A = sigma/sqrt(q), a = (q - r)/(sigma sqrt(q)).
#   f'(tau) = (3/pi) int A exp(-tau^2/(2q)) [ -(tau/q) Phi(-a tau) - a phi(a tau) ] dtheta.
# Unimodality (proved analytically in the write-up): f' = phi (-1 + A0 + 3 B0) with A0, B0 nonincreasing,
# so f'(t1) > 0 > f'(t2) localizes the maximiser in [t1, t2].
from flint import arb, acb, ctx
import mpmath as mp

ctx.prec = 200

R_ = arb(2) * arb(2).sqrt() / 3
S_ = arb(1) / 3
PI = arb.pi()
SQ2 = arb(2).sqrt()
SQ2PI = (2 * PI).sqrt()


def Phi_acb(x):
    return (-(x) / acb(SQ2)).erfc() / 2


def phi_acb(x):
    return (-(x * x) / 2).exp() / acb(SQ2PI)


def make_integrand(tau, deriv=False):
    tau = acb(tau)
    r = acb(R_)
    s = acb(S_)

    def f(th, analytic):
        sig = s * th.cos()
        q = sig * sig + r * r
        sq = q.sqrt(analytic=analytic)
        A = sig / sq
        a = (q - r) / (sig * sq)
        E = (-(tau * tau) / (2 * q)).exp()
        if not deriv:
            return A * E * Phi_acb(-a * tau)
        return A * E * (-(tau / q) * Phi_acb(-a * tau) - a * phi_acb(a * tau))
    return f


def f_val(tau, deriv=False):
    I = acb.integral(make_integrand(tau, deriv), 0, acb(PI) / 3)
    return (I * 3 / acb(PI)).real


# locate tau* approximately with mpmath
mp.mp.dps = 40
r = 2 * mp.sqrt(2) / 3
s = mp.mpf(1) / 3


def fp_mp(t):
    def g(th):
        sig = s * mp.cos(th)
        q = sig ** 2 + r ** 2
        A = sig / mp.sqrt(q)
        a = (q - r) / (sig * mp.sqrt(q))
        E = mp.e ** (-t ** 2 / (2 * q))
        return A * E * (-(t / q) * mp.ncdf(-a * t) - a * mp.npdf(a * t))
    return 3 / mp.pi * mp.quad(g, [0, mp.pi / 3])


def f_mp(t):
    def g(th):
        sig = s * mp.cos(th)
        q = sig ** 2 + r ** 2
        A = sig / mp.sqrt(q)
        a = (q - r) / (sig * mp.sqrt(q))
        E = mp.e ** (-t ** 2 / (2 * q))
        return A * E * mp.ncdf(-a * t)
    return 3 / mp.pi * mp.quad(g, [0, mp.pi / 3])


tstar = mp.findroot(fp_mp, mp.mpf('-0.068'))
print('mpmath tau* =', mp.nstr(tstar, 25), ' F(mp) =', mp.nstr(f_mp(tstar), 25))

c = arb(mp.nstr(tstar, 30))
h = arb('1e-15')
t1 = c - h
t2 = c + h
fp1 = f_val(t1, True)
fp2 = f_val(t2, True)
print("f'(t1) =", fp1)
print("f'(t2) =", fp2)
assert fp1 > 0 and fp2 < 0
fc = f_val(c)
ball = arb(c, h)  # [t1, t2]
fpb = f_val(ball, True)
print('f(c) =', fc)
print("f'([t1,t2]) =", fpb)
bound = abs(fpb).upper() * h
lower = fc.lower()
upper = fc.upper() + bound
print('F lower =', lower.str(30, radius=False))
print('F upper =', upper.str(30, radius=False))
print('width =', (upper - lower).str(5, radius=False))
# value at 0 equals c3 = 3 arccos(5/6)/(4 pi)
f0 = f_val(arb(0))
c3 = 3 * arb(5 / arb(6)).acos() / (4 * PI)
print('f(0) =', f0, '\nc3   =', c3, '\ndiff =', f0 - c3)
fp0 = f_val(arb('-1e-30'), True)
print("f'(0-) ~", fp0)