#!/usr/bin/env python3
"""cert_S3P_K18.py: the certificate of record for the paper's Lemma 16 (Theorem 1 on [X_2', X_3] at K = 18).
Statement certified (from Lemmas 11-15 of the paper; this script is the computation):
  every even N with X2' = exp(2736.86) <= N <= X3 = exp(exp(32.7)) is p + eta, eta >= 2, eta free of primes below N^{1/40}, Omega(eta) <= 18.
Inputs:
  * Lemma 15, with the coefficient c of (I4') taken PIECEWISE in L = log N:
      c >= 9.8036 on [L0, 3000]   (../js_cert/js_cert.json routeB_BJV22.min_lower_case_b = 9.80368416..., case (a) >= 38.417)
      c >= 26.218 on [3000, L3]   (route A at L_s = 3000, tail_case_b = 26.21829386..., case (a) >= 39.28; the independent audit gives > 26.21)
    so on each piece the condition is  Ltop Phi_P/((1-delta) c U_low) + (L0^2 e^{-delta L0} + 5 L3^2 e^{-L0})/(c U_low) < 1, Ltop = 3000 or L3.
  * Lemma 13: z = N^{1/40} >= e^{68.4} >= 8e9, so the a = 1 coefficient c_theta becomes G (1 + 480/(840 log z0)),
    G = 9/32, with + 4 exceptions (folded into the L^2/N term as the factor 5); Pi cancels against U_N >= U_low Pi (good > c U_N N/L^2).
  * Lemma 12 in its crude form at z0 (valid: 1/(z0 - 1) <= eps), and Lemma 11 with J_m(39) <= exact-cell-integral upper sums
    (log(1 + 1/(H+i)) times 1 + 1e-14, lower-corner cover), direct np.convolve inflated (1 + 1e-10) per level;
    U_low = 2 e^{-gamma} TP_lo (rigorous, ../js_cert/tp_bound.out).
  As the paper's computations section says, the computed numbers bound a majorant of the left side as printed (N_j(w) by
  sum_i C(j,i) eps^i J_{j-i}(w), the common L0/L3 remainder, and U_low from the twin-prime lower bound).
Outputs cert_S3P_K18.out. Run with OMP/OPENBLAS/MKL_NUM_THREADS=1."""
import os, sys, math, json
for v in ('OMP_NUM_THREADS', 'OPENBLAS_NUM_THREADS', 'MKL_NUM_THREADS'): assert os.environ.get(v) == '1', v
import numpy as np
from mpmath import iv, mp
iv.dps = 40; mp.dps = 40
def I(x): return iv.mpf(x)
def lo(x): return mp.make_mpf(x._mpi_[0])
def hi(x): return mp.make_mpf(x._mpi_[1])
HERE = os.path.dirname(os.path.abspath(__file__))
jc = json.load(open(os.path.join(HERE, '../js_cert/js_cert.json')))['routeB_BJV22']
C_LOW, C_HIGH = mp.mpf('9.8036'), mp.mpf('26.218')
assert mp.mpf(jc['min_lower_case_b']) > C_LOW and mp.mpf(jc['min_lower_case_a']) > C_LOW, 'route B piece'
assert mp.mpf(jc['tail_case_b']) > C_HIGH and mp.mpf(jc['tail_case_a']) > C_HIGH, 'route A tail piece'
assert jc['L1'] == '3000.0'
W = 39; K = 18
def J_up(H):
    n_max = (W - 1) * H
    w = np.log1p(1.0 / (H + np.arange(n_max + 1, dtype=np.float64))) * (1 + 1e-14)
    J = {0: 1.0}; Dm = w.copy()
    for m in range(1, W):
        if m > 1: Dm = np.convolve(Dm, w)[: n_max + 1] * (1 + 1e-10)
        J[m] = float(np.sum(Dm[: (W - m) * H + 1])) * (1 + 1e-10)
    return J
L0 = I('2736.86'); L3 = iv.exp(I('32.7')); l0 = L0 / 40; z0 = iv.exp(l0)
assert lo(l0) > hi(iv.log(I(8) * I(10) ** 9)), 'Lemma P gate z0 >= 8e9'
D = lambda x: 1 / (10 * x ** 2) + I(4) / (15 * x ** 3)
EPS = 2 * D(l0) + 1 / z0
assert hi(1 / (z0 - 1)) <= lo(EPS), 'crude rep needs 1/(z-1) <= eps'
REP = 39 * iv.exp(2 * EPS) / (z0 - 1)
G = I(9) / 32; C1 = G * (1 + I(480) / (840 * l0))
TPlo = I(open(os.path.join(HERE, '../js_cert/tp_bound.out')).read().split('TP_lo_str ')[1].split()[0])
U = 2 * iv.exp(-iv.euler) * TPlo
DELTA = I('0.01'); assert lo(DELTA * L0) >= 2
def A_up(a): return a * (l0 / z0 ** a + l0 / ((a - 1) * z0 ** (a - 1)) + 1 / (I(a - 1) ** 2 * z0 ** (a - 1)))
with open(os.path.join(HERE, 'cert_S3P_K18.out'), 'w') as out:
    def say(t): print(t); out.write(t + '\n')
    say('# cert_S3P_K18.py: S\'\'\' + Lemma P, K = 18, piecewise c (9.8036 on [L0, 3000], 26.218 on [3000, L3]); delta = 0.01; U_low = 2e^-g TP_lo')
    for H in (200, 400):
        J = J_up(H)
        mu = {0: I(1)}
        for j in range(1, W):
            Nj = sum((math.comb(j, i) * EPS ** i * I(J[j - i]) for i in range(0, j + 1)), I(0))
            mu[j] = Nj / math.factorial(j) + (REP if j >= 2 else I(0))
        PhiP = sum((C1 * mu[k - 1] + sum((A_up(a) * mu[k - a] for a in range(2, k + 1)), I(0)) for k in range(K + 1, W + 1)), I(0))
        second = lambda c: (L0 ** 2 * iv.exp(-DELTA * L0) + 5 * L3 ** 2 * iv.exp(-L0)) / (I(c) * U)
        lowp = I(3000) * PhiP / ((1 - DELTA) * I(C_LOW) * U) + second(C_LOW)
        highp = L3 * PhiP / ((1 - DELTA) * I(C_HIGH) * U) + second(C_HIGH)
        ok = hi(lowp) < 1 and hi(highp) < 1
        say(f"H={H}: Phi_P upper-majorant enclosure {PhiP}; "
            f"low-piece upper-majorant enclosure {lowp}; "
            f"high-piece upper-majorant enclosure {highp}; "
            f"{'CERTIFIED' if ok else 'NOT ESTABLISHED'}")
        assert ok