# Copied 2026-10-08 from an earlier exploratory run.
# (sha256 667e6d8cbf4226cd62170e2c411b44dc904519d904f6e874d0b93024b4c6c6ee); unchanged below this header. Run from the repository root:
#   .venv/bin/python experiments/reproducibility/F/arb_F_global.py 0.1401582524   (one or more target upper bounds)
# Arb (python-flint) interval certificate for F = sup_tau f_*(tau), ideal star r = 2 sqrt2/3 (s = 1/3 exactly).
# Rigorous ingredients:
#  (i)  P(tau) := Pr(Y0>tau, all Yi>tau) is nonincreasing in tau (event shrinks), and Phi(-tau) is decreasing,
#       so for tau in [a,b]:  f_*(tau) = Phi(-tau) - P(tau) <= Phi(-a) - P(b).
#  (ii) P(tau) = Phi(-zl) - (3/pi) int_0^{pi/3} (sig/sqrt q) exp(-tau^2/(2q)) Phi(-(sqrt q/sig)(zl - r tau/q)) dth,
#       sig = s cos th, q = sig^2 + r^2, zl = tau (tau<0) or tau/r (tau>=0)  [re-derived by the referee,
#       and checked against an independent nested inclusion-exclusion formula to 30 digits, fstar_mp.py].
#       The integral is evaluated with acb.integral (rigorous for analytic integrands; sqrt uses analytic flag).
#  (iii) tails: tau <= -T: f <= Pr(min_i Yi <= tau) <= 3 Phi(tau) <= 3 Phi(-T); tau >= T: f <= Phi(-tau) <= Phi(-T).
import sys, time
from flint import arb, acb, ctx
ctx.prec = 80

R = 2*arb(2).sqrt()/3
S = arb(1)/3
PI = arb.pi()

def Phi_arb(x):  # x arb
    return (-x/arb(2).sqrt()).erfc()/2

def P_arb(tau):
    tau = arb(tau)
    zl = tau if tau < 0 else tau/R
    if not (tau < 0) and not (tau >= 0):
        raise ValueError("tau straddles 0")
    ta = acb(tau); zla = acb(zl); Ra = acb(R); Sa = acb(S)
    def integrand(th, analytic):
        sig = Sa*th.cos()
        q = sig*sig + Ra*Ra
        sq = q.sqrt(analytic=analytic)
        arg = -(sq/sig)*(zla - Ra*ta/q)
        Phi_c = (-arg/acb(2).sqrt()).erfc()/2
        return (sig/sq)*(-(ta*ta)/(2*q)).exp()*Phi_c
    I = acb.integral(integrand, 0, PI/3, rel_tol=arb(2)**-60, abs_tol=arb(2)**-70)
    return Phi_arb(-zl) - 3*I.real/PI

def f_arb(tau):
    tau = arb(tau)
    return Phi_arb(-tau) - P_arb(tau)

def certify(target, T=2, maxdepth=60):
    target = arb(target)
    # tails
    tail_lo = 3*Phi_arb(arb(-T)); tail_hi = Phi_arb(arb(-T))
    assert (tail_lo < target) and (tail_hi < target), "tails not certified"
    cache = {}
    def Pc(t):
        key = t.mid().str(40, radius=False)
        if key not in cache: cache[key] = P_arb(t)
        return cache[key]
    # boxes on [-T,0] and [0,T] separately (zl formula changes at 0); endpoints are exact dyadics/rationals
    stack = [(arb(-T), arb(0), 0), (arb(0), arb(T), 0)]
    nbox = 0; worst = arb(0); nfail = 0
    while stack:
        a, b, d = stack.pop()
        ub = Phi_arb(-a) - Pc(b)
        if ub < target:
            nbox += 1
            if ub.upper() > worst.upper(): worst = ub
            continue
        if d >= maxdepth:
            nfail += 1
            print("FAIL box", a.str(20), b.str(20), "ub", ub.str(15)); continue
        m = (a+b)/2
        stack.append((a, m, d+1)); stack.append((m, b, d+1))
    return nbox, worst, nfail, len(cache)

if __name__ == "__main__":
    ts = arb("-0.068137131081692573475")
    v = f_arb(ts)
    print("f_*(tau*) enclosure:", v.str(25), " -> rigorous lower bound F >=", v.lower().str(22))
    print("f_*(0) enclosure:", f_arb(0).str(25), " c3 =", (3*(arb(5)/6).acos()/(4*PI)).str(25))
    for tgt in sys.argv[1:]:
        t0 = time.time()
        nbox, worst, nfail, nev = certify(tgt)
        print("target %s: accepted boxes %d, failures %d, P-evaluations %d, max box upper bound %s, time %.1fs"
              % (tgt, nbox, nfail, nev, worst.upper().str(15), time.time()-t0), flush=True)