"""cert_C3.py: ball-arithmetic certificate (python-flint arb) for the condition of the paper's Lemma 6 (the middle range).

Let 4e18 <= Y_low <= X, L_low = log Y_low, L_X = log X, 0 < delta < 1 with delta*L_low >= 1,
K >= 1, k_max = floor(L_X/log 13). Put S'(X) := log log X + B + d(X) - sigma_11, C_2(z) := sum_{k>=2} (z^k/k) P13(k)
with P13(k) := sum_{p>=13} p^{-k}, and R*_j := [z^j] exp(z S'(X) + C_2(z)) (R*_0 = 1). Then every even N in [Y_low, X] is p + eta with p an odd
prime, eta >= 3 odd, (eta, 1155) = 1 and Omega(eta) <= K, provided
  L_low e^{-delta L_low} + (1 - delta)^{-1} sum_{k=K+1}^{k_max} [ c_theta R*_{k-1} + sum_{a=2}^{k} A_a R*_{k-a} ]
      <  G_0 (1 - 1e-18) - 1e-18,
with c_theta = 1 + 1/36260, A_a := sum_{p>=13} a log p / p^a, G_0 = 135 (1/480 - 1/(840 log(4e18 - 2))). (The paper writes C_11(z) for C_2(z).)

Rigour: arb balls throughout (prec 200); every quantity is a sum of nonnegative terms, so the UPPER end of each ball is used where an upper bound
is needed; prime sums are exact over 13 <= p <= P0 plus explicit integral tails; A_a for a > A_MAX uses the crude bound
A_a <= a [log 13/13^a + log 13/((a-1) 13^{a-1}) + 1/((a-1)^2 13^{a-1})] (sum over all integers n >= 13, integral comparison).
The final print shows the full Arb enclosures and says NOT ESTABLISHED where the condition is not proved.
Usage: cert_C3.py LX_decimal Ylow_decimal delta K [K ...]
"""
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')                       # B = 0.2614972128476427837... rounded up (T2_NOTE uses the same constant)
SIG11 = arb(1) / 2 + arb(1) / 3 + arb(1) / 5 + arb(1) / 7 + arb(1) / 11
CTH = 1 + arb(1) / 36260
G0 = 135 * (arb(1) / 480 - 1 / (840 * (arb('4e18') - 2).log()))
TARGET = G0 * (1 - arb('1e-18')) - arb('1e-18')
A_MAX = 80

def up(x):            # upper end of a ball, as an exact arb point
    return arb(x.upper())

def P13(k, KBIG=40):
    """upper bound for sum_{p >= 13} p^{-k}, k >= 2"""
    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, Ylows, delta_s, Ks):
    LX = arb(LXs); Llow = arb(Ylows).log(); delta = arb(delta_s)
    assert (delta * Llow - 1) > 0, "need delta*L_low >= 1 for the monotonicity of L e^{-delta L}"
    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)))
    # certify kmax: kmax*log13 < LX < (kmax+1)*log13 is not needed for an upper bound if we sum to kmax' >= true kmax; take kmax+1 for safety
    kmax_safe = kmax + 1
    assert (arb(kmax_safe + 1) * arb(13).log() - LX) > 0
    J = kmax_safe + 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)
    # R*_j = sum_i S^{j-i}/(j-i)! g_i
    Spow = [arb(1)] * (J + 1)
    for i in range(1, J + 1):
        Spow[i] = Spow[i - 1] * S / i                   # S^i / 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)]
    out = []
    for K in Ks:
        tot = arb(0)
        for k in range(K + 1, kmax_safe + 1):
            tot += CTH * R[k - 1] + sum((A[a] * R[k - a] for a in range(2, k + 1)), arb(0))
        lhs = Llow * (-delta * Llow).exp() + tot / (1 - delta)
        ok = (TARGET - lhs) > 0
        out.append((K, lhs, ok))
        print(f"C''' log X = {LXs}, Y_low = {Ylows}, delta = {delta_s}, K = {K}: "
              f"majorant enclosure {lhs.str(25)}; target enclosure {TARGET.str(25)} "
              f"-> {'CERTIFIED' if ok else 'NOT ESTABLISHED'}", flush=True)
    return out

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