"""cert_qrhB3p.py: the grid certificate for Theorem 2 on 51046 <= log N <= 159011874857757 (the paper's Lemma 20, with the presieve of
Lemma 13); the 1155-presieve variant of cert_qrhB3.py.

author/cert_qrh_range.py's margin_iv(L1, L2, M, lam, 2, M - 1) is imported read-only. With K = M - 1 it returns [main - rem - e^{-L1}]
with bad = 0; it checks the gates on M, L1, lambda, H, s1, s2, F(s1) and the V factors, and of the gates on z1 only z1 >= 694 and
z1^(M-1) >= 10372. This script applies the gate z1 >= 10372 itself: its search stops below log 10372 + 0.01, and an Arb assertion z1 > 10372
holds on every interval. It then subtracts, per N in the interval,
  bad'''/N <= e^{-delta L1} + (1/((1 - delta) L1)) sum_{k=K+1}^{M-1} [ c_theta mu_{k-1} + sum_{a=2}^{k} A_a(z1) mu_{k-a} ],
with mu_j := N_j(M - 1)/j! + mu_{j-2}/(z1 - 1) (j >= 2; non-squarefree cofactors q^2 n''), N_j(w) <= sum_i C(j,i) eps^i J_{j-i}(w),
eps := 2 D(l1) + 1/z1, z1 := e^{L1/M} >= 10372 (gate). Each term decreases in L (z grows), and 1/L <= 1/L1, e^{-delta L} <= e^{-delta L1}.
J_m(w) <= sum_{n <= (w - m) H} D_m(n): the discrete upper-step convolution (cells of width 1/H, weight 1/(H + i)), with float
inflation 1 + 1e-10 per convolution. L1 enters Arb exactly, from its binary value, and the final sign test is in Arb.
Single-thread BLAS is required.
The gate z1 >= 8e9 is asserted in Arb inside bad_up, on every interval.
PRESIEVE: a bad n = N - p must have N - n = p prime, so p n' avoids N mod q for each q | 1155 with q not dividing N. For each
cofactor n' the allowed classes of p mod 1155 number A(N) = prod_{q | 1155, q not | N} (q - 2) * prod_{q | 1155, q | N} (q - 1), so by
(I2) (q = 1155 <= 10^4: c_theta <= 1/840, x >= 8e9; gate log z1 >= log 8e9) sum_{a allowed} theta(x; 1155, a) <= (A(N)/480) x (1 + 480/(840 log x)).
Since U_N >= U_low prod_{q | N, q | 1155} (q-1)/(q-2) and A(N)/480 = G prod_{q | N} (q-1)/(q-2) with G := prod_{q | 1155} (q-2)/(q-1) = 9/32,
the ratio bad/main improves uniformly in N by G: the a = 1 coefficient c_theta becomes G (1 + 480/(840 log z1)).
Usage: cert_qrhB3p.py K Lstart Lend ratio delta H
"""
import os, sys, math, json
assert os.environ.get("OPENBLAS_NUM_THREADS") == "1" and os.environ.get("OMP_NUM_THREADS") == "1", "set thread limits to 1"
sys.dont_write_bytecode = True
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "author"))  # the author's helper module author/cert_qrh_range.py, imported read-only
from fractions import Fraction
import numpy as np
from mpmath import mpf
import cert_qrh_range as Q
from flint import arb, ctx
ctx.prec = 200

K = int(sys.argv[1]); Lstart = float(sys.argv[2]); Lend = float(sys.argv[3]); ratio = float(sys.argv[4])
delta_s = sys.argv[5]; H = int(sys.argv[6])
delta = arb(delta_s)
WMAX = 60
n_max = (WMAX - 1) * H
cvec = 1.0 / (H + np.arange(n_max + 1, dtype=np.float64)) * (1 + 1e-15)
Dtab = {1: cvec.copy()}
for m in range(2, WMAX):
    Dtab[m] = np.convolve(Dtab[m - 1], cvec)[: n_max + 1] * (1 + 1e-10)
def Jup(m, w):                      # w integer
    if m == 0: return 1.0
    if w < m: return 0.0
    return float(np.sum(Dtab[m][: (w - m) * H + 1])) * (1 + 1e-10)
CTH = 1 + arb(1) / 36260
GPRE = arb(9) / 32
def Dl(l): return 1 / (10 * l ** 2) + arb(4) / (15 * l ** 3)
def exact_mpf_tuple(t):                                   # [Fix 2] exact conversion of an mpmath mpf tuple to Arb
    sign, mantissa, exponent, bitcount = t
    return arb(-mantissa if sign else mantissa) * arb(2) ** exponent

