#!/usr/bin/env python3
"""cert_grh_sieve.py: the sieve range under GRH (the paper's Lemmas 24 and 25).
Claim certified, for K in {9, 8}: under GRH, every even N with 73 <= log N <= L_TOP (L_TOP > e^{15.85}) is p + eta, eta >= 2, eta free of primes
below N^{1/11}, Omega(eta) <= K.  Above exp(exp(15.85)) Bordignon-Starichkova (arXiv:2211.08844v1, Corollary 1, GRH) give Omega <= 2.
Ingredients:
  * good: Johnston's thesis Prop 4.18 (GRH): N >= e^73, z = N^{1/11}: S(A, P, z) > 0.146 U_N N/L^2; U_N >= U_low = 2e^{-gamma} TP_lo.
  * bad (Omega(eta) >= K+1 among z-rough eta < N = z^11, Omega <= 10): <= N^{1-delta} + N Phi/((1-delta) L), Lemma 12
    with M = 11, recursive mu, and Lemma 11 with eps(z) := S(z) + D(log z) + 1/z, where S(z) := sup_{x >= z} (sum_{p<=x} 1/p - loglog x - B)^+:
    Dusart Thm 6.10 gives <= D(log x) for x >= 10372 and the lower bound -D(log x) for all x > 1; below 10372, S(z) is bounded above by
    interval prime sums, a preceding-prime suffix maximum, and the Dusart tail bound, computed from the primes (E(x) is maximal at primes on
    each gap). This covers z >= e^{73/11} ~ 762.
  * optional Lemma 13 when z >= 8e9: a = 1 coefficient G(1 + 480/(840 log z)), G = 9/32, + 4 exceptions (the K = 8 run of Remark 26).
  * condition per N:  L Phi/((1-delta) 0.146 U) + (L^2 N^{-delta} + 5 L^2/N)/(0.146 U) < 1; on [L1, L2]: L Phi <= L2 Phi(z1), L^2 e^{-delta L} <=
    L1^2 e^{-delta L1} (delta L1 >= 2), L^2/N <= L2^2 e^{-L1}; Phi(z) decreases in z (eps, 1/(z-1), A_a decrease).
J_m(10) upper bounds: exact cell integrals log(1 + 1/(H+i)) times 1 + 1e-14, lower-corner cover, direct np.convolve (times 1 + 1e-10 per level)."""
import os, sys, math
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__))
M = 11; W = M - 1; H = 400
B = I('0.26149721284764278375542683860869585905156664826120'); B = iv.mpf([lo(B) - mp.mpf('1e-45'), hi(B) + mp.mpf('1e-45')])
def D(l): return 1 / (10 * l ** 2) + I(4) / (15 * l ** 3)
# --- S(z): sup_{x >= z} E(x)^+, E(x) = sum_{p<=x} 1/p - loglog x - B, computed at primes below 10372 ---
NPR = 10372
sv = np.ones(NPR + 1, dtype=bool); sv[:2] = False
for i in range(2, int(NPR ** 0.5) + 1):
    if sv[i]: sv[i * i::i] = False
PR = [int(p) for p in np.nonzero(sv)[0]]
Ep = {}; acc = I(0)
for p in PR:
    acc += I(1) / p
    Ep[p] = acc - iv.log(iv.log(I(p))) - B
TAIL = hi(D(iv.log(I(NPR))))                     # sup over x >= 10372 (Dusart upper bound), D decreasing
suffix = {}; run = TAIL
for p in reversed(PR):
    run = max(run, hi(Ep[p]), mp.mpf(0)); suffix[p] = run
def S_of(z):                                     # directed selection: sup_{x>=z} E(x)^+; D(log z) for z >= 10372
    zlo = lo(z)
    if zlo >= NPR:
        return hi(D(I(lo(iv.log(z)))))
    pk = max(p for p in PR if p <= zlo)
    return suffix[pk]
def eps_of(l1):                                  # l1 = log z (interval)
    z = iv.exp(l1)
    return I(S_of(z)) + D(l1) + 1 / z
# --- J_m(10) ---
n_max = (W - 1) * H
wv = np.log1p(1.0 / (H + np.arange(n_max + 1, dtype=np.float64))) * (1 + 1e-14)
J = {0: 1.0}; Dm = wv.copy()
for m in range(1, W):
    if m > 1: Dm = np.convolve(Dm, wv)[: n_max + 1] * (1 + 1e-10)
    J[m] = float(np.sum(Dm[: (W - m) * H + 1])) * (1 + 1e-10)
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
GOOD = I('0.146') * U
CTH = 1 + I(1) / 36260
def Phi(l1, K, presieve):
    z = iv.exp(l1); eps = eps_of(l1); iz = 1 / (z - 1)
    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) + (iz * mu[j - 2] if j >= 2 else I(0))
    def Az(a): return a * (l1 / z ** a + l1 / ((a - 1) * z ** (a - 1)) + 1 / (I(a - 1) ** 2 * z ** (a - 1)))
    c1 = (I(9) / 32) * (1 + I(480) / (840 * l1)) if presieve else CTH
    return sum((c1 * mu[k - 1] + sum((Az(a) * mu[k - a] for a in range(2, k + 1)), I(0)) for k in range(K + 1, M)), I(0))
def run(K, ratio, L_TOP, use_P):
    L1 = mp.mpf(73); rows = 0; worst = None; worst_at = None
    while L1 < L_TOP:
        L2 = min(L_TOP, L1 * ratio)
        l1 = I(L1) / M
        presieve = use_P and lo(l1) > hi(iv.log(I(8) * I(10) ** 9))
        Ph = Phi(l1, K, presieve)
        best = None
        for d in ('0.30', '0.25', '0.20', '0.15', '0.10', '0.05', '0.02', '0.01'):
            dl = I(d)
            if lo(dl * I(L1)) < 2: continue
            v = I(L2) * Ph / ((1 - dl) * GOOD) + (I(L1) ** 2 * iv.exp(-dl * I(L1)) + 5 * I(L2) ** 2 * iv.exp(-I(L1))) / GOOD
            if best is None or hi(v) < best: best = hi(v)
        assert best is not None
        rows += 1
        if worst is None or best > worst: worst, worst_at = best, (float(L1), float(L2), presieve)
        if best >= 1: print('  FAIL on [%s, %s] (presieve %s): %s' % (mp.nstr(L1, 8), mp.nstr(L2, 8), presieve, mp.nstr(best, 6)))
        L1 = L2
    return rows, worst, worst_at
if __name__ == '__main__':
    L_TOP = mp.mpf('7.7e6')
    assert lo(I(L_TOP)) > hi(iv.exp(I('15.85'))), 'cover beyond e^{15.85}'
    print('S(z) at z = e^{73/11}: %s; eps there: %s' % (mp.nstr(S_of(iv.exp(I(73) / 11)), 6), mp.nstr(hi(eps_of(I(73) / 11)), 6)))
    for K, P in ((9, False), (8, False), (8, True)):
        rows, worst, at = run(K, mp.mpf('1.05'), L_TOP, P)
        print(f"K={K}, presieve={P}: {rows} intervals; computed upper-bound enclosure {I(worst)} "
              f"at {at}; {'CERTIFIED' if worst < 1 else 'NOT ESTABLISHED'}")