"""cert_C3.py: the exact power-series coefficients R*_j of the paper's Lemma 4, which replace Rankin's bound, with interval enclosures.

For N <= X and j >= 0:  sum_{m <= N, (m,2310)=1, Omega(m)=j} 1/m  <=  [z^j] prod_{13 <= p <= N} (1 - z/p)^{-1}
   <=  R*_j := [z^j] exp( z S'(X) + c_2(z) ),   c_2(z) := sum_{k>=2} P13(k) z^k / k,   P13(k) := sum_{p >= 13} p^{-k},
because log prod (1 - z/p)^{-1} = z sum_{13<=p<=N} 1/p + sum_{k>=2} (z^k/k) sum_{13<=p<=N} p^{-k} has nonnegative coefficients, each <= the
corresponding coefficient of z S'(X) + c_2(z) (S'(X) >= sum_{13<=p<=X} 1/p by Dusart), and exp is coefficientwise monotone
on series with nonnegative coefficients. (Rankin's bound is R*_j <= tau^{-j} exp(tau S' + c_2(tau)) for every 0 < tau < 13.)
(The paper writes C_11(z) for c_2(z).)
R*_j by the recurrence n e_n = S' e_{n-1} + sum_{k=2}^{n} P13(k) e_{n-k}, e_0 = 1 (all terms >= 0, so intervals do not blow up).
P13(k): primes 13 <= p <= 10^5 for k <= 10 and 13 <= p <= 1000 for k > 10, plus the tail sum_{n > P} n^{-k} <= P^{1-k}/(k-1).
Everything else (A_a, the a > 80 tail, k_max, the small-n allowance L_X/(2e9), G_0) is exactly as in cert_C2.py (imported).
The left side is a fixed number for each (X, K) (no free parameter). Every truncation is carried as an ENCLOSURE, not as a point upper
bound: Mertens' B in [0.26149721284764, 0.26149721284765], each A_a's prime tail as a*[0, tail], the a > 80 terms as [0, tailA], P13's
tails as [0, bound]. So the interval encloses the left side itself: its UPPER end proves the condition holds and its LOWER end proves
that it fails.
"""
import math, sys
from mpmath import iv, mpf
import cert_C2 as C
iv.prec = 120

PR_BIG = C.PR                                   # 13 <= p <= 10^5
B_ENC = iv.mpf(['0.26149721284764', '0.26149721284765'])     # Mertens' constant 0.2614972128476427837...
_AENC = None
def A_enclosures():
    """A_a = sum_{p >= 13} a log p / p^a, enclosed: exact prime sum over 13 <= p <= 10^5 plus a * [0, tail]."""
    global _AENC
    if _AENC is not None: return _AENC
    A = {}
    for a in range(2, C.A_MAX + 1):
        s = iv.mpf(0)
        for p in C.PR:
            s += iv.log(iv.mpf(p)) / iv.mpf(p) ** a
        A[a] = a * (s + C.tail_log_over_power(a, C.P0) * iv.mpf([0, 1]))
    _AENC = A
    return A
PR_SMALL = [p for p in PR_BIG if p <= 1000]

def P13(k):
    P, pr = (10 ** 5, PR_BIG) if k <= 10 else (1000, PR_SMALL)
    s = iv.mpf(0)
    for p in pr:
        s += iv.mpf(p) ** (-k)
    return s + iv.mpf(P) ** (1 - k) / (k - 1) * iv.mpf([0, 1])     # tail in [0, P^{1-k}/(k-1)]

def exact_coeffs(Sp, n):
    Pk = [None, None] + [P13(k) for k in range(2, n + 1)]
    e = [iv.mpf(1)]
    for m in range(1, n + 1):
        s = Sp * e[m - 1]
        for k in range(2, m + 1):
            s += Pk[k] * e[m - k]
        e.append(s / m)
    return e

def lhs(LX, K, cache={}):
    LXs = LX
    LX = iv.mpf(LX)
    Sp = iv.log(LX) + B_ENC + 1 / (10 * LX ** 2) + 4 / (15 * LX ** 3) - C.SIG11
    kmax = int(float(mpf(LX.b)) / math.log(13))
    assert (kmax * iv.log(iv.mpf(13))).b <= LX.a and ((kmax + 1) * iv.log(iv.mpf(13))).a > LX.b
    if LXs not in cache: cache[LXs] = exact_coeffs(Sp, kmax)
    Rs = cache[LXs]
    A = A_enclosures()
    tailA = iv.mpf(0)
    for a in range(C.A_MAX + 1, C.A_MAX + 60):
        tailA += a * (iv.log(iv.mpf(13)) / iv.mpf(13) ** a + C.tail_log_over_power(a, 13))
    tailA = 2 * tailA * iv.mpf([0, 1])          # the a > 80 terms are >= 0 and at most this
    tot = LX / (2 * iv.mpf(10) ** 9)
    for k in range(K + 1, kmax + 1):
        tot += 2 * C.CTH * Rs[k - 1]
    for j in range(0, kmax + 1):
        amin = max(2, K + 1 - j); amax = min(C.A_MAX, kmax - j)
        if amin <= amax:
            tot += 2 * Rs[j] * sum((A[a] for a in range(amin, amax + 1)), iv.mpf(0))
        if kmax - j > C.A_MAX:
            tot += 2 * Rs[j] * tailA
    rhs = C.G0 * (1 - iv.mpf('1e-18')) - iv.mpf('1e-18')
    return tot, rhs

if __name__ == '__main__':
    print(f"# python3 {sys.version.split()[0]} + mpmath iv ({iv.prec} bits); Lemma C''' (exact coefficients) certificate, v2 (all truncations enclosed)")
    for (LX, K) in [('2736.86', 20), ('2736.86', 19), ('350', 14), ('350', 13), ('73', 10), ('73', 9)]:
        tot, rhs = lhs(LX, K)
        up = math.ceil(float(mpf(tot.b)) * 1e6) / 1e6; lo_l = math.floor(float(mpf(tot.a)) * 1e6) / 1e6
        rl = math.floor(float(mpf(rhs.a)) * 1e7) / 1e7
        verdict = 'HOLDS' if tot.b < rhs.a else ('FAILS' if tot.a > rhs.b else 'UNDECIDED')
        print(f"Lemma C''' for all even N in [4e18, exp({LX})], K = {K}: left side in [{lo_l:.6f}, {up:.6f}] vs G0' > {rl:.7f}: {verdict}")