"""cert_C2.py: interval certificate (mpmath iv) for the middle-range count in its Rankin form, and the shared ingredients that cert_C3.py
and cert_C3s.py import (the prime sums A_a, k_max, the small-n allowance, G_0). The paper uses this file only through those two scripts.

For even N in [Y_low, X] (Y_low >= 4e18; equality permitted), L = log N, E = {n <= N: (n, 2310) = 1, Omega(n) >= K + 1}:
  #E <= sqrt N + (2/L) sum_{k=K+1}^{k_max} S_k,   S_k := sum_{n <= N, (n,2310)=1, Omega(n)=k} log n,
  S_k <= sum_{a>=1} sum_{p>=13} a log p #{n' <= N/p^a : (n',2310)=1, Omega(n') = k - a}
      <= c_theta N R_{k-1} + N sum_{a=2}^{k} A_a R_{k-a},
  with c_theta = 1 + 1/36260 (Dusart 2010 Prop. 5.1: theta(x) < (1 + 1/36260) x, x > 0), A_a := sum_{p>=13} a log p/p^a, and
  R_j := tau_j^{-j} exp(tau_j S'(X) + c_2(tau_j)) (R_0 := 1) >= sum_{n' <= N, (n',2310)=1, Omega(n')=j} 1/n' (Rankin, 0 < tau_j < 13).
  So the bad (log-weighted) part divided by N is at most
      L_X e^{-L_low/2} + 2 sum_{k=K+1}^{k_max} [ c_theta R_{k-1} + sum_{a=2}^{k} A_a R_{k-a} ],
  which must be < G_0 (1 - 1e-18) - 1e-18.
Rigour: mpmath iv; A_a = exact sum over primes 13 <= p <= P0 plus a tail bound; sum over a > A_MAX bounded by a geometric tail.
"""
import math, sys
from fractions import Fraction
from mpmath import iv, mpf
iv.prec = 120

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]
B_UP = iv.mpf('0.26149721284765')
SIG11 = sum((iv.mpf(1) / p for p in (2, 3, 5, 7, 11)), iv.mpf(0))
G0 = 135 * (iv.mpf(1) / 480 - 1 / (840 * iv.log(iv.mpf('4e18') - 2)))
CTH = 1 + iv.mpf(1) / 36260
A_MAX = 80

def tail_log_over_power(a, P):
    """upper bound for sum_{n > P} log n / n^a (a >= 2, P >= 3): int_P^inf log t/t^a dt + log P/P^a."""
    P = iv.mpf(P)
    return iv.log(P) / ((a - 1) * P ** (a - 1)) + 1 / ((a - 1) ** 2 * P ** (a - 1)) + iv.log(P) / P ** a

_A_CACHE = None
def A_coeffs():
    global _A_CACHE
    if _A_CACHE is not None: return _A_CACHE
    A = {}
    for a in range(2, A_MAX + 1):
        s = iv.mpf(0)
        for p in PR:
            s += iv.log(iv.mpf(p)) / iv.mpf(p) ** a
        A[a] = a * (s + tail_log_over_power(a, P0))
    _A_CACHE = A
    return A

def c2_upper(tau):
    ta = iv.mpf(tau.numerator) / tau.denominator
    s = iv.mpf(0)
    for p in PR:
        if p > 1000: break
        x = ta / p
        s += -iv.log(1 - x) - x
    return s + ta * ta / (2 * (1 - ta / 1009)) * iv.mpf('0.001')

_c2cache = {}
def R(j, Sp):
    if j == 0: return iv.mpf(1)
    tau = Fraction(max(2, min(51, round(4 * j / float(mpf(Sp.a))))), 4)     # quarter steps in [0.5, 12.75], near j/S'
    if tau not in _c2cache: _c2cache[tau] = c2_upper(tau)
    ta = iv.mpf(tau.numerator) / tau.denominator
    return iv.exp(ta * Sp + _c2cache[tau] - j * iv.log(ta))

def lhs(LX, Llow, K):
    LX = iv.mpf(LX); Llow = iv.mpf(Llow)
    Sp = iv.log(LX) + B_UP + 1 / (10 * LX ** 2) + 4 / (15 * LX ** 3) - 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
    A = A_coeffs()
    Rs = [R(j, Sp) for j in range(kmax + 1)]
    Rmax = max(r.b for r in Rs)
    # sum_k sum_{a=2}^{k} A_a R_{k-a} = sum_j R_j sum_{a >= max(2, K+1-j), a <= kmax - j} A_a ; tail a > A_MAX bounded by
    # sum_{a > A_MAX} a * sum_{n>=13} log n/n^a <= sum_{a > A_MAX} a * (log 13/13^a + tail_log_over_power(a, 13)), geometric
    tailA = iv.mpf(0)
    for a in range(A_MAX + 1, A_MAX + 60):
        tailA += a * (iv.log(iv.mpf(13)) / iv.mpf(13) ** a + tail_log_over_power(a, 13))
    tailA *= 2      # covers a >= A_MAX + 60 (ratio < 1/12 per step)
    # [the small-n allowance L e^{-L/2} <= L_X N^{-1/2} <= L_X/(2·10^9) for N >= 4·10^18,
    # taken exactly instead of from a floating log(4e18), which rounded upward by 2.07e-15]
    tot = LX / (2 * iv.mpf(10) ** 9)
    for k in range(K + 1, kmax + 1):
        tot += 2 * CTH * Rs[k - 1]
    for j in range(0, kmax + 1):
        amin = max(2, K + 1 - j); amax = min(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 > A_MAX:
            tot += 2 * Rs[j] * tailA
    rhs = G0 * (1 - iv.mpf('1e-18')) - iv.mpf('1e-18')
    return tot, rhs, tot.b < rhs.a

if __name__ == '__main__':
    print('# python3 3.14.4 + mpmath 1.3.0; Lemma C\'\' certificate')
    for (LX, K) in [('2736.86', 22), ('2736.86', 21), ('350', 16), ('350', 15), ('73', 12), ('73', 11)]:
        tot, rhs, ok = lhs(LX, '0', K)     # Llow is no longer used (RC5-1); Y_low = 4·10^18 inclusive (RC5-2)
        # [RC5-4: directed displays: upper bounds rounded up, lower bounds rounded down]
        up = math.ceil(float(mpf(tot.b)) * 1e6) / 1e6; lo = math.floor(float(mpf(rhs.a)) * 1e7) / 1e7
        print(f"Lemma C'' for all even N in [4e18, exp({LX})], K = {K}: bad/N < {up:.6f} vs G0' > {lo:.7f}: CERTIFIED = {ok}")