def bad_up(L1, M, K):
    """upper bound (arb point) for bad'''/N on an interval starting at L1, parameters M, K (needs z1 >= 10372)"""
    L1a = arb(L1); l1 = L1a / M; z1 = l1.exp()          # [Fix 1] exact binary endpoint, as margin_iv uses
    eps = arb((2 * Dl(l1) + 1 / z1).upper())
    W = M - 1
    mu = {0: arb(1)}
    iz = arb((1 / (z1 - 1)).upper())
    for j in range(1, W):
        Nj = sum((arb(math.comb(j, i)) * eps ** i * arb(Jup(j - i, W)) for i in range(0, j + 1)), arb(0))
        # non-squarefree cofactors n' = q^2 n'' (q >= z): sum 1/n' <= (1/(z-1)) mu_{j-2}  (recursive; n'' z-rough, Omega = j-2, n'' <= z^W)
        mu[j] = arb((Nj / arb(math.factorial(j)) + (iz * mu[j - 2] if j >= 2 else arb(0))).upper())
    lz = l1
    def Az(a):
        return arb((a * (lz / z1 ** a + lz / ((a - 1) * z1 ** (a - 1)) + 1 / ((a - 1) ** 2 * z1 ** (a - 1)))).upper())
    assert (lz - arb(8e9).log()) > 0, "gate: z1 >= 8e9 (BMOR x_theta)"
    cpre = arb((GPRE * (1 + arb(480) / (840 * lz))).upper())
    Phi = arb(0)
    for k in range(K + 1, M):
        Phi += cpre * mu[k - 1] + sum((Az(a) * mu[k - a] for a in range(2, k + 1)), arb(0))
    # + 4 e^{-L1}: the at most four n = N - q with q in {3,5,7,11} (where N - n = q itself is the prime, outside the presieve argument)
    return arb(((-delta * L1a).exp() + 4 * (-L1a).exp() + Phi / ((1 - delta) * L1a)).upper())
_SM = {}
def bad_float(L1, M, K):
    # search-only approximation (the certificate uses bad_up): drops the eps- and A_a-terms, which are < 1e-6 relative here
    if M - 1 >= WMAX: return float('inf')
    key = (M, K)
    if key not in _SM:
        _SM[key] = sum(Jup(k - 1, M - 1) / math.factorial(k - 1) for k in range(K + 1, M))
    d = float(delta_s)
    return math.exp(-d * L1) + (9 / 32) * (1 + 480 / (840 * L1 / M)) * _SM[key] / ((1 - d) * L1) * 1.001
rows = []; fails = 0; L1 = Lstart; minrel = None
while L1 < Lend:
    L2 = min(Lend, L1 * ratio)
    cands = []
    for M in range(16, WMAX):
        if L1 / M < math.log(10372) + 0.01: break
        for lam in (0.1, 1, 10, 100):
            r = Q.margin_float(L1, L2, M, lam, 2.0, M - 1)
            if r is None or r[0] <= 0: continue
            mrg = r[0] - (bad_float(L1, M, K) if K < M - 1 else 0.0)
            if mrg > 0: cands.append((mrg / r[1], M, lam))
    cands.sort(reverse=True)
    best = None; m = None; info = None
    for cand in cands:                     # take the best candidate that passes every Lemma T gate in interval arithmetic
        m_, info_ = Q.margin_iv(L1, L2, cand[1], cand[2], Fraction(2), cand[1] - 1)
        if m_ is not None:
            best, m, info = cand, m_, info_
            break
    if best is None:
        fails += 1; rows.append([L1, L2, None]); print("NO CANDIDATE PASSES (float or gates)", L1, L2, len(cands), flush=True)
    else:
        rel, M, lam = best
        if m is None:
            fails += 1; rows.append([L1, L2, M, lam, "gate", info]); print("GATE FAIL", L1, L2, info, flush=True)
        else:
            assert (arb(L1) / M - arb(10372).log()) > 0, "gate: z1 >= 10372"          # [Fix 1] Arb-only gate
            b = bad_up(L1, M, K) if K < M - 1 else arb(0)
            marg = exact_mpf_tuple(m._mpi_[0]) - b              # [Fix 2] the final sign test stays in interval arithmetic
            ok = marg > 0
            if not ok: fails += 1; print("CERT FAIL", L1, L2, M, lam, flush=True)
            relm = float(marg / exact_mpf_tuple(info[0]._mpi_[1]))   # display/search summary only
            minrel = relm if minrel is None else min(minrel, relm)
            rows.append([L1, L2, M, lam, relm, bool(ok)])
    L1 = L2
assert rows[0][0] == Lstart and all(rows[i][1] == rows[i + 1][0] for i in range(len(rows) - 1)) and rows[-1][1] == Lend
print(f"Theorem B sieve range with R''' + 1155-presieve: K = {K} on [{Lstart}, {Lend}] ratio {ratio}, delta {delta_s}, H {H}: {len(rows)} intervals, "
      f"failures {fails}; coverage asserted; least relative margin {minrel}", flush=True)
json.dump({"K": K, "Lstart": Lstart, "Lend": Lend, "ratio": ratio, "delta": delta_s, "H": H, "rows": rows, "failures": fails},
          open(f"cert_qrhB3p_K{K}_{int(Lstart)}_{int(Lend)}.json", "w"))