"""cert_C4.py: ball-arithmetic certificate for the paper's Lemma 8 (the middle range, class-restricted Hardy-Ramanujan on the squarefree
part), on a sub-range [Y_i, X] of the middle range; it gives the numbers of Lemma 10.
An earlier variant kept the Hardy-Ramanujan terms a >= 2 unrestricted. This one removes them: n = s r (s squarefull, r squarefree, coprime),
with Hardy-Ramanujan applied to r only.

Bad part / N  <=  L_i e^{-delta L_i} + 2 zeta(3/2) L_i e^{-delta2 L_i / 2}
    + (1 - delta - delta2)^{-1} sum_{k=K+1}^{k_max} [ (|A|/480)(1 + eps_Y) R*_{k-1} + c_theta W_k ],
  eps_Y := 480/(840 log Y), Y := 8e9 (the x_0 of Bennett-Martin-O'Bryant-Rechnitzer for q = 1155),
  W_k  <=  [ c_theta (1 + log Y) R*_{k-2} + sum_{a=2}^{k-1} A_a R*_{k-1-a} ] / (L_i - log Y)      (cofactors n' in (N/Y, N]),
good >= |A| (1/480 - 1/(840 log(4e18 - 2))) y - log 2  with |A| >= 135.
Dividing by Lambda := |A|/480 >= 135/480, the condition bad < good follows from
  (480/135) [ L_i e^{-delta L_i} + 2 zeta(3/2) L_i e^{-delta2 L_i/2} + (1-delta-delta2)^{-1} sum_k c_theta W_k ]
     + (1-delta-delta2)^{-1} (1 + eps_Y) sum_k R*_{k-1}   <   (1 - 480/(840 log(4e18 - 2))) (1 - 1e-18) - (480/135) 1e-18.
R*_j = [z^j] exp(z S'(X) + C_2(z)) and A_a exactly as in cert_C3.py (code copied, not imported, so this file stands alone).
The final print shows the full Arb enclosure, displayed outward.
Usage: cert_C4.py LX_decimal Li_decimal delta delta2 K [K ...]   (Li = log Y_i; needs delta*Li >= 1, delta2*Li >= 2, Li > log Y)
"""
import sys, math
from flint import arb, ctx
ctx.prec = 200

def primes_upto(n):
    s = bytearray([1]) * (n + 1); s[0] = s[1] = 0
    for i in range(2, int(n ** 0.5) + 1):
        if s[i]: s[i * i::i] = bytearray(len(s[i * i::i]))
    return [i for i in range(n + 1) if s[i]]

P0 = 10 ** 5
PR = [p for p in primes_upto(P0) if p >= 13]
PR_SMALL = [p for p in PR if p <= 2000]
B_UP = arb('0.26149721284765')
SIG11 = arb(1) / 2 + arb(1) / 3 + arb(1) / 5 + arb(1) / 7 + arb(1) / 11
CTH = 1 + arb(1) / 36260
A_MAX = 80
def up(x): return arb(x.upper())
def P13(k, KBIG=40):
    if k <= KBIG:
        s = sum((arb(p) ** (-k) for p in PR), arb(0)); tail = arb(P0) ** (1 - k) / (k - 1)
    else:
        s = sum((arb(p) ** (-k) for p in PR_SMALL), arb(0)); tail = arb(2000) ** (1 - k) / (k - 1)
    return up(s + tail)
def A_coeff(a):
    if a <= A_MAX:
        s = sum((a * arb(p).log() / arb(p) ** a for p in PR), arb(0)); P = arb(P0)
        tail = a * (P.log() / ((a - 1) * P ** (a - 1)) + 1 / ((a - 1) ** 2 * P ** (a - 1)) + P.log() / P ** a)
        return up(s + tail)
    t = arb(13)
    return up(a * (t.log() / t ** a + t.log() / ((a - 1) * t ** (a - 1)) + 1 / ((a - 1) ** 2 * t ** (a - 1))))

def run(LXs, Lis, delta_s, delta2_s, Ks):
    LX = arb(LXs); Li = arb(Lis); delta = arb(delta_s); delta2 = arb(delta2_s)
    logY = arb('8e9').log()
    # certificate gates
    assert Li >= arb('4e18').log() and Li <= LX
    assert delta > 0 and delta2 > 0 and delta + delta2 < 1
    assert delta * Li >= 1 and delta2 * Li >= 2 and Li > logY
    assert all(isinstance(K, int) and K >= 1 for K in Ks)
    d = 1 / (10 * LX ** 2) + arb(4) / (15 * LX ** 3)
    S = up(LX.log() + B_UP + d - SIG11)
    kmax = int(math.floor(float(LXs) / math.log(13))) + 1          # one extra term: an upper bound
    assert kmax >= 0 and arb(kmax + 1) * arb(13).log() > LX           # the cutoff covers every possible Omega
    J = kmax + 2
    c = [arb(0)] * (J + 1)
    for k in range(2, J + 1): c[k] = up(P13(k) / k)
    g = [arb(0)] * (J + 1); g[0] = arb(1)
    for n in range(1, J + 1):
        g[n] = up(sum((k * c[k] * g[n - k] for k in range(2, n + 1)), arb(0)) / n)
    Spow = [arb(1)] * (J + 1)
    for i in range(1, J + 1): Spow[i] = Spow[i - 1] * S / i
    R = [up(sum((Spow[j - i] * g[i] for i in range(0, j + 1)), arb(0))) for j in range(J + 1)]
    A = [arb(0), arb(0)] + [A_coeff(a) for a in range(2, J + 1)]
    epsY = arb(480) / (840 * logY)
    zeta32 = arb(1.5).zeta()
    G0n = (1 - arb(480) / (840 * (arb('4e18') - 2).log())) * (1 - arb('1e-18')) - arb(480) / 135 * arb('1e-18')
    fac = 1 / (1 - delta - delta2)
    for K in Ks:
        main = arb(0); wsum = arb(0)
        for k in range(K + 1, kmax + 1):
            main += R[k - 1]
            W = (CTH * (1 + logY) * R[k - 2] + sum((A[a] * R[k - 1 - a] for a in range(2, k)), arb(0))) / (Li - logY)
            wsum += CTH * W
        T = Li * (-delta * Li).exp() + 2 * zeta32 * Li * (-delta2 * Li / 2).exp()
        restpart = arb(480) / 135 * (T + fac * wsum)
        mainpart = (1 + epsY) * fac * main
        lhs = restpart + mainpart
        ok = (G0n - lhs) > 0
        # directed display
        print(f"C4 log X={LXs}, log Y_i={Lis}, delta={delta_s}, delta2={delta2_s}, K={K}: "
              f"upper-majorant enclosure {lhs.str(25)}; main enclosure {mainpart.str(25)}; "
              f"rest enclosure {restpart.str(25)}; target enclosure {G0n.str(25)} -> "
              f"{'CERTIFIED' if ok else 'NOT ESTABLISHED'}", flush=True)

if __name__ == "__main__":
    run(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], [int(v) for v in sys.argv[5:]])