"""cert_C3s.py: interval certificate for the condition of the paper's Lemma 6 (the middle range), with the small-n split at N^{1-delta}.
Uses cert_C3.py (all truncations enclosed) for R*_j and A_a.

The bad part: E = {n <= N : (n, 2310) = 1, Omega(n) > K}; the bad part of the good count is at most L |E|, and
  |E| <= #{n <= N^{1-delta}} + #{n in E : n > N^{1-delta}} <= N^{1-delta} + (1/((1-delta) L)) sum_{k > K} S_k,
since log n > (1-delta) L for n > N^{1-delta}. Hence bad/N <= L e^{-delta L} + (1-delta)^{-1} sum_{k>K} [c_theta R*_{k-1} + sum_a A_a R*_{k-a}].
For N in [4e18, X], L e^{-delta L} is decreasing in L once delta L > 1, so it is at most L_low e^{-delta L_low}, L_low = log(4e18).
The condition:  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.
(delta = 1/2 gives the factor 2 of the split at sqrt N, with L_low e^{-L_low/2} in place of the cruder L_X/(2e9).)
"""
import math, sys
from fractions import Fraction
from mpmath import iv, mpf
import cert_C2 as C
import cert_C3 as C3
iv.prec = 120

def lhs_split(LX, K, delta):
    LXs = LX; LX = iv.mpf(LX)
    d = iv.mpf(delta.numerator) / delta.denominator
    Llow = iv.log(iv.mpf(4) * iv.mpf(10) ** 18)
    assert (d * Llow).a > 1                       # L e^{-delta L} decreasing on [L_low, inf)
    Sp = iv.log(LX) + C3.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
    Rs = C3.exact_coeffs(Sp, kmax)
    A = C3.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])
    s = iv.mpf(0)
    for k in range(K + 1, kmax + 1):
        s += 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:
            s += Rs[j] * sum((A[a] for a in range(amin, amax + 1)), iv.mpf(0))
        if kmax - j > C.A_MAX:
            s += Rs[j] * tailA
    tot = Llow * iv.exp(-d * Llow) + s / (1 - d)
    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'''_delta (split at N^(1-delta)) certificate")
    for (LX, K, delta) in [('2736.86', 19, Fraction(1, 5)), ('2736.86', 19, Fraction(21, 100)), ('2736.86', 20, Fraction(1, 5)), ('2736.86', 18, Fraction(1, 5))]:
        tot, rhs = lhs_split(LX, K, delta)
        up = math.ceil(float(mpf(tot.b)) * 1e6) / 1e6; lo = 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'''_delta, delta = {delta}, all even N in [4e18, exp({LX})], K = {K}: left side in [{lo:.6f}, {up:.6f}] vs G0' > {rl:.7f}: {verdict}")