"""cert_S4.py: ball-arithmetic recomputation, with rational cells, of the paper's Lemma 16 (the Johnston-Starichkova range [X_2', X_3]
at K = 18, from Lemma 15), beside the certificate of record third/t2_check/cert_S3P_K18.py; also the K = 17 diagnostic. It uses
  (i) Lemma 13's a = 1 coefficient G(1 + 480/(840 log z0)), G = 9/32, in place of c_theta (z0 = e^{L0/40} >= 8e9, asserted), and the at most
      4 exceptions (N - n in {3,5,7,11}) folded into the L^2/N term as 5 L3^2 e^{-L0};
  (ii) the coefficient of (I4') taken piecewise: c >= 9.8036 on [L0, 3000], with L_top = 3000, and c >= 26.218 on [3000, L3], with L_top = L3.
The J_m(39) upper bounds, N_j, mu and A_a(z) are re-implemented here, not imported.
J_m(w) <= sum_{n <= (w-m)H} D_m(n), D_m = m-fold discrete convolution of 1/(H + i) (cells [1 + i/H, 1 + (i+1)/H), weight >= int du/u);
np.convolve is a direct sum (not FFT); each convolution is inflated by 1 + 1e-10. Single-thread BLAS is required.
The final print shows the full Arb enclosure of an upper majorant ("upper-majorant enclosure").
Usage: cert_S4.py H delta K [K ...]
"""
import os, sys, math
assert os.environ.get("OPENBLAS_NUM_THREADS") == "1" and os.environ.get("OMP_NUM_THREADS") == "1", "set thread limits to 1"
import numpy as np
from flint import arb, ctx
ctx.prec = 200

H = int(sys.argv[1]); delta_s = sys.argv[2]; Ks = [int(v) for v in sys.argv[3:]]
W = 39
L0 = arb('2736.86'); L3 = arb('32.7').exp()
l0 = L0 / 40
D = lambda x: 1 / (10 * x ** 2) + arb(4) / (15 * x ** 3)
z0 = l0.exp()
eps = arb((2 * D(l0) + 1 / z0).upper())
delta = arb(delta_s)
assert (delta * L0 - 2) > 0
assert (l0 - arb('8e9').log()) > 0                      # Lemma P needs z >= 8e9 (BMOR x_0 at q = 1155)
gamma = arb.const_euler()
Ulow = arb((2 * (-gamma).exp() * arb('0.6601618158')).lower())
rep = arb((39 * (2 * eps).exp() / (z0 - 1)).upper())
n_max = (W - 1) * H
cvec = 1.0 / (H + np.arange(n_max + 1, dtype=np.float64)) * (1 + 1e-15)
Jup = {0: 1.0}
Dm = cvec.copy()
for m in range(1, W):
    if m > 1:
        Dm = np.convolve(Dm, cvec)[: n_max + 1] * (1 + 1e-10)
    Jup[m] = float(np.sum(Dm[: (W - m) * H + 1])) * (1 + 1e-10)
def Nj(j):
    return sum((arb(math.comb(j, i)) * eps ** i * arb(Jup[j - i]) for i in range(0, j + 1)), arb(0))
mu = {0: arb(1)}
for j in range(1, W):
    mu[j] = arb((Nj(j) / arb(math.factorial(j)) + (rep if j >= 2 else arb(0))).upper())
def A_z(a, z=z0):
    lz = z.log()
    return arb((a * (lz / z ** a + lz / ((a - 1) * z ** (a - 1)) + 1 / ((a - 1) ** 2 * z ** (a - 1)))).upper())
cpre = arb((arb(9) / 32 * (1 + arb(480) / (840 * l0))).upper())
for K in Ks:
    PhiP = arb(0)
    for k in range(K + 1, W + 1):
        PhiP += cpre * mu[k - 1] + sum((A_z(a) * mu[k - a] for a in range(2, k + 1)), arb(0))
    for Ltop, c_s, lab in ((arb(3000), "9.8036", "[L0, 3000]"), (L3, "26.218", "[3000, L3]")):
        cc = arb(c_s)
        first = Ltop * PhiP / ((1 - delta) * cc * Ulow)
        second = (L0 ** 2 * (-delta * L0).exp() + 5 * L3 ** 2 * (-L0).exp()) / (cc * Ulow)
        lhs = first + second
        ok = (1 - lhs) > 0
        print("S3+P H=%d delta=%s K=%d on %s with c=%s: Phi_P < %s; upper-majorant enclosure %s -> %s" % (
            H, delta_s, K, lab, c_s, PhiP.upper().str(6, radius=False), lhs.str(20),
            "CERTIFIED" if ok else "not established"), flush=True)