"""Star-bound certificate with the arc bound enclosed in Arb (2026-10-10). See Proposition 7.5 of the paper. This is certify_star_bounds_arb.py with one change: the 3-star arc bound for the seven mixed 3-star types is no longer computed with libm exp/arccos and scipy.special.owens_t (lowered by 1e-11), but enclosed in ball arithmetic. Everything else (Lemmas K1 to K3, grids, tables, thresholds, charging, exact check) is unchanged. * Arc bound (Lemma K2, arc version). For thresholds a (doubles, exact), c = -a, rho_a = max(0, max a_i) and a double rho_b > rho_a, the value exp(-rho_a^2/2) - exp(-rho_b^2/2) - (1/pi) sum_i J_i, J_i = int_{r1}^{r2} arccos(c_i/rho) rho exp(-rho^2/2) d rho, r1 = max(rho_a, |c_i|), r2 = rho_b, is a lower bound for P_min(a) for any rho_b. J_i is evaluated by the closed form J_i = [-exp(-rho^2/2) arccos(c_i/rho)]_{r1}^{r2} + sign(c_i) int_{y1}^{y2} exp(-c_i^2 (1+x^2)/2) / (1+x^2) dx, y_j = sqrt(r_j^2 - c_i^2)/|c_i| (the last integral is 2 pi (T(|c_i|,y2) - T(|c_i|,y1)) with Owen's T). Everything is an Arb ball: exp, acos and sqrt are Arb functions, and the Owen integral is a rigorous Arb integration (acb.integral, Gauss-Legendre with certified error bounds) of the defining integrand between the exact midpoints of the balls y_j, widened by the radii of y_j (the integrand is at most 1). For |c_i| < 2^-30 the integral is enclosed instead by [atan y2 - atan y1 - (c_i^2/2)(y2 - y1), atan y2 - atan y1]. The lower end of the final ball, rounded down to a double, is used. * The lower bound used for a mixed type at a z-point is max(Lambda(a), arc), where Lambda is the simplex-truncation bound of gauss_lo (rigorous, as for all other types). The Arb arc bound is computed only at points where a floating-point estimate of the arc bound (the formulas of certify_star_bounds_arb.py, with scipy.special.owens_t) exceeds Lambda(a) + 1e-10; elsewhere Lambda(a) is used. The estimate and the floating-point choice of rho_b only decide which valid lower bound is used; no floating-point value of exp, arccos or Owen's T enters an inequality. * Phi, phi, the simplex tables, threshold rounding, the sums I^lo and the exact rational re-check of every cell are as in certify_star_bounds_arb.py (see its docstring). The linear program for the transfers is solved in floating point (scipy.optimize.linprog); the transfers are used only through the exact check. Outputs go to star-bound-arb-owen/ (cells, transfers, I^lo); the 2026-10-07 files in star-bound-arb/ are kept. Run: .venv/bin/python experiments/arb-certificates/certify_star_bounds_arb_owen.py [problem names] """ import os, sys, itertools, json, time from fractions import Fraction import numpy as np from scipy.special import owens_t from scipy.optimize import linprog import multiprocessing as mp from flint import arb, acb, ctx ctx.prec = 80 HERE = os.path.dirname(os.path.abspath(__file__)) OUT = os.path.join(HERE, "star-bound-arb-owen") DELTA = 1e-10 # Arb arc bound only where the float estimate beats Lambda by more than this NPROC = 14 EPS = 2.0 ** -52 G = 2 ** 17 # grid for Phi, phi: multiples of 2^-17 KMIN, KMAX = -8 * G, 13 * G DZ = 2 ** 14 # z-step 2^-14 = 4 grid steps STRIDE = G // DZ DX = 2.0 ** -19 # grid of the simplex tables XMAX = 12 NG = XMAX * 2 ** 19 MA = 3e-13 SQ2 = arb(2).sqrt(); SQ2PI = (2 * arb.pi()).sqrt() def dn(x): return np.nextafter(x, -np.inf) def up(x): return np.nextafter(x, np.inf) def arb_bounds(vals): """list of arb -> (lo, hi) double arrays enclosing them.""" n = len(vals) lo = np.empty(n); hi = np.empty(n) for i, v in enumerate(vals): m = float(v.mid()); r = float(v.rad()) lo[i] = m - r - 4 * EPS * abs(m) - 1e-300 hi[i] = m + r + 4 * EPS * abs(m) + 1e-300 return np.nextafter(lo, -np.inf), np.nextafter(hi, np.inf) # ---- Arb tables on the grids (built once; workers inherit them by fork) def build_tables(): t0 = time.time() xs = [arb(k) / G for k in range(KMIN, KMAX + 1)] global PHI_LO, PHI_HI, PDF_LO, PDF_HI PHI_LO, PHI_HI = arb_bounds([(1 + (x / SQ2).erf()) / 2 for x in xs]) PDF_LO, PDF_HI = arb_bounds([(-(x * x) / 2).exp() / SQ2PI for x in xs]) PDF_LO = np.maximum(PDF_LO, 0.0) print("grid tables: %.0fs" % (time.time() - t0), flush=True) # simplex tables global FLO xg = [arb(k) / 2 ** 19 for k in range(NG + 1)] pdf_lo, _ = arb_bounds([(-(x * x) / 2).exp() / SQ2PI for x in xg]) f2_lo, _ = arb_bounds([(x / SQ2).erf() for x in xg]) Flo = {2: np.maximum(f2_lo, 0.0)} for d in (3, 4, 5): kap = dn(np.sqrt(d / (d - 2))) k = np.arange(NG, dtype=np.float64) idx = np.minimum(np.floor(k * kap * (1 - 2.0 ** -50)).astype(np.int64), NG) cell = d * pdf_lo[1:] * Flo[d - 1][idx] * DX * (1 - 8 * 2.0 ** -53) ti = np.floor(np.maximum(cell, 0.0) * 2.0 ** 61).astype(np.int64) cs = np.concatenate([[0], np.cumsum(ti)]) assert cs[-1] < 2 ** 62 Flo[d] = np.maximum(dn(cs.astype(np.float64) / 2.0 ** 61 * (1 - 2.0 ** -50)), 0.0) FLO = Flo print("simplex tables: %.0fs" % (time.time() - t0), flush=True) def S_lo(d, a): """lower bound for S_d(a) = G_d(-a); a is an upper bound for the true threshold.""" h = -np.asarray(a, float) idx = np.clip(np.floor(h * 2.0 ** 19), 0, NG).astype(np.int64) return np.where(h <= 0, 0.0, FLO[d][idx]) def _grid(a): ac = np.clip(a, KMIN / G, KMAX / G - 1.0 / G) k = np.floor(ac * G) return k.astype(np.int64) - KMIN, ac - k / G def Phi_ub(a): a = np.asarray(a, float) i, dl = _grid(a) v = PHI_HI[i] + PDF_HI[i] * dl + 0.121 * dl * dl v = np.where(a < KMIN / G, PHI_HI[0], np.where(a >= KMAX / G, 1.0, v)) return np.minimum(up(v * (1 + 8 * EPS) + 1e-300), 1.0) def Phi_lb(a): a = np.asarray(a, float) i, dl = _grid(a) v = PHI_LO[i] + PDF_LO[i] * dl - 0.121 * dl * dl v = np.where(a < KMIN / G, 0.0, np.where(a >= KMAX / G, PHI_LO[-1], v)) return np.maximum(dn(v * (1 - 8 * EPS) - 1e-300), 0.0) def gauss_lo(a, equal): """Union bound and simplex truncation (Lemma 2 of Theorem M5), d = a.shape[1] in 3..5.""" d = a.shape[1] lb = np.maximum(0.0, dn(1 - d * Phi_ub(a[:, 0]) * (1 + 4 * EPS)) if equal else dn(1 - Phi_ub(a).sum(1) * (1 + 8 * EPS))) if equal: return np.maximum(lb, S_lo(d, a[:, 0])) Pu = Phi_ub(a) for j in range(d): lvl = a[:, j] Pl = Phi_lb(lvl) pen = np.sum(np.where(a > lvl[:, None], Pu - Pl[:, None], 0.0), axis=1) * (1 + 8 * EPS) lb = np.maximum(lb, dn(S_lo(d, lvl) - pen)) return lb def arc_estimate(a): """Floating-point estimate of the 3-star arc bound (libm / scipy.special.owens_t), and the outer radius rho_b. Used only to choose rho_b and the points where the Arb arc bound is evaluated (formulas of certify_star_bounds_arb.py).""" c = -a; n = len(c) neg = np.where(c < 0, -c, 0.0); ra = neg.max(1) def gsum(rho): with np.errstate(divide='ignore', invalid='ignore'): q = np.where(rho[:, None] > 0, c / rho[:, None], np.sign(c) * np.inf) return np.arccos(np.clip(q, -1, 1)).sum(1) lo = ra.copy(); hi = np.full(n, 12.0) ok = gsum(lo * (1 + 1e-12) + 1e-300) < np.pi for _ in range(44): mid = 0.5 * (lo + hi); below = gsum(mid) < np.pi lo = np.where(below, mid, lo); hi = np.where(below, hi, mid) rb = lo val = np.exp(-ra ** 2 / 2) - np.exp(-rb ** 2 / 2) for i in range(3): ci = c[:, i]; ac = np.abs(ci) r1 = np.maximum(ra, ac); r2 = rb act = r2 > r1 r1s = np.where(act, r1, 1.0); r2s = np.where(act, r2, 2.0); cs = np.where(ac > 0, ci, 1e-300); acs = np.abs(cs) q1 = np.clip(cs / r1s, -1, 1); q2 = np.clip(cs / r2s, -1, 1) y1 = np.sqrt(np.maximum(r1s ** 2 - acs ** 2, 0)) / acs; y2 = np.sqrt(np.maximum(r2s ** 2 - acs ** 2, 0)) / acs J = -np.exp(-r2s ** 2 / 2) * np.arccos(q2) + np.exp(-r1s ** 2 / 2) * np.arccos(q1) + np.sign(cs) * 2 * np.pi * (owens_t(acs, y2) - owens_t(acs, y1)) J0 = (np.pi / 2) * (np.exp(-r1s ** 2 / 2) - np.exp(-r2s ** 2 / 2)) J = np.where(ac > 0, J, J0) val = val - np.where(act, J, 0.0) / np.pi return np.where(ok, np.maximum(0.0, val), 0.0), rb ARC_PREC = 64 HSMALL = 2.0 ** -30 def owen_integral(h, y1, y2): """Arb ball containing int_{y1}^{y2} exp(-h^2 (1+x^2)/2) / (1+x^2) dx for an exact arb h > 0 and arb balls 0 <= y1 <= y2.""" if float(h.upper().mid()) < HSMALL: at = y2.atan() - y1.atan() lo = at - h * h / 2 * (y2 - y1) return lo.union(at) m1 = arb(y1.mid()); m2 = arb(y2.mid()) k = acb(h * h / 2) def f(x, analytic): u = 1 + x * x return (-(k * u)).exp() / u I = acb.integral(f, acb(m1), acb(m2), rel_tol=2.0 ** -55, abs_tol=2.0 ** -62).real return I + arb(0, y1.rad() + y2.rad()) def arc_arb(a, rb): """Arb ball containing the arc integral over [rho_a, rb] for thresholds a (three doubles), rho_a = max(0, max a_i). It is a lower bound for P_min(a) for every rb; returns None if rb <= rho_a.""" ra = max(0.0, float(max(a))) if not rb > ra: return None Ra = arb(ra); Rb = arb(rb); pi = arb.pi() val = (-(Ra * Ra) / 2).exp() - (-(Rb * Rb) / 2).exp() for ai in a: ci = -float(ai); h = abs(ci) r1 = max(ra, h); r2 = rb if not r2 > r1: continue # empty range: for rho < |c_i| with c_i > 0 the arc is empty R1 = arb(r1); R2 = arb(r2) e1 = (-(R1 * R1) / 2).exp(); e2 = (-(R2 * R2) / 2).exp() if ci == 0: J = pi / 2 * (e1 - e2) else: Cc = arb(ci); H = arb(h) y1 = (R1 * R1 - H * H).sqrt() / H if r1 > h else arb(0) y2 = (R2 * R2 - H * H).sqrt() / H J = -e2 * (Cc / R2).acos() + e1 * (Cc / R1).acos() + (1 if ci > 0 else -1) * owen_integral(H, y1, y2) val = val - J / pi return val def arc_arb_lo(a, rb): """double lower bound from arc_arb (0 if rb <= rho_a or the enclosure is not finite).""" old = ctx.prec; ctx.prec = ARC_PREC try: v = arc_arb(a, rb) finally: ctx.prec = old if v is None or not v.is_finite(): return 0.0, 0.0 lo, _ = arb_bounds([v]) return max(float(lo[0]), 0.0), float(v.rad()) # ---- the correlations (floats checked against Arb) and the vertex types def r_arb(i, j): return 1 / (((((arb(i) - 1).log() + (arb(j) - 1).log()) / 4)).cosh()) RA = {(i, j): r_arb(i, j) for i in (3, 4, 5) for j in (3, 4, 5)} RF = {k: float(v.mid()) for k, v in RA.items()} def types(d): return [t for t in itertools.product(range(d + 1), repeat=3) if sum(t) == d] ALL = [(d, t) for d in (3, 4, 5) for t in types(d)] def rs_of(d, t): return [RF[(d, 3)]] * t[0] + [RF[(d, 4)]] * t[1] + [RF[(d, 5)]] * t[2] def I_lo(d, ty, tk): """Lower bound for P(Y_0 > t, all Y_i > t), t = tk 2^-17, over stars with incident correlations >= r.""" rs = np.array(rs_of(d, ty)); s = np.sqrt(1 - rs ** 2) mode_arc = d == 3 and max(ty) < 3 step = STRIDE * (2 if mode_arc else 1) nz = (10 * DZ) // (2 if mode_arc else 1) k = tk + step * np.arange(nz + 1) z = k / G t = tk / G ph_lo = PHI_LO[k - KMIN]; ph_hi = PHI_HI[k - KMIN] mass = dn(ph_lo[1:] - ph_hi[:-1]) * (1 - 4 * EPS) mass = np.maximum(mass, 0.0) a = (t - rs[None, :] * z[:-1, None]) / s[None, :] + MA equal = len(set(rs.tolist())) == 1 lb = gauss_lo(a, equal) st = [0, 0.0] # Arb arc evaluations, largest ball radius if mode_arc: est, rb = arc_estimate(a) for i in np.nonzero(est > lb + DELTA)[0]: v, rad = arc_arb_lo(a[i], float(rb[i])) lb[i] = max(lb[i], v); st[0] += 1; st[1] = max(st[1], rad) prod = lb * mass * (1 - 8 * 2.0 ** -53) tot = int(np.floor(prod * 2.0 ** 62).astype(np.int64).sum()) return max(float(dn(float(tot) / 2.0 ** 62 * (1 - 2.0 ** -50))), 0.0), st def I_all(tk): out = {}; n = 0; rad = 0.0 for d, ty in ALL: v, st = I_lo(d, ty, tk) out[(d, ty)] = v; n += st[0]; rad = max(rad, st[1]) return tk, out, (n, rad) # ---- cells and charging def cell_rows(Ilo, k0, prob): """Rows (coefficients of the transfers, C) of the charging inequalities on the cell with left end k0 2^-17.""" degs = prob["degs"]; a = prob["a"]; rows = []; Cs = [] top = float(PHI_HI[-k0 - KMIN]) # Phi(-t0) <= top for (d, ty), I in Ilo.items(): 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 prob["cubic_cap"] is not None and prob["cubic_cap"] <= a[0]: continue # known bound for the pure cubic star: no transfers involved C = float(up(top - I)) if (d, ty) == (3, (3, 0, 0)) and prob["cubic_cap"] is not None: C = min(C, float(prob["cubic_cap"])) 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] rows.append((d, row)); Cs.append(C) return rows, Cs def exact_slack(rows, Cs, tr, a): """minimum over rows of a_d - C - row.tr in exact rational arithmetic (a: Fractions, C, tr: doubles).""" trf = [Fraction(x) for x in tr] return min(a[d - 3] - Fraction(C) - sum(c * x for c, x in zip(row, trf)) for (d, row), C in zip(rows, Cs)) def cell_check(Ilo, k0, prob): rows, Cs = cell_rows(Ilo, k0, prob) a = prob["a"] A = np.array([r + [-1.0] for _, r in rows]); b = np.array([float(a[d - 3]) - C for (d, _), C in zip(rows, Cs)]) res = linprog([0, 0, 0, 1], A_ub=A, b_ub=b, bounds=[(-1, 1)] * 3 + [(None, None)], method="highs") tr = [float(v) for v in res.x[:3]] return exact_slack(rows, Cs, tr, a), tr F_UP = Fraction("0.14015825233145217") # >= F (Arb: 0.1401582523314521636939405140175, Proposition 6.9 of the paper) C3 = 3 / (4 * np.pi) * np.arccos(5 / 6) def Fr(*xs): return tuple(Fraction(x) for x in xs) PROBLEMS = { "M": dict(a=Fr(F_UP, "0.2575", "9.9"), degs=(3, 4), cubic_cap=F_UP, lo=-3, hi=3), "M5-P": dict(a=Fr(F_UP, "0.2575", "0.3667"), degs=(3, 4, 5), cubic_cap=F_UP, lo=-3, hi=3), "M5-R": dict(a=Fr("0.1455", "0.2530", "0.3600"), degs=(3, 4, 5), cubic_cap=F_UP, lo=-3, hi=3), "M5-S": dict(a=(Fraction(5, 32), Fraction(1, 4), Fraction(11, 32)), degs=(3, 4, 5), cubic_cap=F_UP, lo=-3, hi=3), # Seed D = (c_3, 0.2627) of the median-switched sweep for maximum degree 4 (Theorem 8.2 of the paper), # added 2026-10-14: thresholds t in [0, 1.09375] only (t < 0 by symmetry; beyond, Phi(-1.09375) < c_3). # The pure cubic star is bounded by c_3 at t >= 0 (manuscript, prop:star-cubic) and needs no row; the other # degree-3 rows use the rational lower bound C3_LO < c_3, which is a stronger requirement than c_3. "S-switched": dict(a=(Fraction("0.139822124008"), Fraction("0.2627"), Fraction("9.9")), degs=(3, 4), cubic_cap=Fraction("0.139822124008"), lo=0, hi=Fraction(35, 32)), } SLACK_MIN = Fraction(1, 10 ** 12) HMIN_K = 1 # smallest cell: 2^-17 def self_test(): rng = np.random.default_rng(1) for name, f in (("Phi_ub", Phi_ub), ("Phi_lb", Phi_lb)): xs = rng.uniform(-9, 14, 400) for x in xs: e = (1 + (arb(float(x)) / SQ2).erf()) / 2 v = f(np.array([x]))[0] assert (v >= min(float(e.upper()), 1.0) - 1e-16) if name == "Phi_ub" else (v <= e.lower()), (name, x, v) # sanity check (not part of the proof): Arb arc bound against the floating-point estimate at random thresholds gap = 0.0; npts = 0 a3 = rng.uniform(-4, 3, (300, 3)) est, rb = arc_estimate(a3) for i in range(300): v, _ = arc_arb_lo(a3[i], float(rb[i])) if est[i] > 0: gap = max(gap, abs(v - est[i])); npts += 1 assert gap < 1e-12, gap print("self-test: Phi bounds hold at 800 random points; Arb arc bound within %.1e of the float estimate at %d points" % (gap, npts), flush=True) for (i, j), v in RA.items(): assert abs(RF[(i, j)] - float(v.mid())) <= 2e-16 and abs(float(v.mid()) - float(v.mid())) == 0 s = (1 - v * v).sqrt() assert abs(np.sqrt(1 - RF[(i, j)] ** 2) - float(s.mid())) < 4e-16 ref = [] for d in (3, 4, 5): for h in (0.5, 1.7, 3.0): e = FLO[d][int(h * 2 ** 19)] ref.append((d, h, e)) print("simplex lower tables (d,h,value):", [(d, h, round(e, 10)) for d, h, e in ref], flush=True) def main(names): build_tables() self_test() ctxm = mp.get_context("fork"); pool = ctxm.Pool(NPROC) cache = {}; stats = [0, 0.0] os.makedirs(OUT, exist_ok=True) def need(ks): ks = sorted(set(k for k in ks if k not in cache)) for k, out, st in pool.imap_unordered(I_all, ks, chunksize=1): cache[k] = out; stats[0] += st[0]; stats[1] = max(stats[1], st[1]) report = {} for name in names: prob = PROBLEMS[name]; t_start = time.time() n0 = int((prob["hi"] - prob["lo"]) * 32); w0 = G // 32; g0 = int(prob["lo"] * G) cells = [(g0 + i * w0, g0 + (i + 1) * w0) for i in range(n0)] done = []; failed = [] while cells: need([c[1] for c in cells]); nxt = [] for k0, k1 in cells: slack, tr = cell_check(cache[k1], k0, prob) if slack >= SLACK_MIN: done.append((k0, k1, float(slack), tr)) elif k1 - k0 > HMIN_K and not (k1 - k0 <= G // 4096 and slack < Fraction(-2, 10 ** 4)): m = (k0 + k1) // 2; nxt += [(k0, m), (m, k1)] else: failed.append((k0, k1, float(slack))) cells = nxt print("%s: %d cells certified, %d pending, %d failed, %d thresholds, %d Arb arc evaluations, %.0fs" % (name, len(done), len(cells), len(failed), len(cache), stats[0], time.time() - t_start), flush=True) done.sort() cover = done[0][0] == g0 and done[-1][1] == int(prob["hi"] * G) and all(done[i][1] == done[i + 1][0] for i in range(len(done) - 1)) report[name] = dict(a=[str(x) for x in prob["a"]], cells=len(done), failed=len(failed), covered=bool(cover), certified=bool(cover and not failed), min_cell_log2=float(np.log2(min(c[1] - c[0] for c in done) / G)) if done else None, min_slack=min(c[2] for c in done) if done else None, seconds=time.time() - t_start) print("RESULT", name, json.dumps(report[name]), flush=True) with open(os.path.join(OUT, name + ".json"), "w") as f: json.dump(dict(problem=dict(a=[str(x) for x in prob["a"]], degs=prob["degs"], cubic_cap=str(prob["cubic_cap"]), grid=G), cells=done, failed=failed), f) ks = sorted(cache) fname = "I_lower_bounds.json" if set(names) >= {"M", "M5-P", "M5-R", "M5-S"} else "I_lower_bounds-%s.json" % "-".join(names) with open(os.path.join(OUT, fname), "w") as f: json.dump({"grid": G, "types": [[d, list(ty)] for d, ty in ALL], "thresholds": {str(k): [cache[k][ty] for ty in ALL] for k in ks}}, f) pool.close() print("Arb arc evaluations: %d in total; largest ball radius %.1e" % tuple(stats), flush=True) return report if __name__ == "__main__": # Default: the four problems of 2026-10-10, whose combined data go to I_lower_bounds.json; S-switched (2026-10-14) # is run on its own (argument S-switched) and writes I_lower_bounds-S-switched.json. main(sys.argv[1:] or ["M", "M5-P", "M5-R", "M5-S"])