"""Independent exact check of the stored certificate (T32): reads star-bound-arb/I_lower_bounds.json (the lower bounds
I^lo(t) at the cell right ends) and star-bound-arb/<problem>.json (cells and transfers), and re-verifies in exact
rational arithmetic that (i) the cells tile [lo, hi], (ii) every charging inequality of Lemma 3 of Theorem M5 holds with
C^ = Phi(-t0) - I^lo(t1), where Phi(-t0) is bounded above by a fresh Arb evaluation, with the stored transfers in [-1,1].
It does not recompute I^lo (that is certify_star_bound_arb.py's job).  Run: .venv/bin/python <this file> [data directory]
The data directory defaults to star-bound-arb/ (run of 2026-10-07); pass star-bound-arb-owen to check the run of
2026-10-10 of certify_star_bounds_arb_owen.py, whose arc bound is enclosed in Arb.

Revision of 2026-10-07 (T33): where a problem's coefficient at degree 3 is the constant F itself (stored as the upper
bound F_UP = 0.14015825233145217 >= F), the mixed degree-3 rows are now checked against the Arb lower bound F_LO <= F
(Proposition 6.9 of the paper), so that they prove C <= F and not only C <= F_UP.  The pure cubic
type is still skipped there: its bound is F by definition.

Revision of 2026-10-14: the problem S-switched (thresholds [0, 35/32], coefficients (c_3, 0.2627) with c_3 stored as the
lower bound 0.139822124008, used for the median-switched sweep of maximum degree 4) is checked as well when its data
exist; its lower bounds I^lo are stored in I_lower_bounds-S-switched.json.  Its pure cubic type is skipped because
cubic_cap <= a_3: at thresholds t >= 0 the pure cubic star is bounded by c_3 (manuscript, prop:star-cubic)."""
import json, os, sys
from fractions import Fraction
from flint import arb, ctx
ctx.prec = 100
D = os.path.join(os.path.dirname(os.path.abspath(__file__)), sys.argv[1] if len(sys.argv) > 1 else "star-bound-arb")
def Phi_up(x):                       # rational upper bound for Phi(x): upper end of the Arb ball, rounded up to 1e-15
    u = ((1 + (arb(x.numerator) / x.denominator / arb(2).sqrt()).erf()) / 2).upper()
    return Fraction(float(u.mid()) + float(u.rad())) + Fraction(1, 10**15)
F_UP = Fraction("0.14015825233145217")
F_LO = Fraction("0.1401582523314521636939405140172")
ok_all = True
names = ["M", "M5-P", "M5-R", "M5-S"] + (["S-switched"] if os.path.exists(os.path.join(D, "S-switched.json")) else [])
for name in names:
    T = json.load(open(os.path.join(D, "I_lower_bounds-%s.json" % name if name == "S-switched" else "I_lower_bounds.json")))
    G = T["grid"]; types = [(d, tuple(ty)) for d, ty in T["types"]]
    P = json.load(open(os.path.join(D, name + ".json")))
    a = [Fraction(x) for x in P["problem"]["a"]]; a_is_F = a[0] == F_UP
    if a_is_F: a[0] = F_LO                 # mixed degree-3 rows must give C <= F; pure cubic row: bound F, skipped
    degs = P["problem"]["degs"]; cap = Fraction(P["problem"]["cubic_cap"])
    cells = sorted(P["cells"]); lo = cells[0][0]; hi = cells[-1][1]
    assert all(c[1] == cells[i + 1][0] for i, c in enumerate(cells[:-1])), "gap"
    worst = None
    for k0, k1, _, tr in cells:
        assert all(abs(x) <= 1 for x in tr) and k0 < k1
        top = Phi_up(Fraction(-k0, G)); I = T["thresholds"][str(k1)]; trf = [Fraction(x) for x in tr]
        for (d, ty), Il in zip(types, I):
            if d not in degs or any(ty[j - 3] > 0 for j in (3, 4, 5) if j not in degs): continue
            if (d, ty) == (3, (3, 0, 0)) and (a_is_F or cap <= a[0]): continue
            C = top - Fraction(Il)
            row = [0, 0, 0]
            if d == 3: row[0] = -ty[1]; row[1] = -ty[2]
            if d == 4: row[0] = ty[0]; row[2] = -ty[2]
            if d == 5: row[1] = ty[0]; row[2] = ty[1]
            s = a[d - 3] - C - sum(c * x for c, x in zip(row, trf))
            if worst is None or s < worst: worst = s
    ok = worst > 0
    ok_all &= ok
    print("%-10s %4d cells tile [%s, %s]; minimum exact slack %.3e  %s" % (name, len(cells), lo / G, hi / G, float(worst), "OK" if ok else "FAIL"))
sys.exit(0 if ok_all else 1)